Galaxy and Mass Assembly (GAMA): Stellar-to-dynamical Mass Relation. I. Constraining the Precision of Stellar Mass Estimates

As outlined in Bellstedt et al. (2020a, 2020b), GAMA provides photometry in 19 bands, from far-UV to the IR, performed with the PRO FOUND package (Robotham et al. 2018). In brief, source detection and deblending are done based on combined images from KiDS r − band and VIKING Z-band, i.e., r + Z. Definitions of initial isophotal outlines (ProFound segments) from the r + Z detection image are then iteratively dilated to obtain a kind of curve-of-growth total flux measurement from the segments, while ensuring that neighboring segments/apertures do not overlap. The resultant FUV-IR Spectral Energy Distributions (SEDs) are corrected for foreground Galactic extinction using the Planck E(B − V) maps (Planck Collaboration et al. 2013) and matched back to spectroscopy and redshifts (DMU: gkvInputCatv02, Bellstedt et al. 2020a), which gives the basis for SPS modeling for stellar masses and ancillary stellar population parameters. To ensure consistency between the SPS-derived stellar-mass estimates and the Sérsic fits, we re-normalize the SEDs to match the Sérsic Z-band magnitude. In other words, we use the multiband ProFound photometry for colors (and hence stellar populations, M_⋆/L, etc.) and use the Sérsic photometry for total magnitude (and hence M_⋆).

The stellar-mass estimates that we use are derived following Taylor et al. (2011), namely using SPS models from Bruzual & Charlot (2003), the IMF of Chabrier (2003), exponentially declining SFHs, and single screen dust with the Calzetti et al. (2000) attenuation law. In the fits, the different bands are weighted to consider an approximately fixed rest-frame wavelength range of 3000–11000 Å. Stellar mass and other parameter estimates are done in a Bayesian way (i.e., marginalizing over a fixed grid in parameter space, with plausible priors), which encapsulates the important covariances and degeneracies in the modeling.

Finally, we add the more sophisticated stellar-mass estimates derived using the SED-fitting code ProSpect (Robotham et al. 2020), which allows for gas-phase metallicity (Z_gas) varying with time and in turn makes it possible to process more complicated SFHs. Notably, as examined in Robotham et al. (2020) and mentioned in Driver et al. (2022), ProSpect masses are systematically 0.06 dex larger than those derived by Taylor et al. (2011) that are used in this study.

To characterize the total flux, effective radii, and Sérsic index for the galaxies in our sample, we use the Sérsic-fit values derived for the VIKING Z-band, which provides a deeper probe in near-infrared with a magnitude limit of ∼23 mag (AB) and 1'' of seeing (see Bellstedt et al. 2020a, Table 1). The process and pipeline for the single-Sérsic fits that we use is presented by Kelvin et al. (2012). Galaxy sizes, i.e., effective radii (θ_e ''), and Sérsic indices (n) are based on 2D Sérsic model fitting to the surface brightness profile, with models truncated at 10 times the angular effective radius to prevent extrapolating flux into the uncertain regime of large radii. While these magnitudes are in agreement with Kron or Petrosian aperture based photometry up to n ∼ 3, they recover an additional ∼0.2 mag in high Sérsic index galaxies (Kelvin et al. 2012), as expected.

An important issue is that the parameters in the Sérsic fits can be strongly covariant: Sérsic index, total magnitude estimate and effective radius are interdependent, as discussed in detail by Kelvin et al. (2012). While the catalogs do not include the correlation coefficients that describe this covariance, it is possible to determine these values empirically using independent measurements in neighboring bands, as described in Appendix A.

The GAMA spectroscopic database incorporates previously obtained spectra and redshifts from multiple sources, including SDSS (Ahn et al. 2014), 2dFGRS (Colless et al. 2001b) and 6dFGS (Jones et al. 2004, 2009), with some targets in these surveys having been observed again by the GAMA team. This heterogeneous spectroscopic data actually presents a challenge in obtaining coherent and consistent velocity dispersion measurements for targets that have spectroscopic data obtained from different instruments and with different data reduction pipelines.

For GAMA DR4, we have used pPXF (Cappellari 2017) with templates from the stellar spectra library MILES (Sánchez-Blázquez et al. 2006) having a resolution of 2.51 Å (Falcón-Barroso et al. 2011) to derive velocity dispersion measurements for all of the spectra in the GAMA database. These measurements, which are given in the DMU VelocityDispersionsv01, are derived by fitting the 1D spectra as a non-negative linear superposition of a set of stellar template spectra, and assuming a Gaussian line-of-sight velocity distribution.

To ensure the coherence and consistency of our measurements, we re-calibrate the results from each independent survey through cross-survey comparison. Specifically, we re-calibrate the velocity dispersions from all other surveys to match those from SDSS, partly because SDSS has the best spectral resolution (see Driver et al. 2022, Table 2) and also provides a connection to SDSS peculiar velocity measurements (Howlett et al. 2022), which forms the basis of our next companion paper. Furthermore, we validate/calibrate the random errors in the measurements through comparisons of repeat observations of the same target within a single survey (i.e., intra-survey comparison). We find that this approach works well for the 6dFGS and GAMA data but not for 2dFGRS, where the process results in larger random and systematic errors. Accordingly, we only use velocity dispersions measured from the 6dFGS, SDSS, and GAMA spectroscopy. We defer further details of cross-survey and intra-survey comparisons and calibration to a future paper.

Finally, we combine these velocity dispersion measurements with the VIKING Z-band effective radii derived through Sérsic fits (Section 2.2), and use the Bertin et al. 2002 prescription for k(n) to calculate our simple dynamical mass estimates. This choice is mainly made because the empirical relation for k(n) of Cappellari et al. (2006) is derived based on only early-type galaxies, while we consider both early and late-type galaxies in this study. However, we also use the M_dyn derived using Equation (3) to study the relation between M_⋆ and M_dyn, and then compare the results to those from when Equation (2) is used.

Our basic sample selection is $\mathrm{log}{M}_{\star }/{M}_{\odot }\gt 10.3$ and 0.01 < z < 0.12. We adopt the following simple quality control cuts for our sample to ensure that our M_⋆ and M_dyn estimates are robust and comparable.

To ensure reliable spectroscopic redshifts, we require that targets have a GAMA redshift quality flag nQ ≥ 3, corresponding to a redshift confidence >90% (see, Equation (15) of Liske et al. 2015). Next, we discard cases where the Sérsic fits have failed or have given extreme values, by requiring that the Sérsic indices are within the range 0.3 < n < 10 to be included in the sample. Additionally, we discard cases where the total flux derived from the Sérsic fits differs from the iteratively dilated ProFound segment flux by more than 0.5 dex, which would imply an overall inconsistency between the SED-aperture and the Sérsic fits. As described in Section 2.3, we only include velocity dispersions measured from 6dFGS, SDSS, and GAMA. Finally, we require the velocity dispersions in the range 60 < σ_e [km s⁻¹] < 450 with estimated measurement uncertainties ε_σ < 0.25σ_e + 25.

The left-hand panel of Figure 1 shows the distribution of stellar masses and redshifts, where the entire sample with reliable redshifts contains 74,594 galaxies. Enforcing the quality controls for the Sérsic and velocity dispersion fits leaves us a parent sample of 32,707 galaxies, shown in orange. Finally, applying the stellar mass and redshift selection, the resulting sample contains 2850 galaxies in VIKING ZYJHK-bands, which is shown in green.

Zoom In Zoom Out Reset image size

Following the arguments of Taylor et al. (2011), Lange et al. (2015), and Baldry et al. (2018), our redshift and mass selection results in a hypothetically volume limited sample.

Due to the significant overlap of red and blue galaxy populations seen in (for example) color-mass, color-Sérsic index diagrams (e.g., Blanton & Moustakas 2009; Taylor et al. 2015, and many others), making an ultimate distinction between early- and late-type galaxies is not a trivial task. A Sérsic cut with n = 2.5 (Shen et al. 2003) can be used where smaller n values correspond to late-type galaxies (LTGs) and greater ones correspond to early-type galaxies (ETGs). A color cut, such as g − i > 0.65 for ETGs, can also be used (Lange et al. 2015).

In this work, we make use of one of the main differences between the two populations, i.e., star formation. While LTGs are star-forming (SF) galaxies, it has diminished for ETGs, thus making them quiescent (Q) galaxies. Star formation can be measured by, e.g., Hα emission lines. Thus, the Hα line will be an absorption feature in the spectra of quiescent galaxies, which means that the equivalent width (EW) will have negative values. As seen in the right-hand panel of Figure 1, Hα EW provides a quite clear bimodality. Following Figure 1 of Taylor et al. (2015), we separate our sample with Hα EW ≥ 1 Å as SF and Hα EW < 1 Å as Q galaxies. Here, we take Hα and Hδ equivalent widths from the DMU SpecLineSFRv05 (Gordon et al. 2017).

We present our final sample in Figure 2, showing the distributions of galaxy properties $r\equiv \mathrm{log}{R}_{e}\,[{h}^{-1}\mathrm{kpc}]$, $s\,\equiv \mathrm{log}{\sigma }_{e}\,[\mathrm{km}\,{{\rm{s}}}^{-1}]$, $i\equiv \mathrm{log}\langle {I}_{e}\rangle [{L}_{\odot }/{\mathrm{pc}}^{2}]$, ${m}^{* }\equiv \mathrm{log}{M}_{\star }[{M}_{\odot }]$, ${\ell }\,\equiv \mathrm{log}L[{L}_{\odot }]$, c = (g − i)_rest rest-frame color, $\nu \equiv \mathrm{log}\,n$ and ${m}^{d}\equiv \mathrm{log}{M}_{\mathrm{dyn}}[{M}_{\odot }]$ among Q and SF galaxies.

Zoom In Zoom Out Reset image size

While M_dyn measurements obtained through simple estimators are usually interpreted as referring to the total mass of the disk and bulge components of a galaxy (see the larger dark matter halo), there is a subtlety depending on how k is calibrated. For example, Ciotti & Lanzoni (1997) prescription oriented k(n) of Bertin et al. (2002) (Equation (2)) strictly applies to a purely stellar bulge; whereas, the Cappellari et al. (2006) calibration (Equation (3)) is derived from projecting Sérsic profiles of 25 E/S0 galaxies to 1R_e. In the latter case, any non-stellar mass component (e.g., molecular or atomic gas) will contribute directly to M_dyn if these components follow a similar density profile. However, if the density profile is different (e.g., the dark matter halo), then the contribution from non-stellar mass changes depending on the amount and the distribution of that mass.

To put this in another way, the dynamical mass estimator in Equation (1) is model dependent, thus it is appropriate to include a correction factor, f_dyn, to account for possible random or systematic errors stemming from the shortcomings of the model. Defining a similar factor for M_⋆ estimates as f_⋆, we have

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\star }}{{M}_{\mathrm{dyn}}}=\displaystyle \frac{{f}_{\star }{M}_{\star }}{{f}_{\mathrm{dyn}}{M}_{\mathrm{dyn}}}=\left({f}_{\star }\displaystyle \frac{{M}_{\star }}{L}\right){\left({f}_{\mathrm{dyn}}k\right)}^{-1}\displaystyle \frac{{GL}}{{\sigma }_{e}^{2}{R}_{e}}\,.\end{eqnarray} \tag{ 5 }$$

Here, the correction factors f_⋆ and f_dyn contain a number of contributing factors, such as corrections to the IMF, dark matter fraction, the contribution of atomic/molecular gas, and the modeling uncertainties in M_⋆/L and k. In general, unambiguously disentangling all possible contributions that are largely if not completely degenerate is an arduous task.

Characterizing the correspondence between simple estimates of stellar and dynamical masses will also lead us to explore potential systematic errors in M_⋆ estimates by studying the residual trends in the stellar-to-dynamical mass ratio (M_⋆/M_dyn) as a function of stellar population parameters, a similar approach adopted by (for example) Taylor et al. (2010), Esdaile et al. (2021), and de Graaff et al. (2021).

For example, if there is some overall deficiency in the SPS modeling, then we might expect to see it in the form of systematic variations in M_⋆/M_dyn as a function of some stellar population parameter, such as (rest-frame) color. In such a scenario, this systematic variation with color could be interpreted as a measure of the correction factor, f_⋆, and might even be considered in terms of a varying IMF. Similarly, if we see a systematic trend with the Sérsic index (n), which is a proxy for structural differences among galaxies, then it might be a sign of structure dependent deficiency in the modeling. Furthermore, we might expect to see the stellar-to-halo mass relation (M_⋆/M_halo) manifesting itself as a trend in M_⋆/M_dyn, which might further be used to derive a refined prescription for k(n). These being said, we have to point out that we will not be able to break the fundamental degeneracy: an observed correlation between M_⋆/M_dyn with, e.g., color, Sérsic index or M_⋆, might imply differential errors in M_⋆ through a varying IMF, or real variation in M_⋆/M_dyn due to variations in gas content and/or dispersion in the M_⋆/M_halo relation (i.e., the dark matter fraction), or any combination thereof.

Considering that there are strong correlations among many galaxy properties, the principal concern here is to uniquely and robustly identify trends with particular observable properties. Many previous studies have adopted an approach in which they studied the residuals from a simple 2D fit to the M_⋆–M_dyn relation (i.e., fitting M_⋆ as a function of M_dyn or vice versa) plotted against other properties, hoping to identify any leading trends or sources of error.

In this work, we instead set out to obtain a complete statistical description of the full high-dimensional distribution of galaxy properties, so that it makes it possible to simultaneously and explicitly account for the elaborate web of real correlations among these properties. This approach enables us to uniquely quantify how much of the random and systematic errors in M_⋆ and/or M_dyn can be directly attributed to different properties.

More specifically, we use a multivariate Gaussian model to describe the observed galaxy distribution in M-dimensional parameter space. ⁸ Saglia et al. (2001) and Colless et al. (2001a) were the first to use a 3D Gaussian model to characterize the Fundamental Plane (FP) of early-type galaxies, showing how the non-negligible (covariant) measurement errors, explicit selection cuts, and large spread of selection weights (implicit selection cuts, e.g., magnitude limits) can be elegantly accounted for within a Gaussian framework. Magoulas et al. (2012) have also shown in their comprehensive work that a 3D Gaussian maximum likelihood (ML) fitting is a statistically rigorous and unbiased method that provides a remarkably successful description of the FP, and they also demonstrated how this scheme is superior to simple regression methods. This approach has become the state of the art for FP studies, which have otherwise been troubled by systematic uncertainties stemming from the choice of fitting formalism/procedure. In this work, we generalize this 3D Gaussian ML method to higher dimensions so that it can simultaneously fit more galaxy parameters and adapt it to a Bayesian framework.

The choice of Gaussian to represent the distribution of galaxies in the observable parameter space is wrong in detail, although that is not necessarily the point—the model does not have to be right to be useful. Even so, we still see (and show in Figure 2) that a Gaussian model does indeed provide a reasonable description of the data, especially of the bivariate linear correlations that are the focus of our work here.

Our motivation for adopting a Gaussian model is completely pragmatic. For the purposes of our study, the most important advantage that comes forward is that once we have the Gaussian description of the full parameter space, it is straightforward to derive any linear relation between the galaxy properties of interest. This advantage stems partly from the facts that (1) any slice through a multivariate Gaussian distribution is itself Gaussian, and (2) any linear combination of independent Gaussian random variables is itself Gaussian. The other aspect of this is that commonly used linear regression methods such as Ordinary Least Squares (OLS) and Orthogonal Distance Regression (ODR) can be expressed in terms of the Gaussian parameters that define the covariance matrix. In other words, our Gaussian model can be seen as just a tool to design a powerful and convenient way to package all of the possible linear fits within the M-dimensional parameter space that we consider.

In Appendix B, we lay out the basic mathematical foundations of our Gaussian model in a Bayesian framework, which includes the details on: handling the error propagation while accounting for the possible covariances among them (Appendix B.1), accounting for the selection effects (Appendix B.2), and implementing a scheme for the identification and rejection of outliers (Appendix B.3). We then describe how we make the modeling problem tractable through using the No U-Turn Sampler (NUTS, Hoffman & Gelman 2014) algorithm implemented in STAN (Carepenter et al. 2017), which makes Markov Chain Monte Carlo (MCMC) sampling some orders of magnitudes faster (Appendix B.4). Finally, we describe the relation between our Bayesian approach to Gaussian forward modeling and the standard frequentist approach to regression, i.e., OLS and ODR (Appendix B.5).

Before moving on, we first present both the data and modeling results in Figure 2. Diagonal panels show the 1D marginal distributions of each parameter in the data set: histograms show the data as observed, and the smooth curves show the Gaussian model description. The off-diagonal panels show the 2D projections of the 8D parameter space, where the upper and lower triangular parts are for SFs and Qs, respectively. Within each of these panels, the points show the data as observed, with median error ellipses given in the top right-hand corner. The large colored ellipses trace the 3σ points of the underlying Gaussian model for each 2D projection of the 8D parameter space. Then, the lines show the forward and inverse linear correlations between each pair of parameters, as derived from the Gaussian fit parameters (see Appendix B.5).

This figure can also be read as a graphical table of our fitting results. In the diagonals, we give the means (μ) and standard deviations (σ). In the off-diagonals, we give the correlation coefficients (ρ). These values give a complete description of the 8D parameter space of galaxies, from which we can derive any other linear correlation.

4.2.1. Linear Regression and Linear Combinations

As previously discussed many times in the literature (e.g., Bernardi et al. 2003; Magoulas et al. 2012), the linear regression coefficients can significantly change with respect to the adopted regression method. OLS regression minimizes the residuals in one of the variables, which makes it useful for predicting one parameter from another (or from a combination of other parameters). The forward fit is the linear relation that best predicts y as a function of x: ${M}_{\star }\propto {M}_{\mathrm{dyn}}^{\alpha };$ and the inverse fit is x as a function of y: ${M}_{\mathrm{dyn}}\propto {M}_{\star }^{\beta }$.

Meanwhile, ODR provides the description of the underlying relation by minimizing the residuals in the orthogonal direction, which is astrophysically the most interesting. This way, we can have insights on the structural differences between the galaxy populations that may cause the difference in these mass estimates.

The relation between M_⋆ and M_dyn can be decomposed in terms of observable galaxy properties as

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\star }}{{M}_{\mathrm{dyn}}}=\displaystyle \frac{2\pi G}{k(n)}\left(\displaystyle \frac{{M}_{\star }}{{L}_{i}}\right)\langle {I}_{e}\rangle {R}_{e}{\sigma }_{e}^{-2}.\end{eqnarray} \tag{ 6 }$$

Since $\mathrm{log}\,{M}_{\star }/{L}_{i}\propto {(g-i)}_{\mathrm{rest}}$ and, as a first order approximation,

$$\begin{eqnarray}\mathrm{log}k(n)\propto \left\{\begin{array}{ll}\mathrm{constant}, & n\lesssim 1.4\\ -\mathrm{log}n, & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 7 }$$

we can study the relation as, e.g., ( m ^*, m ^d ∣ r , s , i , ν , c ). However, considering that our M_dyn estimates are obtained through Equation (1), the inclusion of r , s , and i at the same time will result in over-fitting of this relation. Moreover, under the assumption of homology, we get ${{\boldsymbol{m}}}^{{\boldsymbol{d}}}={\boldsymbol{r}}+2s+\mathrm{constant}$, which causes linear dependence leading to a singular matrix. Therefore, we study a 5D relation by taking a subset of the 8D parameter space as Y = ( m ^*, m ^d , s , ν , c ). This 5D relation can now be evaluated from the corresponding model ${\boldsymbol{Y}}\sim { \mathcal N }({{\boldsymbol{\mu }}}_{{\boldsymbol{Y}}},{{\boldsymbol{\Sigma }}}_{{\boldsymbol{Y}}})$, where the covariance matrix ${{\boldsymbol{\Sigma }}}_{{\boldsymbol{Y}}}={{\boldsymbol{J}}}_{{\boldsymbol{Y}}}{\boldsymbol{\Sigma }}{{\boldsymbol{J}}}_{{\boldsymbol{Y}}}^{\top }$ and the mean vector μ _Y = J _{Y μ} are obtained from the parent model via

$$\begin{eqnarray}{{\boldsymbol{J}}}_{{\boldsymbol{Y}}}=\begin{array}{c}\begin{array}{cccccccc}r & s & i & {m}^{* } & {\ell } & c & \nu & {m}^{d}\end{array}\\ \begin{array}{c}{m}^{* }\\ {m}^{d}\\ s\\ \nu \\ c\end{array}\left(\begin{array}{ccccccccc} & & & & 1 & & & & \\ & & & & & & & & 1\\ & & 1 & & & & & & \\ & & & & & & & 1 & \\ & & & & & & 1 & & \end{array}\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where the omitted elements are zero.

4.2.2. Marginal versus Conditional Relations

The subspace Y can be partitioned as Y = ( Y _a , Y _b ). Similarly, the corresponding mean vector and covariance matrix can also be partitioned in block form as

$$\begin{eqnarray}{{\boldsymbol{\mu }}}_{{\boldsymbol{Y}}}=({{\boldsymbol{\mu }}}_{{\boldsymbol{a}}},{{\boldsymbol{\mu }}}_{{\boldsymbol{b}}})\quad \mathrm{and}\quad {{\boldsymbol{\Sigma }}}_{{\boldsymbol{Y}}}=\left(\begin{array}{cc}{{\boldsymbol{\Sigma }}}_{{\boldsymbol{a}}{\boldsymbol{a}}} & {{\boldsymbol{\Sigma }}}_{{\boldsymbol{ab}}}\\ {{\boldsymbol{\Sigma }}}_{{\boldsymbol{ba}}} & {{\boldsymbol{\Sigma }}}_{{\boldsymbol{b}}{\boldsymbol{b}}}\end{array}\right)\end{eqnarray} \tag{ 9 }$$

where, ${{\boldsymbol{\Sigma }}}_{{\boldsymbol{ba}}}={{\boldsymbol{\Sigma }}}_{{\boldsymbol{ab}}}^{\top }$. The marginal distribution of Y _a is ${{\boldsymbol{Y}}}_{{\boldsymbol{a}}}\sim { \mathcal N }({{\boldsymbol{\mu }}}_{{\boldsymbol{a}}},{{\boldsymbol{\Sigma }}}_{{\boldsymbol{a}}{\boldsymbol{a}}})$ and the conditional distribution of Y _a given Y _b = y ₀ becomes ${{\boldsymbol{Y}}}_{{\boldsymbol{a}}}| {{\boldsymbol{Y}}}_{{\boldsymbol{b}}}\sim { \mathcal N }(\tilde{{\boldsymbol{\mu }}},\tilde{{\boldsymbol{\Sigma }}})$ where,

$$\begin{eqnarray}\begin{array}{rcl}\tilde{{\boldsymbol{\mu }}} & = & {{\boldsymbol{\mu }}}_{{\boldsymbol{a}}}+{{\boldsymbol{\Sigma }}}_{{\boldsymbol{ab}}{{\rm{\Sigma }}}_{{\boldsymbol{bb}}}^{-1}}\left({{\boldsymbol{y}}}_{{\bf{0}}}-{{\boldsymbol{\mu }}}_{{\boldsymbol{b}}}\right)\\ \tilde{{\boldsymbol{\Sigma }}} & = & {{\boldsymbol{\Sigma }}}_{{\boldsymbol{a}}{\boldsymbol{a}}}-{{\boldsymbol{\Sigma }}}_{{\boldsymbol{ab}}}{{\boldsymbol{\Sigma }}}_{{\boldsymbol{b}}{\boldsymbol{b}}}^{-1}{{\boldsymbol{\Sigma }}}_{{\boldsymbol{ab}}}^{\top }.\end{array}\end{eqnarray} \tag{ 10 }$$

It is crucial to emphasize here that fits from conditional distributions are more suitable to examine the dependence of M_⋆–M_dyn relation on other parameters because, as per the name, everything else is marginalized out in marginal distributions. In that case, it is just the projection of the 8D parent model on the 2D M_⋆–M_dyn relation (in other words, fitting the M_⋆–M_dyn pair as a function of one another), in which strong correlations with other observable galaxy properties (e.g., M_⋆ − R_e , M_⋆ − σ_e , R_e − 〈I_e 〉) are indistinguishably included. Using conditional distributions provides an elegant way to perform partial linear regression (e.g., Oh et al. 2022; Barsanti et al. 2022, and references therein), which makes it possible to find the true correlation (or, as we call it, the intrinsic relation) between two properties, while isolating the cross-correlations from other properties.

We now examine the M_⋆–M_dyn relation separately for each galaxy population by inferring the coefficients, slope, and intercept from the conditional and marginal distributions, taking Y _a = ( m ^*, m ^d ) and Y _b = ( s , ν , c ). Here, using the relevant covariance matrices, we calculate both the coefficients that minimize the sum of the orthogonal distances (ODR) and the coefficients that minimize the residuals in the y − direction (parallel or direct fit: OLS). We denote the ODR and OLS slopes of the forward (inverse) relation with α^⊥ and α^∥ (β^⊥ and β^∥), respectively. The intrinsic scatters are calculated via defining a projection vector from the relevant slopes (α^⊥, α^∥, β^⊥, β^∥), as P = (1, − slope). This way, we calculate the parallel intrinsic scatter from ${\sigma }^{\parallel }=\sqrt{P{{\rm{\Sigma }}}^{{\prime} }{P}^{\top }}$, where ${{\boldsymbol{\Sigma }}}^{{\prime} }$ is the relevant covariance matrix; ${{\boldsymbol{\Sigma }}}^{{\prime} }=\tilde{{\boldsymbol{\Sigma }}}$ or ${{\boldsymbol{\Sigma }}}^{{\prime} }={{\boldsymbol{\Sigma }}}_{{\boldsymbol{a}}{\boldsymbol{a}}}$. Then, orthogonal intrinsic scatter can be obtained with, e.g., ${\sigma }^{\perp }\,=\sqrt{P{{\rm{\Sigma }}}^{{\prime} }{P}^{\top }}/\sqrt{1+{{\alpha }^{\perp }}^{2}}$ where P = (1, − α^⊥).

In principle, in the absence of gross systematic errors in the measurements, M_⋆ should strictly be less than M_dyn. While this could be enforced by some choice of prior, to fairly test this expectation, we choose not to implement this prior. For instance, Taylor et al. (2010) have seen that simple dynamical masses calculated under homology assumption (k = constant) have resulted in M_⋆ > M_dyn for M_⋆ < 10^10.5 M_⊙. However, this inconsistency has been removed when dynamical masses were calculated with k(n), which suggests a systematic error in the values of M_dyn estimated using a constant k. Alternatively, we can use the expected (soft) limiting constraint due to (for example) the choice of IMF. As a simple example, adding 0.3 dex to the M_⋆ estimates to approximate the change of IMF from a Chabrier (2003) to a Salpeter (1955) IMF violate the M_⋆ < M_dyn expectation. Because we see that our fits in any case satisfy this expectation, we do not see the need to enforce a M_⋆ > M_dyn prior. (But see 6.3 for some additional nuance around this point.)

In Figure 3, we compare the values of $\mathrm{log}{M}_{\star }$ and $\mathrm{log}{M}_{\mathrm{dyn}}$ for both types of galaxies where the red/blue dots show the data for Qs/SFs. Here we give the inverse relation, M_dyn, as a function of M_⋆. In each panel, the large circles show the median and their associated error bars show the 16/84 percentiles of M_dyn in bins of M_⋆. The solid ellipses in each panel show the 3σ Gaussian ellipse derived from the marginal in orange and for the conditional distribution in green. Dashed and dotted–dashed lines show the best-fitting ODR and OLS relations, respectively, for both distributions. The black lines show the one-to-one relation. Finally, in Table 1, we give the values of the slopes, intrinsic scatters, and the correlation coefficients derived from both conditional and marginal distributions, for both forward and inverse fits, and for both galaxy populations.

Zoom In Zoom Out Reset image size

Table 1. The Relation Between M_⋆ and M_dyn for Both Galaxy Populations Inspected in Different Cases: Forward $({M}_{\star }\propto {M}_{\mathrm{dyn}}^{\alpha })$ and Inverse $({M}_{\mathrm{dyn}}\propto {M}_{\star }^{\beta })$ Fits Derived from Conditional and Marginal Distributions

		${M}_{\star }\propto {M}_{\mathrm{dyn}}^{\alpha }$			${M}_{\mathrm{dyn}}\propto {M}_{\star }^{\beta }$
		Quiescent	Star-forming		Quiescent	Star-forming
Conditional	α⊥	1.008 ± 0.009	1.358 ± 0.027	β⊥	0.992 ± 0.009	0.737 ± 0.015
	α∥	0.877 ± 0.007	0.839 ± 0.016	β∥	0.865 ± 0.007	0.551 ± 0.011
	ρ	0.871 ± 0.003	0.68 ± 0.009	⋯	⋯	⋯
	${\sigma }_{\alpha }^{\perp }$	0.077 ± 0.001	0.101 ± 0.001	${\sigma }_{\beta }^{\perp }$	0.077 ± 0.001	0.101 ± 0.001
	${\sigma }_{\alpha }^{\parallel }$	0.105 ± 0.001	0.147 ± 0.002	${\sigma }_{\beta }^{\parallel }$	0.105 ± 0.001	0.119 ± 0.001
Marginal	α⊥	0.849 ± 0.006	0.889 ± 0.013	β⊥	1.178 ± 0.008	1.125 ± 0.016
	α∥	0.794 ± 0.005	0.748 ± 0.011	β∥	1.074 ± 0.007	0.909 ± 0.012
	ρ	0.923 ± 0.002	0.825 ± 0.006	⋯	⋯	⋯
	${\sigma }_{\alpha }^{\perp }$	0.109 ± 0.001	0.147 ± 0.002	${\sigma }_{\beta }^{\perp }$	0.109 ± 0.001	0.147 ± 0.002
	${\sigma }_{\alpha }^{\parallel }$	0.141 ± 0.002	0.189 ± 0.003	${\sigma }_{\beta }^{\parallel }$	0.164 ± 0.001	0.208 ± 0.002

Note. We give both ODR and OLS slopes and the corresponding intrinsic scatters. We also give the correlation coefficients but only for the forward fits because they will be the same for the inverse relations. This table unambiguously shows that, regardless of galaxy type, for both empirical prediction and astrophysical interpretation purposes, the ideal M_⋆–M_dyn relation with high correlation and low intrinsic scatter is achieved in conditional distributions: at fixed velocity dispersion, Sérsic index, and color.

Download table as:  ASCIITypeset image

Figure 3 and Table 1 show that both orthogonal and parallel (or direct) fits obtained from marginal distributions do not differ much between populations—they are generally consistent within their errors. However, the fits obtained from conditional distributions differ considerably. Furthermore, the correlation between M_⋆ and M_dyn is significantly smaller when fitted using conditionals. As expected, this shows that some part of the higher correlation encountered in the marginal fits are driven by the cross-correlations between ( m ^*, m ^d ) and ( s , ν , c ). Notably, tighter correlations exist in all cases for Qs.

Since the orthogonal slopes of forward and inverse relations are reciprocal of each other (α^⊥ = 1/β^⊥, which also leads to ${\sigma }_{\alpha }^{\perp }={\sigma }_{\beta }^{\perp }$), as seen from Table 1, the orthogonal slope of either relation for Qs is close to unity, albeit only when derived from conditional distributions.

The slopes derived from parallel fits (α^∥ and β^∥) to both forward and inverse relations are consistent with each other for Qs. The small parallel scatters for Qs $({\sigma }_{\alpha }^{\parallel }$ and ${\sigma }_{\beta }^{\parallel }=0.105)$ ⁹ show that both M_⋆ and M_dyn can be predicted well from one another. However, in the case of SFs, it is harder to predict M_⋆ from M_dyn because the scatter in the direction of M_⋆ reaches to 0.147 dex.

Our OLS-derived slope of the M_⋆-M_dyn relation (marginal) for Qs is consistent with previous studies of Taylor et al. (2010), Cappellari et al. (2013), and Zahid & Geller (2017), who find a slope close to unity, with the implication that the central dark matter content is constant and/or a small fraction of the stellar content. The same is not true for SFs, where the ODR-derived slope is considerably smaller. Similarly, after performing the same analysis with the assumption of homology (see Appendix C), we reach the same conclusion as Taylor et al. (2010) and Zahid & Geller (2017) in the aspect that this direct proportionality, M_⋆ ∝ M_dyn, is only possible if (radial) non-homology is taken into account, at least for Qs. Interestingly, the orthogonal slopes in Table 3 differ from unity by only ∼3σ, suggesting that the effect of radial non-homology is still weak. However, this does not seem to apply to SFs, even when radial non-homology is accounted for by k = k(n): we find that regardless of homology, the slope of the relation between M_⋆ and M_dyn for given σ_e , n and (g − i)_rest, largely differs from unity (see Tables 1 and 3).

That said, we emphasize that OLS is not the right method to infer the intrinsic, underlying correlation. Instead, we need to look at the ODR, which shows that the slopes of the M_⋆–M_dyn relation for both Qs and SFs are (1) substantially different from unity, and (2) broadly comparable between the two populations: α = 0.849 ± 0.006 and α = 0.889 ± 0.13 for Q and SF galaxies, respectively. These results imply that f_DM varies as a function of mass in the inner parts of both Q and SF galaxies.

Next, we check the effects of using the structure correction factor of Cappellari et al. (2006) in Equation (3), which has been derived empirically, and thus potentially provides a better prior on structural differences as a function of Sérsic index. By replacing m ^d calculated using Equation (2) with m ^d calculated via Equation (3) and then repeating the analysis, we find that the results given in Table 1 do not undergo any significant changes, as seen in Table 4 which is given in Appendix C. However, the choice of k(n) is expected to affect the dependence of M_⋆/M_dyn relation on other galaxy properties.

Of course, we could have just fitted the M_⋆–M_dyn relation from the very beginning. The advantage of this approach is that it allows us to conveniently account for all of the interrelations within the parameter space and it makes it possible to disentangle the trends stemming from other parameters, which is the topic of Section 5.

To address the difficulties outlined in Section 3.1, our approach is to define a stellar-mass proxy $({\hat{M}}_{\star })$. In the light of our results summarized in Figure 4, which shows that ${M}_{\star }/{M}_{\mathrm{dyn}}\sim {M}_{\star }/{\sigma }_{e}^{2}{R}_{e}$ ratio has strong variations with σ_e , n and (g − i)_rest, we can write ${\hat{M}}_{\star }$ as a hyperplane given by,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{\hat{M}}_{\star } & = & {a}_{0}\mathrm{log}({\sigma }_{e}^{2}{R}_{e})+{a}_{1}\mathrm{log}{\sigma }_{e}+{a}_{2}\mathrm{log}n\\ & & +\,{a}_{3}{\left(g-i\right)}_{\mathrm{rest}}+{a}_{4},\end{array}\end{eqnarray} \tag{ 12 }$$

which we then take as an SED-independent predictor for the stellar mass, and we fit for the a_i coefficients.

From this definition of ${\hat{M}}_{\star }$, one way to understand this approach is thinking of it as fitting M_⋆ vs M_dyn with $\mathrm{log}k(n)\sim {a}_{2}\mathrm{log}n$. When compared to Equations (5) and (11), this can also be interpreted as $\mathrm{log}[({f}_{\star }/{f}_{\mathrm{dyn}})k(n)]\sim {a}_{2}[\mathrm{log}({f}_{\star }/{f}_{\mathrm{dyn}})+\mathrm{log}n]$, where the correction factors f_⋆ and f_dyn are eventually absorbed in the zero-point, a₄. In this sense, by fitting for the values of a_i , we can calibrate a prescription for M_dyn as the one that best predicts M_⋆, and then other coefficients can be understood as describing variation in M_⋆/M_dyn as a function of velocity dispersion, structure, and color.

We set the values of the coefficients by extracting the best-fit relation between the SED-derived stellar-mass estimates, M_⋆, and this proxy variable, ${\hat{M}}_{\star }$. Since our interest is in predicting the value of M_⋆, we use the forward linear OLS regression (see, Section 4.2), which is derived by applying the linear transformation ${\boldsymbol{\omega }}:\hat{{\boldsymbol{y}}}\to {{\boldsymbol{Y}}}^{{\prime} }=({{\boldsymbol{m}}}^{* },2{\boldsymbol{s}}+r,{\boldsymbol{s}},{\boldsymbol{\nu }},{\boldsymbol{c}})$ to our parent model, as shown in Appendix B.5, and then calculating the coefficients from the covariance matrix of ${{\boldsymbol{Y}}}^{{\prime} }$. Note that this is equivalent to calculating the coefficient of each x in the right-hand side of Equation (12) from the conditional distribution $({{\boldsymbol{m}}}^{* },{\boldsymbol{x}}\,| \,{{\boldsymbol{Y}}}^{{\prime} }\setminus \{{m}^{* },x\})$, i.e., OLS regression of ( m ^*, x ) at fixed everything else.

Here we take an extra step and fit our model (Equation (B10)) separately to another 8D parameter space constructed by just replacing the stellar masses from Taylor et al. (2011) with the ones from ProSpect of Robotham et al. (2020). The a_i coefficients and the correlation between M_⋆ and ${\hat{M}}_{\star }$ acquired for both data sets are given in Table 2.

To quantify the mismatch between M_⋆ and ${\hat{M}}_{\star }$, we define a mismatch parameter, ${\delta }_{{m}^{* }}$, as the intrinsic scatter of $\mathrm{log}\,{M}_{\star }/{\hat{M}}_{\star }$. As in Section 4.2, this can be calculated via constructing a projection vector from the best-fitting hyperplane coefficients as P = (1, − a₀, − a₁, − a₂, − a₃). The total observed scatter σ_obs can be evaluated as the rms of $\sqrt{{\boldsymbol{P}}({{\boldsymbol{\Sigma }}}_{{m}^{* }}+{{\bf{E}}}_{j,{m}^{* }}){{\boldsymbol{P}}}^{\top }}$ where ${{\boldsymbol{\Sigma }}}_{{m}^{* }}+{{\bf{E}}}_{{\bf{j}},{m}^{* }}$ is the convolution of covariance and error matrices of parameters included in Equation (12) and can be obtained by applying the Jacobian of the transformation ω to the original convolution matrix Σ + E _j . Since observed scatter is the square root of the quadrature sum of intrinsic scatter and scatter due to measurement errors, ${\sigma }_{\mathrm{obs}}^{2}={\sigma }_{\mathrm{int}}^{2}+{\sigma }_{\mathrm{err}}^{2}$, the intrinsic is $\sqrt{{\boldsymbol{P}}{{\boldsymbol{\Sigma }}}_{{m}^{* }}{{\boldsymbol{P}}}^{\top }}$ and scatter from errors can thus be calculated from the rms of $\sqrt{{\boldsymbol{P}}{{\bf{E}}}_{{\bf{j}},{m}^{* }}{{\boldsymbol{P}}}^{\top }}$. The values for σ_obs, σ_int and σ_err are also given in Table 2.

The basic summary of our results for the Q and SF populations are given in Figure 5. In the first panel (Figure 5), we show M_⋆ plotted against the proxy variable ${\hat{M}}_{\star }$ with the observed, intrinsic, and remaining scatters listed. In each panel of Figure 5, we show the isolated trend, Δ_x (y) ≡ (∂y/∂x)x, which, in addition to Equation (10), can also be calculated as the remaining trend of the parameter in the x − axis after those from the other fitted parameters are subtracted from the quantity in the y − axis. As an example, the isolated trend of $\nu \equiv \mathrm{log}n$ is ${{\rm{\Delta }}}_{\nu }(\mathrm{log}{\hat{M}}_{\star })\equiv \mathrm{log}{\hat{M}}_{* }-[{a}_{0}\mathrm{log}({\sigma }_{e}^{2}{R}_{e})+{a}_{1}\mathrm{log}{\sigma }_{e}+{a}_{3}{\left(g-i\right)}_{\mathrm{rest}}+{a}_{4}]={a}_{2}\mathrm{log}n$, which means that the slope of the trend shown in each panel of Figure 5 is just the slope of the relevant parameter in Equation (12). Finally, this figure also gives the correlation coefficients between $\mathrm{log}{\hat{M}}_{\star }$ and each parameter, which are again calculated via $({{\boldsymbol{m}}}^{* },{\boldsymbol{x}}\,| \,{{\boldsymbol{Y}}}^{{\prime} }\setminus \{{m}^{* },x\})$ and then using Equation (10).

Zoom In Zoom Out Reset image size

Figure 5 shows that the major deterministic factor is the dynamical mass estimator; ${\sigma }_{e}^{2}{R}_{e}$, whereas, structure is only a weak factor for SFs and it essentially has no effect for Qs. This is in fact surprising because a lack of dependence of ${M}_{\star }/{\sigma }_{e}^{2}{R}_{e}$ on structure is not expected considering the results for M_⋆/M_dyn given in Figure 4. As for color, it is significant only for Qs despite their narrow range in (g − i), as seen in Figure 2. In addition, for both populations, the mass proxy ${\hat{M}}_{\star }$ does not strongly depend on s , as the ${M}_{\star }/{\sigma }_{e}^{2}{R}_{e}$ does (Figure 4).

Now that we have calibrated an SED-independent predictor of mass which tightly correlates with the SED-derived stellar masses (ρ ∼ 0.9), we can say that almost all of the observed scatter (σ_obs) in this relation at fixed ${\hat{M}}_{\star }$ (i.e., at fixed σ_e , R_e , n and (g − i)_rest) is due to the uncertainties in the SED-derived masses. By analyzing the residuals as a function of stellar population parameters, we can interpret how much mismatch correlates with these parameters and take them as indicative of systematic errors in M_⋆.

The way that we are going to do this is by first plotting the residuals against rest-frame color (g − i)_rest, equivalent widths of Hα and Hδ lines, 4000 Åbreak strength (D_n 4000), dust extinction [E(B − V)], specific star formation rate (sSFR = SFR/M_⋆), luminosity-weighted mean stellar age $[\mathrm{log}\langle {t}_{\star }{\rangle }_{\mathrm{LW}}/\mathrm{Gyr}]$, and metallicity $[\mathrm{log}\,{Z}_{\star }/{Z}_{\odot }]$. Notice that these are the SP parameters that are not included in our multidimensional fit, except for (g − i)_rest. We then proceed to fitting a linear relation to the residuals as $\mathrm{log}\left({M}_{\star }/{\hat{M}}_{\star }\right)\equiv {\rm{\Delta }}\mathrm{log}{M}_{\star }={Ax}+B$ where x represents each SP parameter. Therefore, the dispersion of x across our sample, σ_x , will propagate through the scatter via A σ_x . Following the notion from Taylor et al. (2020), we will name this quantity as implied dispersion, σ_implied,x ≡ A σ_x , and it can be regarded as a lower limit on the contribution from x to the scatter. We use σ_implied and ρ as metrics to identify which parameters explain the most/least variation, assuming that these residual trends are only due to systematic errors in the SED-derived stellar masses.

Even though what follows might seem to be a simple task of fitting lines to the residual trends, providing estimates of σ_implied,x that are as accurate as possible relies on a rigorous fit that also accounts for the covariant uncertainties in the SPS parameters: measurement errors in stellar mass are strongly correlated to those in color, luminosity-weighted age, and dust extinction. ¹⁰ As such, we implement a 2D-Gaussian model in STAN , which is essentially the 2D and simplified version of our 8D model (e.g., M = 2 in Equation (B4)), to fit ${\rm{\Delta }}\mathrm{log}{M}_{\star }$ as a function of each SPS parameter separately. A, σ_x and the correlation (ρ) between the residuals and each x can be obtained from the corresponding covariance matrix resulting from these 2D fits.

The residuals plotted as a function of SP parameters are given in Figure 6, which also shows the fitted lines. It should be noted here that even though some of the trends are in the same direction as the correlated errors, we can safely say that none of these trends are driven by the correlated errors because they have also been accounted for in our 2D-fitting algorithm.

Zoom In Zoom Out Reset image size

The next question is what fraction of the mismatch between these mass estimates can be explained by contributions from the strong trends with SP parameters? Here, we quantify this by the fraction of the implied scatter from an SP parameter to the mismatch parameter, calculated in quadrature as ${f}_{x,\mathrm{int}}={\sigma }_{\mathrm{implied},x}^{2}/{\delta }_{{m}^{* }}^{2}$. The values inferred from our 2D-Gaussian models for A, σ_x , ρ, σ_x,implied and finally f_x,int are given in Appendix D (Table 5).

6.2.1. Potential Systematics for Star-forming Galaxies

6.2.2. Potential Systematics for Quiescent Galaxies

6.2.3. Projection Effects

The projected axis ratio, denoted as q = b/a where a and b are the semimajor and semiminor axes, respectively, has a significant effect on dynamical mass estimates, as shown by van der Wel et al. (2022). Though, our 8D parameter space does not include q, neither does our stellar-mass proxy in Equation (12). Therefore, it is reasonable to expect some residual trends with q, regardless of whether M_dyn in the 8D parameter space is calculated with Equation (4).

To address this concern, we apply the same procedure of fitting the residuals, this time as a function of projected axis ratio, obtained through Sérsic fits, and also redshift. We give the results in Figure 7, the left-hand panel shows that q is a significant source of scatter in M_⋆ for both Q and SF populations, if it is not accounted for. The effect is particularly pronounced for SFs, corresponding to ∼18% of the intrinsic scatter, which is larger than the total scatter that can be associated with SPS parameters. We also plot the $\mathrm{log}k(q)$ from Equation (4) (van der Wel et al. 2022) with the smooth black curve in the left-hand panel, showing a good match to the trend. We then fit the residuals that are $\mathrm{log}k(q)$ subtracted in the middle panel. This correction successfully removes the trend for SFs, but seems to over-correct the Qs and leaves a trend even slightly stronger in the opposite direction. We note that k(q) correction has no effect whatsoever on our results so far. Finally, as seen in the right-hand panel, the residuals do not show notable variations with redshift for either population.

Zoom In Zoom Out Reset image size

After taking into account the significant scatter stemming from projection effects, the initial intrinsic scatters of 0.108 dex and 0.147 dex drop to 0.103 dex and 0.133 dex for Qs and SFs, respectively. From Figure 6 and Table 5, we can see that SP parameters can explain ∼65% of the intrinsic scatter for quiescent galaxies; however, they can account for only ∼13% for star-forming galaxies. Under these circumstances, at a fixed set of SP parameters, the intrinsic dispersion between the SED-M_⋆ and our dynamical estimates ${\hat{M}}_{\star }$ becomes 0.055 dex and 0.122 dex for Qs and SFs, respectively.

These numbers represent the combination of "true" variations and unmodelled sources of error in the measurements, and therefore provide a soft upper limit on the unmodelled errors. The biggest caveat here is whether the errors in M_⋆ due to the IMF are closely correlated to the errors in M_dyn due to the f_DM. In that case, the correlation between IMF and f_DM can cancel one another in M_⋆/M_dyn. Assuming that most of these unmodelled errors stem from the M_⋆ measurements rather than M_dyn, we can conclude that these values put a conservative limit on the ultimate precision in M_⋆ estimates without (for example) explicit IMF modeling for individual galaxies. By this argument, we come to a fundamental limit to precision of stellar-mass estimates (in the face of IMF variations and any other stochastic factors) being better than 0.05 dex and 0.12 dex for Q and SF galaxies, respectively (see also, e.g., Conroy et al. 2010).

Looking again at Figure 3 with these unmodelled sources of random error in mind, we make the following observations. While it is true that the main (i.e., the mean) relation satisfies the expectation that M_⋆ < M_dyn, it is also true that there is a high M_⋆/M_dyn tail to the population that seems to violate this constraint. This is most naturally explained by unmodelled errors in the stellar masses, dynamical masses, or both. We can subtract off our estimates for the unmodelled errors in the stellar masses, in which case the intrinsic scatters become 0.13 and 0.15 dex for Qs and SFs, respectively. This alleviates, but does not eliminate, the apparent violation of the M_⋆ < M_dyn expectation. Alternatively, we could have tried to enforce M_⋆ < M_dyn constraint by truncating the model (assuming that the mass estimates and their uncertainties are robust), and/or globally scale the mass estimates through changes to the IMF or f_DM. At this point, however, we reach a point where we have exhausted what is possible to infer from single-fiber spectroscopy—to go further will require proper dynamical modeling from spatially resolved spectroscopy (but still subject to the IMF/dark matter degeneracy).

We now apply the same analysis to the ProSpect masses. This analysis is particularly valuable considering that there are not-so-negligible contributions to the scatter from metallicity and star formation rates, as seen in Table 5 and Figure 6.

Figure 8 gives the residuals in M_⋆ from ProSpect as a function of SP parameters. Here, sSFR and Z_gas are provided in ProSpect, but dust and age are not. We also have to assume there is no covariance between the uncertainties of the parameters considered here because this information is currently not available in ProSpect data products.

Zoom In Zoom Out Reset image size

When compared to Figure 6, we see that the apparent effect of SFR on scatter has vanished as expected because ProSpect implements more complicated SFHs through time dependent Z_gas. However, we can see that there is a substantial contribution from Z_gas to the mismatch, particularly in SFs, for which potential uncertainties in Z_gas can explain ∼11% of the scatter. Aside from Z_gas and the minor contribution of Hδ, the SP parameters that are considered here do not seem to have any effect on the overall scatter. In this case, the remaining intrinsic scatter becomes 0.112 dex and 0.12 dex for Qs and SFs, respectively. Under the same aforementioned assumptions, this remaining scatter can be attributed to the joint effects of uncertainties in age, dust, and IMF.

Due to measurement errors, each galaxy in the observed $\hat{{{\boldsymbol{y}}}_{j}}$ has an uncertainty ε_j . As mentioned at the beginning of this section, and just like the ML method used in Magoulas et al. (2012), Springob et al. (2014), and Said et al. (2020), the Gaussian approach enables us to easily account for the errors along with their possible covariances. This is possible by incorporating uncertainties into the error matrix for each galaxy and these matrices will not be diagonal due to covariant errors.

For a data set such as $\hat{{\boldsymbol{y}}}=({\boldsymbol{r}},{\boldsymbol{s}},{\boldsymbol{i}},{{\boldsymbol{m}}}^{* },{\boldsymbol{c}},{\boldsymbol{\nu }})$, the error matrix E _j can be constructed as,

$$\begin{eqnarray}{{\boldsymbol{E}}}_{{\boldsymbol{j}}}=\left(\begin{array}{cccccc}{\varepsilon }_{{r}_{j}}^{2} & 0 & {\varepsilon }_{{r}_{j},{i}_{j}} & 0 & 0 & {\varepsilon }_{{r}_{j},{\nu }_{j}}\\ 0 & {\varepsilon }_{{s}_{j}}^{2} & 0 & 0 & 0 & 0\\ {\varepsilon }_{{r}_{j},{i}_{j}} & 0 & {\varepsilon }_{{i}_{j}}^{2} & 0 & 0 & {\varepsilon }_{{i}_{j},{\nu }_{j}}\\ 0 & 0 & 0 & {\varepsilon }_{{m}_{j}^{* }}^{2} & {\varepsilon }_{{m}_{j}^{* },{c}_{j}} & 0\\ 0 & 0 & 0 & {\varepsilon }_{{m}_{j}^{* },{c}_{j}} & {\varepsilon }_{{c}_{j}}^{2} & 0\\ {\varepsilon }_{{r}_{j},{\nu }_{j}} & 0 & {\varepsilon }_{{i}_{j},{\nu }_{j}} & 0 & 0 & {\varepsilon }_{{\nu }_{j}}^{2}\end{array}\right)\end{eqnarray} \tag{ B5 }$$

where ${\varepsilon }_{{x}_{j}}$ is the uncertainty of the observable parameter x for j − th galaxy. The off-diagonal elements represent the covariances between the errors of parameter pairs as ${\varepsilon }_{{x}_{j},{y}_{j}}={\rho }_{{xy}}^{\varepsilon }{\varepsilon }_{{x}_{j}}{\varepsilon }_{{y}_{j}}$ where ${\rho }_{{xy}}^{\varepsilon }$ is the correlation coefficient of the errors in x and y parameters. Calculation of ${\rho }_{{xy}}^{\varepsilon }$ is detailed in Appendix A.

Due to selection effects, the data $\hat{{\boldsymbol{y}}}$ is a truncated Gaussian, thus, the likelihood $p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}})$ needs to be renormalised by,

$$\begin{eqnarray}&&\begin{array}{l}p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}})\to \displaystyle \frac{p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}})}{{\hat{f}}_{j}},\\ {\hat{f}}_{j}={\displaystyle \int }_{{\hat{{\boldsymbol{y}}}}_{\mathrm{cut}}}^{\infty }p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}}){d}^{M}\hat{{{\boldsymbol{y}}}_{j}},\end{array}\end{eqnarray} \tag{ B6 }$$

where ${\hat{f}}_{j}$ is the normalization integral for the multivariate normal distribution and $\hat{{{\boldsymbol{y}}}_{\mathrm{cut}}}$ is the vector corresponding to the lower limits of the parameters included in the $\hat{{\boldsymbol{y}}}$ vector. ${\hat{f}}_{j}$ ensures that the likelihood $\int p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}}){d}^{M}{\hat{{\boldsymbol{y}}}}_{j}=1$ within the domain defined by the data. Due to the instrumental limits, we expect to see explicit cuts in velocity dispersion (or M_⋆, as discussed in Section 2.4), apparent magnitude and redshift. The redshift and magnitude limits can be simultaneously corrected by selection function, S, defined in Magoulas et al. (2012), which is inversely proportional to ${V}_{\max }$ weighting (Schmidt 1968): $S\propto 1/{V}_{\max }$. When the weight w_j = 1/S_j is applied to the renormalised likelihood as ${[p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}})/{\hat{f}}_{j}]}^{{w}_{j}}$, the natural logarithm of Equation (B3) becomes,

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}p({\boldsymbol{\mu }},{\boldsymbol{\Sigma }}\,| \,\hat{{\boldsymbol{y}}})\propto \mathrm{ln}p({\boldsymbol{\mu }},{\boldsymbol{\Sigma }})\\ +\,\displaystyle \sum _{j=1}^{N}{w}_{j}\left[\mathrm{ln}p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}})-\mathrm{ln}{\hat{f}}_{j}\right].\end{array}\end{eqnarray} \tag{ B7 }$$

An important aspect of weighting here is that the w_j weights are small for our sample: w_j = 1 for 2793 galaxies (i.e., no weighting needed at all) out of 2850 and maximum weight encountered is w_j = 2.54. This maximum value of the weights is minuscule, especially considering that Taylor et al. (2011) and Magoulas et al. (2012) have chosen their samples discarding the galaxies with too large weights, w > 30 and w > 20, respectively, to prevent these galaxies from skewing the fits. The reason why we obtain such low weights is the fact that GAMA provides a deep probe of the massive galaxy population with $\mathrm{log}{M}_{\star }\gtrsim 10.5$ out to z ≈ 0.25, as pointed out in Taylor et al. (2011), which is the bulk of our sample.

The log-likelihood in Equation (B7) enables sampling from a truncated and censored multivariate Gaussian, and thus accounts for the selection effects and provides an unbiased fit.

Conventional weighting of data points with the inverse square of their observational uncertainties $({w}_{j}=1/{\varepsilon }_{j}^{2})$ is prone to be sensitive to outliers (e.g., Hogg et al. 2010; Taylor et al. 2015). A proposed solution to this "bad data" problem comes from Gaussian mixture modeling, in which a secondary component is added to the likelihood function. This enables us to model the outliers as some fraction of the data being drawn from a bad distribution, whereas the "good data" that we are interested in is drawn from the good distribution.

We model this secondary component as a Gaussian because we simply do not have any knowledge of the bad distribution. A key point here is that no a priori assumptions are made as to the goodness and/or badness of the data points, so that the model can treat every point as having an equal probability of being bad; thus, it can objectively identify the outliers (Taylor et al. 2015). If the fraction of good data is denoted as f_good, then the net likelihood becomes,

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\mathrm{net}}(\hat{{\boldsymbol{y}}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }},{f}_{\mathrm{good}})=\displaystyle \prod _{j=1}^{N}[{f}_{\mathrm{good}}p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }})\\ +\,(1-{f}_{\mathrm{good}})p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{{\boldsymbol{\Sigma }}}_{\mathrm{bad}})],\end{array}\end{eqnarray} \tag{ B8 }$$

where Σ_bad is the covariance matrix of the bad distribution. It should also be noted that both good and bad distributions here are assumed to have the same mean μ . Equation (B8) shows that each data point $\hat{{{\boldsymbol{y}}}_{j}}$ may have arisen from either good or bad distribution. Inserting Equation (B8) to the posterior in Equation (B7), we obtain the log-likelihood:

$$\begin{eqnarray}&&\begin{array}{l}{lnp}({\boldsymbol{\mu }},{\boldsymbol{\Sigma }},{f}_{\mathrm{good}}| \hat{{\boldsymbol{y}}})\propto \displaystyle \sum _{j=1}^{N}\mathrm{ln}\left[{f}_{\mathrm{good}}{\left(\displaystyle \frac{p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}})}{{\hat{f}}_{j}}\right)}^{{w}_{j}}\right.\\ \ \left.+\,(1-{f}_{\mathrm{good}})p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{{\boldsymbol{\Sigma }}}_{\mathrm{bad}})\right].\end{array}\end{eqnarray} \tag{ B9 }$$

It is noteworthy that, as seen from Equation (B9), f_good is also declared as a parameter. Furthermore, we can quantify the badness of each data point, i.e., the probability of a data point being bad, by f_bad p_bad,j/p_net,j, where f_bad = 1 − f_good. One final remark is that we also assume the bad distribution is not truncated.

Within this study, we typically use MCMC samples from four Markov chains with each chain consisting of 1000 draws from the posterior in Equation (B9), 500 of which is warm up, thus adding up to 2000 iterations in total, after discarding the warm up. Within this Bayesian MCMC framework, it is crucial to implement convoluted statistical models, like the one we present in this work, in an efficient way to avoid computationally expensive practices. To do this, we make use of optimized and robust built-in functions of STAN .

The critical and possibly the most computationally challenging part of applying Equation (B9) is the normalization factor ${\hat{f}}_{j}$ which accounts for the selection cuts in velocity dispersion and stellar mass. We give the derivation of an expression for ${\hat{f}}_{j}$ in Appendix B.6 in terms of the built-in cumulative and probability distribution functions.

Another challenge arises from the Gaussian mixture. As seen from Equation (B9), the log-posterior is the natural logarithm of the sum of two terms, which can be handled by one of the mainstays of statistical computing, log-sum of exponentials, where $\mathrm{ln}(a+b)=\mathrm{ln}[\exp (\mathrm{ln}a)+\exp (\mathrm{ln}b)]=\mathrm{log}-\mathrm{sum}\,-\exp (\mathrm{ln}a,\mathrm{ln}b)$. This way, Equation (B9) becomes

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}p({\boldsymbol{\mu }},{\boldsymbol{\Sigma }},{f}_{\mathrm{good}}| \hat{{\boldsymbol{y}}})\\ \propto \,\displaystyle \sum _{j=1}^{N}\mathrm{log}-\mathrm{sum}-\exp \left(\mathrm{ln}{f}_{\mathrm{good}}+{w}_{j}[\mathrm{ln}p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}})\right.\\ -\left.\ \mathrm{ln}{\hat{f}}_{j}],\mathrm{ln}(1-{f}_{\mathrm{good}})+\mathrm{ln}p(\hat{{{\boldsymbol{y}}}_{j}}| {\boldsymbol{\mu }},{{\boldsymbol{\Sigma }}}_{\mathrm{bad}}\right)\end{array}\end{eqnarray} \tag{ B10 }$$

which makes efficient sampling possible by using STAN 's optimized functions log_sum_exp and log1m $(x)\,=\mathrm{ln}(1-x)$. Equation (B10) is the final form for the posterior of our model from which the MCMC samples will be drawn in STAN.

All information regarding the relations between the main galaxy parameters in $\hat{{\boldsymbol{y}}}$, and other parameters that can be derived from them, are encoded in the covariance matrix, Σ. The relations between parameters of interest can be extracted from Σ using linear transformations,

$$\begin{eqnarray}\begin{array}{rcl}T:{{\mathbb{R}}}^{M} & \longrightarrow & {{\mathbb{R}}}^{K}\\ {\boldsymbol{y}} & \longrightarrow & {\boldsymbol{\omega }}={\boldsymbol{Ay}}+{\boldsymbol{B}},\end{array}\end{eqnarray} \tag{ B11 }$$

where K is the number of parameters of interest. The Jacobian of this transformation, ${\boldsymbol{J}}={\boldsymbol{A}}\in {{\mathbb{R}}}^{K\times M}$, can then be used in,

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\omega }={\boldsymbol{A}}{\boldsymbol{\Sigma }}{{\boldsymbol{A}}}^{\top },\qquad \bar{{\boldsymbol{\omega }}}={\boldsymbol{A}}\bar{{\boldsymbol{y}}}+{\boldsymbol{B}}\end{eqnarray} \tag{ B12 }$$

to find the covariance matrix $({{\boldsymbol{\Sigma }}}_{\omega }\in {{\mathbb{R}}}^{K\times K})$ and the mean vector $(\bar{{\boldsymbol{\omega }}}\in {{\mathbb{R}}}^{K})$ for those parameters. The OLS coefficients for a linear relation between these K parameters can then be calculated from Σ_ω and $\bar{{\boldsymbol{\omega }}}$. Furthermore, ODR coefficients can also be obtained from the eigenvectors of Σ_ω .

In this section, we show the main steps to derive an expression for the integral normalization factor ${\hat{f}}_{j}$, similar to the works of Magoulas et al. (2012) and Dam (2020). For the data set $\hat{{\boldsymbol{y}}}=({\boldsymbol{r}},{\boldsymbol{s}},{\boldsymbol{i}},{{\boldsymbol{m}}}^{* },{\boldsymbol{\ell }},{\boldsymbol{c}},{\boldsymbol{\nu }},{{\boldsymbol{m}}}^{{\boldsymbol{d}}})$, the selection limits in our final sample are m^* ≥ m_cut = 10.3, 60 < σ_e [km s⁻¹] < 450, 0.01 ≤ z ≤ 0.12 and m_r ≤ 19.65, as discussed in Section 2. Since magnitude and redshift limits can be accounted for by selection probability function $S\propto 1/{V}_{\max }$, we can say that there are no cut-off values for any component of $\hat{{\boldsymbol{y}}}$ except for m ^* and s , which means that the integral over all the other components will just be the marginalization of the 8D probability distribution function (PDF). In this case, the likelihood normalization factor f_j in Equation (B10) can be expressed in terms of a bivariate normal PDF,

$$\begin{eqnarray}&&\begin{array}{l}{f}_{j}={\displaystyle \int }_{{m}_{\mathrm{cut}}}^{\infty }{\displaystyle \int }_{{s}_{\mathrm{low}}}^{{s}_{\mathrm{up}}}p(\hat{{{\boldsymbol{y}}}_{j}^{{\prime} }}| {{\boldsymbol{\mu }}}^{{\prime} },{{\boldsymbol{C}}}_{{\boldsymbol{j}}}^{{\prime} }){d}^{2}\hat{{{\boldsymbol{y}}}_{j}^{{\prime} }}\\ ={\displaystyle \int }_{{m}_{\mathrm{cut}}}^{\infty }{\displaystyle \int }_{{s}_{\mathrm{low}}}^{{s}_{\mathrm{up}}}\displaystyle \frac{\exp \left[-\tfrac{1}{2}{\left(\hat{{{\boldsymbol{y}}}_{j}^{{\prime} }}-{{\boldsymbol{\mu }}}^{{\prime} }\right)}^{\top }{\left({{\boldsymbol{C}}}_{{\boldsymbol{j}}}^{{\prime} }\right)}^{-1}(\hat{{{\boldsymbol{y}}}_{j}^{{\prime} }}-{{\boldsymbol{\mu }}}^{{\prime} })\right]}{2\pi \sqrt{| {{\boldsymbol{C}}}_{{\boldsymbol{j}}}^{{\prime} }}| }\\ \quad \times {{ds}}_{j}{{dm}}_{j}\end{array}\end{eqnarray} \tag{ B13 }$$

where $\hat{{{\boldsymbol{y}}}_{j}^{{\prime} }}=({s}_{j},{m}_{j}^{* })$, mean vector ${{\boldsymbol{\mu }}}^{{\prime} }=(\bar{s},{\bar{m}}^{* })$ and ${{\boldsymbol{C}}}_{{\boldsymbol{j}}}^{{\prime} }$ is constructed with the s, m elements of the covariance matrix convolved with the error matrix, which can be done via

$$\begin{eqnarray}{{\boldsymbol{C}}}_{{\boldsymbol{j}}}^{{\prime} }={{\boldsymbol{J}}}^{{\prime} }({\boldsymbol{\Sigma }}+{{\boldsymbol{E}}}_{{\boldsymbol{j}}}){{\boldsymbol{J}}}^{{\prime} \top },\quad {{\boldsymbol{J}}}^{{\prime} }=\left(\begin{array}{cccccccc}0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\end{array}\right).\end{eqnarray} \tag{ B14 }$$

Using the Cholesky decomposition of ${{\boldsymbol{C}}}_{{\boldsymbol{j}}}^{{\prime} }={L}_{j}{L}_{j}^{\top }$, where L_j is a lower triangular matrix, provides an easier way to reduce Equation (B13) to a more applicable form for efficient MCMC sampling in STAN . This way, Equation (B13) can be expressed as,

$$\begin{eqnarray}&&\begin{array}{l}{f}_{j}={\displaystyle \int }_{{u}_{j,\mathrm{cut}}}^{\infty }\phi ({u}_{j})\left[{\rm{\Phi }}({\alpha }_{\mathrm{up}}+\beta {u}_{j})-{\rm{\Phi }}({\alpha }_{\mathrm{low}}+\beta {u}_{j})\right]{{du}}_{j}\\ \mathrm{where},\ {u}_{j}=\displaystyle \frac{{m}_{j}^{* }-{\bar{m}}^{* }}{{L}_{j,11}},\quad {u}_{j,\mathrm{cut}}=\displaystyle \frac{{m}_{\mathrm{cut}}-{\bar{m}}^{* }}{{L}_{j,11}},\\ {\alpha }_{\mathrm{low}}=\displaystyle \frac{{s}_{\mathrm{low}}-\bar{s}}{{L}_{j,22}},\quad {\alpha }_{\mathrm{up}}=\displaystyle \frac{{s}_{\mathrm{up}}-\bar{s}}{{L}_{j,22}}\ \mathrm{and}\ \beta =-\displaystyle \frac{{L}_{j,21}}{{L}_{j,22}}\end{array}\end{eqnarray} \tag{ B15 }$$

Here, Φ is the cumulative distribution function (CDF) and ϕ is the PDF of the standard normal distribution, having the relation,

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}(x)=\displaystyle \frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{x}{e}^{-{t}^{2}/2}{dt}\\ =\,{\int }_{-\infty }^{x}\phi (t){dt}\ \mathrm{and}\ \phi (x)=\displaystyle \frac{d{\rm{\Phi }}(x)}{{dx}}.\end{array}\end{eqnarray} \tag{ B16 }$$

Using the tables provided by Owen (1980),

$$\begin{eqnarray}\begin{array}{rcl}{f}_{j} & = & {\rm{\Phi }}\left(\displaystyle \frac{\alpha }{\sqrt{1+{\beta }^{2}}}\right)\\ & & -\,{\left.\mathrm{BvN}\left(\displaystyle \frac{\alpha }{\sqrt{1+{\beta }^{2}}},{u}_{j,\mathrm{cut}},\displaystyle \frac{-\beta }{\sqrt{1+{\beta }^{2}}}\right)\right|}_{\alpha ={\alpha }_{\mathrm{low}}}^{\alpha ={\alpha }_{\mathrm{up}}},\end{array}\end{eqnarray} \tag{ B17 }$$

where BvN is the CDF of the standard bivariate normal distribution given as

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{BvN}(h,k,\rho )=\displaystyle \frac{1}{2\pi \sqrt{1-{\rho }^{2}}}\\ \quad \times \,{\displaystyle \int }_{-\infty }^{k}{\displaystyle \int }_{-\infty }^{h}\exp \left[-\displaystyle \frac{1}{2}\left(\displaystyle \frac{{x}^{2}-2\rho {xy}+{y}^{2}}{1-{\rho }^{2}}\right)\right]{dxdy}.\end{array}\end{eqnarray} \tag{ B18 }$$

We use the algorithm of Boys (1989) to evaluate BvN, also making use of STAN 's optimized built-in functions for Φ and ϕ.