Stellar Dynamical Models for 797 z ∼ 0.8 Galaxies from LEGA-C

The formation and evolutionary history of galaxies is encoded in the motions of their stars. The dynamical state of the Milky Way is now being revealed in ever more detail with Gaia (Gaia Collaboration et al. 2018a, 2018b) and stellar kinematics of nearby galaxies has been a very active area of research since the 1970s (e.g., Bertola & Capaccioli 1975; Binney 1978; Kormendy & Illingworth 1982; Davies et al. 1983). In the present-day Universe, much has been learned from resolved kinematics measurements made with single-slit observations (e.g., Faber & Jackson 1976; Tully & Fisher 1977; Djorgovski & Davis 1987; van der Marel 1991). More recently, this has been superseded by integral field observations, revealing kinematic patterns in two dimensions (e.g., Emsellem et al. 2004; Krajnović et al. 2011) and qualifying the dynamics of various types of galaxies (Emsellem et al. 2011; Graham et al. 2018). In particular, early-type galaxies are found to exhibit very diverse rotational patterns and are often split into fast rotators and slow rotators. Fast rotators are consistent with being axisymmetric, oblate rotators, whereas slow rotators are more dispersion dominated and often have complex dynamical structures (Emsellem et al. 2007). These observations have in turn sparked significant theoretical effort (e.g., Jesseit et al. 2009; Bois et al. 2010, 2011; Naab et al. 2014; Lagos et al. 2018a, 2018b; van de Sande et al. 2019), to understand the large variety in observed kinematic signatures and explain the varying degrees of ordered motion by investigating the relative importance of, e.g., (minor) merging and accretion. So far, a clear observational picture has not been produced at higher redshift. Yet it is highly desirable to understand the structure of high-redshift galaxies, when they were in the middle of forming the bulk of their stars and shortly thereafter. Ever following in the footsteps of studies of the local Universe, at higher redshift, we have largely relied on photometry to guide progress. We have learned that the population of galaxies at z ∼ 1 is bimodal in color and can be separated into quiescent and star-forming galaxies based on optical colors (e.g., Labbé et al. 2005; Williams et al. 2009). Like in the local Universe, these groups show distinct structural properties (e.g., Bell et al. 2012) and correlations (e.g., Franx et al. 2008; van der Wel et al. 2014a). In particular, projected ellipticities have been used to study intrinsic shapes of galaxy populations which, for lack of more direct measures, have been associated with a galaxy's state of rotational support (Chang et al. 2013; van der Wel et al. 2014b). Oblate, axisymmetric galaxies are found to form a larger fraction of the total galaxy population with increasing redshift, which is interpreted as the increased rotational support at high redshift.

Kinematic studies of ionized gas have revealed that large, star-forming disks exist at least as early as z ∼ 2 (Förster Schreiber et al. 2006; Genzel et al. 2008), but that those are more pressure-supported than gas disks in the present-day Universe (e.g., Kassin et al. 2012; Wisnioski et al. 2015). Moreover, an increasing fraction of lower mass galaxies have not yet formed regular, rotating structures (e.g., Förster Schreiber et al. 2009). This evolutionary trend has been successfully reproduced by hydrodynamical simulations (Ceverino et al. 2014; Martig et al. 2014; Pillepich et al. 2019). However, there are limitations to the use of ionized gas as a tracer of dynamical structure: while to first-order, the integrated velocity dispersions of ionized gas and stars agree fairly well (Bezanson et al. 2018a), there is large scatter for several reasons. First, it reveals only the state of the youngest parts of galaxies. Second, the spatial distribution of the emission is different from the total mass distribution: it has a patchy structure, which makes it a problematic tracer (Varidel et al. 2019). The stellar component does not suffer from these drawbacks: they are the dominant mass component in most massive galaxies (within the inner 5–10 kpc for most galaxies in our sample), and the stellar light itself is a reasonable proxy for the shape of the mass distribution.

But the necessary data to access the stellar dynamical structure–high-signal-to-noise (S/N) continuum spectroscopyis is challenging to obtain and have so far been limited to small samples (van der Marel & van Dokkum 2007; van der Wel & van der Marel 2008; Newman et al. 2015, 2018; Guérou et al. 2017; Toft et al. 2017) that only include very massive, passive galaxies. Even when these measurements can be made, the interpretation is complicated by the impact of observing conditions, such as the slit geometry and beam smearing. To make progress, detailed dynamical models can help. At the expense of making assumptions about the underlying galaxies (e.g., axisymmetry and regular motions) one can account for the impact of the large slit width, inclination, and seeing and estimate the intrinsic dynamical properties of galaxies.

The LEGA-C survey (van der Wel et al. 2016) is the first survey to provide resolved stellar kinematics for a magnitude-limited sample of galaxies at z ∼ 0.8. Bezanson et al. (2018b) presented a first, qualitative look at the stellar rotation curves and dispersion profiles for a sample of 104 quiescent galaxies in the LEGA-C survey. It was found that quiescent galaxies at z ∼ 0.8, like their local counterparts, show a decrease of rotational support with increasing mass and that, like in the local population, there is a subset of massive galaxies seemingly devoid of significant rotation altogether. But crucially, the LEGA-C galaxies appeared to show significantly more rotation than local galaxies. However, the interpretation of the observed kinematics is complicated by the observational setup: the slit has a considerable width, the seeing is in many cases comparable to the galaxy size, and the slits are not necessarily oriented along the galaxy's major axis. An accurate quantitive description of the dynamical structure of these galaxies therefore requires dynamical modeling.

In this paper we present, for the first time, a concerted effort to infer the dynamical masses and dynamical states for a large, representative sample in the LEGA-C survey of quiescent and star-forming galaxies. We achieve this by using axisymmetric Jeans models. The Jeans method is frequently used to model galaxies in the local Universe, for instance to constrain mass-to-light ratios and dark matter fractions (van der Marel 1991; Williams et al. 2009; Cappellari et al. 2013), or constrain the internal rotation and anisotropy of galaxies (Satoh 1980; Zhu et al. 2016). Similarly, at z > 0 they have been applied to study the dynamical masses (Guérou et al. 2017), initial mass functions (Newman et al. 2017), and rotation rates (van der Wel & van der Marel 2008; Newman et al. 2018).

In this paper, we present the Jeans models for 797 galaxies that are part of the third data release of the LEGA-C survey (van der Wel et al. in preparation). In Section 2, we introduce the data and the sample selection used in this analysis. The models are described in Section 3. In Section 4, we show the results of the fits from the models to the observations, as well as a number of consistency checks. Some basic results from the derived quantities are shown in Section 7, after which Section 8 summarizes and looks forward to anticipated follow-up work. Throughout the paper, a standard ΛCDM cosmology is assumed with H₀ = 70 km s⁻¹, Ω_Λ = 0.7, and Ω_b = 0.3.

2.1. LEGA-C Spectroscopy

2.2. Ancillary Imaging Data

2.3. Sample Selection for Jeans Modeling

The third (and final) LEGA-C release (van der Wel et al. 2021) contains 3028 spectra of K_s -band selected "primary" galaxies in the redshift range 0.6 < z < 1. We construct dynamical models for a subset of these, those (1) with the major axis is not misaligned with the slit by more than 45° (1554 left), (2) without a morphology flag (mergers, irregulars, etc.) or a spectroscopy flag (significant AGN contributions or severe flux calibration problems) as defined in the data-release paper (1381 left), and (3) with S/N > 10 in at least three spatial resolution elements (861 left).

Due to the observing strategy and the way in which galaxies are prioritized, some galaxies are observed in multiple masks. This means we have independent measurements of the same galaxy under slightly different seeing conditions and slit alignments. In the sample of 861, we have 795 unique galaxies: there are 65 "duplicates" (7.6% of the sample), which we use for consistency checks of the method (Section 5.2).

Figure 1 shows the distribution of z_spec and σ_*, the integrated stellar velocity dispersion, and how the selection criteria affect these distributions. Figure 2 shows the bias against edge-on star-forming galaxies: these galaxies have higher dust attenuation and often fall short of our S/N criterion. This is an important consideration when attempting to quantify the relative numbers of flat and round star-forming galaxies, and the evolution thereof. The trend does not exist for the mostly dust-free quiescent galaxies. For those (and for face-on star-forming galaxies) the S/N can be too low due to low stellar mass and/or high redshift.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the single-slit spectra we have one kinematic axis and a spatially resolved profile for v_⋆ and σ_⋆ for typically ∼8–12 spatial pixels, reaching over two effective radii for some of the galaxies in our sample. Figure 3 shows the range probed for the galaxies in the sample, in arcsec as well as in units of the galaxy effective radius R_e.

Zoom In Zoom Out Reset image size

3.1. Axisymmetric Jeans models

The dynamical models used in this work are based on the Jeans equations, a set of stellar hydrodynamic equations relating (average) kinematic properties to the stellar light and the underlying matter distribution. They are derived by taking velocity moments of the Collisionless Boltzmann Equation, describing the evolution of the phase-space density, assuming the system of particles to be collisionless and in steady state (e.g., Binney & Tremaine 2008).

In this work, we implement this using the Jeans Anisotropic MGE models (JAM; Cappellari 2008). ¹² This software uses the computationally convenient Multi-Gaussian Expansion (described below) to model the surface brightness, the gravitational potential and the point-spread function (PSF, see Sections 3.2 and 3.3). The Jeans equations become solvable by making observationally motivated assumptions on the velocity ellipsoid. Via the anisotropy parameter ${\beta }_{z}\equiv 1-\langle {v}_{z}^{2}\rangle /\langle {v}_{R}^{2}\rangle$ (in cylindrical coordinates) we can model the flattening of the velocity ellipsoid in the direction perpendicular to the disk of the galaxy. The velocity ellipsoid is otherwise assumed to be aligned with the main cylindrical axes such that cross-terms (〈v_z v_R 〉, 〈v_ϕ v_R 〉, 〈v_ϕ v_z 〉) all vanish. This leaves us with one independent (possibly spatially varying) second velocity moment at every position. Given the inclination, this can then be integrated along the line of sight as done in Equation (28) of Cappellari (2008). This prediction can then be compared to the observed combination of first- and second-order line-of-sight velocity moments ${v}_{\mathrm{rms}}\equiv \sqrt{{v}_{\star }^{2}+{\sigma }_{\star }^{2}}$, allowing us to constrain the kinematic structure and gravitational potential with only a few parameters.

The parameter v_rms allows us to constrain the total rms motion in the galaxy as a sum of rotation and dispersion, but it gives no information about the relative importance of each of the components. To distinguish between these we make the additional assumption that the dispersion in the azimuthal direction is related to the dispersion in the radial direction such that

$$\begin{eqnarray}&&\overline{{v}_{\phi }}=\kappa \sqrt{\overline{{v}_{\phi }^{2}}-\overline{{v}_{R}^{2}}}\end{eqnarray} \tag{ 1 }$$

where the constant rotation parameter κ is a fitting parameter (van der Marel & van Dokkum 2007; Binney & Tremaine 2008).

3.2. Light and Mass Modeling

An important input of the dynamical models is a detailed parameterization of the tracer light (stellar), from which the observed kinematics originates, as well as the underlying gravitational potential. In our model, the potential is due to two mass components: the stellar body and dark matter. That latter, in actuality, represents any component that breaks the mass-follows-light assumption inherent in the former, but this will be discussed in full below. Both are described with a set of elliptical Gaussians using the Multi-Gaussian Expansion code (Cappellari 2002). ¹³ Extra Gaussian components are iteratively added until the fit does not improve significantly. All Gaussians are assumed to have the same, global PA (no isophotal twists). In addition, we impose a regularization to limit the axis ratios around the global value as long as it does not affect the quality of the fit, which is most relevant if we add residuals or decompose mass profiles (see below). This is important since the minimum axis ratio limits the minimum allowed inclination in the dynamical fit. In some cases (in particular, for Sérsic profiles that fall off faster than a single Gaussian function, n < 0.5), we also allow negative Gaussian components while ensuring that the sum stays positive.

3.2.1. Light Tracer and Stellar-mass Model

We use the ACS/F814W imaging as the light tracer model. The HST/ACS instrumental PSF has a non-negligible size relative to the galaxy sizes, which is why we use the best-fit Sérsic profiles to describe the deconvolved surface brightness distributions, following van der Wel et al. (2012) and van der Wel et al. (2016). These are then subsequently fitted with the MGE, as illustrated in Figure 4, where we show that this procedure accurately follows the true light profile over six magnitudes in surface brightness, out to 2.5 effective radii.

Zoom In Zoom Out Reset image size

We note that, for these relatively low S/N F814W data, direct MGE fits systematically underestimate the fluxes, on average, 20% as inferred from a comparison with ground-based photometry, where we synthesized F814W photometry from the SED fits to the ground-based photometry described in Section 2.2. Low surface brightness light is missed by the MGE, leading to underestimated effective radii. The Sérsic model fluxes, on the other hand, show no such offset, justifying the extrapolation of the profiles out to large radii. The 1σ scatter for both the MGE and Sérsic fluxes is 15% with respect to ground-based fluxes.

3.2.2. Dark Matter Model

In addition to the luminous matter, we also include a dark matter component. Later, in Section 5.2, we will see that the primary role of the dark matter fit component is to allow us to fit mass-to-light variations within the inner 10 kpc or so of the galaxy, having no predictive power outside of that radius. Consequently, while we choose a functional form for this component that is rooted in cosmology, its detailed functional form is not important and does not significantly influence the results of our overall fits.

We adopt a spherically symmetric NFW profile (Navarro et al. 1997):

$$\begin{eqnarray}&&\rho =\displaystyle \frac{{\rm{\Delta }}{\rho }_{\mathrm{crit}}}{(r/{r}_{s}){\left(1+(r/{r}_{s})\right)}^{2}}\end{eqnarray} \tag{ 2 }$$

where ρ_crit is the critical density, r_s a characteristic scale, and Δ the threshold level delimiting the halo, as a matter of convention often chosen to be 200. The halo mass M₂₀₀ then corresponds to the mass inside the radius r₂₀₀ for which the average density is 200ρ_crit. Furthermore, it is convenient to define the concentration parameter c₂₀₀ = r₂₀₀/r_s. It is well known that there is a correlation between halo mass and concentration, such that more massive haloes tend to be less concentrated, though the scatter is large. From the simulations in Dutton & Macciò (2014) we can use the relation between mass, concentration, and redshift to eliminate one of the two halo parameters:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left({c}_{200}\right)=a+b{\mathrm{log}}_{10}\left({M}_{200}/[{10}^{12}{h}^{-1}{M}_{\odot }]\right)\end{eqnarray} \tag{ 3 }$$

where the parameters a and b are given (for Δ = 200) by

$$\begin{eqnarray}&&a=0.520+\left(0.905-0.520\right){e}^{-0.617{z}^{1.21}}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&b=-0.101+0.026z.\end{eqnarray} \tag{ 5 }$$

This allows us to reduce the dark matter component to a one-parameter model: its mass M₂₀₀. These models can be easily included in the dynamical fit by expanding the NFW profile in a series of Gaussian functions and adding them to the MGE terms of the stellar matter.

We note that, in practice, this implementation means that the dark matter component always has the same slope (the inner slope from the NFW profile) and the only variable is the normalization, as dictated by M₂₀₀. The reason for this is that r_s always lies beyond the spatial extent of our observations. We stress again that we interpret the dark matter component not as a constraint on the mass of the entire dark matter halo, but as a constraint on deviations from the mass-follows-light assumption, which physically can be due to dark matter, but also contributions from gas and stellar M/L variations. Global differences in stellar M/L are accounted for by varying the M/L of the light component, as are gas disks that follow the same profile as the stellar light. However, if gas or stellar-mass profiles differ from the stellar light profile, then the "dark matter" component will attempt to account for this. The assumption of a spherical dark matter distribution is then not correct, but for the enclosed mass estimate, this has a minor effect.

We note that we have to allow for a dark component with a negative mass. This may be counterintuitive, since setting a strictly positive prior seems physically self-evident. The two reasons for allowing a negative mass here are (1) to avoid a positive bias in the total mass and (2) to allow for inverted stellar M/L gradients. The positive bias occurs when the dark component is not constrained by the data, such that the median in the posterior is always positive. Allowing for negative masses improves the calculation of the uncertainties by symmetrizing the posterior around zero if there are no constraints on the dark component, avoiding the aforementioned bias. A side effect is that, for a small set of galaxies, we find a negative total mass (see Section 6), which should be interpreted as unavoidable noise in the measurements.

3.3. Effective Seeing

We correct for beam smearing on an object-by-object basis. Thus, we wish to quantify the effective seeing of each galaxy, which reflects a combination of atmospheric seeing as well as observational and instrumental effects that broaden the light distribution. We do this by comparing the wavelength-collapsed LEGA-C spectra (equivalent to a radial light profile of that part of the galaxy covered by the slit) with a model light curve. The model light profile is obtained by convolving the ACS/F814W images successively with a seeing kernel and the slit. Then we re-bin this to the 0205 VIMOS pixel scale, and integrate the resulting image.

For the seeing kernel, we use a Moffat profile (Moffat 1969):

$$\begin{eqnarray}&&M(\alpha ,\gamma )=A{\left(1+\displaystyle \frac{{\left(x-{x}_{0}\right)}^{2}}{{\gamma }^{2}}\right)}^{-\alpha }\end{eqnarray} \tag{ 6 }$$

where A is the normalization and γ and α are the Moffat parameters for which we fit. Given α and γ, it is straightforward to show that the full-width-half-maximum (FWHM) is given by ${FWHM}=2\gamma \sqrt{{2}^{1/\alpha }-1}$. Generally, a Moffat profile has been shown to give a better description of the effect of atmospheric turbulence than a Gaussian function (Trujillo et al. 2001) because of the wings, controlled by the α parameter. We allow for a possible gradient in the background flux (due to, e.g., imperfections in the spectroscopic data reduction) by adding a variable linear component to the Moffat-convolved model light profile. If α is poorly constrained, we fix it at α = 4.765 (Trujillo et al. 2001) and only fit for γ. For objects where we did not find a good fit, we use α = 4.765 and γ derived from the mask-average FWHM and α.

It is encouraging that, despite atmospheric dispersion, drift, flexure, and hysteresis the true PSF (typically, 075) is in fact better than the DIMM seeing reported at the telescope (typically, 095, Straatman et al. 2018). Interestingly, the variation from mask to mask is smaller than the variation from object to object within a given mask. This suggests that technical details of the instrument (focus, positional accuracy, alignment, and location within the mask) are more important in determining the outcome of a spectrum than the variation of atmospheric conditions averaged over 20 hr of integration (good conditions were required to execute observations).

3.4. Systematic Effects in Slit Spectroscopy

In general, there are several things to keep in mind when using slit spectroscopy, in particular at the targeted redshift range where the slit width of 1'' is comparable to the effective radii of galaxies.

3.4.1. Slit Misalignment

During observations, the slits are centered on the galaxies as accurately as possible, but this might still be prone to both random and systematic offsets. This will lead to an offset in the measured velocity profile, as well as asymmetries in the v and σ profiles if the slit is not aligned with the galaxy's major axis, even if the galaxy is perfectly axisymmetric. Although for a small subset of individual sources, this offset can be determined (see Appendix), finding constraints systematically for all sources was not possible. Nevertheless, the sources for which we can constrain the offset suggests that the slits are generally well centered, with offsets of typically 0, 1, or sporadically more pixels, where 1 pixel corresponds to 0205. This motivates us to use a Gaussian prior on the slit center x₀ in the dynamical models, centered at 0, with a standard deviation of 1.5 pixels.

3.4.2. Slit Orientation

In deriving the rotational velocity at any particular spatial position, the light is tacitly assumed to come from a homogeneous light distribution that fills the entire 1 arcsec slit in the wavelength direction (this is how the spectral resolution R in instrument manuals is usually defined). However, this is not true for galaxies that have light disributions that change substantially in wavelength direction of the slit, which is the case for stars or distant galaxies. This results in spurious velocity shifts if the major axis is misalignment with respect to the slit. We derive and subtract this offset as follows. First, we convolve the F814W MGE with the PSF. For every pixel, we then find the light-weighted center, which we convert into a shift in wavelength, using the dispersion of 0.6 Å per pixel. Using the median of the observed spectral coverage of each object, we then convert this into a shift in velocity. This is illustrated in Figure 5. This artificial offset reaches up to 20 km⁻¹, being highest away from the center and for the flattest, most misaligned galaxies.

Zoom In Zoom Out Reset image size

3.4.3. Pixel Convolution

To compare the JAM predictions to the observations, we need to convolve the model grid with a kernel that is a combination of the seeing PSF as well as the slit. The former is an MGE of the Moffat profile derived in Section 3.3. The latter is a top-hat function of 1'' in width (the size of the slit in the wavelength direction).

3.5. Inclination and Anisotropy Priors

The observed kinematics depend on the inclination, but the only a priori constraint we have is the observed axis ratio. If we know the distribution f(p, q) of the intrinsic axis ratios p ≡ a/b and q ≡ b/c of a given class of objects (e.g., quiescent galaxies in a given redshift range), we can construct the likelihood of observing the observed axial ratio ${q}^{{\prime} }={a}^{{\prime} }/{b}^{{\prime} }$ for a given inclination. This will then serve as the prior on the inclination in the fits of the Jeans models to the total rms motion in the galaxy.

The intrinsic axis ratios p and q are related to the observed axis ratio ${q}^{{\prime} }$ through two viewing angles. In the oblate axisymmetric case (with p = 1), this becomes a simple relation between the inclination, the observed axis ratio ${q}^{{\prime} }$ and the intrinsic short-to-intermediate ratio q:

$$\begin{eqnarray}&&{q}^{{\prime} }={\cos }^{2}i+q{\sin }^{2}i.\end{eqnarray} \tag{ 7 }$$

We use the methodology described by Chang et al. (2013) to reconstruct the intrinsic q distribution, assumed to be Gaussian, from an observed set of axis ratios, as shown in Figure 6. For our sample, we find a mean μ_q = 0.41 and a scatter σ_q = 0.18. This distribution allows us to derive a prior for the inclination: given an observed, projected axis ratio, what is the probability distribution for the intrinsic shape and, hence, the inclination? The mathematical tools for this inversion are presented by van de Ven & van der Wel (2021).

While we see in Figure 6 that star-forming galaxies, on average, have flatter intrinsic shapes than quiescent galaxies (Chang et al. 2013; van der Wel et al. 2014b), adopting different shape distributions as priors for the inclination would introduce a systematic difference in rotational properties for quiescent and star-forming galaxies, even if they were kinematically the same. This spurious result would be explained, trivially, by the priors and not by the stellar kinematics. In order to avoid this, we adopt the same intrinsic shape distribution (that is, the same inclination prior) for different types of galaxies such that any difference in dynamical structure can be attributed to differences in the observed kinematics. As we will see, despite not assigning different priors to different types of galaxies, the model successfully distinguishes between rotating and non-rotating galaxies. Given the assumed inclination, the anisotropy β_z can in principle be well constrained with 2D kinematics (Cappellari 2008) but the long-slit data is expected to put weak constraints on the anisotropy. We therefore use a uniform prior of −0.5 < β_z < 0.5.

Following the methods and assumptions described in Section 3, we start with a fit of the line-of-sight second velocity moment of the JAM to the observed ${v}_{\mathrm{rms}}=\sqrt{{v}^{2}+{\sigma }^{2}}$. For axisymmetric galaxies that are centered in the slit, the velocity profile is symmetric. Observing such a galaxy, one would expect observations at + r and − r to be independent observations of the same value. This is not true if the slit is slightly offset (see also Section 3.3) which is why we do not symmetrize the data. However, to ensure that asymmetric data points do not influence the fits too much, the uncertainties are set to the maximum of the formally measured uncertainties in v_rms and half the difference between the v_rms values at + r and − r. In addition, at all radii, a minimum uncertainty of 5 km s⁻¹ is imposed to avoid having data points with spuriously low uncertainties dominate the fits. The fits extend to where we have measurements of σ, even though the velocity measurements tend to extend further out. Finally, we measure and subtract the systematic velocity as a function of the slit center x₀ by fitting an isotropic (β_z = 0), stellar-mass only, model at every value of x₀.

The fit is done using the Markov Chain Monte Carlo software emcee ¹⁴ (Foreman-Mackey et al. 2013), which samples the parameter space with an ensemble of "walkers." The five free parameters (see Table 1 for an overview) in the fit are the stellar-mass-to-light ratio (M/L), the virial mass of the dark matter halo parameterized by the circular velocity (v_c ), the inclination (i), the anisotropy β_z , and the slit centering x₀. We use walkers, with a burn-in run of 200 iterations. Convergence was verified by visually examining the parameter values as the fit progresses, which usually occurs after several thousand iterations.

Table 1. Free Parameters in the Fit, With Assumed Priors

Parameter	Prior
β	Flat prior on −0.5 < β < 0.5
M/L	Unconstrained
vDM	Flat prior on −600 < vdm < 600 km s−1
Inclination	q-dependent prior, assuming an intrinsic thickness distribution of ${ \mathcal N }(0.41,0.18)$ (see Section 3.5)
Slit center, x0	Gaussian, centered on 0 with a 1.5 pixel standard deviation

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

As an example, Figure 7 shows the covariance between the various parameters, for a galaxy that was chosen to be typical in terms of the observed velocity dispersion (192 km s⁻¹) and radial extent of the observed kinematics (five pixels or 1'', one-sided). The inclination is unconstrained by the data and is equal to the prior, as is the slit center x₀. β_z is weakly constrained, as expected. Dark matter mass and stellar mass are degenerate, but the total mass within R_e ($\mathrm{log}{M}_{\mathrm{dyn},\lt \mathrm{Re}}$) is well constrained. Note that this is not an independent parameter but is computed from M₂₀₀ and M/L. The best-fit models have a median reduced χ² of 1.03, so given the simplicity of the models, the fits are more than adequate.

Zoom In Zoom Out Reset image size

The posterior for the rotation parameter κ is shown in Figure 7, but is not fit simultaneously (see Section 3.1). Instead, the best-fit κ is computed at every step of each Markov chain (β_z , M/L, i, M₂₀₀, x₀), without influencing the overall likelihood. The rotation parameter κ has been used in the past as a computationally convenient way of disentangling rotational from random motion, and for an isotropic oblate rotator, κ = 1 by definition. As long as the rotation is not larger than v_rms, κ can be arbitrarily higher than 1, which means that the model contains more rotation than is required to explain the observed flattening, and significant anisotropy (in the vertical direction) contributes to the thickness of the system. A maximum $\kappa ={\kappa }_{\max }$ can be defined for which v_rms = v_rot. If the best-fit is found for $\kappa \gt {\kappa }_{\max }$, we set $\kappa ={\kappa }_{\max }$. These cases can occur for various reasons: if the galaxy has an anisotopy that is not described by β_z (see Section 3.1); if the assumptions of the axisymmetric modeling are violated; if the assumption of a spatially constant κ is a poor one; or if κ is an improper (or too simplistic) parameterization in the first place.

Examples of the observed kinematic profiles and corresponding best-fit models in Figure 8 are shown for 12 galaxies spanning the full dynamic range of mass and kinematic structure. The axisymmetric models are sufficiently flexible to reproduce the full range of observed kinematic and structural properties: non-rotating and rotating galaxies, dispersion profiles that peak and dip in the center, late-, and early-type galaxies, face-on and edge-on galaxies, aligned and misaligned galaxies, and compact and extended galaxies.

Zoom In Zoom Out Reset image size

5.1. Is a Dark Component Required?

If we leave out the dark component, our model reverts to the commonly used mass-follows-light approach, where the total mass profile is identical to the MGE stellar light model, with only the M/L scaling as a free parameter. We show the comparison in Figure 9. Importantly, the formal uncertainties are very small for the mass-follows-light model (typically, just 6%). This is not realistic for total mass estimates of high-redshift galaxies with limited spatial resolution. The uncertainties for the default model are much larger and we verify in Section 5.2 that these are realistic based on duplicate observations. The median differences in mass are modest (<0.10 dex), but for individual galaxies, the default model often produces a mass that is two times higher than the mass-follows-light mass. Moreover, there are strong and systematic differences in the κ parameter that traces the degree of rotational support. For the mass-follows-light models, the round (near face-on) galaxies have the highest κ values. This is counterintuitive and implies a bias: round galaxies should have similar or lower κ values compared to flat galaxies.

Zoom In Zoom Out Reset image size

On the basis of the biased κ estimates and the unrealistically small uncertainties for the mass-follows-light models, we conclude that the dark component is an essential ingredient. It allows for the required flexibility in modeling the observed kinematic data, especially for face-on disks. As discussed in Section 3.2.2, the dark component does not necessarily reflect dark matter alone, but accounts for deviations from the mass-follows-light assumption due to stellar M/L variations and the contribution from gas as well.

5.2. Duplicate Observations

A number of galaxies have been observed multiple times in LEGA-C as a consequence of the observing approach. These provide independent kinematic measurements of the same object, so we can leverage this to test the robustness of the parameters we derive and obtain an independent estimate of all uncertainties associated with the spectroscopic data. In the sample with dynamical models, there are 65 such duplicates. Figure 10 compares the dynamical mass and rotation parameter κ for these galaxies. The two sets of observations scatter around the one-to-one line with a scatter that corresponds reasonably well with the expectation from the formal (random) uncertainties. Theses uncertainties, inferred from the posterior distributions, are on average 50% larger than the scatter for the mass estimate. We attribute this to the flexibility in our model afforded by a mass component that does not follow the light distribution. The comparison of duplicate observations may imply that this conservative approach leads to overestimated random uncertainties. On the other hand, the light tracer, inclination prior, and other elements in the model are identical for the duplicates and any uncertainties associated with those remain hidden in this test: the duplicate results are only independent in terms of kinematic measurements, but the light tracer, stellar-mass model, inclination, and anisotropy priors are identical. The prudent decision, then, is to use the uncertainties as inferred from the posterior distributions.

Zoom In Zoom Out Reset image size

The dynamical modeling results are made available electronically, and Table 2 shows an excerpt. We provide the total mass in five different spherical apertures: 1, 5, 10 kpc, R_e, ¹⁵ and 2R_e; these are spherical apertures, rather than projected, cylindrical apertures. A value is only included in the table if that aperture does not amount to a significant extrapolation given the radial extent of the stellar kinematic data. We allow for a maximum of 20% extrapolation in radial extent. Given this constraint, the number of objects with mass estimates within these five apertures are 861, 809, 206, 730, and 319, respectively. As one would expect, for most objects, the radial extent of the kinematic measurements is similar to the effective radius: the median kinematic radial extent-to-effective radius ratio is 1.45 (see Figure 3).

Table 2. Stellar Dynamical Parameters

ID1	ID2	${r}_{\max }$	Re	logM1	logM5	logM10	$\mathrm{log}{M}_{{R}_{e}}$	$\mathrm{log}{M}_{2{R}_{e}}$	κ	${V}_{{R}_{e}}$	${\sigma }_{{R}_{e}}$
		kpc	kpc	M⊙	M⊙	M⊙	M⊙	M⊙		km s−1	km s−1
5	4792	6.9	5.0	${9.99}_{-0.06}^{+0.08}$	${10.88}_{-0.24}^{+0.27}$	⋯	${10.88}_{-0.24}^{+0.27}$	⋯	${1.10}_{-0.04}^{+0.07}$	${179}_{-34}^{+27}$	${90}_{-4}^{+5}$
15	6556	4.7	3.2	${10.68}_{-0.03}^{+0.03}$	${10.73}_{-0.35}^{+0.53}$	⋯	${10.84}_{-0.12}^{+0.18}$	⋯	${0.34}_{-0.34}^{+0.20}$	${69}_{-90}^{+7}$	${204}_{-16}^{+13}$
16	6859	4.9	7.6	${10.83}_{-0.03}^{+0.04}$	${11.33}_{-0.16}^{+0.17}$	⋯	⋯	⋯	${0.36}_{-0.06}^{+0.15}$	${120}_{-2}^{+4}$	${214}_{-11}^{+9}$
19	7419	13.0	11.5	${9.72}_{-0.31}^{+0.25}$	${11.00}_{-0.20}^{+0.19}$	${11.46}_{-0.20}^{+0.21}$	${11.54}_{-0.21}^{+0.23}$	⋯	${0.78}_{-0.09}^{+0.18}$	${123}_{-7}^{+7}$	${160}_{-13}^{+13}$
26	10462	4.4	3.9	${9.80}_{-0.12}^{+0.14}$	${10.70}_{-0.39}^{+0.52}$	⋯	${10.55}_{-0.36}^{+0.46}$	⋯	${0.57}_{-0.57}^{+4.42}$	${80}_{-19}^{+6}$	${58}_{-12}^{+10}$
27	10902	5.9	12.0	${10.22}_{-0.05}^{+0.07}$	${11.02}_{-0.25}^{+0.25}$	⋯	⋯	⋯	${1.04}_{-1.04}^{+3.96}$	${211}_{-71}^{+34}$	${119}_{-10}^{+11}$
38	14375	6.2	2.8	${10.76}_{-0.02}^{+0.02}$	${11.27}_{-0.16}^{+0.17}$	⋯	${11.11}_{-0.08}^{+0.07}$	${11.29}_{-0.18}^{+0.19}$	${0.70}_{-0.06}^{+0.10}$	${241}_{-6}^{+5}$	${208}_{-31}^{+14}$
39	14729	6.3	10.4	${10.36}_{-0.06}^{+0.06}$	${11.21}_{-0.18}^{+0.16}$	⋯	⋯	⋯	${1.18}_{-0.08}^{+0.30}$	${306}_{-53}^{+43}$	${175}_{-17}^{+16}$
62	25994	4.7	3.1	${10.41}_{-0.05}^{+0.05}$	${11.10}_{-0.23}^{+0.19}$	⋯	${10.91}_{-0.16}^{+0.11}$	⋯	${1.08}_{-0.12}^{+0.05}$	${249}_{-49}^{+22}$	${145}_{-7}^{+7}$
68	26283	7.8	4.2	${10.51}_{-0.03}^{+0.04}$	${10.95}_{-0.35}^{+0.33}$	⋯	${10.91}_{-0.28}^{+0.26}$	${11.01}_{-0.73}^{+0.71}$	${0.13}_{-0.13}^{+0.04}$	${29}_{-34}^{+1}$	${157}_{-9}^{+7}$
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

Note. (1) LEGA-C ID (DR3); (2) UltraVISTA ID (Muzzin et al. 2013b); (3) extent of stellar kinematic data; (4) effective radius of the galaxy; (5–9) model mass enclosed in spherical apertures; (10) rotation parameter; (11) rotational velocity at R_e (edge-on); (12) velocity dispersion within R_e (face-on). Values and uncertainties are based on 16th, 50th, and 84th percentiles of posterior parameter distributions. The machine-readable table has 861 entries. We choose the following notation for a negative mass (see text): logM = 10 means M = −10¹⁰ (not M = 10⁻¹⁰).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Please note the presence of a small subset of galaxies with negative mass estimates. This arises due to the prior on the dark mass component that is centered at zero to avoid a positive bias in the total mass estimates (see Section 3.2.2). All negative mass estimates have large uncertainties, so that a positive mass is never ruled out with more than 95% confidence, as is expected given the derived uncertainties.

The rotation parameter κ is described in Section 3.1. The rotational velocity ${V}_{{R}_{{\rm{e}}}}$ is calculated at R_e for the edge-on, deprojected model. The effective velocity dispersion ${\sigma }_{{R}_{{\rm{e}}}}$ is calculated for the light-weighted, face-on deprojected model within R_e. The catalog contains best values and lower and upper error bars. These are obtained from the posterior distributions, where the best value is the median and the upper (lower) error bars reflect the difference between the 84th (16th) percentile and the median.

A rotation proxy that is popular in the literature is the projected quantity λ_R, which is frequently used by local IFU studies (and first introduced for SAURON in Emsellem et al. 2007). Calculating λ_R directly from our data is not useful given the large seeing convolution kernel and single kinematic axis, so instead we could infer it from the model. But, as shown by Harborne et al. (2019, 2020), even surveys of relatively nearby galaxies suffer from significant seeing convolution and the results would not comparable, and it would be misleading to name such a quantity λ_R. Moreover, given the dependence on the model, it makes little sense to leave in inclination-dependent effects, hence our choice to calculate face- and edge-on projected quantities. A direct comparison with present-day galaxies is, at this point, not possible: the only way to infer evolution in the intrinsic dynamical structure of galaxies is create the equivalent dynamical models for present-day galaxies.

In our own previous work (Bezanson et al. 2018b), we measured (v₅/σ₀) without correcting for seeing or slit losses—instead we used IFU data from CALIFA (Sánchez et al. 2012), downgraded to resemble our high-redshift data, in order to infer redshift evolution. Our models now provide the intrinsic (v₅/σ₀) and we find the effect of downgrading the CALIFA data agrees very well with the change from intrinsic (model) (v₅/σ₀) and directly observed (v₅/σ₀). This validates the conclusions drawn by Bezanson et al. (2018b) that high-redshift quiescent galaxies show more rotation than their present-day counterparts, which we interpreted as direct evidence for merging driving galaxy growth and angular momentum reduction.

It is beyond the scope of this paper to explore rotational support, angular momentum, and the masses of galaxies in great detail, and in forthcoming papers we will explore the connection between star formation history and dynamical structure, the various mass components of galaxies, and the comparison with dynamical modeling based on ionized gas kinematics. But to illustrate the richness of our data we show in Figures 11 and 12 the distribution of galaxy shapes, ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ and dynamical mass ${M}_{{R}_{e}}$.

Zoom In Zoom Out Reset image size

The left-hand panel of Figure 11 is strongly reminiscent of integral field surveys of present-day galaxies (e.g., Emsellem et al. 2011; van de Sande et al. 2017; Falcón-Barroso et al. 2019), with a strong correlation between projected shape and (projected) degree of rotational support. Low Sérsic index (late-type) galaxies generally have higher ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ than high Sérsic index (early-type) galaxies, even at fixed projected ellipticity. In particular, we see a population of rather flat ( ∼ 0.4), but non-rotating high-n galaxies. Such galaxies also exist in the present-day Universe and can be either anisotropically flattened systems or dynamically complex galaxies with, for example, counter-rotating disks (Rix et al. 1992). The spatial resolution of our data is obviously insufficient to distinguish between the different explanations, but we can conclude that the net rotation is low compared to the degree of flattening. There is also a small set of data points at very large ${(V/\sigma )}_{{R}_{{\rm{e}}}}$, away from the general locus: here the best-fit dynamical model has a high anisotropy (with very small vertical dispersion), but with large uncertainties (not shown in the figure, but included in the catalog).

The crowding in the lower-left part of the figure (round, non-rotating galaxies) is largely due to projection effects. Therefore, in the right-hand panel, we show the deprojected ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ (as explained in Section 6). Even for nearly round galaxies with ∼ 0.1 our modeling successfully recovered the degree of rotational support, as indicated by the lack of a correlation between ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ and for low-n galaxies. This population must consist of galaxies with a fairly homogeneous dynamical structure, with ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ ∼ 1 − 2, and for which the scatter in projected ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ is mostly explained by viewing angle (inclination). The high-n galaxies are less homogeneous in dynamical structure. Like in the present-day Universe, there is a population of intrinsically round, slowly rotating galaxies, mixed with a population that is fast rotating, and similar in dynamical structure to the low-n population. We note that our models successfully distinguish between rotation- and dispersion-dominated galaxies; our inclination prior, based on an intrinsically oblate shape, did not preclude this.

It is well documented that Sérsic index correlates well with star-forming activity at all redshifts (Bell et al. 2012), and it is informative to split the population into quiescent and star-forming galaxies on the basis of their rest-frame U−V and V−J colors, computed as explained in our second data-release paper (Straatman et al. 2018). In Figure 12, we show the ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ distribution as a function of dynamical mass ${M}_{{R}_{e}}$ for quiescent and star-forming galaxies. For the quiescent galaxies, we see that the most massive galaxies, as in the present-day Universe, are the most likely to be slowly rotating, while for the star-forming population we see that the most massive galaxies are the most rotation dominated. This immediately implies different evolutionary pathways for the most massive galaxies, depending on their star formation history.

Zoom In Zoom Out Reset image size

The trend for star-forming galaxies may, at first sight, seem counterintuitive, given that the most massive star-forming galaxies are more bulge-dominated than low-mass star-forming galaxies (which we see for our sample by an increase in Sérsic index). Yet these more bulge-dominated galaxies have higher ${(V/\sigma )}_{{R}_{{\rm{e}}}}$. But these trends are not mutually exclusive: indeed, σ increases with mass, signifying the presence of more massive bulges, but V increases faster with mass than σ, signifying the presence of a large, extended disk. Interestingly, however, the galaxies with the highest rotation speeds are quiescent: these are compact objects with high-velocity dispersions (>200 km s⁻¹), yet flattened and rotation dominated, with ${V}_{{R}_{e}}\sim 350$ km s⁻¹, reminiscent of what has been referred to as "fossil" galaxies in the present-day Universe (e.g., Ferré-Mateu et al. 2012; van den Bosch et al. 2012).

If we then consider the interesting set of low-mass star-forming galaxies with little rotation, visual inspection of the HST images and the rotation curves suggest that this can occur for several reasons: some of these galaxies are somewhat irregular systems for which the assumption of axisymmetry is likely invalid, but the assessment of low ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ still accurate; some galaxies are truly non-rotating systems and appear to be young, spheroidal galaxies that just fall short of being classified as quiescent; for other objects, the stellar light distribution obviously does not follow the stellar-mass distribution, and bright, blue-star-forming regions or spiral arms affect the analysis of the HST light profiles, and thereby the geometry of the adopted axisymmetric system in the dynamical model.

We present axisymmetric Jeans models for 797 z = 0.6–1 galaxies with spatially resolved stellar kinematic measurements from the LEGA-C survey (van der Wel et al. 2016; Straatman et al. 2018; van der Wel et al. 2021). The K_s-band-selected parent sample of over 3000 galaxies makes no selection on color or morphology, and we select the subset objects that satisfy additional criteria (S/N, slit alignment, and regular morphology), but is otherwise unbiased except for a bias against edge-on disks due to attenuation. This is by far the largest sample of galaxies with resolved stellar kinematics at this redshift to date, by a factor 30, and the only sample comprising star-forming/late-type galaxies, as previous samples were focused solely on passive galaxies. The application of Jeans models is important because unlike galaxies in the present-day Universe, high-redshift galaxies are at most a few arcseconds in extent: not much lager than both the width of the slits (1'') and the ground-based seeing (∼08). Our primary concern is to correct the observed rotation and velocity dispersion profiles (Bezanson et al. 2018b) for the convolution by slit and seeing, and an unknown inclination angle.

We use HST F814W imaging to describe the tracer light. The underlying gravitational potential is modeled with a mass-follows-light stellar component plus an NFW dark matter halo. Other free parameters are the inclination, the velocity anisotropy, and the positioning of the galaxy within the slit; these three parameters are not well constrained by our data, but we assess their impact on rotation and mass estimates by marginalizing over well-informed prior distributions. We provide a full treatment of the seeing and slit geometry on a galaxy-by-galaxy basis, and the JAM methodology is adapted to accept non-square pixels (Section 3.4). We derive (total) dynamical masses and proxies for the total, intrinsic rotational support of galaxies from these models (Table 2).

This data product provides the community with the tools to study the masses and kinematic properties of a large sample of galaxies, which was not possible before at large look-back time. For the first time, we show the ${(V/\sigma )}_{{R}_{{\rm{e}}}}$ distribution of galaxies at large look-back time (Figure 11) and the complexity in evolutionary pathways revealed by the large variety in dynamical structure seen even for galaxies with the same (dynamical) mass (Figure 12). These results are timely because cosmological hydrodynamical simulations have now reached the resolution to predict the stellar dynamical structure for large, representative samples. Further analysis of the models and the related data products will improve our understanding of these galaxies even further, and this will be published in a number of projects in the near future.

This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement 683184). Based on observations made with ESO Telescopes at the La Silla or Paranal Observatories under programs ID 194.A2005 and ID 1100.A-0949 (the LEGA-C Public Spectroscopic Survey). JvdS acknowledges support of an Australian Research Council Discovery Early Career Research Award (project number DE200100461) funded by the Australian Government.

Software: EAZY (Brammer et al. 2008), JAM (Cappellari 2008), MGE (Cappellari 2002), emcee (Foreman-Mackey et al.2013), corner.py (Foreman-Mackey 2016).

For most purposes, the slit can be assumed to be centered on the center of the galaxy. This is how the masks are designed, with accurate coordinates available from HST astrometry. Furthermore, given the angular sizes of galaxies compared to the slit width, even with a one-pixel offset in either direction, most light would still be captured. Even so, the kinematics will be affected by an offset, leading to (small) asymmetries and offsets in the observed velocity and dispersion profiles.

Here we ask the question of whether it is possible to constrain the slit position by comparing the observed LEGA-C light profile with the expected light profile from a rebinned, convolved, HST image with a slit placed at different offsets. For single, symmetric sources there are too many degeneracies to make this feasible, but galaxies with asymmetric profiles, or galaxies with a close neighbor that happens to fall within the slit, make this a viable option (Figure A1). Owing to their relative position, the flux ratio between the two sources is a function of the assumed offset (perpendicular to the slit, in the wavelength direction), which constrains the observed offset.

Zoom In Zoom Out Reset image size

Unfortunately, such sources are scarce, with on average about five useful examples per mask. This is insufficient for a robust determination, as there is no guarantee that galaxies within the same mask or even the same quadrant will share the same offset. However, given the ensemble of offsets across the entire survey, the majority suggest offsets of one to two pixels around 0, with incidental suggestions of bigger offsets, justifying the assumed prior in Section 3.