THE SINS SURVEY: MODELING THE DYNAMICS OF z ∼ 2 GALAXIES AND THE HIGH-z TULLY–FISHER RELATION^*

Förster Schreiber et al. (2009) describe in detail the selection criteria of the SINS z ∼ 2 sample. Briefly, the SINS galaxies were selected to be massive star-forming galaxies using a variety of complementary techniques in ultraviolet (UV) through IR wavebands (BX/BM criterion, Adelberger et al. 2004; s-BzK, Daddi et al. 2004a; K magnitude, Cimatti et al. 2002; IRAC flux, Kurk et al. 2008). These selection criteria sample fairly luminous (L_Bol ∼ 10¹¹ –10¹² L_☉) and massive galaxies (M ∼ 10⁹–10^11.3 M_☉) with SFRs of 10–300 M_☉ yr⁻¹. These represent the bulk of the cosmic star formation activity and stellar mass density at these redshifts (Reddy et al. 2005; Rudnick et al. 2006; Grazian et al. 2007; Caputi et al. 2007). Here we present the dynamical modeling performed on a subsample of 18 galaxies (see Table 1, Figure 2 below) that are selected for showing ordered rotational signatures, i.e., a regular velocity pattern with a clear velocity gradient and a velocity dispersion that peaks in the center, as well as for their high signal-to-noise ratio (S/N) observations.

In order to distinguish between rotationally dominated galaxies and sources with merger-induced complex dynamics, we used, where possible, the method described in Shapiro et al. (2008). Using kinemetry (Krajnović et al. 2006), we were able to quantify asymmetries in both the velocity and velocity dispersion maps of the warm gas, in order to empirically differentiate between major merging and nonmerging (or minor merging) systems at high redshift. These criteria take full advantage of the wealth of information available with integral field data, which—unlike broadband morphology—give direct information on the dynamical state of the galaxies. The method was first applied to a subsample of galaxies with the best quality data (see Shapiro et al. 2008), and it has been now extended to a larger subset of SINS galaxies (K. Shapiro et al. 2009, in preparation). In total, 11 out of 18 galaxies of the sample discussed here were classified as rotationally dominated by kinemetry. For the remaining seven sources (Q1623-BX663, Q2346-BX416, D3a-7144, D3a-4751, ZC1101592, GK2471, GK1084), we were not able to derive a quantitative classification due to the lower S/N of the data. We include them in this study since they closely resemble at a visual inspection the ones classified as disks in their dynamical properties (regular velocity pattern with clear velocity gradient, velocity dispersion peaking in the center), with no evidence of more complex merger-induced dynamics.

The disklike systems selected for this study include seven BX/BM galaxies selected by their rest-frame UV colors down to R < 25.5, taken from the Hα survey of Erb et al. (2006c), two galaxies selected with IRAC/Spitzer m_AB(4.5 μm) < 23 in the framework of the GMASS survey (Kurk et al. 2008), and 10 rest-frame optical-selected galaxies meeting the s-BzK color criterion introduced by Daddi et al. (2004a). These are taken from the K20 survey (Daddi et al. 2004b), Deep-3a (Kong et al. 2006), and z-Cosmos (Lilly et al. 2007) surveys.

We observed the selected galaxies with SINFONI (Eisenhauer et al. 2003) at the VLT UT4 telescope during several campaigns between July 2004 and November 2007 as part of guaranteed time observations. The H- or K-band grating was used to sample the Hα emission, depending on the redshift of the galaxy. Most of the observations were carried out in seeing-limited mode, using the 0250 × 0125 pixel scale, which provides a total field of view of 8'' × 8''. The average resolution obtained for our 18 galaxies, as sampled by point-spread function (PSF) reference star observations before and after every hour of target integration, is ∼05. The spectral resolution provided is about 80 km s⁻¹ full width at half-maximum (FWHM) in the K band and 100 km s⁻¹ in the H band. For one source, D3a-15504, the presence of a suitable nearby reference star allowed us to use the correction provided by Natural Guide Star (NGS) Adaptive Optics (AO) observations. For two more objects, ZC782941 and Q2346-BX482, the Laser Guide Star (LGS) system PARSEC (Rabien et al. 2004; Bonaccini et al. 2006) was used to provide an average resolution of 02 (corresponding to only 1.25 kpc at z = 2). In the AO observations of these three sources the 0100 × 0050 pixel scale was used, providing a field of view of 32 × 32. The observations and the data reduction were performed as described in Förster Schreiber et al. (2006), Abuter et al. (2006), and Davies (2007).

The final result of the data-reduction procedures is a flux-calibrated datacube, containing the full image of the galaxy observed at each wavelength. We are thus able to extract from the cube both the one-dimensional spectrum at each pixel (or integrated over a larger region) and monochromatic images of the field of view at different wavelengths.

To increase the S/N, we median filtered the reduced cubes spatially using an FWHM = 3 pixels, obtaining an effective resolution of ∼06, or ∼02 for the galaxies observed with the 005 pixel scale and AO. We then extracted from the smoothed cubes the continuum maps, line emission, velocity and velocity dispersion fields of Hα. This was done by fitting a function to the continuum-subtracted spectral profile at each spatial position in the datacube, masking out the wavelength range contaminated with the [N ii] line emission. The function fitted was a convolution of a Gaussian with a spectrally unresolved emission line profile of a suitable sky line, i.e., a high signal and nonblended line in the appropriate band. A minimization was performed in which the parameters of the Gaussian were adjusted until the convolved profile best matched the data. During the minimization, pixels in the data that consistently deviated more than 3σ from the average were rejected, and were not used in the analysis. The uncertainties on the measured velocity and velocity dispersion maps are evaluated through Monte Carlo realizations, perturbing the input data assuming Gaussian uncertainties. All details and issues of the kinematic extraction with the ${\sf IDL}$ routine ${\sf LINEFIT}$ are described in more details in Förster Schreiber et al. (2009) and R. Davies et al. (2009, in preparation).

We create the disk models using the ${\sf IDL}$ code ${\sf DYSMAL}$. The code uses a set of input parameters to derive a datacube with two spatial and one spectral (velocity) axis, from which it is possible to extract two-dimensional morphological and kinematical maps. It begins by creating a face-on model of an axisymmetric galaxy with three spatial dimensions [X₀, Y₀, Z₀], in which Z₀ is normal to the galaxy plane. The radial mass profile is specified by the user and can be based on exponential, Gaussian, Moffat, or power-law functions. The radial luminosity profile can be specified independently, since it can refer to different tracers that do not necessarily follow the mass distribution. The scale height of the disk must also be provided. The three-dimensional model is then rotated to the required inclination and position angles, to create a cube that has axes [X_s, Y_s, Z_s] where [X_s, Y_s] is the projection on the sky. Although a warped disk can be modeled by specifying these two parameters as functions of radius, in all cases here the disk was modeled in a single plane. To make this coordinate transformation more efficient, it is first necessary to identify those pixels in the resulting [X_s, Y_s, Z_s] cube containing flux, and only fluxes for these pixels are calculated. The flux in each pixel is interpolated based on the subpixel location in the original [X₀, Y₀, Z₀] cube. A second matching cube is created in which each pixel is assigned an appropriate line-of-sight velocity. The velocities are interpolated from a rotation curve that is itself generated from the combined mass profile (including a black hole mass if one is specified) and scaled to the total mass at a given radius. The next step is, for each pixel [X_s, Y_s, Z_s], to generate a Gaussian profile to model the line emission, whose scaling depends on the flux, whose center depends on the line-of-sight velocity, and whose width depends on both the instrumental broadening specified by the user and a local isotropic velocity dispersion depending on the z-scale height. This z-velocity dispersion σ₀₁ is estimated in one of the two ways:

• 1.  
  using an approximation appropriate for large thin disks, in which it depends on the scale height and mass surface density (Binney & Tremaine 2008, Equation (4.302c)):

  $$\begin{equation} \sigma _{01} \sim \sqrt{\frac{v^2(R)\ h_z}{R}}, \end{equation} \tag{ 1 }$$

  where h_z is the scale height and v(R) is the rotational velocity at radius R;
• 2.  
  assuming that the vertical motion is restricted to the inner part of an extended mass distribution and does not feel the entire gravity:

  $$\begin{equation} \sigma _{01}\sim \frac{v(R)\ h_z}{R}, \end{equation} \tag{ 2 }$$

  which is more suitable for a thick or compact disk (Genzel et al. 2008); this was used in the modeling discussed below.

The total z-velocity dispersion in the model is then given by $\sigma _0=\sqrt{\vphantom{a^a}\smash{\hbox{$\sigma _{01}^2+\sigma _{02}^2$}}}$, where σ₀₂ is an additional component of isotropic velocity dispersion throughout the disk that can be added by the user. This term was introduced in ${\sf DYSMAL}$ to account for instrumental broadening. It can also be used to include any additional possible contribution to the observed velocity widths from, e.g., rapid mass accretion, strong feedback from star formation, or unresolved noncircular motions (hereafter globally called "random motions" for simplicity). This yields a line-of-sight velocity profile for each pixel in [X_s, Y_s, Z_s]. Then, at each [X_s, Y_s] location, the velocity profiles for all the [Z_s] positions are summed to yield the full line-of-sight velocity distribution at that point. This creates a cube [X_s, Y_s, V] with two spatial and one velocity axes. The final step is to convolve each [X_s, Y_s] plane in this cube with the spatial beam and to rebin to the desired pixel scale.

The reason for beginning with a galaxy in three spatial dimensions is to ensure that there is always enough spatial resolution along the projected minor axis to yield a reliable line-of-sight velocity distribution. This can be an issue when the inclination is high, and in such cases a direct two-dimensional projection is insufficient. It also enables one to apply a mask that modulates the flux distribution, for example, to simulate a knotty or clumpy structure. The code allows one either to generate a mask containing random knots; or to input one if, for example, a more ordered pattern such as spiral arms is required. In either case, the mask is applied to the galaxy while it is face-on, in the [X₀, Y₀] plane. For this high-z study, however, a plain exponential disk was used to keep the number of free parameters as low as possible.

The velocity dispersion across the galaxy disk is an important parameter to consider. When there are significant random motions, as implied by the large observed dispersion, they should also support some of the mass. This is currently not implemented in DYSMAL, which assumes that all the mass is supported by ordered rotation in the disk plane. Thus, if $v/\sigma _0 \lesssim \sqrt{3}$, the total system mass may be underestimated by a factor of 2 or more. The average V/σ₀ for our sample is ∼4.4, producing an effect of ∼12%.

We account for spatial beam smearing from the PSF and velocity broadening from instrumental spectral resolution in our modeling, as we compare this spatially and spectrally convolved model disk to the observations. Hence, the best-fit quantities represent the intrinsic properties derived by properly taking into account these observational effects.

Given the large number of parameters and the potential for many local minima in the likelihood/χ² space, it is important to use a reliable fitting algorithm, i.e., robust and efficient one to find the global extremum. We selected for this task genetic algorithms, which are heuristic search techniques that provide a means to efficiently sample the parameter space, reproducing in a computational setting the biological effect of evolution by natural selection (e.g., Holland 1975). These kinds of algorithms have the advantage of being robust and efficient in finding the global maximum of a given function even in a very complex parameter space, with a reasonable number of iterations and without the need of a good set of initial guesses. It thus represents an ideal tool for our specific problem.

We use a modified version of PIKAIA, a genetic algorithm-based optimization routine by Charbonneau (1995). The algorithm starts with a random population of "individuals" that represent a single location in parameter space with a complete set of the parameters we want to fit. The fitness of each individual is evaluated through their χ², and the probability to survive in the next iteration is set to be proportional to this value. In this way the more suitable parameter sets are selected, but at the same time the parameter space is efficiently sampled by "breeding" between different individuals and allowing random "mutations" and "crossovers" between the parameters. Therefore, in each iteration the whole population becomes more and more suitable, while still sampling the parameter space. When the requested tolerance is reached, or at the end of a given number of iterations, the best-fitting "individual" is chosen as the solution.

The genetic algorithm allows us to efficiently minimize the differences between the exponential disk model and the observed velocity and velocity dispersion fields of the Hα emission line. Simulations on test galaxies have shown that leaving the rotation center as a free parameters in the fitting significantly worsens the final results, even after increasing the number of maximum iterations allowed (see also Shapiro et al. 2008). We therefore locate the center using the continuum light distribution of the galaxies, which corresponds roughly to rest-frame R band as observed in the NIR at z ∼ 2. In the SINFONI data of our galaxies, the continuum is detected, although typically at a lower S/N than the emission lines. The resulting intensity map is generally smoother and distinct than the emission line (Förster Schreiber et al. 2006; Genzel et al. 2008; Shapiro et al. 2008), as it is also the case for local disk galaxies (e.g., Daigle et al. 2006). These continuum intensity maps show a clear peak in the central region, especially in the BzK-selected galaxies, and represent a better tracer of the mass distribution (dominated by stars) in the galaxies than the more clumpy line emission (tracing active sites of star formation). The center is identified by computing the intensity weighted mean of the light profile of the central pixels along both minor and major axes, and it is kept fixed in the disk fitting procedure.

The other parameter that is provided as an input to the fitting routine is the scale length R_d of the disk. This parameter is in fact better constrained by the emission intensity than by the dynamical properties of the gas. In this case we used the higher S/N line emission image of the galaxies to quantify it. The disk Hα distribution has been found to follow the continuum for local spiral galaxies, with consistent derived scale lengths (e.g., Tacconi & Young 1986; Koopmann et al. 2001; Hanish et al. 2006). From an analysis of different SINS and simulated galaxies, we found that the most robust way to measure the scale length of the disks given the S/N and resolution of the data is to fit a linear Gaussian profile along the major axis of the line intensity maps. This method has the advantage of using just the bulk of the brighter line emission, avoiding the problems related to the surface brightness at the edges of the galaxies. Following Bouché et al. (2007), we assumed that the derived HWHM, once corrected for the observed seeing, corresponds to the scale length R_d of the disk. This assumption is supported by model disk fits, incorporating the beam smearing due to the ∼05 FWHM resolution of our data. These simulations show that R_d measured in this way is likely to be overestimated by no more than ≲15%. Using available Hubble Space Telescope (HST)/NICMOS continuum imaging (see Förster Schreiber et al. 2009) and deep ISAAC imaging in GOODS (data release ver. 1.5: http://www.eso.org/science/goods/releases/20050930/) of a subsample of these sources (one GMASS, two K20, and five BX), we find reasonable agreement between our measured scale length and the ones derived by fitting an exponential profile to these K-band images with GALFIT (Peng et al. 2002). The results are shown in Figure 1.

Zoom In Zoom Out Reset image size

According to simulations, a choice of a scale length R_d 15% smaller would provide a dynamical mass ∼10% lower, and a ∼30% smaller one a mass ∼20% lower. We note that the value obtained for the maximum rotational velocity V_max is instead much less sensitive to a variation of the scale length, as the fitting routine has to reproduce the observed velocity pattern. However, we also note that if a flatter mass distribution is a more appropriate description of the structure of these high-z disks, e.g., the ringlike structures discussed in Genzel et al. (2008), the true enclosed masses would then be up to ∼10%–30% lower than the one obtained with the method presented here.

The inclination of the disk with respect to the line of sight i, the position angle on the plane of the sky θ, the total enclosed mass M_dyn, the velocity zero point v₀ at the rotation center, and the constant global isotropic dispersion term σ₀₂ (see Equation (2)) are left as free parameters in the minimization. The model disk obtained is then rescaled as it would appear at the redshift of the source, rebinned to the pixel size of the current observation, and smeared using the FWHM of the PSF observed (see Section 2). In every iteration of the genetic algorithm, the χ² is evaluated for each "individual" from the squared difference between the observed velocity and velocity dispersion maps and the modeled ones. The difference in each pixel is weighted using the uncertainties on the velocity and on σ derived as described in Section 2 at that location, e.g.,

$$\begin{equation} \chi ^2({\rm vel})=\sum _{i} \frac{(V({\rm obs})_i-V({\rm mod})_i)^2}{\Delta V({\rm obs})_i^2} \times \frac{1}{N_{\rm pix}-N_{\rm param}}, \end{equation} \tag{ 3 }$$

where N_pix is the number of pixels with S/N over the threshold (see Section 2.3) used in the analysis. The best-fitting parameter set is therefore taken as the solution. The uncertainties are evaluated following the Avni (1976) method to calculate confidence limits for a multiple number of parameters, studying the variation of χ² around the optimized minimum. We assumed as 90% confidence interval where the reduced χ² is increased by Δχ² = 6.25 for the three "interesting parameters" (inclination, dynamical mass, and σ₀₂) used in the analysis.

Thanks to the multiband coverage of the different fields from which we drew our SINS targets (see Section 2.1), we were able to perform an SED analysis using the available broadband photometry from the observed optical to NIR/MIR. In particular, as described in detail in Förster Schreiber et al. (2009), we used data from published works by Erb et al. (2006b, BX), Daddi et al. (2004b, K20), Kong et al. (2006, D3a), Kurk et al. (2008, GMASS), Capak et al. (2007, COSMOS), as well as photometry kindly provided to us by K. G. McCracken et al. (2009, in preparation, COSMOS). SED modeling has been carried out by these authors but using different model ingredients or assumptions. In order to ensure more uniformly derived properties among our sample, we remodeled all the galaxies following the prescriptions described in Förster Schreiber et al. (2004; see also Förster Schreiber et al. 2009). Our results do not differ significantly from the published ones. We used the solar metallicity Bruzual & Charlot (2003) stellar population synthesis models, with a Chabrier initial mass function (IMF) between 0.1 and 100 M_☉, and the Calzetti et al. (1994) reddening law. We considered three combinations of star formation history and dust content: constant star formation rate and dust (CSF), an exponentially declining SFR with e-folding timescale of t = 300 Myr and dust (τ300), and a dust-free single stellar population formed instantaneously. We adopted the best of those three cases based on the reduced χ² value of the fits (for all the galaxies in the sample, this is either CSF or τ300). The redshift is kept fixed to the spectroscopic redshift measured from the emission lines detected, while the stellar mass M_*, SFR, star formation history, stellar population age, and visual extinction A_V are free parameters. We evaluate the uncertainties through Monte Carlo simulations, where the input photometry is varied randomly assuming the photometric uncertainties are Gaussian. The fitting results for the galaxies of the sample are reported in Table 1.

The reported stellar masses are affected by various uncertainties and systematic errors. In particular, different stellar population synthesis models do not reproduce a consistent picture of evolution in the rest-frame NIR (see Marchesini et al. 2008). For example, Maraston et al. (2006) suggest that the Bruzual & Charlot (2003) stellar population synthesis models used here may overestimate the stellar masses up to a factor of ∼2 when the bulk of the population in the age range between ∼0.5 and ∼2 Gyr. This is due to a different approach in the treatment of the thermally pulsing asymptotic giant branch (TP-AGB) phase. However, findings from much larger samples, where the effects of SED modeling can be better and more reliably investigated on a statistical basis, suggest that the correction between the two different approaches can be negligible or up to ∼30% (see, e.g., Berta et al. 2007; Marchesini et al. 2008; Wuyts et al. 2007). Accordingly, using the different library of synthetic spectra of S. Charlot & G. Bruzual (2009, in preparation), the stellar masses estimated for our small sample change individually by no more than 20%, while the median of the whole sample remains consistent within 2%. In any case, this correction would make the stellar masses lower, and the offset in the zero point of the high-z TFR discussed in Section 5 will be correspondingly larger.

In four objects of the sample, there are independent indications of the presence of an active galactic nucleus (AGN): dynamical and UV spectral evidence suggests the presence of an AGN in D3a-15504 (see Genzel et al. 2006; E. Daddi et al. 2009, in preparation), BX663 and D3a-7144 show AGN features in optical spectroscopy (Erb et al. 2006c; E. Daddi et al. 2009, in preparation, respectively), while for D3a-6004 a 24 μm excess was observed by E. Daddi et al. (2009, in preparation), and a large [N ii]/Hα ∼ 0.63 ratio was measured in our data in the nuclear region of the galaxy (P. Buschkamp et al. 2009, in preparation). Inspection of the optical to K-band available SED for these galaxies suggests that this AGN contribution may be biasing the stellar mass estimate toward higher stellar mass in the latter two cases (see Förster Schreiber et al. 2009). We therefore conservatively consider the stellar mass estimate for D3a-6004 and D3a-7144 as an upper limit.

The stellar masses derived for the sample are between (1–30) × 10¹⁰ M_☉, with an average of 6.3 × 10¹⁰ M_☉ confirming that all the selected galaxies are drawn from the massive star-forming galaxy population at high-z. Thanks to the SED results and our dynamical model, we are able to compare the stellar and dynamical mass of the galaxies. In Figure 3, these are compared for galaxies in which multiband photometry is available to constrain the SED fitting. We observe remarkable correlation between the two quantities: the Spearman correlation test gives a probability P = 2 × 10⁻⁵ that the two masses are uncorrelated.

Zoom In Zoom Out Reset image size

The resulting average dynamical mass at R = 10 kpc is ∼3.3 times larger than the average stellar mass, excluding the sources with a possible AGN contribution. Therefore, given a dark matter contribution of ∼40% at this radius to the total dynamical mass of such disks (see Genzel et al. 2008), we infer a large fraction of gas in the sample galaxies (M_gas ∼ M_*, or ∼30%), although with very large uncertainties. This estimate may be compared with a gas mass derived from the SFR surface density using the Schmidt–Kennicutt relation (see Bouché et al. 2007). In this way, we derive an average gas contribution of 23% to the total dynamical mass of the galaxies, where the difference can be accounted for by the large uncertainties in both our estimate and the extinction correction needed to derive the SFRs from the detected Hα fluxes. The gas fraction estimate is also in agreement with the findings of Daddi et al. (2008), who directly observed molecular gas in two z ∼ 1.5 BzK galaxies, deriving gas fraction as high as M_gas ∼ M_*. An alternative possibility is that the stellar mass is significantly underestimated for most of the sample, presumably because of a faded old stellar population that is not detected in the multicolor SED, but would still contribute to the total dynamical and stellar mass (see the discussion by, e.g., Shapley et al. 2001; Papovich et al. 2001; and for the galaxies with available IRAC photometry Wuyts et al. 2007), or that the dynamical mass is overestimated due to the systematic effects discussed in Section 3.2.

Although the stellar population age is not strongly constrained by the SED fitting due to the uncertainty on the star formation history and the degeneracy with extinction, we note that the objects with the higher M_dyn/M_* generally correspond to younger best-fitting ages, similar to what was discussed by Erb et al. (2006b), and that none of the youngest galaxies has a small mass ratio (see Figure 4). The significance of the correlation according to the Spearman test gives a probability P = 0.08 that the two quantities are uncorrelated. The scatter is in fact high, probably suggesting that these galaxies are not closed systems, and that an important contribution from continuous gas accretion must be considered. This again suggests that the galaxies in the sample are star-forming disks with significant gas fraction and gas accretion, where the gas consumption correlates with the stellar population age.

Zoom In Zoom Out Reset image size

With the available modeled kinematics of z ∼ 2 galaxies from our SINS data, we are now able to push, for the first time with a sizeable sample, the investigation of the evolution of the TFR to even higher redshift. When comparing the evolution of the high-z TFR with that inferred from local (or lower redshift) samples of star-forming galaxies, one must be aware that we are comparing different classes of objects, which are not necessarily linked from an evolutionary point of view. Therefore, any evolution in the relation cannot be interpreted straightforwardly as evolution of any individual galaxy, as from z = 2 to z = 0 part of our galaxies may undergo significant (morphological/dynamical) evolution driven by efficient secular/internal dynamical processes and/or merging (e.g., Genzel et al. 2008; Genel et al. 2008). Rather, evolution of TFR represents the evolution of ensemble properties of disklike populations at different redshifts.

The available full three-dimensional coverage of the dynamical features will not only provide a more complete and detailed view of the kinematics, but will also allow us to avoid the biases and limitations due to slit spectroscopy, such as misalignment of the slit and the major axis, and to identify disturbed rotational patterns avoiding the contamination of our disklike sample with more complex dynamics. Using similar three-dimensional SINFONI data, van Starkenburg et al. (2008) have presented the position of a single galaxy at z = 2.03 on the stellar mass TFR, yielding an offset of (∼2σ) from the local one (see Figure 5). In view of the uncertainties of the quantities involved, which are not insignificant at high z, larger samples are required to constrain more robustly any offset in the TFR at z ∼ 2.

Zoom In Zoom Out Reset image size

We have therefore used the best-fitting V_max from our dynamical modeling, together with the stellar mass from SED fitting discussed in Section 4.1, to build the TFR for our sample, shown in Figure 5. The solid line shows the z = 0 best-fit relation obtained by Bell & de Jong (2001). Both are corrected to convert their "diet" Salpeter and Salpeter IMF to the Chabrier IMF used here, using a factor of 1.19 and 1.7, respectively. We note that Pizagno et al. (2005) have obtained a shallower (M_* ∝ V^3.05) local stellar mass TFR than Bell & de Jong (2001, M_* ∝ V^4.5), but the two relations agree at V = 200 km s⁻¹. However, they used the circular velocity at 2.2R_d instead of V_max to estimate of the rotational velocity, and differences in the sample selection may also account for the different result. McGaugh (2005) obtained an intermediate slope, M_bar ∝ V^4.0, but they explored the total baryonic TFR, including the contribution of the gaseous mass component. In this case, higher gas fraction in lower mass galaxies may explain the flatter slope. Recently, Meyer et al. (2008) presented the local stellar mass TFR for a large H i-selected sample of galaxies, obtaining a slope consistent to the Bell & de Jong (2001) value adopted here for comparison with previous studies.

The circles in Figure 5 show the data for the z ∼ 2.2 galaxies from SINS, while the filled triangles show the z ∼ 1.5 sources. The open circle is the z ∼ 2 galaxy presented in van Starkenburg et al. (2008). The average uncertainty is shown as an error bar in the lower right corner. This does not include possible systematic effects both in the SED fitting (see Section 4.1) and in the dynamical modeling (see Section 3). The z ∼ 2.2 galaxies show an evident evolution in the zero point of the TFR. The dashed line shows the best-fitting relation keeping the slope fixed at the z = 0 value found by Bell & de Jong (2001). The current data are in fact not yet able to fully constrain the slope of the relation due to the limited statistics and mass range. The best-fitting zero point is ZP_2.2 = −0.09 ± 0.11, where the error is the 1σ uncertainty, to be compared with the local ZP₀ = 0.32, giving a significance of ∼3.7σ. The z = 2.2 relation can be written as log(M_*) = −0.09 + 4.5 × log(V_max), finding an offset in log(M_*) of 0.41 for a given rotational velocity. We do not find evidence for significant evolution in our z ∼ 1.5 sample with respect to the local relation, although the data are still too sparse in this redshift range to place a firm constraint. The evolution observed at z ∼ 2.2 is larger than the intermediate 0.36 dex found at z ∼ 0.6 by Puech et al. (2008), although they use a different local reference, suggesting an increasing of the evolution of the TFR with redshift.

Interestingly, the amount of observed evolution in our sample is consistent with the factor ∼1.25 of increased velocity for a given stellar mass predicted by Somerville et al. (2008) in their semianalytic formation model of a disk with M_* = 6 × 10¹⁰ M_☉ at this redshift. They follow the classical paradigm for disk formation within massive extended dark matter halos initially proposed by Fall & Efstathiou (1980) but use a Navarro et al. (1997) mass density profile instead of a singular isothermal sphere as halo mass density profile, account for the expected evolution with time of the profile concentration c_vir (Bullock et al. 2001), and for the effect of self-gravity of the baryons in the central part of the halo (so-called "adiabatic contraction"). They find that with this refined model the disk rotation velocity V_max at fixed disk mass evolves much more gradually than the halo virial velocity at fixed halo mass, remaining nearly constant and showing the largest evolution between z ∼ 1.5–3. If this prediction of the model is correct, our sample therefore probes the best suited cosmic epoch to test the evolution of the TFR.

We also compare the observed TFR with the smoothed particle hydrodynamic (SPH) simulations of disk formation and evolution in a CDM hierarchical clustering scenario by J. Sommer-Larsen et al. (2009, in preparation) and Portinari & Sommer-Larsen (2007). They include in their modeling star formation and stellar feedback effects to delay the collapse of protogalactic gas clouds and allow the infalling gas to preserve a large fraction of its angular momentum. Previous numerical simulations have in fact shown that when only cooling processes are included, the infalling gas loses too much angular momentum, resulting in disks that are much smaller than the observed ones (Navarro & White 1994; Navarro & Steinmetz 2000). Moreover, chemical evolution with noninstantaneous recycling, metallicity-dependent radiative cooling and the effects of a metagalactic UV field, including simplified radiative transfer, were implemented in the simulations. This is an improved version of the models developed by Sommer-Larsen et al. (1999, 2003) to address the angular momentum problem in disk formation models and provide realistic z = 0 galaxies. It has proven successful in reproducing better than previous numerical models the observed sizes of disk galaxies, the luminosity profile of both disks and bulges, bulge to disk ratios, integrated colors and TFR of local galaxies, as well as the observed peak in the cosmic SFR at z ∼ 2 (Sommer-Larsen et al. 2003). The models are characterized by a formation mechanism where a large fraction of the gas is accreted onto the forming disk through both rapid cold gas accretion in filamentary structures (dominant at z ≳ 2) and gradual cooling flows from a surrounding hot phase (dominant for z ≲ 2), without abrupt changes of angular momentum orientation, rather than by rapid mergers of massive cold clumps. These simulated galaxies at z ∼ 2.2 are shown as open diamonds in Figure 5. Again, we obtain very good agreement between our data and the predicted position of z ∼ 2.2 galaxies on the TFR for a range of stellar masses. In the model the predicted zero-point offset is produced by a relevant inside-out mass accretion between z ∼ 2.2 and z ∼ 1, that do not correspond to an equal increase of the rotation velocity of the galaxies. This issue will be investigated in a forthcoming paper (J. Sommer-Larsen et al. 2009, in preparation). However, if the TFR evolution is driven mainly by dark matter as suggested by the models, this may suggest that the galaxies we are observing at high-z are "stable," in contradiction with inferences from our detailed studies of several of our disks which appear to be globally unstable to star formation (see Genzel et al. 2008).

Kassin et al. (2007) have constructed a stellar mass TFR at z ≲ 1 using the $S_{0.5}=\sqrt{\vphantom{a^a}\smash{\hbox{$0.5\cdot V^2_{\rm rot}+\sigma _{0}^2$}}}$ estimator, in order to take into account disordered or noncircular motions through the gas. In this way, they found a tighter relation up to z = 1.2, with no significant evolution, and conclude that the observed scatter in the TFR is partly due to disturbed and peculiar morphologies, corresponding to lower rotation velocities for their masses. In Figure 6, the result obtained for our z ∼ 2 sample using the S_0.5 index is shown, using the fitted σ₀₂ as a measure of the intrinsic dispersion of the gas not due to rotation. Although the σ used here does not exactly match the definition adopted by Kassin et al. (2007), it is the best indicator of turbulent and disturbed motion obtained from our modeling. Using the S_0.5 index we also detect an evolution of the relation between z ∼ 1 and z ∼ 2.2, further supporting the claim of an evolution of the TFR. However, due to the different definition of σ used here, we do not attempt to quantify this evolution. The resulting scatter is slightly smaller than the already remarkably low one obtained for the classical TFR (see Figure 5). We suggest that this is due to the optimal selection of our rotational supported systems at z ∼ 2, as the SINFONI three-dimensional data are ideal to ensure a much lower contamination from complex dynamics and irregular motions than in longslit data sets.

Zoom In Zoom Out Reset image size