The EDGE-CALIFA Survey: Exploring the Star Formation Law through Variable Selection

TL;DR

This study uses multilinear analysis with LASSO to identify key predictors of star formation in galaxies, revealing the influence of molecular gas, stellar densities, metallicities, and galaxy mass on star formation rates.

Contribution

It introduces a novel application of LASSO to galaxy data, identifying significant factors influencing star formation that were previously unquantified.

Findings

Higher molecular and stellar surface densities increase star formation rates.

Older stellar populations and higher metallicities correlate with lower star formation.

Star formation rate increases with galactocentric radius normalized by disk scale length.

Abstract

We present a multilinear analysis to determine the significant predictors of star formation in galaxies using the combined EDGE-CALIFA sample of galaxies. We analyze 1845 kpc-scale lines of sight across 39 galaxies with molecular line emission measurements from EDGE combined with optical IFU data drawn from CALIFA. We use the Least Absolute Shrinkage and Selection Operator (LASSO) to identify significant factors in predicting star formation rates. We find that the local star formation rate surface density is increased by higher molecular gas surface densities and stellar surface densities. In contrast, we see lower star formation rates in systems with older stellar populations, higher gas- and stellar-phase metallicities and larger galaxy masses. We also find a significant increase in star formation rate with galactocentric radius normalized by the disk scale length, which suggests additional parameters regulating star formation rate not explored in this study.

Tables

Table 1. Table 1: List of variables used for the multi-dimensional model

Property	Type1	Units (before transforming)	Description
$\log \left({\Sigma }_{\mathrm{SFR}}\right)$	C, R	$M_{\odot }{\mathrm{pc}}^{-2}{\mathrm{yr}}^{-1}$	Log of star formation rate surface density
$\log \left({\Sigma }_{\mathrm{mol}}\right)$	C, R	$M_{\odot }{\mathrm{pc}}^{-2}$	Log of molecular gas surface density
$\log \left({\sigma }_{\mathrm{CO}}\right)$	C, R	${\text{km s}}^{-1}$	Log of molecular gas velocity dispersion
$\log \left(V_{\mathrm{rot}}\right)$	C, G	$\mathrm{km}{s}^{-1}$	Maximum rotation speed estimated from molecular gas
$12+\log \left(O/H\right)$	C, R	$\cdots$	Gas-phase metallicities from Marino et al.(2013)
$\log \left({\Sigma }_{\star }\right)$	C, R	$M_{\odot }{\mathrm{pc}}^{-2}$	Log of stellar surface density
$\log \left({\sigma }_{\star }\right)$	C, R	${\text{km s}}^{-1}$	Log of stellar velocity dispersion
$Z_{\star }$	C, R	$\log \left(Z/Z_{\odot }\right)$	Stellar metallicity relative to solar
$M_{\mathrm{gal}}$	C, G	$M_{\odot }$	Stellar mass of galaxy
$\log \left({\tau }_{\star ,m}\right)$	C, R	yr	Mass-weighted stellar age from SSP analysis
$\log \left({\tau }_{\star ,l}\right)$	C, R	yr	Luminosity-weighted stellar age from SSP analysis
$A_V$	C, R	mag	Implied $V$-band extinction estimated from across CALIFA bandpass
$\log \left(R/R_s\right)$	C, R	kpc	Distance from the galaxy centre normalized by the scale length of the disk
$T$	D, G	$\cdots$	Numerical Hubble type
Bar	D, G	$\cdots$	Galaxy is classified as having a bar
Ring	D, G	$\cdots$	Galaxy is classified as having a ring
Multiple	D, G	$\cdots$	Galaxy is classified as Multiple
Centre	D, R	$\cdots$	Line-of-sight in centre of galaxy
1 C = continuous, D = discrete, R = resolved, G = global

Table 2. Table 2: Variance inflation factors (VIFs) and LASSO regression coefficients for different factors contributing to the star formation law. The VIFs indicate how correlated a given variable is with the other independent variables in the data set, with VIF > 5 absent 5 >5 taken as a signature of significant correlation. The variable C i subscript 𝐶 𝑖 C_{i} is the coefficient on the multilinear model and the uncertainties are established by resampling the data with the uncertainties and bootstrapping the sample. Values of 0.00 indicate where the LASSO method has removed the variable from the model.

Parameters	$\log \left({\Sigma }_{\mathrm{SFR}}\right)$			$\log \left({\tau }_{\mathrm{dep}}\right)$
	VIF	$C_i$		VIF	$C_i$
$\log \left({\Sigma }_{\mathrm{mol}}\right)$	3.1	$0.43\pm 0.03$		$\cdots$	$\cdots$
$\log \left({\sigma }_{\mathrm{CO}}\right)$	2.2	$-0.21\pm 0.06$		2.0	$0.5\pm 0.1$
$\log \left(V_{\mathrm{rot}}\right)$	4.1	$-0.28\pm 0.08$		4.1	$0.3\pm 0.1$
$12+\log \left(O/H\right)$	1.8	$-2.7\pm 0.3$		1.8	$2.6\pm 0.4$
$\log \left({\Sigma }_{\star }\right)$	5.5	$0.70\pm 0.03$		3.0	$-0.22\pm 0.03$
$\log \left({\sigma }_{\star }\right)$	2.9	$0.19\pm 0.02$		2.7	$-0.33\pm 0.03$
$\log \left(Z_{\star }/Z_{\odot }\right)$	2.6	$-0.46\pm 0.08$		2.6	$0.5\pm 0.1$
$\log \left(M_{\mathrm{gal}}\right)$	2.8	$-0.17\pm 0.02$		2.8	$0.13\pm 0.03$
$\log \left({\tau }_{\star ,m}\right)$	1.8	$-0.56\pm 0.07$		1.7	$0.4\pm 0.1$
$\log \left({\tau }_{\star ,l}\right)$	2.5	$-0.12\pm 0.04$		2.4	$0.00\pm 0.04$
$A_V$	1.8	$0.00\pm 0.04$		1.8	$0.00\pm 0.04$
$\log \left(R/R_s\right)$	2.8	$0.58\pm 0.07$		2.8	$-0.5\pm 0.1$
$T$	1.6	$0.030\pm 0.004$		1.7	$-0.04\pm 0.02$
Bar	2.2	$-0.038\pm 0.005$		2.2	$0.037\pm 0.005$
Ring	1.7	$0.00\pm 0.01$		1.7	$0.00\pm 0.01$
Multiple	1.1	$0.00\pm 0.03$		1.1	$-0.08\pm 0.03$
Centre	1.4	$0.11\pm 0.03$		1.4	$-0.09\pm 0.04$

Equations

$$\frac{\mathrm{SFR}}{M_{\odot }\text{mbox}{yr}^{-1}}=7.9\times 10^{-42}\left(\frac{L_{\mathrm{H}\alpha }}{\text{mbox}{ergs}^{-1}}\right)10^{A_V/2.5}.$$

$${\Sigma }_{\mathrm{SFR}}=A{\Sigma }_{\mathrm{mol}}^N;$$

$$\log ({\Sigma }_{\mathrm{SFR}})=\log A+N\log ({\Sigma }_{\mathrm{mol}}).$$

$${\tau }_{\mathrm{dep}}\equiv \frac{{\Sigma }_{\mathrm{mol}}}{{\Sigma }_{\mathrm{SFR}}}.$$

$$Y_j=C_0+\sum _{i=1}^p{C}_i{X}_{ij}.$$

$$S=\sum _{j=1}^n{\left[Y_j-\left(C_0+\sum _{i=1}^p{C}_i{X}_{ij}\right)\right]}^2.$$

$$\sum _{j=1}^n{\left[Y_j-(C_0+\sum _{i=1}^p{C}_i{X}_{ij})\right]}^2+\alpha \sum _{i=1}^p|C_i|=S+\alpha \sum _{i=1}^p|C_i|$$

$${\Sigma }_{\mathrm{SFR}}=10^{C_0}{\Sigma }_{\mathrm{mol}}^a{\Sigma }_{\star }^b,$$

$$\Delta Y_j=Y_j-C_0-\sum _{i\ne k}{C}_i\left[X_{ij}-⟨X_{ij}⟩\right].$$

$$\Delta \log {\Sigma }_{\mathrm{SFR}}=\log {\Sigma }_{\mathrm{SFR}}-C_0-0.30(\log {\Sigma }_{\star }-⟨\log {\Sigma }_{\star }⟩)-0.22B,$$

$${\Sigma }_{\mathrm{SFR}}\propto {\Sigma }_{\mathrm{mol}}^{0.4}{\Sigma }_{\star }^{0.7}{\left(\frac{R}{R_s}\right)}^{0.6}{\tau }_{\star ,m}^{-0.6}{\sigma }_{\star }^{0.2}{Z}_{\mathrm{gas}}^{-2.7},$$

$$\widetilde{C_i}=C_i\frac{{\sigma }_{X_i}}{{\sigma }_Y},$$

$${\tau }_{\mathrm{dep}}\propto {\sigma }_{\star }^{-0.3}{\sigma }_{\mathrm{CO}}^{0.5}{T}^{-0.04}{\Sigma }_{\star }^{-0.2}{Z}_{\mathrm{gas}}^{2.6}{\left(\frac{R}{R_s}\right)}^{-0.5}$$

Full text

1 Introduction

Star formation is one of the dominant secular evolution processes shaping galaxies over cosmic time. The focus on how gas is converted into stars has guided the exploration of both nearby and distant galaxies (Kennicutt & Evans, 2012; Madau & Dickinson, 2014). Historically, the study of the “star formation law” has been through the conjecture by Schmidt (1959) that the star formation rate is a function of the volume density of gas, establishing the local star formation rate through the free-fall time. Empirically, the work of Kennicutt (1989, 1998b) initiated the study of the star formation law in terms of the surface densities of (neutral) gas and the star formation rate, which represent the observationally tractable quantities in galaxies. Since local observations demonstrated that only the molecular phase of the interstellar medium (ISM) is associated with star formation, the star formation law studies began focusing on how the molecular gas content drives star formation (Wong & Blitz, 2002; Kennicutt et al., 2007). While these works found clear empirical relationships between star formation rate and gas content, though the extremes of the star formation environment, in low-metallicity systems and (U)LIRGs, show significantly different behaviours indicating that relationship defined in galaxy disks is not universal (Daddi et al., 2010).

Using a wealth of of multiwavelength data on nearby galaxies, it has become possible to explore the driving physics behind the star formation law, making quantitative tests of theories (Leroy et al., 2008; Leroy et al., 2013; Saintonge et al., 2011a, b). In particular, two drivers of the star formation rate have appeared: the molecular gas content of the system (e.g., Leroy et al., 2013) and the depth of the stellar potential (e.g., Saintonge et al., 2011b). The contrast between these two studies is notable since they identify separate drivers for the star formation law, but they also employ different methods. In particular, studies that find the local molecular gas content within galaxies drives the star formation rate typically favour line-of-sight based studies where kpc-scale regions of disks are being studied. The studies favouring stellar content as the driving factor are usually drawn from studies of a wide range of galaxy types that focus on the global properties of galaxies.

Since star formation could be driven by multiple physical processes, several studies opt to consider several physical parameters simultaneously (Dopita & Ryder, 1994; Shi et al., 2011; Leroy et al., 2013; Shi et al., 2018). From these works, the focus on a pressure-based formalism (e.g. Ostriker et al., 2010) provides a means to synthesize the stellar potential and gas content contributions to the star formation law. Other models such as dynamical regulation (Silk, 1997; Elmegreen, 1997) have also been explored empirically (Kennicutt, 1998b; Leroy et al., 2013). Recently Colombo et al. (2018) also found that dynamical models could explain the empirical star formation law but that Hubble morphological type was a clearer factor in determining how the galaxies form stars. A similar conclusion was reached based on optical data alone (González Delgado et al., 2014, 2015, 2016).

Thanks to dedicated campaigns, the number of multi-waveband data sets suitable for exploring galaxy evolution continues to grow. In particular, the recent Extragalactic Database for Galaxy Evolution (EDGE, Bolatto et al., 2017) study combined a uniform sample of CO (1-0) emission from 126 galaxies as observed with the Combined Array for Millimeter Astronomy (CARMA) with the integral-field unit spectroscopy made by the Calar Alto Legacy Integral Field Area survey (CALIFA; Sánchez et al., 2012). The EDGE-CALIFA combined sample provides sufficient resolution to extend the line-of-sight based studies of molecular gas properties to a broader set of galaxy types while backing those data with the optical integral field unit (IFU) spectroscopy and modelling of the CALIFA survey. The EDGE-CALIFA data set provides an excellent testbed for variations in the star formation law and has already been applied to explore this relationship, focusing on variations in depletion time in the centres of galaxies (Utomo et al., 2017) and the dynamical influences on star formation (Colombo et al., 2018).

In this work, we explore the star formation law in the EDGE-CALIFA data from a different perspective, taking advantage of the rich database of galactic environments provided through the combined observational data. While the star formation law has been previously explored in the context of specific physical models, here we explore the data using machine learning methods designed to identify significant factors. This provides a data-driven selection of the important physical quantities, potentially offering new insights into the key factors shaping star formation in galaxies. Given the large number of observed quantities, we frame the question as one of variable selection where we seek to identify scalings that are the drivers of the star formation rate. The primary challenge making this selection stems from many of the variables being physically linked. For example, nearly every property (gas content, stellar surface density, metallicity) scales with distance from the centre of the galaxy, leading to large covariance among the physical parameters. Fortunately, the machine learning community has developed several methods for approaching such multi-covariant data sets and we rely on their established methods to identify quantities that are nonetheless significant. This approach, however, suffers from not having an underlying physical basis, but the outcomes can be used to identify the physical effects that could be considered in future explorations that aim to unravel the star formation law. Our overall strategy is to determine how the star formation surface density ($\Sigma_{\mathrm{SFR}}$) and the molecular gas depletion time ($\tau_{\mathrm{dep}}\equiv\Sigma_{\mathrm{mol}}/\Sigma_{\mathrm{SFR}}$) depend on other measurements of the local galactic environment (e.g., stellar surface density, metallicity) and on the galaxy properties as a whole (e.g., morphological type). This approach is similar to the fundamental plane analysis that is applied to early type galaxies, where we seek physical correlations in a high-dimensional space that are not apparent in bivariate analyses alone (e.g., Bernardi et al., 2003).

In Section 2, we outline the data used in this study and present the machine learning methods that guide our approach in Section 3. We present the results of these analyses in Section 4.

2 Data

We use the data acquired as a part of the EDGE-CALIFA survey (Bolatto et al., 2017). The EDGE survey constitutes the largest interferometric CO survey of galaxies in the nearby Universe and provides CO maps of 126 galaxies. The survey data have an angular resolution of $\sim 4.5^{\prime\prime}$, which translates into a typical spatial resolution of $\sim 1.4$ kpc given that the median distance to the galaxies is 64 Mpc. The CO data were obtained using the Combined Array for Millimeter-wave Astronomy (CARMA) interferometric array. These data were combined with the Integral Field Spectroscopy data obtained using the Calar Alto Legacy Integral Field Area (CALIFA) survey, obtained using the Calar Alto 3.5-m telescope, which provides an angular resolution of $\sim 2.5^{\prime\prime}$ (García-Benito et al., 2015). Both the CALIFA and EDGE data are smoothed to a common angular resolution of $7^{\prime\prime}$. The maps are then sampled on a hexagonal grid with $3.5^{\prime\prime}$ spacing so the data are nearly independent. The combined data sets allow us to construct maps of gas, stellar metallicities, extinctions, extinction corrected star formation rates, stellar mass densities and age. In addition, we also have the global parameters of these galaxies obtained from the LEDA catalogue (Makarov et al., 2014).

The joint EDGE-CALIFA data set provides a huge number of possible variables for consideration in this analysis. The CO data provide information on CO emission brightness, line width, and line-of-sight velocity. The CALIFA data include both gas emission line and simple stellar population (SSP) analyses as provided by the Pipe3D pipeline (v2.2, Sánchez et al., 2016b, a). We use global galaxy data as tabulated in Walcher et al. (2014). From the possible set of variables we chose the set of resolved and global variables as shown in Table 1 that could potentially be considered as key factors in controlling the star formation process. The variables under consideration are either resolved or global. A resolved variable is measured on each line of sight through the galaxy and global variables hold for galaxies as a whole. Variables can also be continuous or discrete, where continuous data can take any value and discrete data are fixed to a small number of values. In particular, we include discrete variables to reflect categorical data such as whether the galaxy has a bar or not.

For each variable under consideration, we examine the variance inflation factor (VIF), which measures the multi-collinearity with other variables in the data set. For example, many properties in galaxies change significantly with radius (star formation, stellar surface density, molecular gas surface density, metallicity, etc.). Therefore, these variables are typically correlated with each other and the VIF is a measure of how this multi-collinearity will affect the results of the regression. The variance inflation factor of a variable $X$ is expressed as $\mathrm{VIF}_{X}=(1-R^{2})^{-1}$ where $R^{2}$ is the coefficient of determination for a linear fit predicting $X$ given all the other independent variables. Quantitatively, the VIF measures how much the uncertainty in the regression is increased by the presence of linear correlations between the predictor variables. If the variance inflation factor is high (usually taken to be $>5$; James et al., 2013), the variable is multi-collinear with the other independent data and does not contribute new information to the analysis. The VIF analysis is particularly useful when compared to the standard covariance matrix analysis since it can identify correlations between more than two variables that limit the statistical power of the analysis. We retain all variables, even though some have large VIFs. Instead, we rely on our regression method (see Section 3) to select which of these correlated variables is most relevant to the star formation law.

In framing our analysis, we also considered the mass-to-light ratio of the stellar population, but found it had a VIF significantly larger than all the factors considered here (15 vs $\sim 2$). We removed this variable from our analysis because it was strongly correlated with several other factors linked to the age of the stellar population that had a clearer physical interpretation. A simple model of the stellar population could be also be formulated using the mass-to-light ratio as a proxy for several other physical effects.

For quantities that depend on the surface brightness of emission, we deproject the quantities by multiplying by $\cos i$ where $i$ is the inclination of the galaxy. For the analysis, we transform most variables using a log-transform (base 10). This reduces power-law scalings to a linear relationship, which is required since the variable selection techniques we will use are formulated in terms of linear models. Categorical variables are treated as indicator variables, taking a value of 1 if the galaxy or line-of-sight is in the category and 0 otherwise.

We measure the molecular gas surface density ($\Sigma_{\mathrm{mol}}$) from the EDGE data, using the CO integrated intensity values ($W_{\mathrm{CO}}$) calculated from Bolatto et al. (2017) and converting these into a molecular gas mass using a uniform CO-to-H2 conversion factor $\Sigma_{\mathrm{mol}}=\alpha_{\mathrm{CO}}W_{\mathrm{CO}}$ with $\alpha_{\mathrm{CO}}=4.35~{}M_{\odot}\mbox{ pc}^{-2}/(\mbox{K km s}^{-1})$ and $W_{\mathrm{CO}}$ being the integrated intensity of the CO emission (Bolatto et al., 2013). The velocity dispersion of the molecular gas ($\sigma_{\mathrm{CO}}$) is calculated from the masked, emission-weighted second moment of the CO line. We use the method described in Appendix B of Levy et al. (2018) to correct $\sigma_{\mathrm{CO}}$ for the effects of beam smearing. We further derive the peak rotation speed in the galaxy rotation curve ($V_{\mathrm{rot}}$) using the approach of Lelli et al. (2016) applied to the CO data.

We calculate star formation rates and gas-phase metallicities as per Bolatto et al. (2017) using the H$\alpha$ surface brightness and comparing the observed vs intrinsic H$\alpha$/H$\beta$ line ratio to infer the extinction $A_{V}$ assuming $R_{V}=3.1$. The extinction-corrected H$\alpha$ surface brightness is converted to a star formation surface density using the calibration of Kennicutt (1998a):

[TABLE]

We use a local gas-phase metallicity measurement from the strong-line N2 indicator as calibrated in Marino et al. (2013).

The parameters of the stellar environment are generated from the CALIFA data by Pipe3D (Sánchez et al., 2016a, b). The Pipe3D package fits a set of simple stellar population templates to the stellar continuum of the galaxy, convolving the set of templates by the spectral response of the instrument and a line-of-sight stellar velocity dispersion ($\sigma_{\star}$) to identify an optimal match to the observed spectrum (González Delgado et al., 2015). From the parameters of this template fitting, we determine the stellar surface density ($\Sigma_{\star}$), the stellar metallicity ($Z_{\star}$), the mass-weighted stellar age ($\tau_{\star,m}$), luminosity-weighted stellar age ($\tau_{\star,l}$), and the extinction to the stellar continuum across the CALIFA waveband (3745 Å– 7300 Å), referencing the value to the $V$-band ($A_{V}$).

Finally, we use the CALIFA and LEDA catalogue values to determine the galaxy orientation parameters (inclination and position angle), the mass of the galaxy ($M_{\mathrm{gal}}$), how the galaxy is classified morphologically (numerical Hubble type $T$), and whether it hosts a bar, ring or is part of a multiple-galaxy system. Finally, we use the orientation parameters to infer the mean galactocentric radius ($R$) for a given line of sight to and whether it is in the centre of the galaxy ($R<1$ kpc) or not. The galaxies in the EDGE-CALIFA samples have different physical sizes, so we normalize $R$ by the exponential scale length of the stellar disk, $R_{s}$ derived in Bolatto et al. (2017).

The full sample in the EDGE-CALIFA overlap region contains $7800$ lines of sight at $7^{\prime\prime}$ resolution. We then apply a series of cuts to the data to ensure clean, significant detections throughout the sample. First, we only include galaxies with inclination $i<70^{\circ}$ so that the deprojection correction is small (leaving 7001 data). We also require a joint detection of CO integrated intensity and H$\alpha$ emission at a 3$\sigma$ threshold (leaving 2074 data). We removed a further 229 data using a series of criteria. First, we also rejected all points that would be considered AGN above the Kauffmann et al. (2003) demarcation line based on a BPT diagram analysis (Baldwin et al., 1981) as well as all points with an equivalent width in H$\alpha<6$ Å since this emission is predominantly associated with older stellar populations ($>0.5\mathrm{~{}Gyr}$, Sánchez et al., 2014). After these down selections, we arrive at a final sample of 1845 independent data drawn from 39 different galaxies. We tested whether there were any significant differences between the global galaxy parameters for the entire EDGE sample (126 galaxies) and those hosting lines of sight in this analysis using a two-sided Kolmogorov-Smirnov test. We found no parameters were significantly different though the sample analyzed here is marginally closer (median distance is 56 Mpc vs. 68 Mpc) and less inclined than the overall EDGE population.

3 Methods

Galaxy-scale star formation is typically parameterized through the "Schmidt-Kennicutt law" or the star formation law (Schmidt, 1959; Kennicutt, 1998b). This relationship uses a power-law scaling between the surface densities of star formation rate (SFR), $\Sigma_{\mathrm{SFR}}$ and molecular gas mass, $\Sigma_{\mathrm{mol}}$:

[TABLE]

or with a logarithmic transform,

[TABLE]

Our line-of-sight based analysis recovers the star formation law (Figure 1). We fit the star formation law using a linear regression in the logarithmic space. Our sample data have significant internal scatter and do not have a well defined selection space in either of the axes (since the noise levels vary between the data sets). Thus, we treat the regression using a Theil-Sen estimator (Theil, 1950; Akritas et al., 1995) which determines the best fitting relationship as the median slope between the pairs of data in the sample. This regression method is moderately robust to outliers and is efficient to calculate. To further mitigate the effects of outliers, we perform random sample consensus (RANSAC; Fischler & Bolles, 1987) regression, which carries out repeated fits over subsets of the data to identify outliers and reduce the weight of these data in the fit. RANSAC operates by drawing random subsamples of the data, performing the Theil-Sen regression within that subset, and then identifying outliers with respect to that regression. Outliers for a subsample are identified as data are more than one median absolute deviation away from the trial regression line. The algorithm iterates, drawing new subsets, and identifying new outliers, until a stable set of outliers is identified, which are then ignored in the regression. The combination of these methods act to produce a robust fit in the presence of significant scatter and contamination by outlier data. The star formation law in this analysis is $\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{mol}}^{1.0\pm 0.1}$, consistent with similar approaches in the analysis (Leroy et al., 2013; Bolatto et al., 2017; Utomo et al., 2017). The outlier-rejection process identifies 30% of the data as outliers. There is a scatter of 0.3 dex around the linear relationship. However, at low $\Sigma_{\mathrm{mol}}$ the star formation rates fall below the fit line, suggesting a non-linear form to the star formation law. Using RANSAC mitigates the influence of these data on the overall relationship. The deviation could reflect a real change in the star formation rate in this sample, potentially associated with other variables we will examine. These data are near the significance limit for our sample, and this the behaviour could reflect biases in our star formation rate estimators at these low rates. We proceed with a linear model in this work for consistency with previous work that considers the effects of completeness and censoring more carefully (Leroy et al., 2013) and return to discuss this choice in Section 5.

When the index $N$ is consistent with unity, the data can be well described in terms of a near-constant molecular gas depletion time ($\tau_{\mathrm{dep}}$) with:

[TABLE]

The depletion time can be considered as the amount of time taken to convert all the gas present in the region into stars at the current star formation rate. It is found that the depletion time in galaxies is almost a constant with a typical value of 2.2 Gyr in nearby galaxies (Bigiel et al., 2008; Leroy et al., 2013; Utomo et al., 2017). In our data set, we find a median depletion time of $\tau_{\mathrm{dep}}=1.6\mbox{ Gyr}$.

3.1 A Multilinear Star Formation Law

We now consider the star formation law in terms of a multilinear fit where the star formation rate surface density, or alternatively the depletion time, are fit using a linear model. We use log-transforms as summarized in Table 1 to linearize the data. The formalism of linear models is particularly well developed in the fields of statistics and machine learning, enabling this study to use relatively common methods. This analysis is completed with the Scikit-learn (Pedregosa et al., 2011) package in Python.

We parameterize the star formation law using a linear model to a set of independent variables ($X_{ij}$ for the $i$th variable and $j$th datum) with a functional form:

[TABLE]

In our exploration, the dependent variable $Y$ is either $\Sigma_{\mathrm{SFR}}$ or $\tau_{\mathrm{dep}}$. Typically, for unweighted fits, Ordinary Least Squares (OLS) regression is optimized by finding $C_{i}$ that minimize the residual sum of squares $S$ over $n$ data:

[TABLE]

However, OLS initially assumes that all variables are relevant explanatory variables for the data, but finding $C_{i}\approx 0$ would indicate that that factor is not relevant in predicting the dependent variable. Since not all variables will necessarily be relevant for predicting $Y$, we can perform variable selection by including a regularization term in Equation 6. One approach is the Least Absolute Shrinkage and Selection Operator (LASSO; Tibshirani, 1996). The LASSO regression finds coefficients $C_{i}$ through a linear fit that minimizes the function:

[TABLE]

The extra term added is the LASSO penalty where the magnitude of the penalization is established by penalizing constant, $\alpha$. The regularization term provides a penalty when a given coefficient ($C_{i}$) is non-zero, so the regression prefers a model with the minimum number of explanatory variables that predict $Y$. Setting $\alpha=0$ returns to OLS and in the limit of $\alpha\to\infty$, $C_{i}\to 0$. Intermediate values of $\alpha$ fixes some $C_{i}=0$ depending upon their significance as explanatory variables. The formulation of this as the sum of the absolute values of the coefficients is important since it favours models where coefficients are exactly zero (James et al., 2013). Hence LASSO allows us to perform variable selection.

The parameter $\alpha$ needs to be determined for the regression. Since a single value of $\alpha$ applies to all variables, we standardize the data in the actual regression by subtracting off the mean values and normalizing each variable by its standard deviation. This makes the numerical values of each axis carry uniform weight. We find the optimal value of $\alpha$ using cross-validation (CV) following the methods outlined in Friedman et al. (2010, their Section 2.5). We first generate a grid of $10^{3}$ values of $\alpha$ ranging from the theoretical maximum ($\alpha_{\mathrm{max}}$) down to $10^{-3}\alpha_{\mathrm{max}}$, perform the regression for each of these $\alpha$ values, and assess which $\alpha$ value provides the most robust regression under CV testing. For CV, we use $K$-fold CV with $K=5$, which entails partitioning the data into five random subsets. Four of these subsets are used to establish the regression and the fifth is used as a test to determine how well the regression predicts their behaviour. The test subset is then rotated among each of the groups while the process is repeated. The optimal value of $\alpha$ is determined as the parameter that produces the most consistent fits between the training and validation sets. After the regression is complete, we transform the resulting coefficients back to the scales that would be observed for the original data values instead of the standardized variables.

LASSO is well established for data without uncertainties, but the theoretical framework used to quantify the errors including data uncertainty is not complete (James et al., 2013). We therefore quantify the errors through a bootstrap resampling of the entire data set, repeating the regression including the full CV process on a subset of the data drawn with replacement from the originals. The dependent variables are redistributed within their uncertainties with each bootstrap iteration. We repeat the bootstrapping for 1000 iterations and report the uncertainties in the derived coefficients as the standard deviation of the coefficient values under the bootstrap. We note that we experimented with using a RANSAC approach to wrapping the entire LASSO process to mitigate the effects of outliers but we found that this did not produce a discernible change in our results. By accounting for all the additional effects in the multilinear model, we find that the results are not significantly altered by outliers in our data. Critically, the bootstrapping shows that the variable selection performed by LASSO is robust for nearly all variables: coefficients found to be zero remain zero for the vast majority of bootstrap iterations.

3.2 Analysis of Mock Data

We demonstrate the LASSO approach to a multilinear star formation model by creating a mock data set consisting of 400 lines of sight. We randomly generate molecular gas and stellar surface densities by drawing these values from correlated log-normal distributions. The stellar and gas surface densities are assumed to have a correlation coefficient of $r=0.4$. Given these data, we create a mock star formation law

[TABLE]

where we select $a=0.6$ and $b=0.3$ and $C_{0}$ is a constant. We randomly assign 20% of the the lines of sight into barred galaxies for which the star formation rate is increased by 0.2 dex. Finally, we include 0.3 dex scatter in the SFR representing observational errors. We then process this mock data set through our analysis pipeline, mimicking the treatment that would be used in the real data. In addition to the mock data, we also include an additional factor in the data drawn as random deviates from the standard normal distribution. This additional factor is not included in the SFR model (Equation 8), but it is included in the regression analysis, representing a variable that LASSO should set to zero.

In Figure 2 we show the results of the analysis. The top panel illustrates the molecular gas star formation law. The data have been coloured to illustrate the trend with stellar surface density and the barred galaxy sample is highlighted in square points. Both higher stellar surface densities and being in barred systems are visibly at higher star formation rates, displaying the properties of our mock data. The RANSAC regression to the molecular gas star formation law finds $\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{mol}}^{1.0\pm 0.1}$, showing a shift from the expected value of $a=0.6$. Performing the LASSO analysis on this model provides estimates of $a=0.58\pm 0.02$ and $b=0.30\pm 0.02$, consistent with the expected model values. We also estimate that being in a barred galaxy increases the star formation rate by $0.22\pm 0.03$ dex and the additional random factor included in the model is correctly set to zero.

The correlation between the gas and stellar surface densities affects the results of the regression. When these parameters are correlated, the slope of the gas-only star formation relationship is steeper than the intrinsic scaling with the molecular gas surface density. With these mock data, we find $\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{mol}}^{1.0\pm 0.1}$, like the observed data (Figure 1). The multi-variate analysis is able to separate these effects, provided the correlation coefficient ($r$) is not near 1.0. In exploring different mock data sets, we find that the model cannot distinguish between the influence of stars and molecular gas when $r>0.9$, but this threshold will vary with the sample. When we separate the effects of gas and stars, we identify shallower indices for these two parameters and their intrinsic correlation leads to a steeper index in the gas-only star formation law. The same conclusion has been reached by Shi et al. (2011, 2018).

To highlight the behaviour of the LASSO regression, we display a partial regression plot in the bottom two panels of Figure 2. Here, we plot the $k$th variable against the residual of the star formation rate after correcting for all other variables in the fit:

[TABLE]

We subtract off the median of the data $\langle X_{ij}\rangle$ to account for the observed data not being centred on 0. For example, the data shown in the middle panel of Figure 2 is given by

[TABLE]

where $B$ is an indicator variable equal to 1 if the datum is in the barred sample in our mock data and 0 otherwise. The scalars (0.30, 0.22) are from the regression to the mock data. This visualization subtracts off the effects of all the other factors in the model ($\Sigma_{\star}$, barred galaxy) to highlight the influence of the molecular gas on the star formation rate. The slope of the line is $0.58$, which is the index of the molecular gas scaling alone. The approach performs well on these mock data, but the model is both simpler and has a cleaner data set than the actual observations.

4 Results

We apply our variable selection methods to predict $\Sigma_{\mathrm{SFR}}$ and $\tau_{\mathrm{dep}}$ in terms of the variables given in Table 1. We report the coefficients for these variables and their VIFs in Table 2 for linear models that explain $\Sigma_{\mathrm{SFR}}$ and $\tau_{\mathrm{dep}}$. For modelling the depletion time, we do not include $\Sigma_{\mathrm{mol}}$ as an independent variable since it is in the depletion time expression. Uncertainties in $\Sigma_{\mathrm{SFR}}$ are typically 0.07 dex given the flux errors in the H$\alpha$, H$\beta$ fluxes, but this does not include a potential systematic error due to changes in $R_{V}$. The depletion time error includes uncertainties from both the CO and the Balmer lines and has a typical magnitude of 0.19 dex.

4.1 The Multilinear Star Formation Law

In Figure 3, we show projections of the star formation law generated by the multilinear model for the selected variables we explored. Because a multilinear model is a plane in the high dimensional parameter space, the shaded region shown in each of the fits represents the projection of that plane into the bivariate plot as shown. The slope of each line is the coefficient $C_{i}$ given in Table 2 and the vertical extent of the region shows the possible values taken on by the multilinear regression based on the 5th to 95th percentile range of values seen in all the other dimensions. We only show a few of the predictor variables, but overall the regression spans the general trends in the data, but there remain several outliers. However, the regression is clearly different than the expectations for examining a single variable at a time (compare Figure 1 to the first panel of Figure 3). Showing the data in the original data-space gives a clear connection to the observations, but the regression hyperplane projects into the entire plotted region. The highlighted regime shows the portion of the hyperplane that would span the 90% of the data closest to the plane.

We also plot the partial regressions of the different factors against the incremental change in the star formation rate introduce that those factors in Figure 4. These plots illustrate the significance of the different predictor variables and they are ordered by the significance of the coefficients from the regression. The top, left panel of the Figure illustrates the regression against molecular gas, for which the tail to low star formation rates seen in Figure 1 is reduced but not eliminated. The strong correlation with stellar surface density is also apparent in this analysis (top right). We see clear declines in star formation with metallicity (gas and stellar), galaxy mass, and the age of the stellar population. In particular, the decline in the star formation rate is particularly strong for the mass-weighted relationship, showing this effect is driven by the older stellar populations since these will contribute relatively more to the mass-weighting than to the luminosity-weighting.

The LASSO method is intended to identify factors that are irrelevant to the star formation law and several factors have low significance and are set to zero in a significant fraction of the bootstrap iterations. This analysis finds that star formation rate can be explained without including the effects of the kpc-scale dust extinction for the stellar population ($A_{V}$) or the presence of a ring. The lack of an effect from extinction suggests that the extinction correction for star formation rate based on the Balmer decrement is largely correct. These marginally significant factors include the CO velocity dispersion ($\sigma_{\mathrm{CO}}$), the rotation speed of the galaxy ($V_{\mathrm{rot}}$), the luminosity-weighted age of the stellar population ($\tau_{\star,l}$) and lines of sight through galaxy centres. The lack of a significant change based on the centres of galaxies arises because the sample mixes regions with enhanced star formation rates relative to their molecular gas content in with regions that have suppressed star formation (Utomo et al., 2017).

If we consider only those factors that have coefficients larger that $8\times$ their uncertainty, this analysis argues for an empirical, multi-linear star formation law of the form:

[TABLE]

where we have defined $Z_{\mathrm{gas}}=10^{\log(\mathrm{O/H})}$.

The multi-linear star formation law highlights several important factors at setting the local star formation rate. As expected, the molecular gas surface density is a clear, positive predictor of the SFR, but the model finds a relatively shallow slope for the scaling with $\log(\Sigma_{\mathrm{mol}})$: $0.43\pm 0.03$. This shallow index contrasts with the gas-only index of the star formation law which shows a coefficient of $1.0\pm 0.1$ in our approach and others. This linear model analysis favours a case where the other strong, positive predictor of star formation is stellar surface density ($\Sigma_{\star}$). The stellar surface density traces the local gravitational field in disks and could be a strong controlling factor in star formation. The positive scaling with surface density is expected in models based on equilibrium density of the neutral medium being set by the self gravity of galactic disks (Ostriker et al., 2010), but the stellar velocity dispersion should have a suppressing effect if that model held in isolation since it would imply a higher stellar scale height and lower midplane stellar densities.

The difference between a gas-only star formation law and a multifactor relationship can be reconciled because the molecular gas surface density and stellar surface density are correlated with each other. However, the application of this analysis suggests that the stellar surface density has a strong, independent influence on the local star formation that is more significant than the molecular gas content. All lines of sight in our analysis, by construction, have detectable molecular gas content. We further illustrate these two strong correlations in Figure 6 which shows the scaling between the star formation rate and the local stellar surface density in six equal quantiles of the molecular gas surface density. Within these narrow ranges of the molecular gas surface density, the star formation rate shows a significant positive scaling with the local stellar surface density. This figure illustrates that both of these factors (stellar and gas surface density) separately have positive correlation with the local star formation rate. Such a link with stellar population has been noted before (Bolatto et al., 2017; Shi et al., 2018) in single-factor analyses. The binned trends in Figure 6 also show some curvature, indicating that a linear model is only an approximation to potentially more complex relationships between these data. Of note, the scaling with stellar surface density appears shallower than the average value of 0.7 for galaxies of moderate surface densities.

Several factors are associated with a reduction in the local SFR, all of which are linked to more evolved galaxies. In particular, larger galaxy masses, older stellar populations (as traced by the mass-weighted age), and earlier morphological types all suppress the local SFR. These can be broadly described as linked to older stellar populations, but the exact physical effects that yield the reduced star formation efficiency are unclear (see also Colombo et al., 2018). With all other effects being held constant, molecular gas in a galaxy with an old stellar population forms stars less quickly compared to the same gas in a late-type galaxies with relatively young stars. We also find that molecular gas found in barred galaxies has lower star formation rates by $-0.038\pm 0.005$ dex.

To display the relative importance of different factors in the model, we construct a standardized slope out of each coefficient $C_{i}$ and summarizes these values in Figure 5. The standardized slope is determined from the standardized data used in the fit, i.e., the data scaled to have zero mean and a standard deviation of unity. The standardized slope is related to the coefficients in Table 2 by

[TABLE]

where $\sigma_{X_{i}}$ and $\sigma_{Y}$ are the standard deviations of the corresponding variables. The uncertainties in the coefficients are scaled by the same ratio and plotted as error bars. This visualization shows the fractional change in the $Y$ variable (here, the star formation rate) predicted by the the change in the $i$th independent variable across the range of data contained in our analysis. Points found close to zero indicate that the importance of this factor is relatively weak and the uncertainties indicate the significance of the estimate.

This visualization separates the dominant influences on the local star formation rate from lesser effects that nonetheless have a significant influence on the data. The most significant effects here are encoded into Equation 11, but several factors emerge as less important but nonetheless still significant. The stellar metallicity and the luminosity-weighted mean age of the stellar population are both linked to reduction of SFR. The mass-weighted age shows a more significant negative effect on the SFR, so whatever effect is measured here must be linked to the old stellar population.

Similarly, lines of sight with higher gas-phase and stellar-phase metallicities also show relatively lower star formation rates. Overall, the reduced levels of star formation associated with older stellar populations, earlier Hubble types, and higher rotation speeds may all be proxies for higher metallicities. However, we do include both local stellar and gas phase metallicity as predictor variables and the LASSO method should eliminate these effects in favour of metallicities if this is the driving effect.

If metallicity is a driving effect, the results we see may stem from our assumption of a constant CO-to-H2 conversion factor underestimating the amount of molecular gas in these lower metallicity regions. This apparent suppression may emerge from a variable CO-to-H2 conversion factor since the conversion factor will increase at lower metallicities. This effect would imply more molecular gas than the CO measurements indicate in our model, which would imply a lower overall efficiency of the star formation process, making low- and high-metallicity systems more consistent in their star formation behaviour. However, the conversion factor is thought to be relatively stable down to metallicities of $<0.5Z_{\odot}$ before increasing sharply (Bolatto et al., 2013). Our data are typically found at higher metallicities, so this effect would probably be insufficient to explain the robustness of the anticorrelation of star formation rate with older stellar populations. Even so, we have repeated our analysis using the variable CO-to-H2 conversion factor prescription presented in (Bolatto et al., 2013) using the gas-phase metallicity derived from the N2 indicator as an estimate of the local variation in metallicity. Including these effects has no significant effect on this analysis. The scaling coefficient for $\Sigma_{\mathrm{mol}}$ increases from $0.43\pm 0.03$ to $0.44\pm 0.03$. Including the variable conversion factor marginally decreases the scaling of star formation rate with gas phase metallicity with the index changing rom $-0.27\pm 0.03$ to $-0.25\pm 0.03$ but there are no other significant changes in the derived star formation scalings. With variable conversion factor, the significant anticorrelation between older stellar populations and reduced star formation rates persists. The variable conversion factor prescription also argues for a dependence on the surface densities of giant molecular clouds measured on $<10^{2}$ pc scales. We are unable to assess these effects with our coarse resolution data, but nearly all of our lines of sight are through regions with surface densities $<10^{2}~{}M_{\odot}\mbox{ pc}^{-2}$ so these opacity effects are unlikely to change our results significantly.

The metallicity effect may be more subtle, where gas in high metallicity systems preferentially highlights “sterile” molecular gas that does not form stars. In high metallicity systems, dust shielding of CO formation is more efficient, allowing CO to form at relatively low column densities. This is the opposite effect from what is seen in low metallicity systems (Bolatto et al., 2013). This low density gas could then be found outside of bound molecular clouds, which would manifest as a reduced star formation rate per unit molecular gas.

The significant positive scaling of star formation rate with normalized galactocentric distance likely indicates physical factors that have not been accounted for in our model but nonetheless influence the local star formation rate. For example, several dynamical factors such as the epicyclic frequency or shear are not well represented in our sample data. However, these effects have been hypothesized to play a major role in regulating the local star formation rate in galaxies (e.g., Meidt et al., 2018). Completing the analysis without normalizing by the disk scale length ($R_{s}$) eliminates the significance of the radial scaling but does not significantly alter the derived scalings for other parameters. This suggests that the hidden parameters that are not being explored in this analysis are closely linked to the disk properties.

A multilinear star formation law does reduce the scatter in the relationship compared to a gas-only star formation law. In Figure 7, we show the distribution of residual data for gas-only and the LASSO selected best-fitting relationship. The intrinsic scatter drops from 0.28 dex to 0.18 dex. The clearest systematic effect in the star formation law is the tail of points to lower star formation rates. These are typically associated with the oldest stellar populations and earliest Hubble types present in our sample. Including these parameters of the stellar populations noticeably reduces the tail to lower star formation rates, but does not eliminate its presence altogether. The scatter in the relationship remains in excess of that expected for the dependent variable ($\Sigma_{\mathrm{SFR}}$) but we do not account for uncertainties in the independent variables and several variables are discrete (Hubble type).

4.2 Multilinear Model for Depletion Time

We also consider a model for the star formation in terms of the depletion time and compute a LASSO regression for the predictions of the depletion time in terms of all the parameters except for molecular gas surface density. The results for the regression and the VIFs for the different factors are given in Table 2. In Figure 8, we show the partial regressions for several selected variables in the depletion time relationship and Figure 9 shows the scaled coefficients for different factors.

The regression leads to empirical model for the depletion time, based on scalings with coefficients $>6\times$ their associated uncertainties:

[TABLE]

These results are consistent with the those seen for the star formation law except stellar velocity dispersion appears as a more significant predictor variable here, more so than stellar surface density, which is the opposite of what was found in the star formation law analysis. We also find weak scalings in the depletion times which are longer for older stellar populations (higher $\tau_{\star,m}$), higher metallicities, higher galaxy masses, and earlier morphological types. Even though higher stellar surface densities and velocity dispersions are associated with older stellar populations, this analysis finds that these are separable effects that act in opposite directions.

One notable difference with respect to the star formation law model is that higher values for the velocity dispersion of the molecular gas $\sigma_{\mathrm{CO}}$ predict longer depletion times. There is a clear trend that larger $\sigma_{\mathrm{CO}}$ leads to longer depletion times, best seen in Figure 8. This is generally consistent with expectations given the dynamical state of the gas. Leroy et al. (2017) showed that the depletion times for molecular gas in M51 increased if the gravitational “boundedness” ($\propto\Sigma_{\mathrm{mol}}/\sigma_{\mathrm{CO}}$) decreased. Our results act in the same sense, but the weak scaling in the star formation law approach (Section 4.1) shows that this effect may manifest from how we are framing the analysis. We discuss this further in Section 5.

The morphological indicator variables (bar, ring, multiple) show more significant effects in terms of the depletion time model compared to the star formation law described in Section 4.1. In Figure 10, we show the distributions of depletion times for the sample split into these different categories. These indicator variables are global: all lines of sight in a galaxy with a bar are considered as being in a barred galaxy, irrespective of whether the line of sight is actually located in the bar region. Future work with resolved environmental descriptions will allow a more careful treatment of these effects.

As before, the linear regression indicates that the presence of a bar leads to a marginal reduction in the star formation rate and a corresponding increase in the depletion time but only by 0.04 dex. We also noted the 0.04 dex decrease in the star formation rate found for the star formation law (Section 4.1). The reduced star formation efficiency in barred galaxies initially appears at odds with the analysis carried out in Saintonge et al. (2012), which analyzed entire galaxy star formation efficiencies in the COLD GASS sample. Like our analysis, they found that the effect of bars was not strongly distinguishable from the remainder of their sample. The COLD GASS sample is a more representative selection from the galaxy population, but the EDGE-CALIFA sample is measuring the resolved star formation rate on smaller scales than in the COLD GASS sample. Thus, our measurements are weighted on a by-area basis within galaxies whereas the COLD GASS sample is measuring overall efficiency of galaxies. While these morphological features may lead to significant changes in the local depletion time within galaxies, our analysis finds minimal effects from these factors on the galactic scale, especially when compared to the local environmental effects we are able to assess from this data set.

In both the star formation law and here with the depletion times, morphological classification as a ring galaxy does not lead to significant changes in the depletion time. The multilinear depletion time model finds that galaxies classified as belonging to multiple systems have measurably shorter depletion times than the population as a whole. However, in Figure 10, most of the data flagged as coming from “multiple” systems have longer depletion times than the sample as a whole. However, there are only 85 of the 1845 of the lines of sight labelled as in multiple-galaxy systems drawn from five systems. Most of the statistical weight in this result comes from the data drawn from the single galaxy NGC 3994 and excluding this system eliminates any difference in the depletion time from multiple-galaxy systems. Examining the distribution of points in Figure 10 shows that most lines of sight in multiple systems appear to have longer depletion times. Thus, we do not consider this result to be significant. In contrast, the ring and barred categories consist of 577 and 727 data respectively.

In Figure 11, we show the residuals of the multilinear depletion time model compared to that of the constant depletion time model. Again, there is a significant tail in the constant depletion time model and the multilinear model toward high depletion times. The tail is somewhat reduced in the multilinear model, but the reduction in scatter is not as large for this model compared the star formation law. The expected uncertainties in the observed depletion time estimate are 0.19 dex, and the reduction of 0.28 dex to 0.20 dex in scatter is approaching the limits of expected from noise.

5 Discussion

This work focuses on using the tools of linear modelling and variable selection drawn from the field of machine learning and applied to the questions of the star formation law. The initial results of this study are promising and connect many different lines of research seen individually through the literature. However, there are clear areas where the work can be extended in future efforts.

5.1 Context for Interpretation of Empirical Results

Two key points emerge from this data-driven analysis of how galaxies form stars. First, the data clearly show the gas-only star formation law has a linear scaling with $\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{mol}}^{1.0}$; however, once additional factors are included, the scaling relationship with molecular gas becomes significantly shallower. The star formation is usually considered to emerge from the properties of the cold molecular medium, as traced by the CO emission. Here, we can consider star formation in the context of factors deriving from the galactic environment (i.e., everything but the properties of the cold ISM). Our analysis of mock data and the actual system argues that the intrinsic scaling with molecular gas surface density can be shallower than the nominal index of 1.0, with the correlations between different environmental factors, in particular stellar surface density, driving the resolved gas-only star formation law to the standard index of 1.0.

Taking our empirical models at face value (Equation 11 and 13), this would imply a strong environmental dependence setting the efficiency of the star formation process with the properties of molecular gas playing a secondary role. Deeper explorations of the single-variable star formation law have also pointed to many of these same environmental effects to address the trends seen in earlier studies (Leroy et al., 2008; Leroy et al., 2013; Meidt et al., 2013, 2018; Colombo et al., 2018). Instead, the analysis constantly identifies the role of the stellar environment as a key driver. Treated as a linear model, the stellar surface density is a positive predictor for the star formation rate.

The second key result that emerges is that we see evidence that the age of the stellar population anti-correlates with the star formation rate. Moreover, the anticorrelation with the age of the stellar population is stronger for the mass-weighted age instead of the luminosity-weighted age. In Section 4.1, we argue that the effects of changing metallicity on the CO-to-H2 conversion factor are unlikely to explain these effects as observational effects. We are faced with a significant and real effect: all else being equal, molecular gas found in environments with old stellar populations is less efficient at forming stars.

The significant dependence on age is likely a proxy for the dynamical state of galaxy: younger systems are dynamically cold (i.e., small random motions in the vertical and radial directions), and this enhances the star formation rate per unit molecular mass. Such effects could change the local stability of the stellar dynamical systems, but both the stellar surface density and velocity dispersion show up as positive predictors for star formation. The velocity dispersions of the molecular gas and the stars are analyzed as the directly observed quantities. They have not been corrected into the components of the velocity ellipsoids for either of these species, though such analyses are becoming available (Zhu et al., 2018). Further, this analysis considers kpc-scale resolution elements, so averaging over these large regions will also introduce significant beam smearing across coherent velocity structures. As higher resolution or more careful analyses become available, a natural extension will be to include refined estimates of the dynamical environment in this type of analysis.

The influence of the local stellar population on the star formation rate also emerges from studies focusing on the optical IFU data alone (e.g., González Delgado et al., 2016; Ibarra-Medel et al., 2016; García-Benito et al., 2017; Belfiore et al., 2017; López Fernández et al., 2018; Sánchez et al., 2018). The more detailed stellar analysis of these works make the case that the age of the stellar population is a good proxy for its current dynamical state, with older populations being dynamically hotter. Moreover, these works also establish the morphological type correlations with the stellar population: earlier Hubble types have older, dynamically hotter stellar populations. Optical IFU analyses can also infer the molecular gas content from extinction measurements (Sánchez et al., 2018), all of which find consitent analysis with the results presented here: molecular gas in older stellar populations is relatively less efficient at forming stars. Star formation history analysis also suggests why systems with higher stellar metallicity have less efficient star formation: old systems had early bursts of star formaton, higher enrichment, and would be dynamically hotter today (Vale Asari et al., 2009). Physically, the likely origin of reduced star formation rate in dynamically hot systems arises from the increased stellar scale height for a given stellar surface density. The molecular gas then follows the shallower stellar potential, reducing its average volume density and increasing its mean free-fall time (Ostriker et al., 2010).

While consistent with the link to older stellar populations, this conjecture is undermined by our finding of a positive scaling with stellar velocity dispersion ($\sigma_{\star}$). In isothermal stellar disks, the stellar volume density $\rho_{\star}\propto\Sigma_{\star}^{2}/\sigma_{\star,z}^{2}$, where $\sigma_{\star,z}$ refers to the vertical velocity dispersion. Thus, we would expect the scaling from the stellar velocity dispersion to have a negative coefficient where high vertical stellar velocity dispersions are associated with lower star formation rates. Instead, we see the opposite: a weak positive correlation with $\sigma_{\star}$. Testing this conjecture will require a dynamical analysis to measure the vertical component of the stellar velocity ellipsoid, separate from the line-of-sight velocity dispersion on which our analysis is based.

The lower star formation rates seen in older stellar populations could also be a proxy for those systems where feedback from active galactic nuclei (AGN) quenches the star formation in these systems (Springel et al., 2005). However, most quenching through feedback requires the AGN winds to drive material out of the galaxy. This analysis highlights, instead, that the molecular gas already residing in such systems has a lower star formation rate. Hence, our analysis shows that the effects of morphological quenching are visible in this sample (Martig et al., 2009), though the analysis cannot quantify the relative importance of AGN-driven and morphological quenching. Colombo et al. (2018) showed the clear scaling of depletion time with the morphological type of galaxies for systems with molecular gas, with earlier types have lower star formation rates. However, Sánchez et al. (2018) also show that AGN are also preferentially found in morphologically earlier types of galaxies, so the effects may be difficult to distinguish.

We also see that the LASSO method identifies different factors as being significant when the analysis is recast in terms of the depletion time. The choice to explore depletion time is generally motivated by the linear relationship between molecular gas and star formation rate surface density. Since our multilinear star formation law shows a non-linear scaling, framing the analysis in terms of depletion time asks a different statistical question of the data. Functionally, working with depletion times asserts that the amount of molecular gas is the primary factor at setting the star formation rate, so other effects are modifying a fundamental $\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{mol}}^{1.0}$ relationship. Asserting this scaling with index of 1.0 highlights other effects as important, notably the velocity dispersion of the molecular gas and the morphological (Hubble) types and shape classifications (bar, ring, multiple). Like the analysis of Colombo et al. (2018), we find that the morphological type is a significant factor at regulating depletion time, of comparable significance and strength as the galaxy mass. While the results for these two types of analyses are broadly consistent, the two approaches remain conceptually distinct when the scaling of star formation rate with molecular gas content is considered in a multi-factor analysis. New factors emerge in the depletion time relationship because the multilinear star formation law finds the index on the molecular gas scaling to be 0.4 not 1.0 (Equation 11). In a study of 15 nearby galaxies, Sun et al. (2018) find a strong correlation between $\Sigma_{\mathrm{mol}}$ and $\sigma_{\mathrm{CO}}$ on 100 pc scales. The stronger scaling with $\sigma_{\mathrm{CO}}$ manifests only the depletion time relationship because index of 0.4 in the star formation law is sufficient to account for the internal correlation between the molecular gas properties.

5.2 Caveats and Future Work

The clearest forthcoming improvement for this analysis is that we will be able to directly test the connection between the dynamical state of the stellar system and the star formation rate of the molecular gas at a given location. This step will require full orbit decomposition of the galaxies, which is becoming possible for nearby systems (Zhu et al., 2018). In particular, we will directly be able to test whether the stellar scale height explains all of the variation that manifests in our empirical star formation law in terms of Hubble type, stellar population age, and stellar velocity dispersion.

In carrying out the analysis, we tested the importance of the log of the distance to the target in the analysis, which does show up as a moderately significant effect in the star formation law ($C_{D}=-0.2\pm 0.5$) but only results in small changes to the coefficients of the star formation law (e.g., $<1\sigma$ so the scaling for molecular gas shifts from 0.42 to 0.41). However, this is an indicator that the EDGE sample, so analyzed, suffers some effects from Malmquist bias. Specifically, the more distant galaxies are also those that are more luminous, with older stellar populations. Ideally, the analysis would be executed at a common linear resolution. However, the analysis already being carried out at kpc scales and reaching a common linear resolution for the whole sample would reduce many galaxies to single points, so future work would require higher resolution CO observations and deeper observations of more distant galaxies.

As the number of galaxies available for studies of the star formation law grows, it will become possible to examine the star formation law in terms of the galaxy cluster environment. The initial studies of cluster effects on the star formation law have largely focused on binary tests (cluster vs. field galaxies, Vollmer et al., 2012; Mok et al., 2016). These studies find that cluster galaxies have longer depletion times compared to the field systems. However, there are a range of different cluster environments, which readily could be integrated into a multilinear star formation model to assess their significance.

One other area for improvement in this work lies in the statistical formulation. The main advantage of this approach is that we can simultaneously assess all the possible factors contributing to the star formation rate. However, to apply the LASSO method, we reduced the problem of star formation to a simple linear model. While the regularization of the fitting through the $\alpha$ parameter (Equation 7) is established through cross-validation, the linear model does not include either uncertainties in the independent variables, which are significant, nor does it include any contributions from the intrinsic variance of the data. We attempt to include the the latter effect through the bootstrapping of the population. We also have not fully modelled the censoring of our data based on the different selection effects present in the work. Thus, the uncertainties for coefficients derived through the analysis (Table 2) may be underestimated. Finally, we adopted a linear model but there appears to be some non-linear behaviour in the full and partial regressions of star formation rate against molecular gas surface density. For consistency with previous work, we used adopted a linear model to describe the star formation law. Without a careful treatment of censoring and possible biases in CO or star formation tracers at low sensitivity, we cannot make strong claims about the shape of the star formation law. Given these limitations, a full Bayesian approach to the problem is merited in future work (Tibshirani, 1996; Gelman et al., 2014). However, the primary purpose of this work is to demonstrate the efficacy of variable selection approaches and identify the strongest positive and negative predictors of the star formation rate in galaxies through an empirical analysis.

6 Summary

Understanding what factors shape the star formation rate in galaxies drives a wealth of theoretical and observational work. In this work, we have used the combination of CO(1-0) line observations of molecular gas with optical IFU data in the EDGE-CALIFA sample to explore a new data-driven approach to understanding what factors shape the star formation law. Specifically, we have used the Least Absolute Shrinkage and Selection Operator (LASSO) modification to a linear model for the star formation rate to identify what variables most strongly predict the star formation rate surface density ($\Sigma_{\mathrm{SFR}}$) and local depletion time of the molecular gas ($\tau_{\mathrm{dep}}$).

Using the combined EDGE-CALIFA sample, we have identified a set of 1845 lines of sight spanning 39 galaxies that have well determined gas and stellar properties. Performing a simple linear model analysis on how the star formation rate scales with local molecular gas surface density recovers a similar star formation law as is found in other samples: $\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{mol}}^{1.0\pm 0.1}$ (cf., Leroy et al., 2013; Bolatto et al., 2017; Utomo et al., 2017). Given the combined EDGE-CALIFA data sets, we identified a set of 18 possible factors (Table 1) that could contribute to the star formation rate. Given the LASSO model analysis, we find:

The star formation rate increases with increasing molecular gas surface density, and stellar surface density. The scaling with stellar properties is statistically more significant and stronger than the scaling with molecular gas surface density. 2. 2.

While the star formation rate increases at higher stellar surface density, it also decreases for several factors associated with older stellar populations including the mass-weighted age of the stellar population, the gas- and stellar-phase metallicity and the mass of the galaxy. Earlier Hubble types are also linked to lower star formation rates. 3. 3.

The star formation rate increases with galactocentric radius normalized by the disk scale length, suggesting additional paramters regulating the star formation rate not explored in this study. One potential candidate for such a regulation mechanism is dynamical regulation of star formation through shear. 4. 4.

Resolved star formation rates in barred galaxies are 0.04 dex lower than in non-barred galaxies after controlling for other factors, but this result is of marginal significance. 5. 5.

The analysis of the molecular gas depletion time finds the same factors that the star formation law analysis does. However, we also find depletion times are significantly longer in regions with higher molecular gas velocity dispersion. 6. 6.

Developing a multilinear model for the star formation rate and the depletion times reduces the scatter in both of these relations, but neither model reaches the limits imposed by observational uncertainties in the dependent variables.

With the advent of resolved, multiwaveband surveys, using a machine learning approach to analyses of the star formation rate will enable the identification of the key factors in galaxy evolution. This approach is complementary to theoretical studies and provides a new way of assessing whether specific factors are actually important in different models.

Acknowledgments

We are grateful to an anonymous referee whose review improved the clarity and robustness of the results. BD is supported by a MITACS Globalink fellowship. ER is supported by Discovery Grant from NSERC of Canada (RGPIN-2017-03987). DU is supported by the National Science Foundation under grants No. 1615105, 1615109, and 1653300. SFS thanks the support of the CONACYT grant CB-285080 and funding from the PAPIIT-DGAPA-IA101217 (UNAM) project. In addition to the software cited in the text, we gratefully acknowledge the use of matplotlib (Hunter, 2007), pandas (McKinney, 2010), and astropy (The Astropy Collaboration et al., 2018; Astropy Collaboration et al., 2013). RGB acknowledges financial support from the Spanish Ministry of Economy and Competitiveness through grant AYA2016-77846-P. We gratefully acknowledge conversations with Eric Koch and John Braun who suggested the LASSO approach to this data set. Support for CARMA construction was derived from the Gordon and Betty Moore Foundation, the Eileen and Kenneth Norris Foundation, the Caltech Associates, the states of California, Illinois, and Maryland, and the NSF. Funding for CARMA development and operations were supported by NSF and the CARMA partner universities. This study uses data provided by the Calar Alto Legacy Integral Field Area (CALIFA) survey (http://califa.caha.es/), based on observations collected at the Centro Astronómico Hispano Alemán (CAHA) at Calar Alto, operated jointly by the Max-Planck-Institut fűr Astronomie and the Instituto de Astrofśica de Andalucá (CSIC).