The evolution of the Tully-Fisher relation between z~2.3 and z~0.9 with KMOS^3D

TL;DR

This study examines how the Tully-Fisher relation for massive star-forming disk galaxies evolves between redshifts 2.3 and 0.9, revealing complex changes in baryonic and stellar mass relations over cosmic time.

Contribution

It provides new insights into the evolution of the Tully-Fisher relation by analyzing spatially resolved data and modeling the effects of pressure support and dark matter contributions.

Findings

Higher baryonic masses at z~2.3 compared to z~0.9 at fixed velocity

No evolution of stellar TFR zero-point from z~0.9 to z~2.3

Negative evolution of TFR zero-points from z=0 to z~0.9

Abstract

We investigate the stellar mass and baryonic mass Tully-Fisher relations (TFRs) of massive star-forming disk galaxies at redshift z~2.3 and z~0.9 as part of the KMOS^3D integral field spectroscopy survey. Our spatially resolved data allow reliable modelling of individual galaxies, including the effect of pressure support on the inferred gravitational potential. At fixed circular velocity, we find higher baryonic masses and similar stellar masses at z~2.3 as compared to z~0.9. Together with the decreasing gas-to-stellar mass ratios with decreasing redshift, this implies that the contribution of dark matter to the dynamical mass at the galaxy scale increases towards lower redshift. A comparison to local relations reveals a negative evolution of the stellar and baryonic TFR zero-points from z=0 to z~0.9, no evolution of the stellar TFR zero-point from z~0.9 to z~2.3, and a positive evolution of the baryonic TFR zero-point from z~0.9 to z~2.3. We discuss a toy model of disk galaxy evolution to explain the observed, non-monotonic TFR evolution, taking into account the empirically motivated redshift dependencies of galactic gas fractions, and of the relative amount of baryons to dark matter on galaxy and halo scales.

Tables

Table 1. Table 1: Median physical properties of our TFR subsamples at z ∼ 0.9 similar-to 𝑧 0.9 z\sim 0.9 ( Y ​ J 𝑌 𝐽 YJ ), z ∼ 1.5 similar-to 𝑧 1.5 z\sim 1.5 ( H 𝐻 H ), and z ∼ 2.3 similar-to 𝑧 2.3 z\sim 2.3 ( K 𝐾 K ), together with the associated central 68 th percentile ranges in brackets.

	$z\sim 0.9$	$z\sim 1.5$	$z\sim 2.3$
	(65 galaxies)	(24 galaxies)	(46 galaxies)
log(M∗ [M⊙])	10.49 [10.03; 10.83]	10.72 [10.08; 11.07]	10.51 [10.18; 11.00]
log(Mbar [M⊙])	10.62 [10.29; 10.98]	10.97 [10.42; 11.31]	10.89 [10.59; 11.33]
SFR [M⊙/yr]	21.1 [7.1; 39.6]	53.4 [15.5; 134.5]	72.9 [38.9; 179.1]
log($\Delta$ MS)aaMS offset with respect to the broken power law relations derived by Whitaker et al. (2014), using the redshift-interpolated parametrization by W15, $\Delta$ MS=$\mathrm{SFR}-{\mathrm{SFR}}_{\mathrm{MS}(z,M_{\ast })\left[\mathrm{W14}\right]}$.	0.20 [-0.21; 0.42]	0.10 [-0.21; 0.45]	-0.01 [-0.29; 0.13]
$R_e^{5000}$ [kpc]	4.8 [3.0; 7.6]	4.9 [3.0; 7.0]	4.0 [2.5; 5.2]
log($\Delta$ M-R)bbOffset from the mass-size relation of SFGs with respect to the relation derived by van der Wel et al. (2014), $\Delta$ M-R=$R_e^{5000}-R_{e,M-R(z,M_{\ast })\left[\mathrm{vdW14}\right]}^{5000}$, after correcting the $H-$band $R_e$ to the rest-frame $5000\AA$.	-0.02 [-0.17; 0.16]	0.08 [-0.10; 0.17]	0.06 [-0.14; 0.17]
$n_S$	1.3 [0.8; 3.1]	0.9 [0.4; 2.2]	1.0 [0.4; 1.6]
$B/T$ccBulge-to-total mass ratio if available, namely for 78, 92, and 89 per cent of our galaxies in $YJ-$, $H-$, and $K-$band, respectively. Values of $B/T=0$ usually occur when the galaxy’s S rsic index $n_S$ is smaller than 1 (cf. Lang et al., 2014).	0.11 [0.00; 0.39]	0.00 [0.00; 0.23]	0.10 [0.00; 0.25]
$v_{\mathrm{rot},\max }$ [km/s]	233 [141; 302]	245 [164; 337]	239 [160; 284]
${\sigma }_0$ [km/s]	30 [9; 52]	47 [29; 59]	49 [32; 68]
$v_{\mathrm{rot},\max }/{\sigma }_0$	6.7 [3.2; 25.3]	5.5 [3.4; 65.6]	4.3 [3.4; 9.1]
$v_{\mathrm{circ},\max }$ [km/s]	239 [167; 314]	263 [181; 348]	260 [175; 315]

Table 2. Table 2: Results from the inverse linear regression fits to Equation ( 2 ) using the least-squares method, including bootstrapped errors of the zero-point. The reference velocity is v ref = 242 subscript 𝑣 ref 242 v_{\rm ref}=242  km/s.

TFR	redshift range	number of galaxies	slope $a$ (local relation)	zero-point $b$ (error)	intrinsic scatter ${\zeta }_{\mathrm{int}}$
			$\left[\frac{\log \left(M\left[M_{\odot }\right]\right)}{\log \left(v_{\mathrm{circ}}\left[\mathrm{km}/s\right]\right)}\right]$	[log($M[M_{\odot }$])]	[dex of $M_{\odot }$]
sTFR		135	3.60 (Reyes et al., 2011)	10.50 ($\pm 0.03$)	0.22
	$z\sim 0.9$	65	3.60 (Reyes et al., 2011)	10.49 ($\pm 0.04$)	0.21
	$z\sim 2.3$	46	3.60 (Reyes et al., 2011)	10.51 ($\pm 0.05$)	0.26
bTFR		135	3.75 (Lelli et al., 2016)	10.75 ($\pm 0.03$)	0.23
	$z\sim 0.9$	65	3.75 (Lelli et al., 2016)	10.68 ($\pm 0.04$)	0.22
	$z\sim 2.3$	46	3.75 (Lelli et al., 2016)	10.85 ($\pm 0.05$)	0.26

Table 3. Table 3: Physical properties of galaxies in our TFR sample in terms of redshift z 𝑧 z , stellar mass M ∗ subscript 𝑀 M_{*} , baryonic mass M bar subscript 𝑀 bar M_{\rm bar} , maximum modelled circular velocity v circ , max subscript 𝑣 circ max v_{\rm circ,max} , and modelled intrinsic velocity dispersion σ 0 subscript 𝜎 0 \sigma_{0} .

#	$z$	log($M_{\ast }\left[M_{\odot }\right]$)	log($M_{\mathrm{bar}}\left[M_{\odot }\right]$)	$v_{\mathrm{circ},\max }$ [km/s]	${\sigma }_0$ [km/s]
1	0.602	10.85	10.93	274.9	30.9
2	0.626	11.00	11.07	314.3	25.8
3	0.669	10.76	10.82	267.5	49.8
4	0.678	10.49	10.58	273.4	38.5
5	0.758	10.66	10.77	313.8	24.3
⋮	⋮	⋮	⋮	⋮	⋮

Equations

$$v_{\mathrm{circ}}(r{)}^2=v_{\mathrm{rot}}(r{)}^2+2{\sigma }_0^2\frac{r}{R_d},$$

$$\log (M[M_{\odot }])=a\cdot \log (v_{\mathrm{circ}}/v_{\mathrm{ref}})+b.$$

$$M_h=\frac{V_h^3}{10G\cdot H(z)};R_h=\frac{V_h}{10H(z)}$$

$$M_{\mathrm{bar}}=\frac{v_{\mathrm{circ}}^3(R_e)}{H(z)}\cdot \frac{[1-f_{\mathrm{DM}}(R_e,z){]}^{3/2}}{m_d^{1/2}(z)}\cdot C$$

$$M_{\ast }=\frac{v_{\mathrm{circ}}^3(R_e)}{H(z)}\cdot \frac{[1-f_{\mathrm{DM}}(R_e,z){]}^{3/2}[1-f_{\mathrm{gas}}(z)]}{m_d^{1/2}(z)}\cdot C^{\prime },$$

$$\Delta{\rm log}(v_{\rm circ})={\rm log}(v_{\rm circ})-\Bigg{[}\frac{-b}{a}+\frac{{\rm log}(M/M_{\odot})}{a}+{\rm log}(v_{\rm ref})\Bigg{]}.$$

$$\Sigma (r)={\Sigma }_0{e}^{-r/R_d},$$

$$M_{\mathrm{bar}}=m_d\cdot M_h;R_{\mathrm{bar}}=r_f\cdot R_h.$$

$$v_{\mathrm{circ}}(r)=\sqrt{v_{\mathrm{bar}}^2(r)+v_{\mathrm{DM}}^2(r)}.$$

$$v_{\mathrm{bar}}^2(r)=4\pi G{\Sigma }_0{R}_d{y}^2[I_0(y)K_0(y)-I_1(y)K_1(y)],$$

$$v_{\mathrm{bar}}^2(R_e)=\frac{M_{\mathrm{bar}}}{R_d}\cdot C^{\prime \prime },$$

$$M_h=R_h^{3}H(z{)}^{2}10^2{G}^{-1}.$$

$$\log \left(\frac{M_{\mathrm{gas},\mathrm{mol}}}{M_{\ast }}\right)\approx$$

$$-0.33\cdot \left[\log (M_{\ast }[M_{\odot }])-10.7\right].$$

Full text

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

††thanks: Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programs 092.A-0091, 093.A-0079, 094.A-0217, 095.A-0047, and 096.A-0025.

1 Introduction

State-of-the-art cosmological simulations in a $\Lambda$CDM framework indicate that three main mechanisms regulate the growth of galaxies, namely the accretion of baryons, the conversion of gas into stars, and feedback. While gas settles down at the centers of growing dark matter (DM) haloes, cools and forms stars, it keeps in its angular momentum an imprint of the dark halo. Conservation of the net specific angular momentum, as suggested by analytical models of disk galaxy formation (e.g. Fall & Efstathiou, 1980; Dalcanton et al., 1997; Mo et al., 1998; Dutton et al., 2007; Somerville et al., 2008), should result in a significant fraction of disk-like systems. In fact, they make up a substantial fraction of the observed galaxy population at high redshift ($1\lesssim z\lesssim 3$; Labbé et al., 2003; Förster Schreiber et al., 2006, 2009; Genzel et al., 2006, 2014b; Law et al., 2009; Epinat et al., 2009, 2012; Jones et al., 2010; Miller et al., 2012; Wisnioski et al., 2015; Stott et al., 2016) and in the local Universe (e.g. Blanton & Moustakas, 2009, and references therein). The detailed physical processes during baryon accretion from the halo scales to the galactic scales are, however, complex, and angular momentum conservation might not be straightforward to achieve (e.g. Danovich et al., 2015). To produce disk-like systems in numerical simulations, feedback from massive stars and/or active galactic nuclei is needed to prevent excessive star formation and to balance the angular momentum distribution of the star-forming gas phase (e.g. Governato et al., 2007; Scannapieco et al., 2009, 2012; Agertz et al., 2011; Brook et al., 2012; Aumer et al., 2013; Hopkins et al., 2014; Marinacci et al., 2014; Übler et al., 2014; Genel et al., 2015). Despite the physical complexity and the diverse formation histories of individual galaxies, local disk galaxies exhibit on average a tight relationship between their rotation velocity $V$ and their luminosity $L$ or mass $M$, namely the Tully-Fisher relation (TFR; Tully & Fisher, 1977). In its mass-based form, the TFR is commonly expressed as $M\propto V^{a}$, or ${\rm log}(M)=a\cdot{\rm log}(V)+b$, where $a$ is the slope, and $b$ is the zero-point offset.

In the local Universe, rotation curves of disk galaxies are apparently generally dominated by DM already at a few times the disc scale length, and continue to be flat or rising out to several tens of kpc (see e.g. reviews by Faber & Gallagher, 1979; Sofue & Rubin, 2001; and Catinella et al., 2006). Therefore, the local TFR enables a unique approach to relate the baryonic galaxy mass, which is an observable once a mass-to-light conversion is assumed, to the potential of the dark halo. Although the luminosity-based TFR is more directly accessible, relations based on mass constitute a physically more fundamental approach since the amount of light measured from the underlying stellar population is a function of passband, systematically affecting the slope of the TFR (e.g. Verheijen, 1997, 2001; Bell & de Jong, 2001; Pizagno et al., 2007; Courteau et al., 2007; McGaugh & Schombert, 2015). The most fundamental relation is given by the baryonic mass TFR (bTFR). It places galaxies over several decades in mass onto a single relation, whereas there appears to be a break in the slope of the stellar mass TFR (sTFR) for low-mass galaxies (McGaugh et al., 2000; McGaugh, 2005).

Observed slopes vary mostly between $3\lesssim a\lesssim 4.5$ for the local sTFR (e.g. Bell & de Jong, 2001; Pizagno et al., 2005; Avila-Reese et al., 2008; Williams et al., 2010; Gurovich et al., 2010; Torres-Flores et al., 2011; Reyes et al., 2011) and between $3\lesssim a\lesssim 4$ for the local bTFR (e.g. McGaugh et al., 2000; McGaugh, 2005; Trachternach et al., 2009; Stark et al., 2009; Zaritsky et al., 2014; McGaugh & Schombert, 2015; Lelli et al., 2016; Bradford et al., 2016; Papastergis et al., 2016). It should be noted that the small scatter of local TFRs can be partly associated to the very efficient selection of undisturbed spiral galaxies (e.g. Kannappan et al., 2002; see also Courteau et al., 2007; Lelli et al., 2016, for discussions of local TFR scatter). Variations in the observational results of low-$z$ studies can be attributed to different sample sizes, selection bias, varying data quality, statistical methods, conversions from $L$ to $M$, or to the adopted measure of $V$ (Courteau et al., 2014; for a detailed discussion regarding the bTFR see Bradford et al., 2016).

Any such discrepancy becomes more substantial when going to higher redshift where measurements are more challenging and the observed scatter of the TFR increases with respect to local relations (e.g. Conselice et al., 2005; Miller et al., 2012). The latter is partly attributed to ongoing kinematic and morphological transitions (Flores et al., 2006; Kassin et al., 2007, 2012; Puech et al., 2008, 2010; Covington et al., 2010; Miller et al., 2013; Simons et al., 2016), possibly indicating non-equilibrium states. Another complication for comparing high-$z$ studies to local TFRs arises from the inherently different nature of the so-called disk galaxies at high redshift: although of disk-like structure and rotationally supported, they are significantly more “turbulent”, geometrically thicker, and clumpier than local disk galaxies (Förster Schreiber et al., 2006, 2009, 2011a, 2011b; Genzel et al., 2006, 2011; Elmegreen & Elmegreen, 2006; Elmegreen et al., 2007; Kassin et al., 2007, 2012; Epinat et al., 2009, 2012; Law et al., 2009, 2012; Jones et al., 2010; Nelson et al., 2012; Newman et al., 2013; Wisnioski et al., 2015; Tacchella et al., 2015b, a).

Despite the advent of novel instrumentation and multiplexing capabilities, there is considerable tension in the literature regarding the empirical evolution of the TFR zero-points with cosmic time. Several authors find no or only weak zero-point evolution of the sTFR up to redshifts of $z\sim 1.7$ (Conselice et al., 2005; Kassin et al., 2007; Miller et al., 2011, 2012; Contini et al., 2016; Di Teodoro et al., 2016; Molina et al., 2017; Pelliccia et al., 2017), while others find a negative zero-point evolution up to redshifts of $z\sim 3$ (Puech et al., 2008, 2010; Cresci et al., 2009; Gnerucci et al., 2011; Swinbank et al., 2012; Price et al., 2016; Tiley et al., 2016; Straatman et al., 2017). Similarly for the less-studied high$-z$ bTFR, Puech et al. (2010) find no indication of zero-point evolution since $z\sim 0.6$, while Price et al. (2016) find a positive evolution between lower-$z$ galaxies and their $z\sim 2$ sample. There are indications that varying strictness in morphological or kinematic selections can explain these conflicting results (Miller et al., 2013; Tiley et al., 2016). The work by Vergani et al. (2012) demonstrates that also the assumed slope of the relation, which is usually adopted from a local TFR in high-$z$ studies, can become relevant for the debate of zero-point evolution (see also Straatman et al., 2017).

A common derivation of the measured quantities as well as similar statistical methods and sample selection are crucial to any study which aims at comparing different results and studying the TFR evolution with cosmic time (e.g. Courteau et al., 2014; Bradford et al., 2016). Ideally, spatially well resolved rotation curves should be used which display a peak or flattening. Such a sample would provide an important reference frame for studying the effects of baryonic mass assembly on the morphology and rotational support of disk-like systems, for investigating the evolution of rotationally supported galaxies as a response to the structural growth of the parent DM halo, and for comparisons with cosmological models of galaxy evolution.

In this paper, we exploit spatially resolved integral field spectroscopic (IFS) observations of 240 rotation-dominated disk galaxies from the KMOS3D survey (Wisnioski et al., 2015, hereafter W15) to study the evolution of the sTFR and bTFR between redshifts $z=2.6$ and $z=0.6$. The wide redshift coverage of the survey, together with its high quality data, allow for a unique investigation of the evolution of the TFR during the peak epoch of cosmic star formation rate density, where coherent data processing and analysis are ensured. In Section 2 we describe our data and sample selection. We present the KMOS3D TFR in Section 3, together with a discussion of other selected high$-z$ TFRs. In Section 4 we discuss the observed TFR evolution, we set it in the context to local observations, and we discuss possible sources of uncertainties. In Section 5 we constrain a theoretical toy model to place our observations in a cosmological context. Section 6 summarizes our work.

Throughout, we adopt a Chabrier (2003) initial mass function (IMF) 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 Data and sample selection

The contradictory findings about the evolution of the mass-based TFR in the literature motivate a careful sample selection at high redshift. In this work we concentrate on the evolution of the TFR for undisturbed disk galaxies. Galaxies are eligible for such a study if the observed kinematics trace the central potential of the parent halo. To ensure a suitable sample we perform several selection steps which are described in the following paragraphs.

2.1 The KMOS3D survey

This work is based on the first three years of observations of KMOS3D, a multi-year near-infrared (near-IR) IFS survey of more than 600 mass-selected star-forming galaxies (SFGs) at $0.6\lesssim z\lesssim 2.6$ with the $K-$band Multi Object Spectrograph (KMOS; Sharples et al., 2013) on the Very Large Telescope. The 24 integral field units of KMOS allow for efficient spatially resolved observations in the near-IR passbands $YJ$, $H$, and $K$, facilitating high-$z$ rest-frame emission line surveys of unprecedented sample size. The KMOS3D survey and data reduction are described in detail by W15, and we here summarize the key features. The KMOS3D galaxies are selected from the 3D-HST survey, a Hubble Space Telescope Treasury Program (Brammer et al., 2012; Skelton et al., 2014; Momcheva et al., 2016). 3D-HST provides $R\sim 100$ near-IR grism spectra, optical to 8 $\mu$m photometric catalogues, and spectroscopic, grism, and/or photometric redshifts for all sources. The redshift information is complemented by high-resolution Wide Field Camera 3 (WFC3) near-IR imaging from the CANDELS survey (Grogin et al., 2011; Koekemoer et al., 2011; van der Wel et al., 2012), as well as by further multi-wavelength coverage of our target fields GOODS-S, COSMOS, and UDS, through Spitzer/MIPS and Herschel/PACS photometry (e.g. Lutz et al., 2011; Magnelli et al., 2013; Whitaker et al., 2014, and references therein). Since we do not apply selection cuts other than a magnitude cut of $Ks\lesssim 23$ and a stellar mass cut of log$(M_{*}~{}[M_{\odot}])\gtrsim 9.2$, together with OH-avoidance around the survey’s main target emission lines H$\alpha$+[Nii], the KMOS3D sample will provide a reference for galaxy kinematics and H$\alpha$ properties of high$-z$ SFGs over a wide range in stellar mass and star formation rate (SFR). The emphasis of the first observing periods has been on the more massive galaxies, as well as on $YJ-$ and $K-$band targets, i.e. galaxies at $z\sim 0.9$ and $z\sim 2.3$, respectively. Deep average integration times of $5.5,7.0,10.5$ h in $YJ,H,K$, respectively, ensure a detection rate of more than 75 per cent, including also quiescent galaxies.

The results presented in the remainder of this paper build on the KMOS3D sample as of January 2016, with 536 observed galaxies. Of these, 316 are detected in, and have spatially resolved, H$\alpha$ emission free from skyline contamination, from which two-dimensional velocity and dispersion maps are produced. Examples of those are shown in the work by W15 and Wuyts et al. (2016, hereafter W16).

2.2 Masses and star-formation rates

The derivation of stellar masses $M_{*}$ uses stellar population synthesis models by Bruzual & Charlot (2003) to model the spectral energy distribution of each galaxy. Extinction, star formation histories (SFHs), and a fixed solar metallicity are incorporated into the models as described by Wuyts et al. (2011).

SFRs are obtained from the ladder of SFR indicators introduced by Wuyts et al. (2011): if Herschel/PACS $60-160\mu$m and/or Spitzer/MIPS $24\mu$m observations were available, the SFRs were computed from the observed UV and IR luminosities. Otherwise, SFRs were derived from stellar population synthesis modelling of the observed broadband optical to IR spectral energy distributions.

Gas masses are obtained from the scaling relations by Tacconi et al. (2017), which use the combined data of molecular gas ($M_{\rm gas,mol}$) and dust-inferred gas masses of SFGs between $0<z<4$ to derive a relation for the depletion time $t_{\rm depl}\equiv M_{\rm gas,mol}/{\rm SFR}$. It is expressed as a function of redshift, main sequence offset, stellar mass, and size. Although the contribution of atomic gas to the baryonic mass within $1-3\,R_{e}$ is assumed to be negligible at $z\sim 1-3$, the inferred gas masses correspond to lower limits (Genzel et al., 2015).

Following Burkert et al. (2016), we adopt uncertainties of $0.15$ dex for stellar masses, and $0.20$ dex for gas masses. This translates into an average uncertainty of $\sim 0.15$ dex for baryonic masses (see § 4.3.1 for a discussion).

2.3 Dynamical modelling

W16 use the two-dimensional velocity and velocity dispersion fields as observed in H$\alpha$ to construct dynamical models for selected galaxies. The modelling procedure is described in detail by W16, where examples of velocity fields, velocity and dispersion profiles, and 1D fits can also be found (see also Figure 1). In brief, radial velocity and dispersion profiles are constructed from $0\farcs 8$ diameter circular apertures every other $0\farcs 2$ along the kinematic major axis using linefit (Davies et al., 2009), where spectral resolution is taken into account. On average, the outermost apertures reach 2.5 times the effective $H$-band radius, corresponding to $\sim$15 and $\sim$12 extracted data points for galaxies at $z\sim 0.9$ and $z\sim 2.3$, respectively. A dynamical mass modelling is performed by fitting the extracted kinematic profiles simultaneously in observed space using an updated version of dysmal (Cresci et al., 2009; Davies et al., 2011).

The free model parameters are the dynamical mass $M_{\rm dyn}$ and the intrinsic velocity dispersion $\sigma_{0}$. The inclination $i$ and effective radius $R_{e}$ are independently constrained from galfit (Peng et al., 2010) models to the CANDELS $H$-band imaging by HST presented by van der Wel et al. (2012). The inclination is computed as ${\rm cos}(i)=[(q^{2}-q_{0}^{2})/(1-q_{0}^{2})]^{1/2}$. Here, $q=b/a$ is the axial ratio, and $q_{0}=0.25$ is the assumed ratio of scale height to scale length, representing the intrinsic thickness of the disk. The width of the point spread function (PSF) is determined from the average PSF during observations for each galaxy. The mass model used in the fitting procedure is a thick exponential disk, following Noordermeer (2008), with a S rsic index of $n_{\rm S}=1$. We note that the peak rotation velocity of a thick exponential disk is about 3 to 8 per cent lower than that of a Freeman disk (Freeman, 1970). For a general comparison of observed and modelled rotation velocities and dispersions, we refer the reader to W16. Another key product of the modelling is the baryonic (or DM) mass fraction on galactic scales, as presented in W16.

The merit of the W16 modelling procedure includes the coupled treatment of velocity and velocity dispersion in terms of beam-smearing effects and pressure support. The latter is of particular importance for our study since high-$z$ galaxies have a non-negligible contribution to their dynamical support from turbulent motions (Förster Schreiber et al., 2006, 2009; Genzel et al., 2006, 2008, 2014a; Kassin et al., 2007, 2012; Cresci et al., 2009; Law et al., 2009; Gnerucci et al., 2011; Epinat et al., 2012; Swinbank et al., 2012; Wisnioski et al., 2012, 2015; Jones et al., 2013; Newman et al., 2013). The resulting pressure compensates part of the gravitational force, leading to a circular velocity which is larger than the rotation velocity $v_{\rm rot}$ alone:

[TABLE]

where $R_{d}$ is the disk scale length (Burkert et al., 2010; see also Burkert et al., 2016; Wuyts et al., 2016; Genzel et al., 2017; Lang et al., 2017).

If not stated otherwise, we adopt the maximum of the modelled circular velocity, $v_{\rm{circ,max}}\equiv v_{\rm{circ}}$, as the rotation velocity measure for our Tully-Fisher analysis. For associated uncertainties, see § 4.3.2. We use an expression for the peak velocity because there is strong evidence that high-$z$ rotation curves of massive star forming disk galaxies exhibit on average an outer fall-off, i.e. do not posses a ‘flat’ part (van Dokkum et al., 2015; Genzel et al., 2017; Lang et al., 2017). This is partly due to the contribution from turbulent motions to the dynamical support of the disk, and partly due to baryons dominating the mass budget on the galaxy scale at high redshift (see also van Dokkum et al., 2015; Stott et al., 2016; Wuyts et al., 2016; Price et al., 2016; Alcorn et al., 2016; Pelliccia et al., 2017). A disk model with a flattening or rising rotation curve as the ‘arctan model’, which is known to be an adequate model for local disk galaxies (e.g. Courteau, 1997), might therefore be a less appropriate choice for high-$z$ galaxies.

2.4 Sample selection

We start our investigation with a parent sample of 240 KMOS3D galaxies selected and modelled by W16. The sample definition is described in detail by W16, and we briefly summarize the main selection criteria here: (i) galaxies exhibit a continuous velocity gradient along a single axis, the ‘kinematic major axis’; (ii) their photometric major axis as determined from the CANDELS WFC3 $H$-band imaging and kinematic major axis are in agreement within 40 degrees; (iii) they have a signal-to-noise ratio within each $0\farcs 8$ diameter aperture along the kinematic major axis of $S/N\gtrsim 5$, with up to $S/N\sim 10-100$ within the central apertures. The galaxies sample a parameter space along the main sequence of star forming galaxies (MS) with stellar masses of $M_{*}\gtrsim 6.3\times 10^{9}~{}M_{\odot}$, specific star formation rates of sSFR $\gtrsim 0.7/t_{\rm Hubble}$, and effective radii of $R_{e}\gtrsim 2$ kpc. The W16 sample further excludes galaxies with signs of major merger activity based on their morphology and/or kinematics.

For our Tully-Fisher analysis we undertake a further detailed examination of the W16 parent sample. The primary selection step is based on the position-velocity diagrams and on the observed and modelled one-dimensional kinematic profiles of the galaxies. Through inspection of the diagrams and profiles we ensure that the peak rotation velocity is well constrained, based on the observed flattening or turnover in the rotation curve and the coincidence of the dispersion peak within $\lesssim 2$ pixels ($\lesssim 0\farcs 4$) with the position of the steepest velocity gradient. The requirement of detecting the maximum velocity is the selection step with the largest effect on sample size, leaving us with 149 targets. The galaxy shown in the fourth row of Figure 1 is excluded from the TFR sample based on this latter requirement.

To single out rotation-dominated systems for our purpose, we next perform a cut of $v_{\rm rot,max}/\sigma_{0}>\sqrt{4.4}$, based on the properties of the modelled galaxy (see also e.g. Tiley et al., 2016). Our cut removes ten more galaxies where the contribution of turbulent motions at the radius of maximum rotation velocity, which is approximately at $r=2.2\,R_{d}$, to the dynamical support is higher than the contribution from ordered rotation (cf. Equation (1)).

We exclude four more galaxies with close neighbours because their kinematics might be influenced by the neighbouring objects. These objects have projected distances of $<20$ kpc, spectroscopic redshift separations of $<300$ km/s, and mass ratios of $>1:5$, based on the 3D-HST catalogue. One of the dismissed galaxies is shown in the third row of Figure 1.

After applying the above cuts, our refined TFR sample contains 135 galaxies, with 65, 24, 46 targets in the $YJ,H,K$ passbands with mean redshifts of $z\sim 0.9,1.5,2.3$, respectively. The median and central 68th percentile ranges of offsets between the morphological and kinematic position angle (PA) are $6.4^{\circ}$ $[0.1^{\circ};18.4^{\circ}]$. This should minimize the possible impact of non-axisymmetric morphological features on the fixed model parameters ($R_{e}$, sin($i$), PA) that are based on single-component S rsic model fits to the observed $H$-band images (see Rodrigues et al., 2017, and also the discussion by W16). The median physical properties of redshift subsamples are listed in Table 1. Individual properties of galaxies in the TFR sample in terms of $z$, $M_{*}$, $M_{\rm bar}$, $v_{\rm circ,max}$, and $\sigma_{0}$, are listed in Table LABEL:tab:ind-properties.

To visualize the impact of our sample selection we show in Figure 2 a ‘first order’ sTFR of all detected and resolved KMOS3D galaxies. Here, $v_{\rm circ}$ is computed from the observed maximal velocity difference and from the intrinsic velocity dispersion as measured from the outer disk region, after corrections for beam-smearing and inclination, as detailed in Appendix A.2 of Burkert et al. (2016). For simplicity, we assume in computing $v_{\rm circ}$ for this figure that the observed maximal velocity difference is measured at $r=2.2R_{d}$, but we emphasize that, in contrast to the modelled circular velocity, this is not necessarily the case. We indicate our parent sample of modelled galaxies by W16 in black, and our final TFR sample in blue. For reference, we also show in orange a subsample of the selection by W16 which is only based on cuts in MS offset ($\pm 0.6$ dex), mass-radius relation offset ($\pm 0.3$ dex), and inclination ($0.5\leq{\rm sin}(i)\leq 0.98$). We emphasize that the assessment of recovering the true maximum rotation velocity is not taken into account for such an objectively selected sample. We discuss in Appendix A in more detail the effects of sample selection, and contrast them to the impact of correcting for e.g. beam-smearing.

The distribution of the TFR sample with respect to the full KMOS3D sample (as of January 2016) and to the corresponding 3D-HST sample in terms of star formation rate and effective radius as a function of stellar mass is shown in Figure 3 (for a detailed comparison of the W16 sample, we refer the reader to W16). We select 3D-HST galaxies with $0.6<z<2.7$, log$(M_{*}~{}[M_{\odot}])>9.2$, $Ks<23$, and for the ‘SFGs only’ subset we apply sSFR $>0.7/t_{\rm Hubble}$, for a total of 9193 and 7185 galaxies, respectively. Focussing on the ‘SFGs only’ subset, the median and corresponding 68th percentiles with respect to the MS relations for the $z\sim 0.9$ and the $z\sim 2.3$ populations are log($\Delta$ MS)=$0.00^{+0.34}_{-0.39}$ and log($\Delta$ MS)=$-0.05^{+0.26}_{-0.35}$, and with respect to the mass-size (M-R) relation log($\Delta$ M-R)=$-0.04^{+0.16}_{-0.28}$ and log($\Delta$ M-R)=$-0.02^{+0.17}_{-0.31}$, respectively. At $z\sim 0.9$, the TFR galaxies lie on average a factor of $\sim 1.6$ above the MS, with log($\Delta$ MS)=$0.20^{+0.42}_{-0.21}$, and have sizes corresponding to log($\Delta$ M-R)=$-0.02^{+0.16}_{-0.17}$. At $z\sim 2.3$, the TFR galaxies lie on average on the MS and M-R relations (log($\Delta$ MS)=$-0.01^{+0.13}_{-0.29}$, log($\Delta$ M-R)=$0.06^{+0.17}_{-0.14}$), but their scatter with respect to higher SFRs and to smaller radii is not as pronounced as for the star-forming 3D-HST sample.

In summary, our analysis accounts for the following effects: (i) beam-smearing, through a full forward modelling of the observed velocity and velocity dispersion profiles with the known instrumental PSF; (ii) the intrinsic thickness of high$-z$ disks, following Noordermeer (2008); (iii) pressure support through turbulent gas motions, following Burkert et al. (2010), under the assumption of a disk of constant velocity dispersion and scale height. The former steps are all included in the dynamical modelling by W16. On top of that, we retain in our TFR sample only non-interacting SFGs which are rotationally supported based on the $v_{\rm rot,max}/\sigma_{0}>\sqrt{4.4}$ criterion, and for which the data have sufficient $S/N$ and spatial coverage to robustly map, and model, the observed rotation curve to or beyond the peak rotation velocity.

3 The TFR with KMOS3D

3.1 Fitting

In general, there are two free parameters for TFR fits in log-log space: the slope $a$ and the zero-point offset $b$. It is standard procedure to adopt a local slope for high$-z$ TFR fits111While the slope might in principle vary with cosmic time, a redshift evolution is not expected from the toy model introduced in Section 1.. This is due to the typically limited dynamical range probed by the samples at high redshift which makes it challenging to robustly constrain $a$. The TFR evolution is then measured as the relative difference in zero-point offsets (e.g. Puech et al., 2008; Cresci et al., 2009; Gnerucci et al., 2011; Miller et al., 2011, 2012; Tiley et al., 2016). In Appendix B we briefly investigate a method to measure TFR evolution which is independent of the slope. For clarity and consistency with TFR investigations in the literature, however, we present our main results based on the functional form of the TFR as given in Equation (2) below. For our fiducial fits, we adopt the local slopes by Reyes et al. (2011) and Lelli et al. (2016) for the sTFR and the bTFR, respectively.222The sTFR zero-point by Reyes et al. (2011) is corrected by $-0.034$ dex to convert their Kroupa (2001) IMF to the Chabrier IMF which is used in this work, following the conversions given in Madau & Dickinson (2014).

To fit the TFR we adopt an inverse linear regression model of the form

[TABLE]

Here, $M$ is the stellar or baryonic mass, and a reference value of $v_{\rm ref}=\overline{v_{\rm circ}}$ is chosen to minimize the uncertainty in the determination of the zero-point $b$ (Tremaine et al., 2002). If we refer in the remainder of the paper to $b$ as the zero-point offset, this is for our sample in reference to $v_{\rm circ}=v_{\rm ref}$, and not to log($v_{\rm circ}$ [km/s])=0. When comparing to other data sets in §§ 3.4 and 4.2 we convert their zero-points accordingly.

For the fitting we use a Bayesian approach to linear regression, as well as a least-squares approximation. The Bayesian approach to linear regression takes uncertainties in ordinate and abscissa into account.333 We use the IDL routine linmix_err which is described and provided by Kelly (2007). A modified version of this code which allows for fixing of the slope was kindly provided to us by Brandon Kelly and Marianne Vestergaard. The least-squares approximation also takes uncertainties in ordinate and abscissa into account, and allows for an adjustment of the intrinsic scatter to ensure for a goodness of fit of $\chi^{2}_{\rm reduced}\approx 1$.444 We use the IDL routine mpfitexy which is described and provided by Williams et al. (2010). It depends on the mpfit package (Markwardt, 2009). To evaluate the uncertainties of the zero-point offset $b$ of the fixed-slope fits, a bootstrap analysis is performed for the fits using the least-squares approximation. The resulting errors agree with the error estimates from the Bayesian approach within $0.005$ dex of mass. We find that the intrinsic scatter obtained from the Bayesian technique is similar or larger by up to $0.03$ dex of mass as compared to the least-squares method. Both methods give the same results for the zero-point $b$ (see also the recent comparison by Bradford et al., 2016).

We perform fits to our full TFR sample, as well as to the subsets at $z\sim 0.9$ and $z\sim 2.3$. The latter allows us to probe the maximum separation in redshift possible within the KMOS3D survey. Due to the low number of TFR galaxies in our $H-$band bin we do not attempt to fit a zero-point at $z\sim 1.5$.

3.2 The TFR at $0.6<z<2.6$

In this paragraph, we investigate the Tully-Fisher properties of our full TFR sample at $0.6<z<2.6$. The sTFR as well as the bTFR are clearly in place and well defined at $0.6<z<2.6$, confirming previous studies (e.g. Cresci et al., 2009; Miller et al., 2011, 2012; Tiley et al., 2016, and other high$-z$ work cited in Section 1). In Figure 4 we show the best fits for the sTFR and the bTFR using the local slopes by Reyes et al. (2011) ($a=1/0.278=3.60$) and Lelli et al. (2016) ($a=3.75$), respectively. The best-fit parameters are given in Table 2.

The intrinsic scatter as determined from the fits is with $\zeta_{\rm int,sTFR}\approx 0.22$ and $\zeta_{\rm int,bTFR}\approx 0.23$ larger by up to a factor of two in dex of mass than in the local Universe (typical values for the observed intrinsic scatter of the local relations used in this study are $\zeta_{\rm int}=0.1-0.13$ in dex of mass; see Reyes et al., 2011; Lelli et al., 2016). A larger scatter in the high$-z$ TFR is expected simply due to the larger measurement uncertainties. It might further be due to disk galaxies being less “settled” (Kassin et al., 2007, 2012; Simons et al., 2016; see also Flores et al., 2006; Puech et al., 2008, 2010; Covington et al., 2010; Miller et al., 2013). This can become manifest through actual displacement of galaxies from the TFR due to a non-equilibrium state (see e.g. simulations by Covington et al., 2010).

Miller et al. (2013) studied the connection between TFR scatter and bulge-to-total ratio, and found that above $z\approx 1$ the TFR scatter is increased due to an offset of bulge-less galaxies from the $B/T>0.1$ galaxy population. $B/T$ measurements for our galaxies come from bulge-disk decompositions based on two-component fits to the two-dimensional CANDELS $H$-band light distribution (Lang et al., 2014). If we select only galaxies with $B/T>0.1$ (57 galaxies), we do not find a decrease in scatter for our sample ($\zeta_{{\rm int,sTFR,}B/T>0.1}=0.22$ and $\zeta_{{\rm int,bTFR,}B/T>0.1}=0.24$). The same is true if we select for galaxies with $B/T<0.1$ (78 galaxies), leading to $\zeta_{{\rm int,sTFR,}B/T<0.1}=0.23$ and $\zeta_{{\rm int,bTFR,}B/T<0.1}=0.22$.

However, the scatter is affected by the sample selection: if we create ‘first order’ TFRs (§ 2.4, Figure 2), i.e. using all detected and resolved KMOS3D galaxies without skyline contamination (316 SFGs), but also without selecting against dispersion-dominated systems, low $S/N$ galaxies, or mergers, we find an intrinsic scatter of $\zeta_{\rm int,sTFR}=0.60$ and $\zeta_{\rm int,bTFR}=0.64$ for these ‘first order’ TFRs (for the parent W16 sample we find $\zeta_{\rm int,sTFR}=0.27$ and $\zeta_{\rm int,bTFR}=0.29$). We caution that this test sample includes galaxies where the maximum rotation velocity is not reached, thus introducing artificial scatter in these ‘first order’ TFRs. In contrast to the properties of our TFR sample, this scatter is asymmetric around the best fit, with larger scatter towards lower velocities, but also towards lower masses where more of the dispersion-dominated galaxies reside (cf. Figures 2 and 8). This underlines the importance of a careful sample selection.

Also the zero-points are affected by the sample selection (see also Figure 8). For our TFR sample, we find $b_{\rm sTFR}=10.50\pm 0.03$ and $b_{\rm bTFR}=10.75\pm 0.03$. If we consider the ‘first order’ samples we find an increase of the zero-points of $\Delta b_{\rm sTFR}=0.37$ dex and $\Delta b_{\rm bTFR}=0.39$ dex (for the parent W16 sample we find $\Delta b_{\rm sTFR}=0.03$ dex and $\Delta b_{\rm bTFR}=0.04$ dex).

It is common, and motivated by the scatter of the TFR, to investigate the existence of hidden parameters in the relation. For example, a measure of the galactic radius (effective, or exponential scale length) has been investigated by some authors to test for correlations with TFR residuals (e.g. McGaugh, 2005; Pizagno et al., 2005; Gnedin et al., 2007; Zaritsky et al., 2014; Lelli et al., 2016). The radius, together with mass, determines the rotation curve (e.g. Equation (D4)). Adopting the local slopes, we do not find significant correlations (based on Spearman tests) of the TFR residuals with $R_{e}$, $B/T$, $n_{\rm S}$, stellar or baryonic mass surface density, offset from the main sequence or the mass-radius relation, SFR surface density $\Sigma_{\rm SFR}$, or inclination. In Appendix C we investigate how the uncertainties in stellar and baryonic mass affect second-order parameter dependencies for TFR fits with free slopes, by example of $R_{e}$ and $\Sigma_{\rm SFR}$.

In summary, we find well defined mass-based TFRs at $0.6<z<2.6$ for our sample. If galaxies with underestimated peak velocity, dispersion-dominated and disturbed galaxies are included, the TFR zero-points are increasing, and also the scatter increases, especially towards lower velocities and masses. Adopting the local slopes, we find no correlation of TFR residuals with independent galaxy properties.

3.3 TFR evolution from $z\sim 2.3$ to $z\sim 0.9$

We now turn to the TFR subsamples at $z\sim 0.9$ and $z\sim 2.3$. We adopt the local slopes by Reyes et al. (2011) and Lelli et al. (2016) to investigate the zero-point evolution. Our redshift subsamples are shown in Figure 5 for the sTFR (left) and bTFR (right), together with the corresponding local relations and the respective fixed-slope fits. The parameters of each fit are given in Table 2.

For the sTFR we find no indication for a significant change in zero-point between $z\sim 0.9$ and $z\sim 2.3$ within the best fit uncertainties. Using the local slope of $a=3.60$ and the reference value $v_{\rm ref}=242$ km/s, we find a zero-point of $b=10.49\pm 0.04$ for the subsample at $z\sim 0.9$, and of $b=10.51\pm 0.05$ for the subsample at $z\sim 2.3$, translating into a zero-point evolution of $\Delta b=0.02$ dex between $z\sim 0.9$ and $z\sim 2.3$.

For the bTFR, however, using the local slope of $a=3.75$, and again the reference value $v_{\rm ref}=242$ km/s, we find a positive zero-point evolution between $z\sim 0.9$ and $z\sim 2.3$, with $b=10.68\pm 0.04$ and $b=10.85\pm 0.05$, respectively, translating into a zero-point evolution of $\Delta b=0.17$ dex between $z\sim 0.9$ and $z\sim 2.3$.

If we consider the ‘first order’ TFR subsamples at $z\sim 0.9$ and $z\sim 2.3$, we find significantly different zero-point evolutions of $\Delta b_{\rm sTFR}=0.23$ dex and $\Delta b_{\rm bTFR}=0.28$ dex between $z\sim 0.9$ and $z\sim 2.3$. Again, this highlights the importance of a careful sample selection for TFR studies. Figure 9 shows that if instead we extend our data set to the sample from W16, we find qualitatively the same trends as for the adopted TFR sample, namely an evolution of $\Delta b_{\rm sTFR}=0.05$ dex and $\Delta b_{\rm bTFR}=0.20$ dex for the zero-point between $z\sim 0.9$ and $z\sim 2.3$ (see Appendix A). Also, if we consider only TFR galaxies with $B/T>0.1(<0.1)$, our qualitative results remain the same.

In summary, we find no evolution for the sTFR, but a positive evolution of the bTFR between $z\sim 0.9$ and $z\sim 2.3$. If galaxies with underestimated peak velocity, dispersion-dominated and disturbed galaxies are included, we find positive evolution of both the sTFR and the bTFR.

3.4 Comparison to other high$-z$ studies

At $z\sim 0.9$ we compare our sTFR (65 KMOS3D galaxies) to the work by Tiley et al. (2016) and Miller et al. (2011). Tiley et al. (2016) have investigated the sTFR at $z\sim 0.9$ using 56 galaxies from the KROSS survey with KMOS (Stott et al., 2016). Miller et al. (2011, 2012) have presented an extensive slit-based sTFR study at $0.2<z<1.7$ with 37 galaxies at $z\sim 1$. From Tiley et al. (2016), we use their best fixed-slope fit to their disky subsample ($a=3.68$). From Miller et al. (2011), we use the $z\sim 1$ fit corresponding to total stellar mass and $v_{\rm rot,3.2}$ ($a=3.78$). For a sTFR comparison at $z\sim 2.3$ (46 KMOS3D galaxies), we consider the work by Cresci et al. (2009). The authors have studied the sTFR at $z\sim 2.2$ for 14 galaxies from the SINS survey ($a=4.5$). Despite the small sample size, the high-quality data based on the 2D modelling of velocity and velocity dispersion maps qualify the sample for comparison with our findings in the highest redshift bin.

In the following, we use $v_{\rm rot,max}$ to ensure a consistent comparison with the measurements presented in these studies. For a comparison with the literature data, we make the simplifying assumption that $v_{\rm rot,max}$ is comparable to $v_{\rm rot,80}$ and $v_{\rm rot,3.2}$ (see § 4.3.3 for a discussion). We adopt the slopes reported in the selected studies to guarantee consistency in the determination of zero-point offsets. The results are shown in Figure 6 as dashed lines, while the original relations from the literature are shown as solid lines. The difference in zero-points, $\Delta b$, is then computed as the zero-point from the KMOS3D fixed-slope fit minus the zero-point from the literature. Given the typical zero-point uncertainty of our fits of $\delta b\approx 0.05$ dex, our results are in agreement with Tiley et al. (2016) ($\Delta b=0.06$) and Cresci et al. (2009) ($\Delta b=0.07$), but in disagreement with Miller et al. (2011) at $z\sim 1$ ($\Delta b=-0.31$). We further note that our findings are in disagreement with the recent study by Di Teodoro et al. (2016) who employed a tilted ring model on a small subset of galaxies from the KMOS3D and KROSS surveys at $z\sim 1$ ($\Delta b=-0.34$; see also Tiley et al., 2016).

A number of complications might give rise to conflicting results of different TFR studies, such as the use of various kinematic models, velocity tracers, mass estimates, or statistical methods. Tiley et al. (2016), who present an extensive comparison of several sTFR studies from the literature, argue that conflicting results regarding the zero-point evolution with redshift depend on the ability of the studies to select for rotationally supported systems. The two-dimensional information on the velocity and velocity dispersion fields is a major advantage of IFS observations as it allows for the robust determination of the kinematic center and major axis.

We test the case of selecting against dispersion-dominated or disturbed systems for our TFR samples. For the full sample of 240 SFGs by W16, which includes some dispersion-dominated systems and cases where the peak rotation velocity might be underestimated by the model, we indeed find that the difference in zero-point, $\Delta b$, with Miller et al. (2011) shrinks by $\sim 30$ per cent. If we now even turn to the purely observational ‘first order sTFR’, this time using only the $z<1.3$ galaxies (122 SFGs) and the $v_{\rm rot,max}$ tracer, we find agreement to Miller et al. (2011) ($\Delta b=0.02$). Again, we caution that this ‘first order’ sample contains not only dispersion-dominated and merging galaxies, but also galaxies for which the maximum velocity is underestimated. This exercise supports the interpretation that the disagreement with Miller et al. (2011) is partly due to our selection of rotation-dominated systems. Beam-smearing corrections could lead to effects of comparable order, as is discussed in more detail in Appendix A and explicitly shown in Figure 8.

The high$-z$ evolution of the bTFR has received less attention in the literature. At intermediate redshift ($z\sim 1.2$), Vergani et al. (2012) found no evolution of the bTFR when comparing to the local relation by McGaugh (2005). We compare our results to the slit-based relation at $z\sim 2$ by Price et al. (2016) using galaxies from the MOSDEF survey (Kriek et al., 2015). Price et al. (2016) use the $S_{0.5}=(0.5\cdot v_{\rm rot}^{2}+\sigma_{g}^{2})^{1/2}$ velocity tracer, which also incorporates dynamical support from disordered motions based on the assumption of isotropic (or constant) gas velocity dispersion $\sigma_{g}$ (Weiner et al., 2006; Kassin et al., 2007). Price et al. (2016) show a plot of the $S_{0.5}-$bTFR of 178 SFGs, of which 35 (15) have detected (resolved) rotation measurements. For resolved galaxies, $S_{0.5}$ is obtained through combining a constant intrinsic velocity dispersion, and $v_{\rm rot,2.2}$. For unresolved galaxies, Price et al. (2016) estimate $S_{0.5}$ through an rms velocity (see their Appendix B for details). We use their fixed-slope fit ($a=1/0.39$) to compare their results to our 46 KMOS3D galaxies at $z\sim 2.3$ in the right panel of Figure 6. Our fixed-slope fit is in agreement with the result by Price et al. (2016) ($\Delta b=-0.03$). This is surprising at first, given the above discussion of IFS vs. slit-based rotation curve measurements, and the fact that the Price et al. (2016) sample contains a large fraction of objects without detected rotation. However, Price et al. (2016) state that their findings regarding the S0.5-bTFR do not change if they consider only the galaxies with detected rotation measurements. This is likely due to the detailed modelling and well-calibrated translation of line width to rotation velocity by the authors. In general, any combination of velocity dispersion and velocity into a joined measure is expected to bring turbulent and even dispersion-dominated galaxies closer together in TFR space, which might further serve as an explanation for this good agreement (see also Covington et al., 2010).555 Partly, this is also the case for the measurements by Miller et al. (2011, 2012), if a correction for turbulent pressure support is performed. Since their velocity dispersions are not available to us, however, only an approximative comparison is feasible. From this, we found agreement of their highest redshift bin ($z\sim 1.5$) with our $0.6<z<2.6$ data in the $v_{\rm circ}$-sTFR plane, but still a significant offset at $z\sim 1$.

In summary, our inferred $v_{\rm rot}$-sTFR zero-points (i.e., not corrected for pressure support) agree with the work by Cresci et al. (2009) and Tiley et al. (2016), but disagree with the work by Miller et al. (2011). Our $S_{0.5}$-bTFR zero-point agrees with the result by Price et al. (2016). We emphasize that the negligence of turbulent motions in the balance of forces leads to a relation which has lost its virtue to directly connect the baryonic kinematics to the central potential of the halo.

4 TFR evolution in context

4.1 Dynamical support of SFGs from $z\sim 2.3$ to $z\sim 0.9$

At fixed $v_{\rm circ}$, our sample shows higher $M_{\rm bar}$ and similar $M_{*}$ at $z\sim 2.3$ as compared to $z\sim 0.9$ (Figure 5). Galactic gas fractions are strongly increasing with redshift, as it has become clear in the last few years (Tacconi et al., 2010; Daddi et al., 2010; Combes et al., 2011; Genzel et al., 2015; Tacconi et al., 2017). In our TFR sample, the baryonic mass of the $z\sim 2.3$ galaxies is on average a factor of two larger as compared to $z\sim 0.9$, while stellar masses are comparable. The relative offset at fixed $v_{\rm circ}$ of our redshift subsamples in the bTFR plane, which is not visible in the sTFR plane, confirms the relevance of gas at high redshift.

Building on the recent work by W16 on the mass budgets of high$-z$ SFGs, we can identify through our Tully-Fisher analysis another redshift-dependent ingredient to the dynamical support of high$-z$ SFGs. The sTFR zero-point does not evolve significantly between $z\sim 2.3$ and $z\sim 0.9$. Since we know that there is less gas in the lower$-z$ SFGs, the ‘missing’ baryonic contribution to the dynamical support of these galaxies as compared to $z\sim 2.3$ has to be compensated by DM. We therefore confirm with our study the increasing importance of DM to the dynamical support of SFGs (within $\sim 1.3\,R_{e}$) through cosmic time. This might be partly due to the redshift dependence of the halo concentration parameters, which decrease with increasing redshift. In the context of the toy model mentioned in Section 1, it is indeed the case that a decrease of the DM fraction as probed by the central galaxy with increasing redshift can flatten out or even reverse the naively expected, negative evolution of the TFR offset with increasing redshift. This will be discussed in more detail in Section 5.

The increase of baryon fractions with redshift is supported by other recent work: W16 find that the baryon fractions of SFGs within $R_{e}$ increase from $z\sim 1$ to $z\gtrsim 2$, with galaxies at higher redshift being clearly baryon-dominated (see also Förster Schreiber et al., 2009; Alcorn et al., 2016; Price et al., 2016; Burkert et al., 2016; Stott et al., 2016; Contini et al., 2016). W16 also find that the baryonic mass fractions are correlated with the baryonic surface density within $R_{e}$, suggesting that the lower surface density systems at lower redshift are more diffuse and therefore probe further into the halo (consequently increasing their DM fraction). Most recently, Genzel et al. (2017) find in a detailed study based on the outer rotation curves of six massive SFGs at $z=0.9-2.4$ that the three $z>2$ galaxies are most strongly baryon-dominated. On a statistical basis, this is confirmed through stacked rotation curves of more than 100 high$-z$ SFGs by Lang et al. (2017).

Given the average masses of our galaxies in the $YJ$ and $K$ subsamples, we emphasize that we are generally not tracing a progenitor-descendant population in our sample, since the average stellar and baryonic masses of the $z\sim 2.3$ galaxies are already higher than for those at $z\sim 0.9$ (Table 1). It is very likely that a large fraction of the massive star-forming disk galaxies we observe at $z\gtrsim 1$ have evolved into early-type galaxies (ETGs) by $z=0$, as discussed in the recent work by Genzel et al. (2017). Locally, there is evidence that ETGs have high SFRs at early times, with the most massive ETGs forming most of their stars at $z\gtrsim 2$ (e.g. Thomas, 2010; McDermid et al., 2015). This view is supported by co-moving number density studies (e.g. Brammer et al., 2011), which also highlight that the mass growth of today’s ETGs after their early and intense SF activity is mainly by the integration of (stellar) satellites into the outer galactic regions (van Dokkum et al., 2010). The observed low DM fractions of the massive, highest$-z$ SFGs seem to be consistent with the early assembly of local ETGs, with rapid incorporation of their baryon content. In future work, we will compare our observations to semi-analytical models and cosmological zoom-in simulations to investigate in greater detail the possible evolutionary scenarios of our observed galaxies in the context of TFR evolution.

4.2 Comparison to the local Universe

In Figure 5 we show the TFR zero-point evolution in context with the recent local studies by Reyes et al. (2011) for the sTFR, and by Lelli et al. (2016) for the bTFR. Reyes et al. (2011) study the sTFR for a large sample of 189 disk galaxies, using resolved H$\alpha$ rotation curves. Lelli et al. (2016) use resolved Hi rotation curves and derive a bTFR for 118 disk galaxies. To compare these local measurements to our high$-z$ KMOS3D data, we assume that at $z\approx 0$ the contribution from turbulent motions to the dynamical support of the galaxy is negligible, and therefore $v_{\rm circ}\equiv v_{\rm rot}$. We make the simplifying assumption that $v_{\rm circ}$ is comparable to $v_{80}$ and $v_{\rm flat}$ used by Reyes et al. (2011) and Lelli et al. (2016), respectively (see § 4.3.3 for a discussion). From Lelli et al. (2016), we use the fit to their subsample of 58 galaxies with the most accurate distances (see their classification).

For the sTFR as well as the bTFR we find significant offsets of the high$-z$ relations as compared to the local ones, namely $\Delta b_{\rm sTFR,z\sim 0.9}=-0.44$, $\Delta b_{\rm sTFR,z\sim 2.3}=-0.42$, $\Delta b_{\rm bTFR,z\sim 0.9}=-0.44$ and $\Delta b_{\rm bTFR,z\sim 2.3}=-0.27$. We have discussed in §§ 3.2 and 3.3 the zero-points of the ‘first order’ TFRs as compared to our fiducial TFRs: while there is significant offset for both the ‘first order’ sTFR and bTFR when comparing the $z\sim 0.9$ and the $z\sim 2.3$ subsamples, the overall offset to the local relations is reduced. The difference between the local relations and the full ‘first order’ samples is only $\Delta b_{\rm sTFR}=-0.06$ and $\Delta b_{\rm bTFR}=0.02$, which would be consistent with no or only marginal evolution of the TFRs between $z=0$ and $0.6<z<2.3$.

For the interpretation of the offsets to the local relations, it is important to keep in mind that we measure the TFR evolution at the typical fixed circular velocity of galaxies in our high$-z$ sample. This traces the evolution of the TFR itself through cosmic time, not the evolution of individual galaxies. Our subsamples at $z\sim 0.9$ and $z\sim 2.3$ are representative of the population of massive MS galaxies observed at those epochs, with the limitations as discussed in § 2.4. Locally, however, the typical disk galaxy has lower circular velocity than our adopted reference velocity, and consequently lower mass (cf. e.g. Figure 1 by Courteau & Dutton, 2015). Figure 5 does therefore not indicate how our galaxies will evolve on the TFR from $z\sim 2$ to $z\sim 0$, but rather shows how the relation itself evolves, as defined through the population of disk galaxies at the explored redshifts and mass ranges. This is also apparent if actual data points of low- and high-redshift disk galaxies are shown together. We show a corresponding plot for the bTFR in Appendix B.

In summary, our results suggest an evolution of the TFR with redshift, with zero-point offsets as compared to the local relations of $\Delta b_{\rm sTFR,z\sim 0.9}=-0.44$, $\Delta b_{\rm sTFR,z\sim 2.3}=-0.42$, $\Delta b_{\rm bTFR,z\sim 0.9}=-0.44$ and $\Delta b_{\rm bTFR,z\sim 2.3}=-0.27$. If galaxies with underestimated peak velocity, dispersion-dominated and disturbed galaxies are included, the overall evolution between the $z=0$ and $0.6<z<2.6$ samples is insignificant.

4.3 The impact of uncertainties and model assumptions on the observed TFR evolution

Before we interpret our observed TFR evolution in a cosmological context in Section 5, we discuss in the following uncertainties and modelling effects related to our data and methods. We find that uncertainties of mass estimates and velocities cannot explain the observed TFR evolution. Neglecting the impact of turbulent motions, however, could explain some of the tension with other work.

4.3.1 Uncertainties of stellar and baryonic masses

A number of approximations go into the determination of stellar and baryonic masses at high redshift. Simplifying assumptions like a uniform metallicity, a single IMF, or an exponentially declining SFH introduce significant uncertainties to the stellar age, stellar mass, and SFR estimates of high$-z$ galaxies. While the stellar mass estimates appear to be more robust against variations in the model assumptions, the SFRs, which are used for the molecular gas mass calculation, are affected more strongly (see e.g. Förster Schreiber et al., 2004; Shapley et al., 2005; Wuyts et al., 2007, 2009, 2016; Maraston et al., 2010; Mancini et al., 2011, for detailed discussions about uncertainties and their dependencies). Most systematic uncertainties affecting stellar masses tend to lead to underestimates; if this were the case for our high$-z$ samples, the zero-point evolution with respect to local samples would be overestimated. However, the dynamical analysis by W16 suggests that this should only be a minor effect, given the already high baryonic mass fractions at high redshift.

An uncertainty in the assessment of gas masses at high redshift is the unknown contribution of atomic gas. In the local Universe, the gas mass of massive galaxies is dominated by atomic gas: for stellar masses of log$(M_{*}~{}[M_{\odot}])\approx 10.5$, the ratio of atomic to molecular hydrogen is roughly $M_{\rm H{\textsc{i}}}/M_{\rm H_{2}}\sim 3$ (e.g. Saintonge et al., 2011). While there are currently no direct galactic Hi measurements available at high redshift,666 But see e.g. Wolfe et al. (2005); Werk et al. (2014) for measurements of Hi column densities of the circum- and intergalactic medium using quasar absorption lines. From these techniques, a more or less constant cosmological mass density of neutral gas since at least $z\sim 3$ is inferred (e.g. Péroux et al., 2005; Noterdaeme et al., 2009). Recently, the need for a significant amount of non-molecular gas in the haloes of high$-z$ galaxies has also been invoked by the environmental study of the 3D-HST fields by Fossati et al. (2017). a saturation threshold of the Hi column density of only $\lesssim 10\,M_{\odot}/{\rm pc}^{2}$ has been determined empirically for the local Universe (Bigiel & Blitz, 2012). The much higher gas surface densities of our high$-z$ SFGs therefore suggest a negligible contribution from atomic gas within $r\lesssim R_{e}$ (see also W16, ). Consequently, the contribution of atomic gas to the maximum rotation velocity and to the mass budget within this radius should be negligible. However, there is evidence that locally Hi disks are much more extended than optical disks (e.g. Broeils & Rhee, 1997). If this is also true at high redshift, the total galactic Hi mass fractions could still be significant at $z\sim 1$, as is predicted by theoretical models (e.g. Lagos et al., 2011; Fu et al., 2012; Popping et al., 2015). Due to the lack of empirical confirmation, however, these models yet remain uncertain, especially given that they under-predict the observed high$-z$ molecular gas masses by factors of $2-5$. Within these limitations, we perform a correction for missing atomic gas mass at high$-z$ in our toy model discussion in Section 5.

Following Burkert et al. (2016), we have adopted uncertainties of $0.15$ dex for stellar masses, and $0.20$ dex for gas masses. This translates into an average uncertainty of $\sim 0.15$ dex for baryonic masses. These choices likely underestimate the systematic uncertainties in the error budget which can have a substantial impact on some of our results, because the slope as well as the scatter of the TFR are sensitive to the uncertainties. For the presentation of our main results, we adopt local TFR slopes, thus mitigating these effects. In Appendix C, we explore the effect of varying mass uncertainties on free-slope fits of the TFR, together with implications on TFR residuals and evolution. We find that measurements of the zero-point are little affected by the uncertainties on mass, to an extent much smaller than the observed bTFR evolution between $z\sim 2.3$ and $z\sim 0.9$.

4.3.2 Uncertainties of circular velocities

We compute the uncertainties of the maximum circular velocity as the propagated errors on the observed velocity and $\sigma_{0}$, including an uncertainty on $q$ of $\sim 20$ per cent. The latter is a conservative choice in the light of the current KMOS3D magnitude cut of $Ks<23$ (cf. van der Wel et al., 2012). For details about the observed quantities, see W15, and W16 for a comparison between observed and modelled velocities and velocity dispersions. The resulting median of the propagated circular velocity uncertainty is 20 km/s.

Maximum circular velocities can be systematically underestimated: although the effective radius enters the modelling procedure as an independent constraint, the correction for pressure support can lead to an underestimated turn-over radius if the true turn-over radius is not covered by observations. For our TFR sample we selected only galaxies where modelled and observed velocity and dispersion profiles are in good agreement, and where the maximum or flattening of the rotation curve is covered by observations. It is therefore unlikely that our results based on the TFR sample are affected by systematic uncertainties of the maximum circular velocity.

4.3.3 Effects related to different velocity measures and models

The different rotation velocity models and measures used in the literature might affect comparisons between different studies. Some TFR studies adopt the rotation velocity at 2.2 times $R_{d}$, $v_{2.2}$, as their fiducial velocity to measure the TFR. We verified that for the dynamical modelling as described above, $v_{\rm{circ,2.2}}$ equals $v_{\rm{circ,max}}$, and $v_{\rm{rot,2.2}}$ equals $v_{\rm{rot,max}}$ with an average accuracy of $\lesssim 1$ km/s. Other commonly used velocity measures are $v_{\rm flat}$, $v_{3.2}$, and $v_{80}$, the rotation velocity at the radius which contains 80 per cent of the stellar light. For a pure exponential disk, this corresponds to roughly $v_{3.0}$ (Reyes et al., 2011). It has been shown by Hammer et al. (2007) that $v_{\rm flat}$ and $v_{80}$ are comparable in local galaxies. For the exponential disk model including pressure support which we use in our analysis, $v_{\rm{rot(circ),max}}$ is on average $\lesssim 15(10)$ km/s larger than $v_{\rm{rot(circ),3.2}}$. Since $v_{3.2}$ and $v_{80}$ are, however, usually measured from an ‘arctan model’ with an asymptotic maximum velocity (Courteau, 1997), reported values in the literature generally do not correspond to the respective values at these radii from the thick exponential disk model with pressure support. Miller et al. (2011) show that for their sample of SFGs at $0.2<z<1.3$, the typical difference between $v_{2.2}$ and $v_{3.2}$, as computed from the arctan model, is on the order of a few per cent (see also Reyes et al., 2011). This can also be assessed from Figure 6 by Epinat et al. (2010), who show examples of velocity fields and rotation curves for different disk models (exponential disk, isothermal sphere, ‘flat’, arctan). By construction, the peak velocity of the exponential disk is higher than the arctan model rotation velocity at the corresponding radius.

We conclude that our TFR ‘velocity’ values derived from the peak rotation velocity of a thick exponential disk model are comparable to $v_{\rm flat}$, and close to $v_{3.2}$ and $v_{80}$ from an arctan model, with the limitations outlined above. The possible systematic differences of $<20$ km/s between the various velocity models and measures cannot explain the observed evolution between $z=0$ and $0.6<z<2.6$.

Another effect on the shape of the velocity and velocity dispersion profiles is expected if contributions by central bulges are taken into account. We have tested for a sample of more than 70 galaxies that the effect of including a bulge on our adopted velocity tracer, $v_{\rm circ,max}$ is on average no larger than 5 per cent. From our tests, we do not expect the qualitative results regarding the TFR evolution between $z\sim 2.3$ and $z\sim 0.9$ presented in this paper to change if we include bulges into the modelling of the mass distribution.

4.3.4 The impact of turbulent motions

The dynamical support of star-forming disk galaxies can be quantified through the relative contributions from ordered rotation and turbulent motions (see also e.g. Tiley et al., 2016). We consider only rotation-dominated systems in our TFR analysis, namely galaxies with $v_{\rm rot,max}/\sigma_{0}>\sqrt{4.4}$. Because of this selection, the effect of $\sigma_{0}$ on the velocity measure is already limited, with median values of $v_{\rm rot,max}=233$ km/s at $z\sim 0.9$, and $239$ km/s at $z\sim 2.3$, *vs. *median values of $v_{\rm circ,max}=239$ and $v_{\rm circ,max}=260$ km/s at $z\sim 0.9$ and $z\sim 2.3$, respectively (Table 1).

However, this difference translates into changes regarding e.g. the TFR scatter: for the $v_{\rm rot,max}$-TFR, we find a scatter of $\zeta_{\rm int,sTFR}=0.28$ and $\zeta_{\rm int,bTFR}=0.31$ at $z\sim 0.9$, and at $z\sim 2.3$ we find $\zeta_{\rm int,sTFR}=0.33$ and $\zeta_{\rm int,bTFR}=0.33$, with those values being consistently higher than the values reported for the $v_{\rm circ,max}$-TFR sample in Table 2. More significantly, neglecting the contributions from turbulent motions affects the zero-point evolution: without correcting $v_{\rm rot,max}$ for the effect of pressure support, we would find $\Delta b_{\rm sTFR,z\sim 0.9}=-0.34$, $\Delta b_{\rm sTFR,z\sim 2.3}=-0.26$, $\Delta b_{\rm bTFR,z\sim 0.9}=-0.33$ and $\Delta b_{\rm bTFR,z\sim 2.3}=-0.09$. The inferred zero-points at higher redshift are affected more strongly by the necessary correction for pressure support (cf. Figure 5).

These results emphasize the increasing role of pressure support with increasing redshift, confirming previous findings by e.g. Förster Schreiber et al. (2009); Epinat et al. (2009); Kassin et al. (2012); W15. It is therefore clear that turbulent motions must not be neglected in kinematic analyses of high$-z$ galaxies. If the contribution from pressure support to the galaxy dynamics is dismissed, this will lead to misleading conclusions about TFR evolution in the context of high$-z$ and local measurements.

5 A toy model interpretation

The relative comparison of our $z\sim 2.3$ and $z\sim 0.9$ data and local relations indicates a non-monotonic evolution of the bTFR zero-point with cosmic time (Figure 5). In this section, we present a toy model interpretation of our results, aiming to explain the redshift evolution of both the sTFR and the bTFR, in particular the relative zero-point offsets at $z\sim 2.3$, $z\sim 0.9$, and $z\sim 0$.

The basic premise is that galaxies form at the centers of DM haloes. A simple model for a DM halo in approximate equilibrium is a truncated isothermal sphere, limited by the radius $R_{h}$ where the mean density equals 200 times the critical density of the Universe. The corresponding redshift-dependent relations between halo radius, mass $M_{h}$, and circular velocity $V_{h}$ are

[TABLE]

(Mo et al., 1998), where $H(z)$ is the Hubble parameter, and $G$ is the gravitational constant. The first equation shows that the relation between $M_{h}$ and $V_{h}$ is a smooth function of redshift.

In theory, the relation between these halo properties and corresponding galactic properties can be complex due to the response of the halo to the formation of the central galaxy (see e.g. the discussions on halo contraction *vs. *expansion by Duffy et al., 2010; Dutton et al., 2016; Velliscig et al., 2014). However, recent studies and modelling of high$-z$ SFGs now provide a number of empirical constraints that implicitly contain information on the DM halo profile on galactic scales.

Relations corresponding to Equations (3) for the central baryonic galaxy can then be derived by assuming a direct mapping between the halo and galaxy mass and radius. Information on the inner halo profile is contained in parameters such as the disk mass fraction $m_{d}=M_{\rm bar}/M_{h}$, or the central DM fraction $f_{\rm DM}(r)=v^{2}_{\rm DM}(r)/v^{2}_{\rm circ}(r)$. For our galaxies, we know their stellar mass $M_{*}$ and effective radius $R_{e}$, their baryonic mass $M_{\rm bar}$ and gas mass fraction $f_{\rm gas}=M_{\rm gas}/M_{\rm bar}$ from empirical scaling relations, and their circular velocity $v_{\rm circ}(r)$ and related central DM fraction $f_{\rm DM}(r)$ from dynamical modelling, as detailed in §§ 2.2 and 2.3 and in the references given there. We further have an estimate of their average baryonic disk mass fraction $m_{d}$ (Burkert et al., 2016). We can combine this information to construct a toy model of the TFR zero-point evolution, where we take the redshift dependencies of these various parameters into account (see Appendix D.1 for a detailed derivation):

[TABLE]

[TABLE]

where $C$ and $C^{\prime}$ are constants. Here, we have assumed that, in contrast to the disk mass fraction, the proportionality factor between DM halo radius and galactic radius is independent of redshift (see e.g. Burkert et al., 2016).

Equations (4) and (5) reveal that the TFR evolution can be strongly affected by changes of $f_{\rm DM}(R_{e})$, $m_{d}$, or $f_{\rm gas}$ with redshift, and does not necessarily follow the smooth evolution of the halo parameters given in Equation (3). There have been indications for deviations from a simple smooth TFR evolution scenario in the theoretical work by Somerville et al. (2008). Also the recent observational compilation by Swinbank et al. (2012) showed a deviating evolution (although qualified as consistent with the smooth evolution scenario).

Evaluating Equations (4) and (5) at fixed $v_{\rm circ}(R_{e})$, we learn the following: (i) if $f_{\rm DM}(R_{e})$ decreases with increasing redshift, the baryonic and stellar mass will increase and consequently the TFR zero-point will increase; (ii) if $m_{d}$ increases with increasing redshift, the baryonic and stellar mass will decrease and consequently the TFR zero-point will decrease; (iii) if $f_{\rm gas}$ increases with increasing redshift, the stellar mass will decrease and consequently the sTFR zero-point will decrease. These effects are illustrated individually in Figure 14 in Appendix D.

We constrain our toy model at redshifts $z=0$, $z\sim 0.9$, and $z\sim 2.3$ as follows: the redshift evolution of $f_{\rm gas}$ is obtained through the empirical atomic and molecular gas mass scaling relations by Saintonge et al. (2011) and Tacconi et al. (2017). At fixed circular velocity, $f_{\rm gas}$ evolves significantly with redshift, where $z\sim 2$ galaxies have gas fractions which are about a factor of eight higher than in the local Universe. The redshift evolution of $f_{\rm DM}(R_{e})$ is constrained through the observational results by Martinsson et al. (2013b, a) in the local Universe, and by W16 at $z\sim 0.9$ and $z\sim 2.3$. We tune the redshift evolution of $f_{\rm DM}(R_{e})$ within the ranges allowed by these observations to optimize the match between the toy model and the observed TFR evolution presented in this paper. $f_{\rm DM}(R_{e})$ evolves significantly with redshift, with $z\sim 2$ DM fractions which are about a factor of five lower than at $z=0$. $m_{d}$ is constrained by the abundance matching results by Moster et al. (2013) in the local Universe, whereas at $0.8<z<2.6$ we adopt the value deduced by Burkert et al. (2016). Details on the parametrization of the above parameters are given in Appendix D.2.

In Figure 7 we show how these empirically motivated, redshift-dependent DM fractions, disk mass fractions, and gas fractions interplay in our toy model framework to approximately explain our observed TFR evolution, specifically the TFR zero-point offsets at fixed circular velocity as a function of cosmic time. In particular, this is valid at $z=0$, $z=0.9$, and $z=2.3$, while we have partially interpolated in between. Our observed KMOS3D TFR zero-points of the bTFR (blue squares) and the sTFR (yellow stars) at $z\sim 0.9$ and $z\sim 2.3$ are shown in relation to the local TFRs by Lelli et al. (2016) and Reyes et al. (2011). The horizontal error bars of the KMOS3D data points indicate the spanned range in redshift, while the vertical error bars show fit uncertainties. For this plot, we also perform a correction for atomic gas at high redshift:777Lelli et al. (2016) neglect molecular gas for their bTFR, but state that it has generally a minor dynamical contribution. we follow the theoretical prediction that, at fixed $M_{*}$, the ratio of atomic gas mass to stellar mass does not change significantly with redshift (e.g. Fu et al., 2012). We use the fitting functions by Saintonge et al. (2011) to determine the atomic gas mass for galaxies with log$(M_{*}~{}[M_{\odot}])=10.50$, which corresponds to the average stellar mass of our TFR galaxies at $v_{\rm ref}=242$ km/s in both redshift bins. We find an increase of the zero-point of $+0.04$ dex at $z\sim 0.9$ and $+0.02$ dex at $z\sim 2.3$. This is included in the figure.

We show as green lines our empirically constrained toy model governed by Equations (4) and (5). This model assumes a redshift evolution of $f_{\rm gas}$, $f_{\rm DM}(R_{e})$, and $m_{d}$ as shown by the blue, purple, and black lines, respectively, in inset (a) in Figure 7 (details are given in Appendix D.2). In this model, the increase in $f_{\rm gas}$ is responsible for the deviating (and stronger) evolution of the sTFR as compared to the bTFR. The decrease of $f_{\rm DM}(R_{e})$ is responsible for the upturn/flattening of the bTFR/sTFR evolution. The increase of $m_{d}$ leads to a TFR evolution which is steeper than what would be expected from a model governed only by $H(z)$ (see also Fig. 14). Our toy model evolution is particularly sensitive to changes of $f_{\rm DM}(R_{e})$ with redshift. We illustrate this by showing as cyan shaded areas in Figure 7 how the toy model evolution would vary if we would change only $f_{\rm DM}(R_{e})$ by $\pm 0.1$ at $z=0$, $z=0.9$, and $z=2.3$.

We note that the toy model zero-point offset at $R_{e}$ as derived from Equations (4) and (5), and based on a thin exponential baryon distribution, is comparable to our empirical TFR offset for a thick exponential disk and using $v_{\rm circ,max}$, since the correction factors for the circular velocity measure from thin to thick exponential disk, and from $v_{\rm circ}(R_{e})$ to $v_{\rm circ,max}\approx v_{\rm circ}(r_{2.2})$, are both of the order of $\sim 5$ per cent and approximately compensate one another. The toy model slope ($a=3$) is shallower than our adopted local slopes. In Appendix C we show that the usage of a reference velocity leads to negligible zero-point differences of TFR fits with different slopes.

Although our toy model is not a perfect match to the observed TFR evolution, it reproduces the observed trends reasonably well: for the sTFR, the zero-point decreases from $z=0$ to $z\sim 1$, but there is no or only marginal evolution between $z\sim 1$ and $z\sim 2$. In contrast, there is a significantly non-monotonic evolution of the bTFR zero-point, such that the zero-point first decreases from $z=0$ to $z\sim 1$, and then increases again up to $z\sim 2$. We note that although we show the TFR evolution up to $z=3$, the constraints on $f_{\rm DM}(R_{e})$ and $m_{d}$ are valid only up to $z\approx 2.6$, as indicated in the figure by the grey shading. Also in the redshift range $0\lesssim z\lesssim 0.8$ the model is poorly constrained because we assume a simplistic evolution of $m_{d}$ (cf. Appendix D.2).

A more complete interpretation of our findings also at intermediate redshift has to await further progress in observational work. With the extension of the KMOS3D survey towards lower mass galaxies and towards a more complete redshift coverage in the upcoming observing periods, we might already be able to add in precision and redshift range to our model interpretation. Our current data and models, however, already show the potential of state-of-the-art high$-z$ studies of galaxies to constrain parameters which are important also for theoretical work.

We would like to caution that our proposed model certainly draws a simplified picture. For instance, the assumption of a common scale length of the atomic gas as well as the molecular gas plus stars, as we did for this exercise, can only be taken as approximate, given the high central surface mass densities of our typical high$-z$ galaxies (see § 4.3.1, and W16). Also, the effective radii predicted by our “best fit” toy model are 10-30 per cent larger than what is observed. Other factors not addressed in our approach might also come into play: we did not explore in detail the possible effects of varying halo spin parameter $\lambda$ or of the ratio between baryonic and DM specific angular momenta $j_{\rm bar}/j_{\rm DM}$, which commonly relate $R_{h}$ to $R_{d}$. We also note that possible conclusions on the NFW halo concentration parameter $c$ are in tension with current models (cf. Appendix D.2). We therefore caution that our proposed toy model perspective can only reflect general trends, in particular the relative TFR zero-point offsets at $z=0$, $z=0.9$, and $z=2.3$, and likely misses other relevant ingredients.

Having in mind the limitations outlined above, we conclude that the observed evolution of the mass-based TFRs can be explained in the framework of virialized haloes in an expanding $\Lambda$CDM universe, with galactic DM fractions, disk mass fractions, and gas fractions that are evolving with cosmic time. Adopting the proposed evolution of the model parameters in Equations (4) and (5) as described above and shown in inset (a) in Figure 7, namely at fixed $v_{\rm circ}$ increasing $f_{\rm gas}$ and $m_{d}$, and decreasing $f_{\rm DM}(R_{e})$ with redshift, leads to a redshift evolution of the TFR which is non-monotonic, in particular for the bTFR.

6 Summary

We have investigated the mass-based Tully-Fisher relations (TFRs) of massive star-forming disk galaxies between redshift $z\sim 2.3$ and $z\sim 0.9$ as part of the KMOS3D survey. All our data are reduced and analyzed in a self-consistent way. The spatially resolved nature of our observations enables reliable modelling of individual galaxies, and allows for a careful selection of objects based on kinematic properties and data quality. We have taken into account inclination, beam-smearing, and instrumental broadening, and we have incorporated the significant effects of pressure support to the gravitational potential at these redshifts in our derivation of the circular velocities.

We find that the TFR is clearly in place already at $0.6<z<2.6$ (§ 3.2). Its scatter increases with redshift, but we did not find any second-order parameter dependencies when adopting a local slope. At fixed $v_{\rm circ,max}$, we find higher $M_{\rm bar}$ but similar $M_{*}$ at $z\sim 2.3$ as compared to $z\sim 0.9$ (§ 3.3). This highlights the important effects of the evolution of $f_{\rm gas}$, where, at the same stellar mass, high$-z$ star-forming galaxies (SFGs) have significantly higher gas fractions than lower$-z$ SFGs. This strengthens earlier conclusions by Cresci et al. (2009) in the context of the interpretation of TFR evolution. Since we do not find a significant evolution of the sTFR between $z\sim 2.3$ and $z\sim 0.9$, our observed TFR evolution together with the decrease of $f_{\rm gas}$ with decreasing redshift, implies that the contribution of dark matter (DM) to the dynamical mass on the galaxy scale has to increase with decreasing redshift to maintain the dynamical support of the galaxy as measured through $v_{\rm circ,max}$. Our results complement the findings in other recent work that higher$-z$ SFGs are more baryon-dominated (§ 4.1).

Comparing to other selected high$-z$ TFR studies, we find agreement with the work by Cresci et al. (2009); Price et al. (2016); Tiley et al. (2016), but disagreement with the work by Miller et al. (2011) (§ 3.4). The significant differences in zero-point offsets of our high$-z$ TFRs as compared to the local relations by Reyes et al. (2011) and Lelli et al. (2016) indicate an evolution of the TFR with cosmic time (§ 4.2). From the local Universe to $z\sim 0.9$ and further to $z\sim 2.3$, we find a non-monotonic TFR zero-point evolution which is particularly pronounced for the bTFR.

To explain our observed TFR evolution, we present a toy model interpretation guided by an analytic model of disk galaxy evolution (Section 5). This model takes into account empirically motivated gas fractions, disk mass fractions, and central DM fractions with redshift. We find that the increasing gas fractions with redshift are responsible for the increasingly deviating evolution between the sTFR and the bTFR with redshift. The decreasing central DM fractions with redshift result in the flattening/upturn of the sTFR/bTFR zero-point evolution at $0.9<z<2.3$. This simple model matches our observed TFR evolution reasonably well.

It will be interesting to make more detailed comparisons between the growing amount of observations that can constrain the TFR at high redshift, and the newest generation of simulations and semi-analytical models. Further investigations of galaxies at lower ($z\lesssim 0.7$) and higher ($z\gtrsim 2.5$) redshifts using consistent reduction and analysis techniques will help to unveil the detailed evolution of the mass-based TFR, and to reconcile current tensions in observational work. Another important quest is to provide data which cover wider ranges in velocity and mass at these high redshifts to minimize uncertainties in the fitting of the data, and to investigate if the TFR slope changes with redshift.

We are grateful to the anonymous referee for a constructive report which helped to improve the manuscript. We thank the ESO Paranal staff for their helpful support with the KMOS observations for this work. We are grateful to Jonathan Freundlich, Susan Kassin, Federico Lelli, Raymond Simons, Jakob Walcher, and in particular to Amiel Sternberg and Simon White, for fruitful discussions, and to Rachel Somerville, Mike Williams, and Dennis Zaritsky for valuable insight into various aspects of this work. We thank Sedona Price for providing us with details on the fits by Price et al. (2016), and we are grateful to Brandon Kelly and Marianne Vestergaard for providing us with a modified version of the Bayesian approach to linear regression code (Kelly, 2007) which allows for fixing the slope. JC acknowledges the support of the Deutsche Zentrum für Luft- und Raumfahrt (DLR) via Project ID 50OR1513. MF and DJW acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG) via Projects WI 3871/1-1, and WI 3871/1-2.

Appendix A The effects of sample selection

For the discussion of the TFR at high redshift it is important to be aware not only of the location of the subsample of ‘TFR galaxies’ within a larger parent sample, but also of the effect of the necessary corrections to the observed velocity which ultimately lead to the high-$z$ TFR. Figure 8 illustrates for three stellar mass bins (log$(M_{*}~{}[M_{\odot}])$<10.3; 10.3<log$(M_{*}~{}[M_{\odot}])$<10.8; 10.8<log$(M_{*}~{}[M_{\odot}])$) how the mean maximum rotation velocity changes through corrections for beam-smearing and pressure support, when selecting for rotating disks, and when eventually selecting for ‘TFR galaxies’ following the steps outlined in § 2.4.

The effect of beam-smearing on the rotation velocity is with differences of $\gtrsim 0.1$ dex significant for our galaxies, translating into an offset in stellar mass of $\gtrsim 0.4$ dex. Considering next the impact of turbulent motions, one can clearly see how this is larger for lower-mass (and lower-velocity) galaxies.888Taking turbulent motions into account also has a larger effect at higher redshift due to the increase of intrinsic velocity dispersion with redshift. This is not explicitly shown in Figure 8. This reflects the larger proportion of dispersion-dominated systems at masses of log$(M_{*}~{}[M_{\odot}])\lesssim 10$. Correcting the observed rotation velocity for these two effects does not involve a reduction of the galaxy sample, and the corresponding data points in Figure 8 include all 316 resolved KMOS3D galaxies. The procedure of selecting galaxies suitable for a kinematic disk modelling (W16; § 2.4) has a noticeable effect in the full mass range explored here. It becomes clear that the further, careful selection of galaxies best eligible for a Tully-Fisher study has an appreciable effect on the mean velocity of about $0.02-0.03$ dex, but is minor as compared to the other effects discussed.

While we consider the selection of the ‘TFR sample’ important due to the $v_{\rm rot,max}/\sigma_{0}$ cut and the reliable recovery of the true maximum rotation velocity, we note that it only leads to a small change in TFR parameters as compared to the W16 sample (Figure 9).

Appendix B An alternative method to investigate TFR evolution

It is standard procedure in investigations of the TFR to adopt a local slope for galaxy subsamples in different redshift bins, and to quantify its evolution in terms of zero-point variations, since high$-z$ samples often span too limited a range in mass and velocity to reliably constrain a slope. This method has two shortcomings: first, potential changes in slope with cosmic time are not taken into account. Second, every investigation of TFR evolution is tied to the adopted slope which sometimes complicates comparative studies.

We consider an alternative, non-parametric approach. In Figure 10 we show our TFR galaxies at $z\sim 2.3$ (red) and $z\sim 0.9$ (blue) together with the local sample by Lelli et al. (2016) (black) in the bTFR plane. In the mass bins labeled ‘A’, ‘B’, and ‘C’, we compute the weighted mean velocity of each redshift and mass subsample. We then compare the weighted mean velocities at different redshifts, as indicated in the figure, and determine an average velocity difference from combining the results from individual mass bins.

Although this approach is strongly limited by the number of galaxies per mass bin, and by the common mass range which is spanned by low- as well as high$-z$ galaxies, its advantage becomes clear: not only is the resulting offset in velocity independent of any functional form usually given by a TFR, but the method would also be sensitive to changes of the TFR slope with redshift if the covered mass range would be large enough.

For our TFR samples, we find an average difference in velocity as measured from the average local velocity minus the average high$-z$ velocity, $\Delta\overline{{\rm log}(v_{\rm circ}~{}[{\rm km/s}])}$, of $-0.119$ between $z=0$ and $z\sim 0.9$, and of $-0.083$ between $z=0$ and $z\sim 2.3$. This confirms our result presented in § 4.2, that the bTFR evolution is not a monotonic function of redshift.

Appendix C The impact of mass uncertainties on slope and residuals of the TFR

The slope and scatter of the TFR are affected by the adopted uncertainties in mass. In Figure 11 we show fit examples to the bTFR of the full sample with varying assumptions for the mass uncertainties, namely $0.05\leq\delta{\rm log}(M_{\rm bar}~{}[M_{\odot}])\leq 0.4$. The corresponding changes in slope (from $a=2.11$ to $a=3.74$) are well beyond the already large fit uncertainties on the individual slopes, confirming that a proper assessment of the mass uncertainties is essential. For simple linear regression, the effect of finding progressively flatter slopes for samples with larger uncertainties is known as ‘loss of power’, or ‘attenuation to the null’ (e.g. Carroll et al., 2006). The relevant quantity for our study, however, is the change in zero-point offset, which is for the explored range only $0.02$ dex. This is due to the use of $v_{\rm ref}$ in Equation (2) which ensures only little dependence of the zero-point $b$ on the slope $a$.

Variations of the TFR slope naturally affect the TFR residuals to the best-fit relation (see also Zaritsky et al., 2014). We define the TFR residuals as follows:

[TABLE]

To demonstrate the effect of changing the slope, we show in Figure 12 the bTFR residuals as a function of $R_{e}$. In the upper panel, we show the residuals to a fit with baryonic mass uncertainties of 0.05 dex, leading to a slope which approximately corresponds to the local slope by Lelli et al. (2016). In the lower panel, we show the same for a fit adopting 0.4 dex uncertainties for $M_{\rm bar}$. While there is no correlation found for the former case (Spearman correlation coefficient $\rho=0.02$ with a significance of $\sigma=0.8059$), we find a weak correlation when adopting $\delta M_{\rm bar}=0.4$ dex ($\rho=-0.19$, $\sigma=0.0295$).

We find a similar behaviour for baryonic (and stellar) mass surface density, with no significant correlation between TFR offset and mass surface density for the $\delta M_{\rm bar}=0.05$ dex fit, but a strong correlation for the $\delta M_{\rm bar}=0.4$ dex fit (not shown). No correlation for the $\delta M_{\rm bar}=0.05$ dex fit residuals is found for SFR surface density ($\rho=-0.08$, $\sigma=0.3557$), but a significant correlation with $\rho=-0.37$ and $\sigma=1.1\times 10^{-5}$ for the $\delta M_{\rm bar}=0.4$ dex fit (Figure 13).

From this exercise it becomes clear that the high$-z$ slope, and with it the TFR residuals, are strongly dependent on the accuracy of the mass and SFR measurements.

Appendix D Derivation of the toy model for TFR evolution

D.1 The theoretical framework

In the following, we give details on the theoretical toy model derivation of the TFR and its evolution. The relationship between the DM halo mass, radius, and circular velocity are given by Equations (3), describing a truncated isothermal sphere. A plausible model for a SFG which has formed inside the dark halo is a self-gravitating thin baryonic disk with an exponential surface density profile

[TABLE]

where $\Sigma_{0}$ is the central surface density, related to the baryonic disk mass as $M_{\rm bar}\propto\Sigma_{0}\,R_{d}^{2}$. In reality, disk galaxies feature a finite thickness. This does not affect the scalings presented here (see e.g. Courteau & Rix, 1999; Binney & Tremaine, 2008, and references therein). To associate the baryonic disk to the dark halo, one can assume a simple model where the corresponding masses and radii are related through a proportionality factor:

[TABLE]

$R_{\rm bar}$ can be expressed through the disk scale length $R_{d}$, or the effective radius $R_{e}$, which for rotation-dominated disks are related through $R_{e}\approx R_{d}\cdot 1.68$. As noted in Section 5, we take $r_{f}$ to be independent of redshift. In standard models of disk galaxy evolution, $r_{f}$ combines information on the halo spin parameter, on the halo concentration parameter, and on the ratios of the angular momenta and masses of baryons and DM (cf. Equation (28) of Mo et al. (1998), accounting for adiabatic contraction). It has however been shown that the ratio between $R_{h}$ and $R_{d}$ is approximately constant for massive SFGs in the redshift range $0.8<z<2.6$ (Burkert et al., 2016). This does also hold for our TFR sample and the average values at $z\sim 0.9$ and $z\sim 2.3$, even though there is substantial scatter for individual objects.

To quantify the contributions of baryons and DM to the circular velocity at a given radius we write

[TABLE]

The baryonic contribution can be computed, for instance, using the expression for an infinitely thin exponential disk (Freeman, 1970),

[TABLE]

where $y=r/(2R_{d})$, and $I_{i}(y)$ and $K_{i}(y)$ are the modified Bessel functions of the first and second kind. At $r=R_{e}$, this equation becomes

[TABLE]

where $C^{\prime\prime}$ is a constant. The DM component can be derived simply through a DM fraction at the radius of interest, $f_{\rm DM}(r)=v_{\rm DM}^{2}(r)/v_{\rm circ}^{2}(r)$, or via adopting a full mass profile (e.g. NFW or Einasto, Navarro et al., 1996; Einasto, 1965).

Equations (3) can be combined to

[TABLE]

By inserting Equations (D2) into Equation (D6), and by substituting $R_{d}$ through a re-arranged Equation (D5), one arrives at Equation (4) given in Section 5. After introducing the gas fraction $f_{\rm gas}=M_{\rm gas}/M_{\rm bar}$, one arrives at Equation (5). These equations predict a TFR evolution with a constant slope, but evolving zero-point with cosmic time, depending not only on $H(z)$, but also on changes in $m_{d}$, $f_{\rm DM}(R_{e})$, and $f_{\rm gas}$ with cosmic time.

We note that deviations from the proposed slope ($a=3$) can be related to additional dependencies on $v_{\rm bar}$, e.g. of the surface density $\Sigma$ (Courteau et al., 2007).

D.2 Observational constraints on the redshift evolution of $f_{\rm gas}$, $m_{d}$, and $f_{\rm DM}(R_{e})$

In the following paragraphs, we discuss the motivation for the adopted redshift evolution of $f_{\rm gas}$, $m_{d}$, and $f_{\rm DM}(R_{e})$ in the toy model context. Figure 14 summarizes the individual and combined effects of adopting the respective redshift evolutions of $f_{\rm gas}$, $m_{d}$, and $f_{\rm DM}(R_{e})$ for the bTFR and sTFR evolution.

D.2.1 The redshift evolution of $f_{\rm gas}$

For our toy model approach, we consider the gas fraction $f_{\rm gas}$ to be the sum of molecular and atomic gas mass divided by the total baryonic mass, $f_{\rm gas}=(M_{\rm gas,mol}+M_{\rm gas,at})/(M_{\rm gas,mol}+M_{\rm gas,at}+M_{*})$. The evolution of the molecular gas mass-to-stellar mass ratio is given through the scaling relation by Tacconi et al. (2017):

[TABLE]

Here, we do not take into account the additional dependencies given in the full parametrization by Tacconi et al. (2017) on MS offset, and offset from the M-R relation, but assume that the model galaxies lie on these relations.

Locally, the galactic gas mass is dominated by atomic gas. To account for atomic gas mass at $z=0$, we use the fitting functions presented by Saintonge et al. (2011). We use a local reference stellar mass of log$(M_{*}~{}[M_{\odot}])=10.94$, i.e. the stellar mass corresponding to our reference velocity $v_{\rm ref}=242$ km/s in the context of the sTFR fit by Reyes et al. (2011).

To account for atomic gas masses at $z>0$, we follow the theoretical prediction that, at fixed $M_{*}$, the ratio of atomic gas mass to stellar mass does not change significantly with redshift (e.g. Fu et al., 2012). We use again the fitting functions by Saintonge et al. (2011) to now determine the atomic gas mass for galaxies with log$(M_{*}~{}[M_{\odot}])=10.50$, which corresponds to the average stellar mass of our TFR galaxies at $v_{\rm ref}=242$ km/s in both redshift bins.

Between $z=0$ and $z=0.9$, we assume a smooth TFR evolution, meaning that at fixed circular velocity, galaxies have decreasing $M_{*}$ with increasing redshift, in order to compute the gas fractions. Although we cannot quantify this assumption with our observations, we note that in comparing to our data, only the relative offset in $f_{\rm gas}$ (or any other parameter discussed below) between $z=0$, $z=0.9$, and $z=2.3$ is relevant. Our assumption therefore serves mainly to avoid sudden (unphysical) offsets in the redshift evolution of $f_{\rm gas}$.

Corresponding values of the gas mass fraction at $z=\{0.0;0.9;2.3\}$ are $f_{\rm gas}\approx\{0.07;0.36;0.58\}$.

D.2.2 The redshift evolution of $m_{d}$

The baryonic disk mass fraction, $m_{d}=M_{\rm bar}/M_{h}$, is not a direct observable, since it depends on the usually unknown DM halo mass. For the local Universe, we use the fitting function by Moster et al. (2013) from abundance matching to determine a stellar disk mass fraction, $m_{d,*}=M_{*}/M_{h}$. For a stellar mass of log$(M_{*}~{}[M_{\odot}])=10.94$, this gives $m_{d,*}\approx 0.012$. Again, we use the fitting functions by Saintonge et al. (2011) to determine the corresponding gas mass, taking into account contributions from helium via $M_{\rm He}\approx 0.33~{}M_{\rm H{\textsc{i}}}$. This results in a baryonic disk mass fraction at $z=0$ of $m_{d}\approx 0.013$.

The recent study by Burkert et al. (2016) finds a typical value of $m_{d}=0.05$ for SFGs at $0.8<z<2.6$ based on a Monte-Carlo NFW modelling of data from the KMOS3D and SINS/zC-SINF (Förster Schreiber et al., 2009; Mancini et al., 2011) surveys. These galaxies have masses similar to the galaxies in our TFR sample. We adopt their value of $m_{d}=0.05$ for $0.8<z<2.6$.

Between $z=0$ and $z=0.8$ we assume a linear increase of $m_{d}$. Clearly, this is a simplifying conjecture. As for the atomic gas masses, we emphasize that this assumption has primarily cosmetic effects, while the crucial quantity is the relative difference in $m_{d}$ between $z=0$, $z\sim 0.9$, and $z\sim 2.3$.

D.2.3 The redshift evolution of $f_{\rm DM}(R_{e})$

For the DM fraction of local disk galaxies, we follow Figure 1 by Courteau & Dutton (2015) which, among others, shows galaxies from the DiskMass survey (Martinsson et al., 2013b, a). At $v_{\rm circ}=242$ km/s, DM fractions of local disk galaxies lie roughly between $f_{\rm DM}(r_{2.2})=0.55$ and $f_{\rm DM}(r_{2.2})=0.75$, with large scatter and uncertainties.

At higher redshift, W16 derived DM fractions from the difference between dynamical and baryonic masses of the KMOS3D subsample of 240 SFGs, which represents our parent sample. Corresponding values, also corrected for mass completeness, are given in their Table 1.

For convenience, we parametrize the evolution of the DM fraction with redshift as follows: $f_{\rm DM}(R_{e})=0.7\cdot{\rm exp}[-(0.5\cdot z)^{2.5}]$. This gives an evolution which is somewhat stronger than what is suggested by just taking the average values provided by Courteau & Dutton (2015) and W16, but easily within the uncertainties presented in both papers. We adopt this marginally stronger evolution to better match our observed TFR offsets with the toy model.

Corresponding values of the DM fraction at $z=\{0.0;0.9;2.3\}$ are $f_{\rm DM}(R_{e})\approx\{0.70;0.61;0.17\}$.

We note that our toy model evolution is particularly sensitive to the parametrization of $f_{\rm DM}(R_{e},z)$ which is in our implementation with the simplistic description for $m_{d}(z)$ responsible for the flattening/upturn of the sTFR/bTFR (see Figure 14). The high value for the local DM fraction (which would at $r=R_{e}$ rather be lower than at $r=r_{2.2}$) as well as the comparably strong evolution at $z>1$ can certainly be challenged.

D.2.4 Comments on the evolution of the halo concentration parameter

The predicted evolution of the halo concentration parameter $c$ between $z=2$ and $z=0$ for haloes of masses that are relevant to this study (i.e. central stellar masses of log$(M_{*}~{}[M_{\odot}])\approx 10.5$ at $z\sim 2$, and log$(M_{*}~{}[M_{\odot}])\approx 10.9$ at $z\sim 0$) goes from $c\approx 4$ at $z=2$ to $c\approx 6$ at $z=1$ and to $c\approx 7$ at $z=0$ (Dutton & Macciò, 2014). This alone would increase the DM fraction at $R_{e}$ by roughly $0.1$.

Starting from the central DM fractions as determined by W16, abundance-matched haloes (Moster et al., 2013) would require concentrations of $c\approx 3$ and $c\approx 12$ at $z\sim 2.3$ and $z\sim 0.9$, respectively (cf. Eq. 19 by Mo et al., 1998). Extending this to $z=0$ is not straight-forward since local late-type galaxies have typically lower circular velocity as required for the extrapolation of the local TFR to our $v_{\rm ref}=242$ km/s (see discussion in § 4.2). However, using the stellar mass-radius relation presented by van der Wel et al. (2014), inferred concentrations of these hypothetical haloes would have to be $c\approx 13$.

This points towards a potential issue in the observational constraints to our toy model because the $m_{d}$ values inferred by Burkert et al. (2016) are based on Monte-Carlo modelling involving standard NFW haloes. One could consider fitting $m_{d}$ to better match the observed TFR zero-point evolution.

In general, the possible effects of adiabatic contraction or expansion of the halo as a response to baryonic disk formation make theoretical predictions of the central DM fractions uncertain (see e.g. the discussions by Duffy et al., 2010; Velliscig et al., 2014; Dutton et al., 2016; and also Dutton & Macciò, 2014 for an overview of predictions of concentration-mass relations from analytical models).

Appendix E Physical properties of galaxies in the TFR sample

In Table LABEL:tab:ind-properties we list redshift $z$, stellar mass $M_{*}$, baryonic mass $M_{\rm bar}$, maximum modelled circular velocity $v_{\rm circ,max}$, and modelled intrinsic velocity dispersion $\sigma_{0}$ of our TFR galaxies. The full table is available in machine readable form.