The Evolution and Origin of Ionized Gas Velocity Dispersion from $z\sim2.6$ to $z\sim0.6$ with KMOS$^{\rm 3D}$

TL;DR

This study analyzes the evolution of ionized gas velocity dispersion in star-forming disk galaxies from redshift 2.6 to 0.6, revealing a decrease over cosmic time and insights into the physical mechanisms driving turbulence.

Contribution

It provides the first comprehensive analysis of ionized gas velocity dispersion evolution over this redshift range using the KMOS$^{3D}$ survey, incorporating advanced modeling to account for observational effects.

Findings

Velocity dispersion decreases from ~45 km/s at z~2.3 to ~30 km/s at z~0.9.

Significant intrinsic scatter suggests dynamic processes like minor mergers influence turbulence.

Ionized and atomic+molecular dispersions evolve similarly, with ionized being slightly higher.

Abstract

We present the $0.6<z<2.6$ evolution of the ionized gas velocity dispersion in 175 star-forming disk galaxies based on data from the full KMOS$^{\rm 3D}$ integral field spectroscopic survey. In a forward-modelling Bayesian framework including instrumental effects and beam-smearing, we fit simultaneously the observed galaxy velocity and velocity dispersion along the kinematic major axis to derive the intrinsic velocity dispersion $\sigma_0$. We find a reduction of the average intrinsic velocity dispersion of disk galaxies as a function of cosmic time, from $\sigma_0\sim45$ km s$^{-1}$ at $z\sim2.3$ to $\sigma_0\sim30$ km s$^{-1}$ at $z\sim0.9$. There is substantial intrinsic scatter ($\sigma_{\sigma_0, {\rm int}}\approx10$ km s$^{-1}$) around the best-fit $\sigma_0-z$-relation beyond what can be accounted for from the typical measurement uncertainties ($\delta\sigma_0\approx12$ km s$^{-1}$), independent of other identifiable galaxy parameters. This potentially suggests a dynamic mechanism such as minor mergers or variation in accretion being responsible for the scatter. Putting our data into the broader literature context, we find that ionized and atomic+molecular velocity dispersions evolve similarly with redshift, with the ionized gas dispersion being $\sim10-15$ km s$^{-1}$ higher on average. We investigate the physical driver of the on average elevated velocity dispersions at higher redshift, and find that our galaxies are at most marginally Toomre-stable, suggesting that their turbulent velocities are powered by gravitational instabilities, while stellar feedback as a driver alone is insufficient. This picture is supported through comparison with a state-of-the-art analytical model of galaxy evolution.

Tables

Table 1. Table 1: Comparison of modelling results from the three setups (S1, S2, S3) described in Section  3.2 .

comparison	quantity	mean	std. dev.
S1 – S2	$\Delta {\sigma }_0$ [km s-1]	0.9	6.0
	$\Delta$log($M_{\mathrm{tot}})$ [dex of $M_{\odot }$]	-0.01	0.03
S1 – S3	$\Delta {\sigma }_0$ [km s-1]	0.5	7.4
	$\Delta$log($M_{\mathrm{tot}})$ [dex of $M_{\odot }$]	-0.06	0.11
S2 – S3	$\Delta {\sigma }_0$ [km s-1]	-1.4	5.3
	$\Delta$log($M_{\mathrm{tot}})$ [dex of $M_{\odot }$]	-0.04	0.10
S3: $H-$band – H$\alpha$aaComparison of the fiducial $R_e$ and $B/T$ as derived from the stellar light $H-$band images (see Section 2) to the modelling results from setup 3, where we fit for $R_e$ and $B/T$ as detailed in Section 3.2.	$\Delta R_e$ [kpc]	-0.6	1.0
	$\Delta B/T$	0.03	0.14

Table 2. Table 2: Spearman rank correlation coefficients, ρ S subscript 𝜌 S \rho_{\rm S} , and their significance σ ρ subscript 𝜎 𝜌 \sigma_{\rho} , between σ 0 subscript 𝜎 0 \sigma_{0}  and respectively R e subscript 𝑅 𝑒 R_{e} , R max / R e subscript 𝑅 max subscript 𝑅 𝑒 R_{\rm max}/R_{e} , σ instrumental subscript 𝜎 instrumental \sigma_{\rm instrumental} , and b / a 𝑏 𝑎 b/a .

quantity	${\rho }_S$	${\sigma }_{\rho }$
$R_e$	0.01	1.2
$R_{\max }/R_e$	-0.05	0.7
${\sigma }_{\mathrm{instrumental}}$	-0.07	0.9
$b/a$	-0.04	0.5

Table 3. Table 3: Results from the linear regression fits of the form σ 0 / km ​ s = a + b ⋅ z + c subscript 𝜎 0 km s 𝑎 ⋅ 𝑏 𝑧 𝑐 \sigma_{0}/{\rm km\,s}=a+b\cdot z+c for our kinematic sample, where a 𝑎 a and b 𝑏 b are the regression coefficients, and c 𝑐 c is the intrinsic random scatter about the regression (see Kelly, 2007 ) . For each parameter a 𝑎 a , b 𝑏 b , the standard deviation of c 𝑐 c , and the derived linear correlation coefficient l corr subscript 𝑙 corr l_{\rm corr} between σ 0 subscript 𝜎 0 \sigma_{0} and z 𝑧 z , we list the median together with the standard deviation of the posterior distribution. For each redshift slice we list the best-fit σ 0 subscript 𝜎 0 \sigma_{0} value corresponding to these medians.

sample	N	$a$	$b$	${\sigma }_c$	$l_{\mathrm{corr}}$	${\sigma }_0$ at $z\sim 0.9$	${\sigma }_0$ at $z\sim 1.5$	${\sigma }_0$ at $z\sim 2.3$
						[km s-1]	[km s-1]	[km s-1]
including upper	175	$21.1\pm 3.0$	$11.3\pm 2.0$	$11.3\pm 1.1$	$0.51\pm 0.08$	31.1	38.3	46.7
limits
excluding upper	147	$26.2\pm 3.1$	$9.2\pm 2.1$	$10.4\pm 1.1$	$0.46\pm 0.09$	34.3	40.3	47.1
limits (robust)
using formal	175	$17.2\pm 3.2$	$12.7\pm 2.2$	$13.2\pm 1.1$	$0.49\pm 0.07$	28.4	36.7	46.1
best-fit ${\sigma }_0$

Table 4. Table 4: Variances of σ 0 subscript 𝜎 0 \sigma_{0}  around the robust best-fit relation: observed variance VAR obs , variance due to measurement uncertainties VAR δ ​ σ 0 𝛿 subscript 𝜎 0 {}_{\delta\sigma_{0}} , and intrinsic variance, VAR int .

measure	$z\sim 0.9$	$z\sim 1.5$	$z\sim 2.3$
VARobs [km2 s-2]	171	208	357	237
VAR${}_{\delta {\sigma }_0}$ [km2 s-2]	87	130	194	133
VARint [km2 s-2]	85	78	163	104
VARint/VARobs	0.50	0.38	0.46	0.44

Table 5. Table 5: Literature data of the 0 < z < 4 0 𝑧 4 0<z<4 velocity dispersion measurements shown in Figure  6 .

Sample/References	$z$	Primary Target Line	N included	Instrument/Method	BS correction	Comment on Selection
KDS	$3.8-3.1$	[Oiii]	14	KMOS/IFU	(1)	their ‘RD’
MOSDEF	$3.8-1.4$	H$\alpha$, H$\beta$, [Oiii]	108	MOSFIRE/slit	(2)	their ‘resolved/aligned’
AMAZE-LSD	$3.7-3.1$	[Oiii]	11	SINFONI/IFU	(1)	their ‘rotating’
Livermore et al. (2015)	$3.7-1.3$	H$\alpha$, H$\beta$	8	SINFONI+NIFS+OSIRIS/IFU+lensing	(3)	their ‘Disk’ with $v_{\mathrm{rot}}/{\sigma }_0>1$
KMOS3D	$2.6-0.6$	H$\alpha$	175	KMOS/IFU	(2)	see Section 3.3
SINS/zC-SINF	$2.5-1.4$	H$\alpha$	25	SINFONI/IFU+AO	(2)	see Section 3.3
SIGMA	$2.5-1.3$	H$\alpha$, [Oiii]	49	MOSFIRE/slit	(4)
Genzel et al. (2017)	$2.4-0.9$	H$\alpha$	6	KMOS+SINFONI/IFU+AO	(2)	see Section 3.3
MASSIV	$1.6-0.9$	H$\alpha$, [Oiii]	53	SINFONI/IFU	(1)
DEEP2	$1.2-0.1$	H$\alpha$, H$\beta$, [Oii], [Oiii]	544	DEIMOS/slit	(4)
KROSS	$1.0-0.8$	H$\alpha$	171	KMOS/IFU	(4)	their ‘sigma0_flag=O’ with $v_{2.2}/{\sigma }_0>1$
DYNAMO	$\sim 0.1$	H$\alpha$	25	SPIRAL+WiFeS/IFU	(1)	their ‘RD’
GHASP	local	H$\alpha$	153	scanning Fabry-Perot	(1)
Übler et al. (2018)	1.4	H$\alpha$, CO(3–2)	1	LUCI/slit+NOEMA/interferometry	(2)	see Section 3.3
Swinbank et al. (2011)	2.4	CO(6–5), CO(1–0)	1	IRAM+EVLA/interferometry	(3)
PHIBSS	$1.5-0.7$	CO(3–2)	7	IRAM+NOEMA/interferometry	(4)	see Section 3.3
HERACLES	local	CO(2–1)	13	IRAM/single dish	(5)
EDGE-CALIFA	local	CO(1–0)	17	CARMA/interferometry	(1)
THINGS	local	Hi	35	VLA/interferometry	(5)
Dib et al. (2006)	local	Hi	13	literature compilation

Table 6. Table 6: Results and standard deviations from the robust least-squares linear regression fits of the form σ 0 / km ​ s = a + b ⋅ z subscript 𝜎 0 km s 𝑎 ⋅ 𝑏 𝑧 \sigma_{0}/{\rm km\,s}=a+b\cdot z to the data sets shown in the bottom panel of Figure  6 .

sample	$a$ [km s-1]	$b$ [km s-1]
ionized gas (best averages)	23.3$\pm$4.9	9.8$\pm$3.5
… fixing slope to atomic+molecular	22.8	11.0 (fixed)
KMOS3D incl. upper limits (Table 3)	21.1$\pm 3.0$	11.3$\pm$2.0
atomic+molecular gas	10.9$\pm$0.6	11.0$\pm$2.0

Table 7. Table 7: Spearman rank correlation coefficients ρ S subscript 𝜌 S \rho_{\rm S} , and their significance σ ρ subscript 𝜎 𝜌 \sigma_{\rho} , between σ 0 subscript 𝜎 0 \sigma_{0} and different galaxy properties for our robust sample before and after accounting for the redshift dependence of σ 0 subscript 𝜎 0 \sigma_{0}  .

	${\sigma }_0\left(z\right)$		${\sigma }_{0,\mathrm{norm}}$
quantity	${\rho }_S$	${\sigma }_{\rho }$	${\rho }_S$	${\sigma }_{\rho }$
z	0.33	4.0	–	–
SFR	0.38	4.6	0.18	2.1
SFRHα	0.36	4.4	0.14	1.7
${\Sigma }_{\mathrm{SFR}}$	0.32	3.9	0.06	1.0
${\Sigma }_{\mathrm{SFR},H\alpha }$	0.30	3.7	0.08	0.9
$M_{\mathrm{gas}}$	0.38	4.6	0.19	2.3
${\Sigma }_{\mathrm{gas}}$	0.31	3.8	0.07	0.9
$M_{\ast }$	0.26	3.1	0.20	2.4
${\Sigma }_{\ast }$	0.26	3.1	0.14	1.6
$M_{\mathrm{bar}}$	0.32	3.9	0.20	2.4
${\Sigma }_{\mathrm{bar}}$	0.30	3.6	0.12	1.5
$\Delta$MS	–	–	0.15	1.8
$\Delta$MR	–	–	-0.05	0.6

Table 8. Table 8: Results from the linear regression fits of the form σ 0 / km ​ s = a + b ⋅ log ⁡ ( 1 + z ) + c subscript 𝜎 0 km s 𝑎 ⋅ 𝑏 1 𝑧 𝑐 \sigma_{0}/{\rm km\,s}=a+b\cdot\log(1+z)+c for our kinematic sample, where a 𝑎 a and b 𝑏 b are the regression coefficients, and c 𝑐 c is the intrinsic random scatter about the regression (see Kelly, 2007 ) . For each parameter a 𝑎 a , b 𝑏 b , the standard deviation of c 𝑐 c , and the derived linear correlation coefficient l corr subscript 𝑙 corr l_{\rm corr} between σ 0 subscript 𝜎 0 \sigma_{0} and z 𝑧 z , we list the median together with the standard deviation of the posterior distribution. For each redshift slice we list the best-fit σ 0 subscript 𝜎 0 \sigma_{0} value corresponding to these medians.

sample	N	$a$	$b$	${\sigma }_c$	$l_{\mathrm{corr}}$	${\sigma }_0$ at $z\sim 0.9$	${\sigma }_0$ at $z\sim 1.5$	${\sigma }_0$ at $z\sim 2.3$
						[km s-1]	[km s-1]	[km s-1]
including upper	175	$19.4\pm 3.8$	$52.6\pm 10.1$	$9.8\pm 1.1$	$0.48\pm 0.08$	33.8	40.6	46.5
limits
excluding upper	147	$19.5\pm 4.5$	$53.1\pm 11.4$	$10.4\pm 1.1$	$0.46\pm 0.09$	34.1	40.9	46.9
limits (robust)
using formal	175	$19.5\pm 4.5$	$53.1\pm 11.3$	$10.4\pm 1.1$	$0.46\pm 0.09$	34.1	40.9	46.8
best-fit ${\sigma }_0$

Table 9. Table 9: Results from the linear regression fits of the form log ⁡ ( σ 0 / km ​ s ) = a + b ⋅ log ⁡ ( 1 + z ) + c subscript 𝜎 0 km s 𝑎 ⋅ 𝑏 1 𝑧 𝑐 \log(\sigma_{0}/{\rm km\,s})=a+b\cdot\log(1+z)+c for our kinematic sample, where a 𝑎 a and b 𝑏 b are the regression coefficients, and c 𝑐 c is the intrinsic random scatter about the regression (see Kelly, 2007 ) . For each parameter a 𝑎 a , b 𝑏 b , the standard deviation of c 𝑐 c , and the derived linear correlation coefficient l corr subscript 𝑙 corr l_{\rm corr} between σ 0 subscript 𝜎 0 \sigma_{0} and z 𝑧 z , we list the median together with the standard deviation of the posterior distribution. For each redshift slice we list the best-fit σ 0 subscript 𝜎 0 \sigma_{0} value corresponding to these medians.

sample	N	$a$	$b$	${\sigma }_c$	$l_{\mathrm{corr}}$	${\sigma }_0$ at $z\sim 0.9$	${\sigma }_0$ at $z\sim 1.5$	${\sigma }_0$ at $z\sim 2.3$
						[km s-1]	[km s-1]	[km s-1]
including upper	175	$1.29\pm 0.05$	$0.77\pm 0.13$	$0.12\pm 0.01$	$0.54\pm 0.08$	31.8	39.9	48.6
limits
excluding upper	147	$1.38\pm 0.05$	$0.61\pm 0.13$	$0.11\pm 0.01$	$0.49\pm 0.09$	35.3	42.2	49.4
limits (robust)
using formal	175	$1.30\pm 0.05$	$0.75\pm 0.13$	$0.12\pm 0.01$	$0.53\pm 0.07$	32.3	40.3	49.0
best-fit ${\sigma }_0$

Equations

$${\sigma }_0[\mathrm{km}{\mathrm{s}}^{-1}]=(21.1\pm 3.0)+(11.3\pm 2.0)\cdot z.$$

$${\sigma }_{0,\mathrm{norm}}={\sigma }_0-(a+b\cdot z),$$

$${\mathrm{VAR}}_{\mathrm{int}}(z)={\mathrm{VAR}}_{\mathrm{obs}}(z)-{\mathrm{VAR}}_{\delta {\sigma }_0}(z),$$

$$\sigma_{z}=5.5~{}{\rm kms^{-1}}\cdot\frac{f_{p}}{(1+\chi)^{1/2}}\bigg{(}\frac{\epsilon_{\rm ff}(\rho_{0})}{0.005}\bigg{)}\bigg{(}\frac{p_{*}/m_{*}}{3000~{}{\rm kms^{-1}}}\bigg{)}.$$

$$Q_{\mathrm{gas}}$$

$$\mathrm{SFR}$$

Full text

galaxies: evolution – galaxies: high-redshift – galaxies: ISM – galaxies: kinematics and dynamics

1 Introduction

Extragalactic surveys over the last decades have produced thousands of spectrally and spatially resolved observations of galaxies from the present day out to $z\sim 4$. For massive galaxies on the star-forming main sequence, these efforts resulted in two main findings regarding their kinematic evolution: (i) already by $z\sim 2$, the majority of star-forming galaxies (SFGs) show ordered rotation, and (ii) their velocity dispersions are higher by factors of 2-5 compared to local SFGs (Labbé et al., 2003; Förster Schreiber et al., 2006, 2009, 2018; Genzel et al., 2006, 2008, 2014; Cresci et al., 2009; Epinat et al., 2009, 2012; Law et al., 2009; Jones et al., 2010; Gnerucci et al., 2011; Wisnioski et al., 2011, 2015, in prep.; Miller et al., 2012; Swinbank et al., 2012b; Stott et al., 2016; Simons et al., 2017). The redshift evolution of the ionized gas velocity dispersion has captured a lot of attention through its potential to constrain feedback and star formation models (Förster Schreiber et al., 2006; Genzel et al., 2006, 2008, 2011; Weiner et al., 2006; Kassin et al., 2007, 2012; Epinat et al., 2009, 2012; Law et al., 2009; Lehnert et al., 2009, 2013; Gnerucci et al., 2011; Wisnioski et al., 2012, 2015; Swinbank et al., 2012a; Newman et al., 2013; Simons et al., 2016, 2017; Turner et al., 2017; Mason et al., 2017; Zhou et al., 2017; Johnson et al., 2018; Girard et al., 2018).

Starting from small scales in the Milky Way, the velocity dispersion in molecular clouds is proportional to cloud size and mass, in a way that suggests molecular clouds are turbulent, with kinetic and gravitational energy being in near equipartition (Larson, 1981; McKee & Ostriker, 2007; Heyer & Dame, 2015, and references therein). However, the lack of dependence of the turbulence level on factors such as environment or local star formation activity points towards larger scale drivers (Heyer & Brunt, 2004; Brunt et al., 2009; but see Heyer & Dame, 2015 for extreme environments).

In nearby galaxies, velocity dispersions of atomic gas are $\sigma_{\rm HI}\approx 10-12$ km s*-1* on scales of $\sim 100$ pc (Dib et al., 2006; Tamburro et al., 2009; Ianjamasimanana et al., 2012; Fukui et al., 2009; Caldú-Primo et al., 2013; Mogotsi et al., 2016; Koch et al., 2019). Molecular gas velocity dispersions are typically lower, with reported ratios in the range $\sigma_{\rm CO}/\sigma_{\rm HI}\approx 0.3-1$ (Tamburro et al., 2009; Ianjamasimanana et al., 2012; Fukui et al., 2009; Wong et al., 2009; Caldú-Primo et al., 2013; Druard et al., 2014; Mogotsi et al., 2016; Levy et al., 2018; Koch et al., 2019). Ionized gas velocity dispersions are substantially higher, with $\sigma_{\rm H\alpha}\approx 24$ km s*-1* (Epinat et al., 2010).

At high redshift, most measurements of gas velocity dispersion are based on ionized gas, which is accessible from the ground in the near-infrared through strong rest-frame optical lines. Typical values are $\sigma=25-100$ km s*-1* for disk galaxies. It is smore challenging to measure accurate velocity dispersions at high redshift because of the combined effects of beam-smearing and limited instrumental spectral resolution (see Davies et al., 2011). The former can be corrected for instance by using the velocity field and the spatial resolution of the observations to create a beam-smearing map (e.g. Green et al., 2010; Gnerucci et al., 2011; Epinat et al., 2012), through model-based look-up tables (e.g. Burkert et al., 2016; Johnson et al., 2018), or through forward-modelling (e.g. Cresci et al., 2009; Genzel et al., 2011; Di Teodoro et al., 2016; Wuyts et al., 2016; Varidel et al., 2019). Typical spectral resolutions of near-infrared spectroscopic observations at $z\sim 1-3$ correspond to velocity dispersions of $\sigma_{\rm instrumental}\approx 30-40$ km s*-1*. However, depending on the signal-to-noise ratio (S/N), it is possible to recover velocity dispersions through forward-modelling down to 1/3 of the instrumental resolution.

It is well established that the galactic gas velocity dispersion is correlated with redshift (e.g. review by Glazebrook, 2013), but the physical processes responsible for driving and maintaining the dispersions are still debated. It has been shown theoretically that constant energy input is necessary to maintain turbulence in the interstellar medium (ISM) because it will otherwise decay within a few Myr (e.g. Mac Low et al., 1998; Stone et al., 1998). A number of potential drivers has been identified, with two main classes: (i) the conversion of kinetic energy through stellar feedback in the form of winds, expanding Hii regions, and supernovae, and (ii) the release of gravitational energy through clump formation, radial flows within the disk, and accretion from the cosmic web. Other possible sources include effects of galactic rotation, fluid instabilities, and galaxy interactions (see Elmegreen & Scalo, 2004, for a review). Generally, the different scales on which the proposed mechanisms operate present a challenge to simulations (see Naab & Ostriker, 2017, for a review).

In this paper, we investigate the intrinsic velocity dispersion of the ionized gas phase in rotation-dimonated, star-forming galaxies from our KMOS3D survey at $0.6<z<2.6$. In Section 2 we briefly present the KMOS3D data set. Our modelling and sample selection is discussed in Section 3. In Section 4 we investigate the evolution of the intrinsic velocity dispersion with redshift and put it into the broader context of multi-phase literature values from $z=4$ to $z=0$. In Section 5 we discuss possible drivers of turbulence, particularly gravitational instabilities and stellar feedback, and compare our data to a state-of-the-art analytical model by Krumholz et al. (2018). We conclude our study in Section 6.

Throughout, we adopt a Chabrier (2003) initial mass function and a flat $\Lambda$CDM cosmology with $H_{0}=70{\rm\,km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\Lambda}=0.7$, and $\Omega_{m}=0.3$.

2 The KMOS3D survey

Our study is based on data from the KMOS3D survey, targeting the H$\alpha$ line emission of primarily main sequence galaxies in three redshift bins centered at $z\sim 0.9$, $z\sim 1.5$, and $z\sim 2.3$. The survey is presented by Wisnioski et al. (2015) and Wisnioski et al., in prep., to which we refer the reader for details. Below, we summarize its main characteristics.

The KMOS3D galaxies were selected from the 3D-HST survey (Brammer et al., 2012; Skelton et al., 2014; Momcheva et al., 2016), providing optical-to-8$\mu$m photometry and, importantly, secure spectroscopic or grism redshifts, so that bright OH skylines at the location of the H$\alpha$ line emission could be avoided. In addition, high-resolution imaging for all galaxies is available through CANDELS (Grogin et al., 2011; Koekemoer et al., 2011; van der Wel et al., 2012), and further multi-wavelength coverage through photometry from Spitzer/MIPS and Herschel/PACS (Lutz et al., 2011; Magnelli et al., 2013; Whitaker et al., 2014, and references therein).

For the KMOS3D survey, we selected galaxies with stellar masses log$(M_{*}/M_{\odot})>9$ and $K_{\rm AB}\lesssim 23$. The selection aimed to provide a homogeneous coverage of the star formation main sequence across stellar mass in the three redshift slices, thus ensuring near equal statistical coverage up to the highest masses. In addition, KMOS3D also extends below the main sequence regime where galaxies are ‘quiescent’, and it contains starburst outliers above the main sequence.

Stellar masses were derived following Wuyts et al. (2011), by fitting the broad- and medium-band optical-to-mid-infrared spectral energy distribution with Bruzual & Charlot (2003) stellar population synthesis models, adopting a Calzetti et al. (2000) extinction law, solar metallicity, and a range of star formation histories. Gas mass measurements are not available for most of our galaxies. We exploit the scaling relation by Tacconi et al. (2018) which depends on redshift, offset from the main sequence, and stellar mass, with the main sequence prescription by Whitaker et al. (2014), to estimate the molecular gas masses ($M_{\rm gas}$) of our galaxies. We don’t account for atomic gas in this study. The derivation of star-formation rates (SFRs) followed the ladder of SFR indicators as described by Wuyts et al. (2011).

Structural properties such as the axis ratio $q=b/a$, the disk effective radius $R_{e}$, and the bulge-to-total stellar mass fraction $B/T$ are based on two-dimensional S rsic models to the stellar light distribution high-resolution $H-$band images from HST observations (van der Wel et al., 2012; Lang et al., 2014). For the effective radius we apply a color correction following van der Wel et al. (2014a).

The survey was conducted during the years 2013 to 2018 with the multiplexing near-infrared integral field spectrograph KMOS (Sharples et al., 2004, 2013) at the Very Large Telescope. The full KMOS3D sample consists of 740 targeted galaxies (Wisnioski et al., in prep.).

3 Dynamical Modelling and Sample Selection

We constrain the intrinsic velocity dispersions by forward-modelling the observed one-dimensional velocity and velocity dispersion profiles extracted from the data cubes. For this work, we use the two-dimensional kinematic information to determine the kinematic major axis, and to distinguish rotation-dominated, dispersion-dominated, and disturbed systems. The full kinematic information on the motions of stars or gas in the plane of a rotating disk can be extracted along its kinematic major axis. Modelling the one-dimensional kinematics instead of the two- or three-dimensional data increases the S/N of our measurements, and thus allows us to study a larger sample of galaxies with reliable modelling. We have verified that this has only a minor impact on the derived dynamical parameters, with an average, non-systematic difference of one-dimensional *vs. *two-dimensional intrinsic velocity dispersion of $\sim 5-10$%.

3.1 One-Dimensional Kinematic Profiles

We derive two-dimensional projected H$\alpha$ velocity and velocity dispersion fields for all KMOS3D galaxies using linefit (Davies et al., 2009, 2011; Förster Schreiber et al., 2009), a code that takes into account the instrument line spread function and fits a Gaussian model for each spaxel of the reduced data cube after continuum subtraction. From these maps we exclude spaxels with $S/N\leq 2$, uncertainties on the velocity or velocity dispersion of $\geq 100$ km s*-1*, as well as off-source fits to noise features. We determine the maximum and minimum of the velocity map through a weighted average of either the 5 % of spaxels of both the highest and lowest velocity values for galaxies with $\geq 50$ suitable spaxels, or otherwise of the five spaxels with highest and lowest velocities. The kinematic major axis is defined as the line going through the maximum and minimum of the velocity map. The kinematic center is defined as the midpoint on the kinematics major axis connecting the maximum and minimum of the velocity map. This method follows the procedures outlined by Wisnioski et al. (2015), and in the KMOS3D data release and final survey paper by Wisnioski et al., in prep.

Along the kinematic major axis, we then extract spectra in circular apertures of diameter $2\times$FWHM of the model-independent point spread function (PSF) associated with each individual galaxy. Here, the flux from all spaxels within an aperture is integrated to create a single spectrum. For the dynamical modelling of our galaxies (see Section 3.2), we repeat this same procedure for each iteration of the model fitting to properly account for any effects related to this integration process. We consider a galaxy to be spatially resolved if we can measure its kinematics over a total of at least $3\times$PSFFWHM along the kinematic major axis. We fit the H$\alpha$ velocity and velocity dispersion from the resulting spectra, providing us with the one-dimensional rotation curve $v_{\rm rot}(r)\cdot\sin(i)$ and dispersion profile $\sigma(r)$, uncorrected for beam-smearing. Uncertainties for each data point are derived using Monte Carlo analysis and have typical values of 6 km s*-1* and 10 km s*-1* for the velocity and dispersion values, respectively.

With this methodology we have successfully extracted kinematic profiles for all 535 KMOS3D H$\alpha-$detected galaxies with secure redshifts.

3.2 Dynamical Modelling with dysmal

As a next step, we consider all galaxies with kinematic profile extractions that are resolved, a total of 456 SFGs. We further exclude targets for which we detect multiple systems within the IFU, and we eliminate merging or potentially interacting systems with larger separations based on projected distances, redshift separations, and mass ratios, as informed through the 3D-HST catalog (Mendel et al., in prep.). Galaxies that are strongly contaminated by sky features, have prominent broad line regions, or have very strong outflows affecting the recovery of the galaxies’ velocity and dispersion, are also excluded. This results in 356 galaxies.

We exploit the dynamical fitting code dysmal (Cresci et al., 2009; Davies et al., 2011; Wuyts et al., 2016; Übler et al., 2018) to model our galaxies. dysmal is a forward-modeling code that allows for a flexible number of components (disk, bulge, halo, etc.) and free parameters. It accounts consistently for finite scale heights and flattened spheroidal potentials (Noordermeer, 2008), and it includes the effects of pressure support on the rotation velocity. It also accounts for the instrument line spread function, and for beam-smearing effects by convolving with the two-dimensional PSF of each galaxy.

For our modelling, we assume a velocity dispersion that is isotropic and constant throughout the disk, motivated by deep adaptive optics imaging spectroscopy on kpc scales of 35 $z=1-2.6$ SFGs in the SINS/zC-SINF sample (Genzel et al., 2006, 2008, 2011, 2017; Cresci et al., 2009; Förster Schreiber et al., 2018, see also Section 5.2). We note that for nearby galaxies radially declining velocity dispersions have been observed for atomic and molecular gas (van der Kruit & Freeman, 1984; Dickey et al., 1990; Boulanger & Viallefond, 1992; Kamphuis & Sancisi, 1993; Meurer et al., 1996; Petric & Rupen, 2007; Tamburro et al., 2009; Wilson et al., 2011; Caldú-Primo et al., 2013; Mogotsi et al., 2016; Sun et al., 2018; Koch et al., 2019), where the velocity dispersion usually reaches a constant level only in the disk outskirts. The observed radial changes in velocity dispersion are however rarely larger than $10-20$ km s*-1*, and such variations on small scales are likely washed out through the coarser spatial resolution of typical high$-z$ observations (but see Section 5.2 for a high-resolution example).

We create a three-dimensional mass model of each galaxy consisting of an exponential disk with the effective radius $R_{e}$ adopted from the $H-$band measurements, with ratio of scale height to scale length $q_{0}=0.25$, and with a central bulge ($R_{e,{\rm bulge}}=1$ kpc, Sérsic index $n_{\rm S,bulge}=4$, e.g. Lang et al., 2014; Tacchella et al., 2015b). The value of $q_{0}=0.25$ is motivated by the fall-off in the $q=b/a$ distribution of SFGs at the mass and redshift of our sample (van der Wel et al., 2014b). For galaxies without an $H-$band based measurement of the bulge mass (see Section 2; ca. 30%) we use average values of $B/T=[0.25;0.35;0.45;0.5]$ for total stellar masses of log$(M_{\star}/M_{\odot})=[<10.8;10.8-11;11-11.4;>11.4]$, following Lang et al. (2014). We fix the physical size of the bulge because individual measurements of $R_{e,{\rm bulge}}$ are very uncertain, in contrast to measurements of $B/T$ (see Tacchella et al., 2015b). In a population-averaged sense, however, $R_{e,{\rm bulge}}=1$ kpc is a robust choice (see Lang et al., 2014). We calculate the galaxy inclination $i$ as $\cos(i)=[(q^{2}-q_{0}^{2})/(1-q_{0}^{2})]^{1/2}$. The mass model is then rotated to match the galaxy’s observed orientation in space, convolved with the line spread function and the PSF of the observation to take into account beam-smearing, and subsequently pixelated to resemble the spatial sampling of the observation. We approximate the PSF as a two-dimensional Moffat function that has been fitted to the standard star observations associated with each KMOS detector and pointing. For our modelling, we assume that light traces mass.

Using dysmal, we simultaneously fit the one-dimensional velocity and velocity dispersion profiles of our galaxies in observed space. The best-fitting intrinsic rotation velocity, $v_{\rm rot}$, is constrained both through the mass model and the intrinsic velocity dispersion via pressure support. We apply Markov chain Monte Carlo (MCMC) sampling to determine the model likelihood based on comparison to the observed one-dimensional kinematic profiles, and assuming Gaussian measurement noise. To ensure convergence of the MCMC chains, we model each galaxy with 400 walkers, a burn-in phase of 50-100 steps, followed by a running phase of another 50-100 steps (>10 times the maximum autocorrelation time of the individual parameters). For each free parameter, we adopt the median of all model realizations as our best fit value, with asymmetric uncertainties corresponding to the 1$\sigma$ confidence ranges of the one-dimensional marginalized posterior distributions.

In order to recover the intrinsic velocity dispersion as best as possible, we consider a total of three setups with varying free parameters and treatment of the kinematic profiles:

In our first setup, we feed the kinematic profiles obtained as described in Section 3.1, with free parameters being the total dynamical mass in the range log($M_{\rm tot}/M_{\odot})=[9;13]$, and the intrinsic velocity dispersion in the range $\sigma_{0}=[5;300]$ km s*-1*. $M_{\rm tot}$ is the total mass distributed in the three-dimensional disk plus bulge structure necessary to reproduce the observed kinematics. Other parameters are fixed, specifically $i$, $R_{e}$, and $B/T$. 2. 2.

Due to extinction, skyline contamination, and noise limitations, some galaxies display asymmetric kinematic profiles. Therefore, we employ a symmetrization technique in a second setup, where the one-dimensional profiles are folded (for the dispersion profile) or rotated (for the rotation curve) around the kinematic centre, interpolated onto a common grid, and averaged by calculating the mean at each radial grid point to obtain symmetric profiles, with uncertainties added in quadrature. Again, free parameters are $M_{\rm tot}$ and $\sigma_{0}$, allowed to vary within the same ranges as for setup 1. 3. 3.

As noted in Section 2, $R_{e}$ and $B/T$ of our galaxies are derived from $H-$band imaging. It is known that the mass distribution derived from the $H-$band light might differ from the corresponding H$\alpha$ flux profiles (Wuyts et al., 2012; Tacchella et al., 2015a; Nelson et al., 2016, e.g.; Wilman et al., in prep.). In particular the dispersion profiles can be sensitive to the central mass concentration. In the third setup we therefore proceed as in setup 2, but additionally leave the disk effective radius $R_{e}$ and the bulge-to-total fraction $B/T$ as free parameters. For $R_{e}$ we use a truncated Gaussian prior centered on the fiducial value with a standard deviation of 1 kpc, and truncated at $\pm 2.5$ kpc from the peak value, with hard bounds of $R_{e}=[0.1;20]$ kpc. For $B/T$ we use a Gaussian prior centered on the fiducial value with a standard deviation of 0.2 and hard bounds of $B/T=[0;1]$.

Comparing results from the three setups, we generally find good agreement for both the derived intrinsic dispersions and the dynamical masses, as listed in Table 1. For setup 3, the model-derived (mass/H$\alpha$) effective radii are systematically higher compared to the $H-$band measurements by $\sim 0.6$ kpc. For the range of $R_{e}\approx 2-10$ kpc and log$(M_{*}/M_{\odot})\approx 9.2-11.5$ in our kinematic sample, this is agreement with the results by Nelson et al. (2016) and Wilman et al., in prep., who find $R_{e,{\rm H\alpha}}/R_{e,H}\approx 1.1-1.2$ from high-resolution HST observations and from our full KMOS3D sample, respectively. The average agreement between the $H-$band-derived $B/T$ and the model-derived $B/T$ is better, however the model-derived value is likely more realistic for cases with only a grid-based $B/T$.

We tested a fourth setup for a subset of our sample, including not only the bulge and disk components but in addition an NFW halo (Navarro et al., 1996), with a prior on the expected dark matter halo mass (Moster et al., 2018) and the concentration parameter fixed to the theoretically expected value (Dutton & Macciò, 2014). The resulting best-fit velocity dispersions are robust in that they agree within the uncertainties with the results from the other three setups with a standard deviation of 5.9 km s*-1*, and there are no systematic effects. However, the limited field-of-view of KMOS (compared to e.g. SINFONI) together with our typical integration times of $5-9$ h per target constrain our ability to map the faint outskirts of galaxies where the kinematics are most sensitive to additional dynamical components with a different mass distribution. Therefore, we do not include fits from this fourth setup in our final sample.

3.3 The Kinematic sample

We inspect the fits from all three model setups to create our best-fit sample. By default, we choose the fit to setup 1, but if it is bad or poorly constrained, we consider setups 2 and 3 in this order. Galaxies with poor fits in all setups are excluded. With this strategy we stay as closely as possible to the original data, but at the same time do not need to disregard galaxies with one-sided extinction or skyline contamination that otherwise show good data quality, and we can choose fits from setup 3 with a more appropriate mass distribution, if necessary.

Finally, we impose a $v_{\rm rot}/\sigma_{0}\geq 1$ cut to focus on rotation-dominated systems. Here, $v_{\rm rot}$ is the model intrinsic rotation velocity at $1.38~{}R_{e}$, which is the location of the peak of the rotation curve for a Noordermeer disk with $n_{\rm S}=1$. Our final sample consists of 175 galaxies, with 80, 47, and 48 galaxies in the redshift slices $z=0.6-1.1$, $z=1.2-1.7$, and $z=1.9-2.6$. Of those galaxies, 56 % are from setup 1, 31 % from setup 2, and 13 % from setup 3. We show examples of galaxies and their fits from different setup in Appendix A. The averaged uncertainties on our derived $\sigma_{0}$ values cover the range $\delta\sigma_{0}=2-29$ km s*-1*, with $68^{\rm th}$ percentiles of $\delta\sigma_{0}=5-15$ km s*-1*, and mean values in the three redshift slices $z\sim 0.9;1.5;2.3$ of $\delta\sigma_{0}=8;10;13$ km s*-1*. Asymmetric uncertainties can be as high as $\delta\sigma_{0}=37$ km s*-1*.

In Figure 1, we compare physical properties of our final sample (blue shading) to the underlying representative population of star-forming galaxies from the 3D-HST survey (grey shading) and to the full KMOS3D sample (pink lines). Compared to our full KMOS3D sample, we have not selected preferentially in redshift. In terms of stellar mass, both our full KMOS3D sample and our kinematic sample include fewer lower mass systems compared to the 3D-HST galaxies, such that our sample is not mass-complete. This is mainly a consequence of the $K_{\rm AB}\lesssim 23$ cut. With respect to the main sequence of star-forming galaxies, however, our kinematic sample follows the distribution of both the full KMOS3D and the 3D-HST sample. The fraction of systems with effective radii below the population average is smaller for our kinematic sample compared to the 3D-HST and KMOS3D samples. This is due to our conservative definition of resolved kinematics, where we request measurements over at least $3\times{\rm PSF}_{\rm FWHM}$, with the primary effect of reducing the number of galaxies with $R_{e}<2$ kpc. Generally, for very small systems it is more challenging to recover the intrinsic velocity dispersion, because the kinematics are often unresolved (but see Wisnioski et al., 2018, for a detailed study of the kinematics of compact galaxies in the KMOS3D survey). Axis ratios of our galaxies are homogeneously distributed, following the KMOS3D and 3D-HST parent samples (see also Section 3.6).

In Figure 2, we show SFR (top) and size (bottom) both as a function of stellar mass for the 3D-HST parent sample (grey density histogram), the full KMOS3D sample (purple diamonds), and our final kinematic sample at redshifts $z\sim 0.9$ (blue circles), $z\sim 1.5$ (green pentagons), and $z\sim 2.3$ (red hexagons). The figure illustrates the homogeneous coverage of the KMOS3D survey of typical main sequence galaxies over more than two orders of magnitude in stellar mass. Similarly, the galaxies from our final sample are distributed along the main sequence and have typical sizes for their redshifts, with a tendency towards higher-than-average sizes particularly at $z\sim 2.3$. This bias at the highest redshifts is introduced through our conservative definition of resolved galaxies, and ensures robust $\sigma_{0}$ measurements even at these high redshifts.

3.4 Upper Limit Cases

Our final sample contains 28 galaxies for which the best-fit $\sigma_{0}$ value within the 1$\sigma$ uncertainties is lower than 10 km s*-1*. In using the H$\alpha$ line we are supposedly tracing emission originating from ionized Hii regions. Due to thermal broadening ($\sigma_{\rm th}\approx 10$ km s*-1*) as well as the expansion of Hii regions ($v_{\rm ex}\gtrsim 10$ km s*-1*), we expect some minimum velocity dispersion for the average galaxy of $\sigma_{0}\approx 10-15$ km s*-1* (Shields, 1990).

This minimum value is lower than the typical spectral resolution of KMOS: the effective FWHM spectral resolution at the H$\alpha$ line measured from the reduced data of galaxies in our KMOS3D survey is $\Delta R=\lambda/\Delta\lambda\sim 3515;3975;3860$ in the $YJ$, $H$, and $K$ bands, respectively (Wisnioski et al., in prep.). For our kinematic sample, the corresponding mean spectral resolutions are $\sigma_{\rm instrumental}\sim 37;32;34$ km s*-1*. However, as discussed in more detail in Wisnioski et al., in prep., within the bands there are variations of the spectral resolution of up to $\Delta R=1000$ for individual IFUs. It is therefore crucial to measure the associated spectral resolution at H$\alpha$ for each individual galaxy from sky or arc lines in order to reliably recover the velocity dispersion, as it is done for KMOS3D.

Our line fitting procedure can recover intrinsic velocity dispersions that are a fraction of the instrumental resolution. However, these measurements get increasingly uncertain for decreasing intrinsic velocity dispersions. For galaxies for which the best-fit $\sigma_{0}$ value within the 1$\sigma$ uncertainties is lower than 10 km s*-1*, we adopt as a conservative upper limit the upper 2$\sigma$ boundary of the marginalized posterior distribution derived from the MCMC chain. The resulting upper limits lie between 18 and 53 km s*-1*.

3.5 Validation of Point Spread Function and Line Spread Function Corrections

Before we investigate in detail the redshift evolution of $\sigma_{0}$ and its potential drivers, we want to exclude any residual effects of beam-smearing. Therefore we consider $\sigma_{0}$ as a function of the effective radius, $R_{e}$, and of the ratio of the outermost measured data point to the effective radius, $R_{\rm max}/R_{e}$.

We do not find significant correlations with $R_{e}$ or $R_{\rm{max}}/R_{e}$, as listed in Table 2 (for $R_{e}$ see also Figure 16). We would expect correlations with these parameters if unresolved rotation enters our measure of velocity dispersion. As mentioned in Sections 3.1 and 3.3, we only consider galaxies for our final sample for which we can extract kinematics over a distance of at least $3\times{\rm PSF}_{\rm FWHM}$, with a mean value of $4\times{\rm PSF}_{\rm FWHM}$. However, the extracted kinematics can still be affected by beam-smearing even in the outer parts of the galaxies. The fact that we do not find correlations with size implies that our forward-modelling procedure properly accounts for beam-smearing even for the smaller systems we include.

Similarly, we test for correlations of $\sigma_{0}$ with instrumental resolution and again we do not find a significant correlation, indicating that both our kinematic fitting code and forward-modelling procedure properly account for the instrumental line-spread function (see Table 2).

3.6 Vertical *vs. *Radial Velocity Dispersion

For local galaxies there exists a correlation between galaxy inclination and line-of-sight velocity dispersion. This is due to the transition from measuring predominantly vertical velocity dispersion in face-on systems to measuring predominantly radial velocity dispersion in edge-on systems, with a typical ratio of $\sigma_{z}/\sigma_{r}\sim 0.6$ (van der Kruit & Freeman, 2011; Glazebrook, 2013). For instance, Leroy et al. (2008) find for the THINGS sample that the Hi line-of-sight velocity dispersion increases for galaxies with $i>60^{\circ}$ ($b/a<0.5$), as does the variation of velocity dispersion in individual galaxies. Intriguingly, evidence for higher velocity dispersions in more edge-on systems has been found in the higher resolution $z\sim 1-2$ data from the SINS survey (Genzel et al., 2011). For our KMOS3D kinematic sample, and in agreement with the earlier results by Wisnioski et al. (2015), we do not find a correlation between $\sigma_{0}$ and $b/a$, as shown in Figure 3 and Table 2, possibly due to the coarser spatial resolution of our data.

4 Velocity Dispersion Increases with Redshift

Previous studies have shown that the velocity dispersion of star-forming galaxies increases with redshift (Förster Schreiber et al., 2006; Genzel et al., 2006, 2011; Weiner et al., 2006; Kassin et al., 2007, 2012; Wisnioski et al., 2012, 2015; Newman et al., 2013; Simons et al., 2016, 2017; Mason et al., 2017; Turner et al., 2017; Johnson et al., 2018), albeit with large uncertainties and scatter. In the following, we confirm and increase the robustness of this conclusion with the highest quality IFU data now available with KMOS3D on sub-galactic scales, over a wider redshift and mass range than previously, and using a sample purely selected on the basis of disk galaxies near the main sequence at each redshift. We further put our results into the broader literature context, including multi-phase gas velocity dispersion and expanding the redshift range to $0<z<4$.

4.1 The KMOS3D Velocity Dispersions from $z=2.6$ to $z=0.6$

In Figure 4 we show the intrinsic velocity dispersion of our KMOS3D galaxies in the kinematic sample as a function of redshift, where upper limits are indicated as arrows (Section 3.4). Our data reflect the known trend of increasing average velocity dispersions with increasing redshift.

To quantify the evolution, we fit a linear relation in $\sigma_{0}-z$ space to our best-fit data.222Our results do not depend on this particular functional form, and we list fits in $\sigma_{0}-\log(1+z)$ and $\log(\sigma_{0})-\log(1+z)$ space in Appendix B. We use the Bayesian approach to linear regression by Kelly (2007) which allows for the inclusion of censored data (i.e., upper limits). The routine requires symmetric uncertainties, which we calculate as the mean of the asymmetric uncertainties on $\sigma_{0}$ from our MCMC modelling.333 We assume an uncertainty on $z$ of five times the spectral resolution in each redshift bin, translating into (negligible) uncertainties of $\delta z\sim 0.001-0.002$. Figure 4 shows the derived fit as a solid line, with average values of $\sigma_{0}\sim 31.1;38.3;46.7$ km s*-1* at $z\sim 0.9;1.5;2.3$. The corresponding fit is described by the equation

[TABLE]

We also perform a ‘robust’ fit where the upper limit cases are not included, but entirely left out. We find a slightly shallower evolution indicated by the dashed line. If, instead, for these galaxies we do not assign upper limits but take the formal median of the posterior distribution at face value, we find a slightly steeper evolution indicated by the dash-dotted line. In Table 3 we list our fit parameters and uncertainties.

The $\sigma_{0}-$evolution we derive between $z\sim 0.9$ and $z\sim 2.3$ is slightly shallower than what has been reported for the first year of data from the KMOS3D survey by Wisnioski et al. (2015). In particular, the authors cite $\sigma_{0}\sim 24.9$ km s*-1* at $z\sim 0.9$ and $\sigma_{0}\sim 47.5$ km s*-1* at $z\sim 2.3$, i.e. a difference of $6-7$ km s*-1* for the lowest redshift bin. We partly attribute this difference to our treatment of upper limits. Indeed, if we take the formal best-fit values of the upper limit cases at face value, we find through linear regression a value of $\sigma_{0}\sim 28.4$ km s*-1* at $z\sim 0.9$ (see Table 3), reducing the difference to $\sim 4$ km s*-1*. This difference is smaller than the uncertainty on the average $\sigma_{0}$ value we derive through our fitting based on the standard deviation of the posterior distribution of the zero-point and slope, which is $\delta\sigma_{0}=4.8$ km s*-1* for the $z\sim 0.9$ bin.

4.2 Quantification of Observational Uncertainties and the Scatter in $\sigma_{0}$

Figure 4 shows substantial scatter in $\sigma_{0}$ at fixed redshift with values from $\sigma_{0}\approx 20$ km s*-1* to $\sigma_{0}\approx 100$ km s*-1*. The question is whether this scatter is physical, or purely driven by observational uncertainties.

As listed in Table 3, our robust best-fit relation has an intrinsic scatter around the regression with a standard deviation of 10.4 km s*-1*, suggesting that part of the scatter is indeed due to real variations of the intrinsic dispersion values, and not just due to measurement uncertainties. To quantify the intrinsic variance in each redshift slice, we first calculate the observed variance around the robust best-fit relation, i.e. the variance of the redshift-normalized dispersion values excluding upper limits. The redshift-normalized values are defined as

[TABLE]

with coefficients $a$ and $b$ as listed in Table 3. Then, we perform a Monte Carlo analysis of the scatter due to uncertainties: for each measurement $i$, we draw 1000 times from a normal distribution $\mathcal{N}(0,\delta\sigma_{0,i})$, where $\delta\sigma_{0,i}$ is the symmetric uncertainty of $\sigma_{0,i}$ derived from our dysmal MCMC modelling, and calculate the corresponding sample variance per redshift slice.

We calculate this intrinsic variance as

[TABLE]

and list the corresponding values in Table 4. From this analysis we conclude that at least $\sim$40–50% of the observed variance, i.e. $\sim$60–70% of the observed scatter, is due to real variations of the intrinsic dispersion values, mostly independent of redshift. We also show a histogram of the redshift-normalized dispersion values in Figure 5, $\sigma_{\rm 0,norm}$, in black, together with a histogram of the Monte Carlo draws from the uncertainty distribution in red. Again, this clearly shows that, even though the uncertainties are substantial, there is residual scatter in our $\sigma_{0}$ distribution beyond what can be accounted for by uncertainties. Further, if we focus on the absolute values listed in Table 4, the intrinsic variance increases above $z\sim 1.5$ such that at $z\sim 2.3$ it has doubled compare to $z\sim 0.9$ and $z\sim 1.5$. This suggests that the population of galaxies in our highest redshift bin is more diverse in ISM conditions compared to the lower-redshift samples.

However, no significant residual trends with $\sigma_{0,{\rm norm}}$ and physical properties related to SFR, mass, size, or rotation velocity remain, as we show in detail in Figure 15 in Appendix C. That means that we cannot identify a physical source for the scatter in $\sigma_{0}$ at fixed redshift. This might be due to the limited dynamical range of our data, or it could imply that the intrinsic scatter is driven through the interplay of more than one parameter. Alternatively, the scatter could be due to real variations of the velocity dispersion on short timescales, for instance caused by a dynamic driver such as minor mergers or variations in gas accretion from the cosmic web. This has recently been proposed by Hung et al. (2019) based on results from the FIRE simulations, where variations of intrinsic dispersion are connected to variations of the gas inflow rate on time scales $\lesssim 100$ Myr.

4.3 Comment on the Effect of Sample Selection

The results presented above and in the remainder of the paper are based on our kinematic sample as defined in Section 3.3, i.e. 175 resolved and rotation-dominated disk galaxies that are well-fit by our dynamical model, without strong contamination from OH lines or outflows, and without close neighbours. If we instead consider all modelled galaxies from setup 1, about twice as many compared to the kinematic sample (see Section 3.2), we find a similar median evolution of $\sigma_{0}\approx 31;40;49$ km s*-1* at $z\sim 0.9;1.5;2.3$, however the mean values in the three redshift slices are systematically higher with $\sigma_{0}\approx 34;45;58$ km s*-1*. While at all redshifts the scatter is substantially increased due to galaxies with higher observational uncertainties or poor fits (VAR${}_{\rm obs}\approx 730;850;1560$ km2 s*-2*), the systematic increase of the mean values is mostly due to the inclusion of dispersion-dominated systems (see e.g. Newman et al., 2013, for a discussion of such galaxies).

4.4 Multi-phase Velocity Dispersions from $z=4$ to $z=0$

To put the evolution of velocity dispersion from $z=2.6$ to $z=0.6$ based on our KMOS3D sample into a broader context, we collect measurements reported in the literature from $z\sim 4$ to $z=0$, covering 12 Gyr of cosmic history (Table 5).

In the top panel of Figure 6 we show again our KMOS3D kinematic sample as clouds of grey circles, including upper limits as arrows, in the $\sigma_{0}-z$ space. The median values at $z\sim 0.9;1.5;2.3$, shown as large circles in blue, green, and red, are based on the best fit plotted in Figure 4 and its uncertainties (see Table 3). We include other individual intrinsic dispersion measurements or averages from ionized gas as colored symbols, and atomic and molecular data as black symbols, which are listed in Table 5. Error bars show the mean uncertainty of individual systems in those samples. In our comparison, we do not apply any corrections or normalizations in mass (cf. Wisnioski et al., 2015) which are expected to be small for main sequence galaxies (Simons et al., 2017).

In Table 5 we also list the different techniques used to correct for beam-smearing effects. As explained in Section 3.2 and in the references listed there, we account for beam-smearing effects through a full forward-modelling of both the velocity and velocity dispersion fields with a unique PSF model for each galaxy. Techniques based on only the velocity information, or on grid-based models or look-up tables, might perform less well in their beam-smearing corrections generally resulting in overestimated intrinsic velocity dispersions. For slit surveys, systematic offsets towards higher values might be expected due to the sometimes uncertain galaxy position angle and the resulting difficulties in disentangling rotational and turbulent motions (see Price et al., 2016, 2019, for a discussion and solution approach). Similarly, the methods chosen to calculate or model the intrinsic velocity dispersion might further introduce systematic differences. We note that recent work by Varidel et al. (2019) on a sample of 20 local SFGs suggests that complex structure in the gas distribution may further impact the derived dispersion values.

Figure 6 shows generally good agreement of the various $\sigma_{0}$ measurements reported in the literature. Comparing slit *vs. *IFU techniques, the slit measurements shown here, i.e. data from DEEP2, SIGMA, and MOSDEF, have a tendency towards higher values as compared to the averages derived from our KMOS3D and SINS/zC-SINF surveys, likely for the reasons discussed above, but agree within their uncertainty with the IFU data where available. Interestingly, the deep measurements obtained for individual targets by Genzel et al. (2017), and particularly for the lensed systems by Livermore et al. (2015); Jones et al. (2010) at $1.5<z<3$ also tend towards higher $\sigma_{0}$ values, but have moderate values at $z>3$ in agreement with the averages obtained from seeing-limited IFU and slit spectroscopy by Gnerucci et al. (2011); Turner et al. (2017); Price et al. (2019). Generally, the statistical power of these time-intensive and challenging individual measurements is still very limited. Systematic differences in $\sigma_{0}$ may arise through selection effects: for instance, the nearby galaxies from the DYNAMO sample are selected to be $z\sim 2$ analogues and have many physical properties, including dispersions, similar to high$-z$ SFGs (see Green et al., 2014; White et al., 2017; Fisher et al., 2019).

In contrast, the molecular and atomic data indicated by black points suggest somewhat lower values on average, particularly at $z\approx 0$. Levy et al. (2018) study 17 nearby, rotation-dominated SFGs in CO and ionized gas. They find consistently higher rotation velocities (<$v_{\rm CO}-v_{{\rm H}\alpha}$>$\approx 14$ km s*-1*) and lower velocity dispersions (<$\sigma_{\rm CO}-\sigma_{{\rm H}\gamma}$>$\approx-17$ km s*-1*) for the molecular gas as compared to the ionized gas. At high redshift, there exist only few multi-phase measurements of the intrinsic gas velocity dispersion. Detailed observations reveal comparable values for ionized and molecular gas (Genzel et al., 2013; Übler et al., 2018), however the uncertainties are larger such that differences like those found locally could be washed out.

\movetabledown

=5cm {rotatetable*}

4.5 Multi-phase Gas Velocity Dispersions Evolve Similarly with Redshift

We quantify the difference between the atomic+molecular and the ionized gas velocity dispersions over cosmic time in the bottom panel of Figure 6. Fitting a robust linear relation444 We use the least trimmed squares method by Cappellari et al. (2013). to the average local and individual high$-z$ measurements of atomic+molecular gas, we find a zero point of $a=10.9\pm 0.6$ km s*-1* and a slope of $b=11.0\pm 2.0$ km s*-1* (grey dashed line). For the ionized velocity dispersion, we choose in addition to our own averages from the KMOS3D and SINS/zC-SINF surveys the other large KMOS surveys, KROSS and KDS, and the local average from the GHASP survey. This choice maximizes the redshift range and avoids systematic effects at $z>0$ through different instrumentation. We find a higher zero-point offset of $a=23.3\pm 4.9$ and a somewhat shallower slope of $b=9.8\pm 3.5$, while the extrapolation of our best fit to the KMOS3D data only gives $a=21.1\pm 3.0$ and $b=11.3\pm 2.0$ (Table 3). Fixing the slope to that of the atomic+molecular fit, the zero point shifts in between these measurements, with $a=22.8$. In Table 6 we list our fit parameters and uncertainties.

This suggests that the redshift evolution of the intrinsic velocity dispersion in all gas phases is quite comparable, but their normalization differs. Typical thermal broadening of the atomic/molecular and the ionized gas due to their characteristic temperatures are $\sim 5$ km s*-1* and $\sim 10$ km s*-1*, respectively, meaning the measured velocity dispersions are super-thermal even in the local Universe. Part of the difference between atomic+molecular and ionized gas velocity dispersions can be explained through the expansion of Hii regions from which the ionized emission originates, with typical values of $10-25$ km s*-1* (Shields, 1990), accounting for another $\sim 5-15$ km s*-1* when added in quadrature. In combination, these effects can explain the difference in the local normalizations of the gas phase velocity dispersions, as well as their average offset of $\sim 10-15$ km s*-1* at fixed redshift.

Clearly, more studies of high$-z$ molecular kinematics are warranted to corroborate our result, which potentially has important implications for work on ionized gas kinematics.

4.6 Comments on Thin *vs. *Thick Disk Evolution

Figure 6 shows a smooth evolution of velocity dispersion with redshift over the past $\sim 12$ Gyr, likely connected to decreasing accretion rates and gas fractions with cosmic time (see Sections 5.3 and 5.4). This evolution suggests that also the typical thickness of the young, star-forming gas disk is lower for lower redshift SFGs, as has also been found in state-of-the-art cosmological simulations (Pillepich et al., 2019).

This is potentially interesting in the context of Galactic archeology: early research of the vertical structure of our Milky Way found evidence for two main, distinct exponential disks with scale heights of $\sim 300$ pc and $\sim 1450$ pc (Gilmore & Reid, 1983). This was confirmed through later work on the Milky Way as well as nearby edge-on galaxies (e.g. Dalcanton & Bernstein, 2002; Yoachim & Dalcanton, 2006; Jurić et al., 2008). The thick disk components have been found to be generally older ($>6$ Gyr) than the thin disks, raising the question of distinct formation periods. Naturally, observations of the typically thick high$-z$ disks also prompted the question of the connection between these early thick disks and modern disk structure (e.g. Elmegreen & Elmegreen, 2006).

To explicitly address the question of distinct formation periods of thin vs. thick disks, we make the simple assumption of a step function describing $\sigma_{0}$ of the ionized gas below and above $z=1$. Unsurprisingly, the resulting fit with $\sigma_{0}=28$ km s-1 at $z<1$ and $\sigma_{0}=42$ km s*-1* at $z>1$ is not a good description of the compiled data, with a goodness of fit that is a factor $\sim 20$ worse compared to the linear fit shown in the bottom panel of Figure 6.

Our results suggest that in the absence of recent major mergers it should depend primarily on the star-formation history (connected to gas accretion) if present-day galaxies have distinct disks of different age and scale height, or if there is rather one component with a vertical age gradient (see also Leaman et al., 2017). This interpretation is in agreement with the recent work by Bovy et al. (2012, 2016); Rix & Bovy (2013) who argue based on elemental abundances that the Milky Way has a continuous range of different scale heights, with no sign of a thin-thick disk bimodality. Simulations by e.g. Burkert et al. (1992); Aumer et al. (2016); Aumer & Binney (2017); Grand et al. (2016) support this picture.

However, in this context it is important to remember that based on stellar and gas masses of our galaxies and results from co-moving number density studies (e.g. Brammer et al., 2011), only the lower mass, lower redshift systems in our sample may evolve into present-day disk galaxies, while the galaxies that have high baryonic masses already at high redshift will most likely evolve into present-day’s early-type galaxies. With our data, we do therefore not necessarily track the change in star-forming scale height over time for progenitor-descendant populations, but rather the change in average star-forming scale height of main sequence galaxies at different epochs.

5 What Drives the Gas Velocity Dispersion?

5.1 Galaxy-scale Velocity Dispersion Correlates with Gas Mass and SFR Properties

The redshift dependence of $\sigma_{0}$ suggests that one or more physical galaxy properties that are themselves redshift-dependent drive velocity dispersion. Consistent with previous findings in the literature (e.g. Johnson et al., 2018), we find several properties positively correlating with $\sigma_{0}$, particularly, SFR, SFR surface density $\Sigma_{\rm SFR}$, gas and stellar mass, and their surface densities. We list direct and residual (after correcting for redshift dependence) Spearman rank correlations in Table 7 and show plots for several quantities in Figure 16 in Appendix C. In Table 7, we also list SFRHα and $\Sigma_{\rm SFR,H\alpha}$ derived from the H$\alpha$ fluxes (see Wisnioski et al., in prep.), tracing the more recent star formation history, but find no appreciable difference in correlations compared to our fiducial SFR properties (see Section 2).

We emphasize that due to the limited dynamical range covered by the individual redshift slices, we do not find significant correlations of $\sigma_{0}$ within one redshift slice with any of the above properties, such that we cannot readily connect the scatter in $\sigma_{0}$ at fixed redshift to a physical driving source. Similarly, if we remove the redshift dependence of $\sigma_{0}$ by normalizing with our best-fit relation, we do not find any significant correlations of the normalized $\sigma_{0}$ with physical properties (see Section 4.2 and Figure 15).

Over the full redshift range covered by our KMOS3D survey, SFR and gas mass correlate most strongly and significantly with intrinsic velocity dispersion. In order to identify which of these two physical quantities is most directly tied to the elevated velocity dispersions at high redshift, we discuss in the following sections the physical mechanisms through which quantities such as SFR and gas mass may affect velocity dispersion, namely stellar feedback and gravitational instabilities, and we comment on the tentative connection to AGN feedback for individual galaxies.

5.2 Stellar Feedback

Turbulence-driving can be provided through thermal and momentum feedback from massive stars. Correlations between intrinsic velocity dispersion and SFR properties have previously been reported in the literature (e.g. Dib et al., 2006; Lehnert et al., 2009, 2013; Green et al., 2010, 2014; Moiseev et al., 2015; Johnson et al., 2018), and often invoked the argument for stellar feedback-driven turbulence.

From a theoretical point of view, feedback-driven turbulence is mainly generated through momentum injection from supernovae into the ISM (contributions to the momentum injection from e.g. expanding Hii regions or stellar winds are minor, see Mac Low & Klessen, 2004; Ostriker & Shetty, 2011). Feedback-driven turbulence should therefore primarily depend on the decay rate of turbulence, the momentum injected per supernova, and the supernova rate, where the latter is the quantity connecting turbulence to SFR and $\Sigma_{\rm SFR}$. Ostriker & Shetty (2011) and Shetty & Ostriker (2012) derive a weak dependence of $\sigma_{0}$ on star formation rate surface density. Even considering the case where feedback-driven turbulence vertically stabilizes the disk, the resulting velocity dispersions are low (Equation (22) by Ostriker & Shetty, 2011):

[TABLE]

Here, $f_{p}$ is a factor characterizing the evolution of turbulence, with $f_{p}=1$ for strong dissipation, and $f_{p}=2$ for weak dissipation. $\chi$ is a measure of the importance of the gas disk’s self-gravitational weight, and is below 0.5 for marginally stable disks, such that the first factor is in the range $\sim 0.8-2$. The mean star formation efficiency $\epsilon_{\rm ff}(\rho_{0})$ is assumed to be approximately constant with a fiducial value of $\epsilon_{\rm ff}(\rho_{0})=0.005$. $p_{*}/m_{*}=3000$ km s*-1* is the fiducial value of momentum injection per supernova (but see e.g. Fisher et al., 2019, for arguments for a $z$-dependent $p_{*}/m_{*}$). As a result, the gas velocity dispersion is expected to vary only mildly due to supernova explosions.

Similar results are obtained by other theoretical models investigating stellar feedback as the sole driver of the turbulence in the ISM, for instance the models discussed by Dib et al. (2006); Joung et al. (2009); Kim et al. (2013). In fact, the resulting velocity dispersions in the ISM do not even seem to depend much on the supernova rate. Rather, very high supernova rates might create super-bubble structures that, instead of stirring the ambient medium, will eventually blow out of the galactic disk, thus transferring energy and metals into the circum-galactic medium (Mac Low et al., 1989; Joung et al., 2009; Baumgartner & Breitschwerdt, 2013; Kim & Ostriker, 2018). This is an important result because at higher redshift also the supernova rate is higher. However, other work indicates that not only the rate but also the location of supernovae is crucial for the efficiency of stellar feedback turbulence driving: considering peak driving, where supernovae go off in the densest ISM regions (e.g. their birth clouds), Gatto et al. (2015) find local H$\alpha$ velocity dispersions of up to 60 km s*-1* for gas mass surface densities of $\Sigma_{\rm gas}\sim 100~{}M_{\odot}~{}$pc*-2*. This is similar to high$-z$ conditions and therefore suggests that stellar feedback can more easily maintain elevated velocity dispersions at higher redshift. Also, some idealized simulations of isolated galaxies are able to produce velocity dispersion of $\sim 50$ km s*-1* from strong stellar feedback (Hopkins et al., 2011).

If stellar feedback is an important factor in powering turbulence, then not only would the (observed) global scaling of velocity dispersion with SFR or $\Sigma_{\rm SFR}$ be expected, but particularly locally elevated velocity dispersion in regions of high star formation rate density (cf. Gatto et al., 2015). We exploit the high-resolution data from the SINS-zC/SINF AO survey (Förster Schreiber et al., 2018) to study local correlations between $\Sigma_{\rm SFR}$ and $\sigma_{0}$. In Figure 7 we show the local intrinsic velocity dispersion per spaxel of galaxy Q2346-BX482 as a function of radius, color-coded by $\Sigma_{\rm SFR}$ (adopted from Figure A1 by Genzel et al., 2011). The local intrinsic velocity dispersion is derived from the observed dispersion map, after correcting all instrumental and beam-smearing effects through modelling. In the vicinity of the giant star-forming clump $\sim 6.5$ kpc South-East from the center (inset), no elevated velocity dispersion can be registered.

In Figure 8 we show the local intrinsic velocity dispersion per spaxel as a function of local $\Sigma_{\rm SFR}$ for ten SINS/zC-SINF galaxies. The velocity dispersions of these galaxies with a mean redshift of $z\sim 2.2$ have somewhat higher values compared to our KMOS3D sample, consistent with their higher average SFR and $\Sigma_{\rm SFR}$. Only two of these galaxies show an intrinsic scaling of $\sigma_{0}$ with $\Sigma_{\rm SFR}$. The best-fit power-law relation derived from this sub-galactic, high-quality data shows a very weak dependence of local $\sigma_{0}$ on $\Sigma_{\rm SFR}$,555This finding does not extend to nuclear regions, since more complex circum-nuclear kinematic structure caused by a combination of nuclear outflows, radial inflow and bulge induced rotation in a number of cases generates unresolved velocity fields that appear as an increased velocity dispersion. To explore its true nature will require $<0.1\arcsec$ IFU spectroscopy on 30 m class telescopes. confirming the earlier findings by (Genzel et al., 2011, see also Patrício et al., 2018; Tadaki et al., 2018; but Swinbank et al., 2012a). Similar results are found for both ionized gas (Varidel et al., 2016; Zhou et al., 2017) and molecular gas (Caldú-Primo & Schruba, 2016) in local galaxies. For atomic gas, several studies of local galaxies find correlations with SFR or $\Sigma_{\rm SFR}$ that are too weak to explain the turbulent velocities in the galaxy outskirts (e.g. Tamburro et al., 2009; Ianjamasimanana et al., 2015; Utomo et al., 2019).

In summary, while global $\sigma_{0}$ correlates with SFR properties, we do not find a direct connection between high, local star-formation activity and elevated $\sigma_{0}$, as suggested by some simulations. Generally, however, simulations and models agree that stellar feedback is able to maintain galaxy-wide turbulence on scales of 10-20 km s*-1*.

5.3 Marginally Toomre-stable Disks

Turning to gravity-driven turbulence, an empirical model to describe the redshift evolution of velocity dispersion is that of marginally stable disks, where (non-interacting) galaxies are subject to gas replenishment from the halo or the cosmic web, and to gas loss through either outflows or star formation (Noguchi, 1999; Silk, 2001; Immeli et al., 2004a, b; Förster Schreiber et al., 2006; Elmegreen et al., 2007; Genzel et al., 2008; Dekel et al., 2009b; Dekel & Burkert, 2014; Bouché et al., 2010; Krumholz & Burkert, 2010; Krumholz & Burkhart, 2016; Cacciato et al., 2012; Davé et al., 2012; Lilly et al., 2013; Saintonge et al., 2013; Rathaus & Sternberg, 2016; Leaman et al., 2017). In this framework, the (in)stability of the disk directly corresponds to the level of turbulence in the interstellar medium, where turbulence is fed through external (accretion) and internal (radial flows, clump formation) events via the release of gravitational energy, creating a self-regulation cycle to maintain marginal stability (Dekel et al., 2009b; Genel et al., 2012a; but see Elmegreen & Burkert, 2010).

For a snapshot in time that represents the observation of a high$-z$ galaxy, this equilibrium situation is captured through the Toomre $Q$ parameter (Toomre, 1964), where generally $Q<Q_{\rm crit}\approx 1$ indicates gravitational instability. Considering the one-component approximation for a gas disk, we can write (Binney & Tremaine, 2008; Escala & Larson, 2008; Dekel et al., 2009a)

[TABLE]

Here, $\kappa$ is the epicyclic frequency, $a$ is a constant taking values of $1$ and $\sqrt{2}$ for Keplerian and constant rotation velocity, respectively, and $v_{c}$ is the circular velocity tracing the dynamical mass.

The framework of Toomre-(in)stability generally refers to the linear regime, where perturbations are assumed to be axisymmetric. The galaxies studied here, however, are in the non-linear limit where the ISM is turbulent and many stars have formed (Mandelker et al., 2014). Inoue et al. (2016) investigated the stability of simulated high$-z$ disks, finding that large parts of the disks are in the non-linear regime with $Q>1-3$. This result however depends on gas fraction, which is generally too low in the simulations. Indeed, for those simulated galaxies with the highest gas fractions ($f_{\rm gas}\sim 0.4$, still lower than for typical $z\sim 2$ SFGs), Inoue et al. (2016) find values more compatible with observational findings. Meng et al. (2018) argue in recent work that the Toomre$-Q$ linear stability analysis is still applicable to simulated high$-z$ galaxies, with values of $Q\sim 0.5-1$ in gas-rich regions (see also Behrendt et al., 2015, for simulations of isolated gas-rich disks).

Generally, for a multi-component system an effective $Q$ parameter has to be computed, $Q_{\rm eff}^{-1}=\sum\nolimits_{i}Q_{i}^{-1}$, where $i$ refers to e.g. stars or different gas phases (e.g. Wang & Silk, 1994; Escala & Larson, 2008; Genzel et al., 2011; Obreschkow et al., 2015, and references therein). Simulations of galaxy formation support a picture where $Q_{\rm eff}\sim 1$ for high$-z$ galaxies, and $Q_{\rm eff}\sim 2-3$ for low$-z$ galaxies where the increasing impact of a stellar disk increases $Q_{\rm crit}$ (Hohl, 1971; Athanassoula & Sellwood, 1986; Bottema, 2003; Immeli et al., 2004a; Kim & Ostriker, 2007; Agertz et al., 2009a, b; Aumer et al., 2010; Ceverino et al., 2010; Hopkins et al., 2011; Genel et al., 2012b; Danovich et al., 2015). For gas-rich, thick disks instead $Q_{\rm crit}$ decreases, such that for $z\gtrsim 1$ galaxies values $Q_{\rm crit}\approx 0.7$ are expected (e.g. Goldreich & Lynden-Bell, 1965; Kim & Ostriker, 2007; Wang et al., 2010; Behrendt et al., 2015).

It has been shown that the gas-rich, star-forming disks observed at high redshift are at most marginally stable to gravitational fragmentation (Genzel et al., 2011; see also Swinbank et al., 2017; Johnson et al., 2018; Tadaki et al., 2018; and Fisher et al., 2017 for local high$-z$ analogues), and Wisnioski et al. (2015) have shown that the redshift evolution predicted by Equation (5) for $Q\sim 1$ gas disks is in remarkable agreement with observations (see also e.g. Green et al., 2014; Turner et al., 2017; White et al., 2017; Johnson et al., 2018). In addition, Genzel et al. (2011) have shown that on sub-kpc spatially resolved scales, values of $Q\sim 0.2$ can be reached in regions of star-forming clumps, possibly demonstrating gravitational fragmentation at work.

We calculate $Q_{\rm gas}$ for our galaxies following Eq. (5) by evaluating the circular velocity at $v_{c}(r=1.38R_{e})$. As mentioned in Section 3.3, this radius corresponds to the theoretical peak of a Noordermeer disk with $n_{\rm S}=1$, such that the local gradient of the rotation curve is flat, leading to $a=\sqrt{2}$. The circular velocity $v_{c}$ is computed from the model rotation velocity corrected for pressure support from the turbulent motions (Burkert et al., 2010, 2016; Wuyts et al., 2016). In Figure 9 we show $f_{\rm gas}=M_{\rm gas}/M_{\rm bar}$, with $M_{\rm bar}=M_{*}+M_{\rm gas}$, as a function of $Q_{\rm gas}$, color-coded by redshift as in Figure 4. Albeit large scatter, an anti-correlation between $f_{\rm gas}$ and $Q_{\rm gas}$ is evident, such that galaxies with higher gas fractions have lower $Q$ (Spearman rank correlation coefficient $\rho_{S}=-0.30$ with significance $\sigma_{\rho}=3.6$). This is in agreement with the theoretical prediction that SFGs that are more gas rich have lower $Q$ values. The average $Q_{\rm gas}$ for our galaxies in the redshift bins $z\sim 0.9;1.5;2.3$ is $Q_{\rm gas}=1.2;0.7;0.5$ (arrows in Figure 9). Our results on the average offset of ionized *vs. *atomic+molecular gas from Section 4.5 suggest that the cold gas tracing the bulk of the gas mass might have lower velocity dispersion by $10-15$ km -1. This would lower the $Q_{\rm gas}$ values by a factor $\sim 1.2-2$. While our calculation of the Toomre$-Q$ parameter is simplified through the omission of the stellar component, this suggests that thick high$-z$ disks with high gas fractions of $\gtrsim 50$ % can be marginally stable even down to $Q_{\rm gas}<0.7$.

5.4 Combining Feedback and Gravity

While gravitational instabilities are likely important drivers of the elevated velocity dispersions at $z>1$, the contribution from stellar feedback-driven turbulence of the order of $10-20$ km s*-1* could become comparable or even dominant for lower$-z$, low$-\sigma_{0}$ galaxies. Therefore, one must consider both processes to get a complete picture.

The combination of stellar feedback and gravitational processes for turbulence driving has recently been investigated through the analytic model for structure and evolution of gas in galactic disks by Krumholz et al. (2018), who combine prescriptions for star formation, stellar feedback, and gravitational instabilities into a unified ‘transport+feedback’ model to explain the range of observed dispersions from $z=3$ to the present day. In their model, gas is in vertical hydrostatic equilibrium and energy equilibrium. This model assumes (isolated) rotating galactic disks built of gas and stars within a quasi-spherical dark matter halo over a wide redshift range. Disks are stable or marginally stable to gravitational collapse, regulated by mass transport through the disk. The gas is in vertical hydrostatic equilibrium, and in energy equilibrium such that losses through the decay of turbulence are balanced by energy input into the system via stellar feedback and the release of gravitational energy via mass transport through the disk.

Consistent with the discussion above, Krumholz et al. (2018) show in their model that stellar feedback may maintain velocity dispersions of $\sim 10$ km s*-1*, creating a dispersion floor, while gravitational instabilities, for instance created through radial mass transport through the disk, are necessary to constantly drive velocity dispersions beyond $\sigma_{0}\sim 20$ km s*-1* for moderate star-formation rates (cf. also Figure 4 by Krumholz et al., 2018). They make a prediction for galactic gas velocity dispersion and its correlation with SFR. Particularly, they show that (see their Equation (60))

[TABLE]

where we have substituted appropriate constants for high$-z$ galaxies following Krumholz et al. (2018). Specifically, we adopt a rotation curve slope of $\beta=0$, an offset between resolved and unresolved star formation law normalizations of $\phi_{a}=3$, a fraction of ISM in the star-forming phase $f_{\rm sf}=1$, a ratio of total pressure to turbulent pressure at the midplane of $\phi_{\rm mp}=1.4$, a star-formation efficiency per free-fall time of $\epsilon_{\rm ff}=0.0015$, an orbital period of $t_{\rm orb,out}=200$ Myr, and a maximum star-formation time-scale of $t_{\rm sf,max}=2$ Gyr.

We make two adjustments to our data to properly compare to the model: here, and for all of Section 5.4, we subtracted 15 km s*-1* in quadrature from our intrinsic dispersion values, denoted by $\sigma_{0,15}$, to ensure consistency with the theoretical model (see Krumholz & Burkhart, 2016 and Krumholz et al., 2018, Appendix B). These 15 km s*-1* represent the average combination of thermal motions and expansion of Hii regions that enter our ionized gas velocity dispersion measurement (see also Sections 3.4 and 4.5). We also modify our gas mass fractions: the corresponding parameter used by Krumholz et al. (2018) describes an effective gas fraction at the mid plane. This has typically higher values than our gas fraction $f_{\rm gas}$ because of the larger stellar scale heights compared to the gas scale heights. For the comparison here we adopt a scaling factor of 1.5 for our gas mass fractions, motivated by measurements in the Solar neighbourhood (McKee et al., 2015; Krumholz et al., 2018; M. Krumholz, private communication).

To compare the model prediction from Equation (6) to our data, we group correlated quantities and separate the star formation properties SFR and $f_{\rm gas}$ from the kinematic tracers $v_{\rm circ}$ and $\sigma_{0}$. We show the result for our kinematic sample in Figure 10, specifically SFR divided by gas fraction as a function of circular velocity squared times intrinsic velocity dispersion. Figure 10 reveals a clear trend between the displayed quantities, with a Spearman rank correlation of $\rho_{\rm S}=0.57$ with significance $\sigma_{\rho}=6.8$. We also show model predictions from Krumholz et al. (2018) as quoted in Equation (6) for three values of $Q$. There is a tendency for higher$-z$ galaxies to have a predicted $Q\lesssim 1$, consistent with our results presented in Figure 9. Generally, however, our galaxies scatter around $Q=1$ at all redshifts. This suggests that SFGs self-regulate at all times such that the population of SFGs evolves roughly along lines of constant $Q$. This result is largely independent from the specific choices of parameters such as $\phi_{a}$ or $f$, which will only affect the average $Q$ value.

It is important to realize that the above correlation is predicted for both the combined ‘transport+feedback’ model and a model without feedback, but not for models lacking the ‘transport’ component accounting for gravitational instabilities. Therefore the correlation between SFR/$f_{\rm gas}$ and $v_{\rm circ}^{2}\sigma_{0}$ illustrates the importance of gravitational instabilities as drivers of turbulence for our kinematic sample.

In the following, we now investigate separately changes of circular velocity and gas fraction in the $\sigma_{0}-$SFR parameter space. In Figure 11 we show for our kinematic sample the intrinsic velocity dispersion as a function of SFR, color-coded by circular velocity. As expected from the main sequence and Tully Fisher relation (Tully & Fisher, 1977), which is in place for our data set at all redshifts (Übler et al., 2017), our data display a gradient such that circular velocity on average increases with increasing SFR. As lines we plot the high$-z$ model by Krumholz et al. (2018), but we modify it such that we vary the galaxy circular velocity from $v_{\rm circ}=50$ km s*-1* to $v_{\rm circ}=450$ km s*-1* in order to appropriately cover the range of observed velocities in our kinematic sample. In the model framework, stellar feedback creates sustains a dispersion floor, represented through the horizontal regime of the model lines. The predicted rapid increase of velocity dispersion with SFRs, the exact location here dependent on circular velocity, requires the release of gravitational energy through radial transport through the disk (see Krumholz et al., 2018, for details). The agreement between the theoretical model and our data is remarkably good: $\sim 60$% of our data are matched by the model for this simple variation of only the circular velocity, with all other parameters being fixed to the fiducial ‘transport+feedback high-z’ parameters.

For the gas fraction we can make only an approximate comparison. As mentioned in Section 2, gas masses for our galaxies are calculated applying the scaling relation by Tacconi et al. (2018), since direct gas mass measurements are not available for most of our galaxies. With this, we get the total gas mass over the total baryonic mass per galaxy. Again, we use a scaling factor of 1.5 for our gas mass fractions. In Figure 12 we show the same parameter space as in Figure 11 but now color coded by gas fraction. While galaxies with SFR $\lesssim 10~{}M_{\odot}$ yr*-1* have on average lower gas fractions, no strong trend is apparent at higher SFRs. We show again lines based on the ‘transport+feedback high$-z$’ model by Krumholz et al. (2018), but now we vary the gas fraction (and with it $f_{g,P}$) from $f_{g,Q}=0.2$ to $f_{g,Q}=1.0$ in order to explore the range of scaled gas fractions of galaxies in our kinematic sample. With solid lines we show models with $v_{\phi}=400$ km s*-1*, and dashed lines show $v_{\phi}=200$ km s*-1*. It becomes clear that in the model framework galaxies at fixed SFR and $\sigma_{0}$ can have higher $f_{g,Q}$ and lower $v_{\phi}$, or lower $f_{g,Q}$ and higher $v_{\phi}$, but rotation velocity has to be varied to cover the full range of SFRs in our observations.

Horizontal variations of the model predictions can be reached through changing the fraction of gas assumed to be in the star-forming ISM, and through changes in the outer rotation curve slope. For instance, assuming only 20$\%$ of the gas to be in the star-forming phase pushes the horizontal floor of the model below 10 km s*-1*, and lowers the predicted SFR by almost an order of magnitude. Assuming a dropping rotation curve, on the other hand, lifts the horizontal floor and increases the predicted SFR. Assuming an outer rotation curve slope of $\beta=-0.5$ increases the horizontal saturation of the model to $\sim 32$ km s*-1*, whilst increasing the star-formation rate only marginally. Lang et al. (2017) have shown that the typical outer rotation curve slope of galaxies in our sample is negative. This is more pronounced at higher redshift, possibly offering an additional reason for the elevated velocity dispersions at $z\gtrsim 2$ in this model framework.

Considering these analytic model prescriptions, and the typical uncertainty of the intrinsic dispersion measurements of $\delta\sigma_{0}\sim 10$ km s*-1* in our kinematic sample, we conclude that galaxies with $\sigma_{0}\gtrsim 35$ km s*-1* are dominated by gravitational instability-driven turbulence. This encompasses more than 60$\%$ of galaxies in our sample, underlining the importance of gravity-driven turbulence in star forming galaxies at $z\sim 1-3$.

5.5 AGN feedback

As a final remark, we briefly want to comment on AGN feedback as a potential additional source for elevated velocity dispersions in the SFGs in our kinematic sample. While we excluded galaxies, or regions of galaxies, that are so strongly affected by the AGN and associated outflows that the disk kinematics cannot be recovered, we do not entirely exclude AGN. This ensures that we can explore the full mass range covered by the KMOS3D survey, including the high-mass end where at log$(M_{*}/M_{\odot})>11$, above the Schechter mass, the fraction of AGN increases rapidly (Förster Schreiber et al., 2014, 2018; Genzel et al., 2014).

While we do not find significant correlations between $z-$normalized $\sigma_{0}$ and mass properties (Table 7), we do note a cloud of galaxies from all redshifts with dispersions above average for the highest stellar (log$(M_{*}/M_{\odot})>11$) and baryonic masses (log($M_{\rm bar}/M_{\odot})\gtrsim 11.3$) as shown in Figure 13. About half of the log($M_{\rm bar}/M_{\odot})\gtrsim 11.3$ above-average dispersion galaxies are known to host an AGN (stars in Figure 13). We speculate that the energy deposited by strong AGN feedback in the form of nuclear outflows could induce turbulence in the disk via the re-accretion of material at larger radii.

It is important to keep in mind that outflow components with velocities similar to the galaxy rotation velocity can broaden the line width but may not be distinguishable from the star-forming regions due to S/N limitations. Comparing to the deep AO data from the SINS/zC-SINF survey that we show in Figure 8, one of the three identified, log($M_{*}/M_{\odot})\gtrsim 11$ AGN (Q2343-BX610) shows above-average velocity dispersions (after excluding the regions clearly affected by the nuclear outflow), while the other two (D3a-6004, D3a-15504) have average dispersions.

6 Conclusions

We have investigated the evolution of the ionized gas intrinsic velocity dispersions, $\sigma_{0}$, between $0.6<z<2.6$ based on data from our full KMOS3D survey. We have selected a high-quality, rotation-dominated ($v_{\rm rot}/\sigma_{0}\geq 1$) sample for which we forward-modelled in a Bayesian framework the one-dimensional galaxy kinematics extracted from the H$\alpha$ velocity and velocity dispersion maps, taking into account instrumental effects, beam smearing, and pressure support. Our main conclusions are as follows:

• •

Assuming an isotropic and radially constant intrinsic velocity dispersion, we find an average decrease of the H$\alpha$ intrinsic dispersion for our kinematic sample from $\sigma_{0}\sim 46$ km s*-1* at $z\sim 2.3$ to $\sigma_{0}\sim 31$ km s*-1* at $z\sim 0.9$, solidifying trends previously reported in the literature (Section 4.1). Putting our sample into the broader context of literature measurements from $z=4$ to $z=0$, tracing ionized, molecular, and atomic gas phases, confirms the general increase of intrinsic galaxy velocity dispersion with redshift (Section 4.4).

• •

Comparing the redshift evolution of ionized and molecular plus atomic gas velocity dispersion, we find that the ionized gas dispersion is on average higher by $\sim 12$ km s*-1* (Section 4.5). This offset can in principle be accounted for through the different gas temperatures together with the line broadening through expansion of Hii regions the ionized gas emission typically originates from.

• •

For our KMOS3D kinematic sample, we find that there is intrinsic scatter in the $\sigma_{0}$ distribution at fixed redshift after accounting for measurement and modelling uncertainties, and it increases for our highest redshift slice (Section 4.2). However, we cannot single out a physical mechanism behind this scatter. This could imply that the velocity dispersion is highly variable in time, due to a dynamic mechanism such as minor mergers or variation in accretion (see Hung et al., 2019, for evidence from simulations). Alternatively, the scatter could be caused by the interplay of different physical properties responsible to maintain marginal stability (see Section 5.4).

• •

Investigating the physical driver of the elevated velocity dispersions at higher redshift, we find that galaxies in our kinematic sample are at most marginally Toomre-stable, i.e. they are consistent with their turbulence being powered through gravitational instabilities in a self-regulated environment (Section 5.3).

• •

We find no evidence from our high-resolution SINS/zC-SINF AO data that stellar feedback as traced through $\Sigma_{\rm SFR}$ typically increases the velocity dispersion on sub-galactic scales beyond the average level, or that the local velocity dispersion correlates strongly with $\Sigma_{\rm SFR}$, suggesting that contributions from stellar feedback to turbulence driving are minor for our $z>1$ SFGs (Section 5.2).

• •

We find good agreement between data from our KMOS3D kinematic sample and predictions from the state-of-the-art analytical model of galaxy formation and evolution by Krumholz et al. (2018), further strengthening the evidence that the majority of our galaxies ($\gtrsim 60$%) are dominated by gravity-driven turbulence (Section 5.4).

The measurement of intrinsic gas velocity dispersion at $z>0$ is challenging. Next-generation instruments such as ERIS+AO at the VLT or HARMONI at the ELT will expand current samples on spatial scales that are currently only achievable for strongly lensed objects, and push spectral scales down to $\sim 15$ km s*-1*. The statistics from these observations will facilitate further investigation of the scatter of the intrinsic velocity dispersion at fixed redshift, and tests of theoretical predictions such as the transition regime from gravity-driven turbulence to feedback-driven turbulence as a function of redshift and mass (Krumholz et al., 2018). Deep, high-S/N observations of particularly molecular gas reaching 1–2 kpc resolution at $z>1$ with NOEMA or ALMA are necessary to test if the redshift evolution of molecular and ionized gas velocity dispersion is indeed comparable.

We are grateful to the anonymous referee for a very constructive report that helped to improve the quality of this manuscript. We thank the ESO Paranal staff for their helpful support with the KMOS observations for this work. We are grateful to Raymond Simons and Susan Kassin for providing us with data from the DEEP2 and SIGMA surveys. We thank Amiel Sternberg, Vadim Semenov, and Chia-Yu Hu for fruitful discussions, and we are particularly grateful to Blakesley Burkhart, Mark Krumholz, and Andreas Schruba for valuable comments on various aspects of this work. DW acknowledges the support of the Deutsche Forschungsgemeinschaft via projects WI 3871/1-1 and WI 3871/1-2. EW acknowledges support by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. MF acknowledges the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 757535). PL acknowledges funding from the ERC under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343).

Appendix A Example Galaxies and Fits

We show examples of galaxies in our kinematic sample together with their best-fit kinematic models in Figure 14. See the figure caption for details.

Appendix B Alternative Fits to our KMOS3D Velocity Dispersions

We list fits to our KMOS3D velocity dispersion data from $z=2.6$ to $z=0.6$ in $\sigma_{0}-\log(1+z)$ space and $\log(\sigma_{0})-\log(1+z)$ space in Tables 8 and 9, respectively. These results agree with our fiducial fits in $\sigma_{0}-z$ space listed in Table 3, and do not change our conclusions.

Appendix C Correlations of Physical Properties with Velocity Dispersion and Redshift-normalized Velocity Dispersion

We show correlations of various physical properties with velocity dispersion after (see Equation (2)) and before correcting for the redshift dependence of $\sigma_{0}$ in Figures 15 and 16 (see also Table 7). While several properties positively correlate with $\sigma_{0}$, particularly SFR and $M_{\rm gas}$, we do not find any significant correlation after correcting for the redshift-dependence of $\sigma_{0}$. This means that we cannot readily identify a single physical driving source behind the intrinsic scatter in $\sigma_{0}$ (see discussions in Sections 4.2 and 5).