The Variability of the Star Formation Rate in Galaxies. I. Star Formation Histories Traced by EW(Hα) and EW(Hδ_A)

The relatively long timescales for the formation of individual stars means that measurements of the rate of star formation must necessarily represent averages over some even longer timescale, say 10⁷ yr. Ideally, we would like to have a change parameter that reflects the change of the SFR, as measured within some fixed time interval, e.g., 10⁷ yr, over some other, longer, time interval, say 10⁹ yr. Unfortunately, this is not possible with current observational material, and it is in fact hard to see how it ever will be. Practicalities, therefore, force us to instead compare the SFRs that are obtained by averaging over different periods of time prior to the epoch of observation, e.g., to compare the SFR averaged over the previous 5 Myr with that averaged over the previous 800 Myr. As noted above, we adopt a shorthand of SFR7 and SFR9 for these quantities, with the ratio denoted by SFR79.

The ratio SFR79 therefore mixes information both on short-term (10⁷ yr) variations in the SFR, i.e., on the "burstiness" of star formation, with longer-term drifts in the SFR of the galaxy taking place on longer timescales (10⁹ yr). For this reason, we prefer to think of the ratio SFR7/SFR9 as a star formation "change parameter" rather than simply as a measure of the "burstiness" of the star formation. To think of "bursts" of star formation implies values of SFR79 greater than unity. This may be appropriate for some subset of the galaxy population, but within the overall population, we would expect to find some values of SFR79 below unity, and indeed the average SFR79 should be roughly unity. To be more precise on this point, we would expect the ratio of the average SFR7 divided by the average SFR9 (which will not be precisely the same as the average SFR79) to be unity, modulo any long-term evolution of the SFR of galaxies with cosmic time.

A galaxy with a constant SFR will have an SFR79 of exactly unity. A galaxy with constant sSFR will have an SFR79 that is greater than unity by an amount that depends on that constant sSFR, because the increase in mass during the last Gyr will have produced an (exponentially) increasing SFR. However, this effect is small if the mass-doubling timescale sSFR⁻¹ is long compared with 1 Gyr, as will generally be the case for galaxies at the present epoch. This effect will be discussed further below.

SFR79 will also give information on the movement of a galaxy in the SFR–mass plane. If we neglect the changes in the stellar mass over the timescales of interest, i.e., if sSFR⁻¹ ≫ 1 Gyr, then the SFR79 will tell us the present location of a galaxy on the SFR7–mass plane compared to the average position it has occupied over the last 10⁹ yr. In this sense, it tells us whether the individual galaxy is broadly moving up or down relative to its SFMS.

The basis of the calibration is that the three chosen observational parameters contain information about the (specific) SFR averaged within 5 Myr and roughly 1 Gyr, and therefore, the change parameter can be derived from a combination of these three diagnostic parameters.

Here, we briefly describe our overall approach to derive the change parameter. The details will then be presented in the following subsections. We first construct millions of mock SFHs of galaxies. These mock SFH span the whole of cosmic time and should cover as much as possible the range of SFH encountered in the real universe. We then generate synthetic spectra of these mock galaxies at the present epoch based on stellar population models for a range of different metallicities. We then measure the three diagnostic parameters of interest. We also compute the actual SFR79 from the mock SFHs. Finally, we search for the solution of SFR79 in terms of the three diagnostic observational parameters.

The three diagnostic parameters have long been used to indicate the recent SFHs on different timescales (e.g., Worthey & Ottaviani 1997; Balogh et al. 1999; Kauffmann et al. 2003; Li et al. 2015; Wang et al. 2017, 2018). In SF galaxies, the Hα emission mainly comes from the recombination of gas ionized by photons from extremely massive stars (>15M_⊙), which is therefore expected to trace the SFHs within the lifetime of these massive stars (∼5 Myr). However, the EW(Hδ_A) traces the recent star formation within a longer timescale of more like 1 Gyr. The Balmer absorption lines arise from intermediate mass main-sequence stars with lifetimes of ∼1 Gyr. They are relatively insensitive to the metal abundance because they depend mostly on the behavior of the main-sequence turn-off temperature rather than the behavior of the red giant branch temperature (Worthey et al. 1994; Worthey & Ottaviani 1997). Finally, the 4000 Å break is determined by the SFHs on still longer timescales with respect to EW(Hδ_A) and is found to be sensitive to the light-weighted stellar age (Balogh et al. 1999).

The evolution of these three diagnostic parameters for single stellar population (SSP) models of different metallicities can be seen by using the Flexible Stellar Population Synthesis code (FSPS; Conroy et al. 2009). FSPS² is a powerful code that can generate spectra and absolute magnitudes of arbitrary stellar populations, with a series of flexible settings, such as metallicity, choice of stellar library, different IMF, and different evolutionary isochrones. Throughout this work, we will adopt the MILES stellar library (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011), a Chabrier (2003) IMF, and the Padova isochrones (e.g., Bertelli et al. 1994, 2008), unless specified otherwise.

Nebular emission is produced by using the FSPS implementation of the photoionization code, CLOUDY (Byler et al. 2017). By simulating physical conditions within a gas cloud, CLOUDY predicts the thermal, ionization, and chemical structure of the cloud, and further produces the resultant spectrum of the diffuse emission (Ferland et al. 2013). In the FSPS model, the ionizing radiation is produced by a point source at the central of a spherical shell of cloud, with assuming a constant gas density of n_H = 100 cm⁻³. The fraction of the ionizing luminosity to escape from the H ii region is assumed to be zero (Byler et al. 2017).

For each of our mock SFH, we produce the current-epoch spectrum for six metallicities (log₁₀ Z/Z_⊙ = 0.0, −0.2, −0.4, −0.6, −0.8, and −1.0), without implementing any dust attenuation. This means that the diagnostic parameters in the calibration are assumed to be dust-free. Being equivalent widths, the observed diagnostic parameters should, ideally, be independent of dust attenuation, but this will only be true if the nebular emission and stellar continuum have the same attenuation, which is unlikely to be the case. There are likely be second-order effects if different stellar populations have different dust obscuration. The correction of the observed diagnostic parameters for these second-order dust effects will be presented in Section 2.6. Although we will only calibrate the SFR79 for the six discrete metallicities, we can obtain the SFR79 calibration for galaxies of other metallicities by linear interpolation in log₁₀ Z/Z_⊙ (see details in Section 3.1).

Figure 1 shows the evolution of EW(Hα), EW(Hδ_A), and D_n(4000) as a function of stellar age for SSP models at the six different metallicities. We present the evolution of EW(Hα) (solid lines) and EW(Hδ_A) (dashed lines) on the left panel of Figure 1. For all of the different metallicities, the EW(Hα) is scaled to 10000 Å, and the EW(Hδ_A) is scaled to 10 Å. The two vertical dotted lines represent the ages of 5 Myr and 800 Myr, respectively. As expected, after a single burst of star formation, EW(Hα) is large at first, but then quickly decays and becomes only 10% of its maximum value after 5 Myr, for all six different metallicities. In addition, EW(Hα) is higher at lower metallicities. The dependence of EW(Hα) on metallicity comes from the nearly equal contribution of the variation in stellar continuum and the variation in Hα emission.

Zoom In Zoom Out Reset image size

For an SSP model with given metallicity, EW(Hδ_A) instead shows a peak at an age of a few hundred Myr. The stellar population age corresponding to the peak EW(Hδ_A) decreases with metallicity. Rather than having a different timescale for each metallicity, we choose instead a standard timescale of 800 Myr for all metallicities that enables the defined change parameter SFR79 to be reasonably well calibrated at all six metallicities. This however makes the SFR79 calibration dependent on metallicity. The D_n(4000) increases with increasing stellar age at a given metallicity, and increases with increasing metallicity at given stellar age.

Our definition of SFR79 is based on the SFR across two orders of magnitude in timescale, which is much larger than that of the widely used "burstiness" based on Hα-to-UV ratio.

Furthermore, as already noted, the diagnostic parameters based on equivalent widths are insensitive to the dust attenuation and can be readily measured from the existing large body of optical spectra of galaxies, and are not strongly model-dependent. These conditions make the three diagnostic parameters to be an ideal choice to study the variability of SFR in galaxies.

In parametric SED modeling, strong priors are usually imposed on the SFHs. One of the widely used ones is the exponentially declining SFH (e.g., Bruzual 1983; Papovich et al. 2001; Shapley et al. 2005; Pozzetti et al. 2010; Carnall et al. 2019), i.e., the SFR is assumed to decline exponentially with some e-fold timescale τ: $\mathrm{SFH}(t)\propto \exp (-t/\tau )$. However, it is clear in the real universe, that the SFH of galaxies may be much more complicated than any assumed analytic formula. Motivated by the fact that the global SFD is well fit by a log-normal in time, Gladders et al. (2013) proposed that the log-normal form might also characterize the SFHs of individual galaxies (Dressler et al. 2013; Oemler et al. 2013; Abramson et al. 2015). Using the Illustris simulation, Diemer et al. (2017) investigated the SFHs for individual galaxies, and found that the log-normal form fits the overall shape of the majority of SFHs very well: 85% of cumulative SFHs are fitted to within a maximum error of 5% of the total stellar mass formed. The log-normal works systematically better than the commonly used exponentially declining model, and appears to be a reasonably good description for the global shape of SFHs for individual galaxies. Therefore, we adopt the log-normal fits of the SFHs for Illustris galaxies (Diemer et al. 2017) as being representative of the global shape for the long-term variation of SFHs in the real universe.

On top of these, smoothly varying underlying SFH must be added short-term stochastic variations. We describe the stochastic variations in SFR in the frequency (or time) domain using the power-spectrum distribution (PSD) of variations in the SFR. To construct the mock SFH, we use, for simplicity and following the work of Caplar & Tacchella (2019), a broken power-law PSD to characterize the possible variations in SFR that are superposed on the broad underlying log-normal SFH. The PSD can be written as

$$\begin{eqnarray}&&\mathrm{PSD}(\nu )=\displaystyle \frac{{\sigma }^{2}}{1+{\left({\tau }_{\mathrm{break}}\nu \right)}^{\alpha }},\end{eqnarray} \tag{ 1 }$$

where ν is the frequency, σ defines the amplitude of the PSD, α is the slope of PSD at the high frequency end, and τ_break defines the break point where the PSD becomes flat toward lower frequency. We refer the reader to Caplar & Tacchella (2019) for more details of the properties of this kind of PSDs.

We then use a public IDL code³ to generate the random time series of variation in SFR with a given power-spectrum distribution. Note that the variations of SFR are generated in logarithmic space. Figure 2 shows examples of the variations of stochastic components with different α and τ_break. Here, we only show the variations with a time range of 2 Gyr, while in practice we generate the variation in time series over the full lifetime of the universe. Following the work of Caplar & Tacchella (2019), the 1σ scatter of the variations are normalized to 0.4 dex, which is comparable to the maximum scatter of SFMS in the observations (e.g., Whitaker et al. 2012; Speagle et al. 2014; Schreiber et al. 2015; Davies et al. 2019). The time resolution is set to be 1 Myr, which is much smaller than the SFH timescale traced by Hα, and smaller than the freefall timescale of molecular clouds (Murray et al. 2010; Hollyhead et al. 2015; Freeman et al. 2017). As shown in Figure 2, larger α and longer τ_break result in slower oscillations and stronger correlation in time domain.

Zoom In Zoom Out Reset image size

In total, there are 29,203 galaxies in the Illustris simulation with stellar mass greater than 10⁹M_⊙ (Diemer et al. 2017). We do not yet exclude the quenched galaxies but will do so later in Section 2.4 based on the EW(Hα) and EW(Hδ_A) of the mock spectra. For each of these 29,203 galaxies, with a given underlying log-normal SFH, we then construct 100 stochastic variations of the SFH by varying both α and τ_break in logarithmic space, Δ_SFH (see examples in Figure 2). The mock SFHs are then constructed by multiplying the log-normal SFH from Illustris galaxies with a factor ${10}^{{{\rm{\Delta }}}_{\mathrm{SFH}}}$. We therefore have 2,920,300 mock SFHs in total. In practice, we make a 10 × 10 grid for α and τ_break, where the two parameters are evenly spaced with α in the range of 1–3, and τ_break in the range of 100–1000 Myr. The ranges of the two parameters are chosen according to the result of Caplar & Tacchella (2019). Further, we find that the constrained slope of PSD is within this range in the second paper of this series. For each point on this grid, we generate a stochastic variation of SFH based on its α and τ_break according to the approach above. We stress that our purpose is not to try to model or reproduce the stochastic variation in SFHs in the real universe but is instead to simply generate a huge range of SFHs, which should cover the range of SFR79 that will be encountered in normal galaxies in the real universe.⁴

For each of the 2.9 million mock SFHs constructed above, we calculate the change parameter SFR79 at the present epoch, i.e., the simple ratio of the SFR averaged over the last 5 Myr to that averaged over the last 800 Myr. Using the mock SFH and the SSP models, we can obtain the mock spectrum of the composite stellar population produced by each mock SFH by convolving the time-varying spectrum of the SSP (at a given metallicity) with the detailed age distribution of each mock SFH. We can further compute the three diagnostic spectral parameters for each mock spectrum at the present epoch, for each of the six metallicities.

In practice, we do not of course need to produce an entire high-resolution composite spectrum but simply measure the relevant input fluxes (or flux deficits) of the SSP once as a function of age from its evolving spectrum and then produce the diagnostic parameters for each mock SFH at the present epoch through a straight convolution of these functions with the age distribution of the SFH in question. For instance, the D_n(4000) is defined as the ratio of the flux density between the 4000 and 4100 Å (f_red) and that between 3850 and 3950 Å (f_blue) (Balogh et al. 1999). We first compute the evolution of f_red(t) and f_blue(t) with time for an SSP at a given metallicity. Then, the f_red (or f_blue) for a given mock SFH is the convolution of the corresponding age distribution n(t) with the f_red(t) (or f_blue(t)) evolution curve. In the similar way, the EW(Hα) and EW(Hδ_A) can also easily be obtained. The bandpasses for calculating D_n(4000) and EW(Hδ_A) are defined in Balogh et al. (1999). Specifically, the blue and red bandpass of wavelength (in Å) in calculating the D_n(4000) are [3850, 3950] and [4000, 4100]. The three bandpasses for the index of Hδ absorption are [4083.50, 4122.25], [4041.60, 4079.75], and [4128.50, 4161.00]. In calculating emission (or absorption) line flux of Hα (or Hδ), the contamination of the absorption (or emission) in Hα (or Hδ) is corrected. We note that the approach in calculating the three diagnostic parameters for the mock spectra is exactly the same as that used in analysis of the observations in Section 3.1 (also see Wang et al. 2018).

2.4.1. SFR79 as a Function of EW(Hα) and EW(Hδ_A)

Now that we have the measurements of SFR79 as well as the diagnostic parameters for millions of mock SFHs, it is straightforward to search for the solution of SFR79 as a function of the three diagnostic parameters. Figure 3 shows the ${\mathrm{log}}_{10}$ EW(Hα) versus EW(Hδ_A) diagram with the color-coding of SFR79 for the six different metallicities. The SFR79 shows clear gradients on this diagram, confirming that the EW(Hδ_A) and ${\mathrm{log}}_{10}$ EW(Hα) indeed contain information on the change parameter. For all metallicities, at fixed EW(Hδ_A), the SFR79 increases with increasing ${\mathrm{log}}_{10}$ EW(Hα); at fixed ${\mathrm{log}}_{10}$ EW(Hα), the SFR79 decreases with increasing EW(Hδ_A). This is as expected, since the two parameters indicate the strength of star formation at two different timescales. Another interesting feature is that the range of EW(Hδ_A) becomes smaller toward smaller metallicity. This is due to the different evolution curves of EW(Hδ_A) for the SSP models at the different metallicities (see Figure 1). For the lowest metallicity (${\mathrm{log}}_{10}\,Z/{Z}_{\odot }=-1.0$), the EW(Hδ_A) of the SSP is greater than zero over the entire age range of 10 Gyr after a single starburst. This is the reason why there are no data points with EW(Hδ_A) below zero at the lowest metallicity. Note that EW(Hδ_A) is here defined to be positive for absorption.

Zoom In Zoom Out Reset image size

After exploring many kinds of combinations of the three parameters, we find that a combination of polynomials to the third-order can well reproduce the values of SFR79 to within a scatter of 0.06–0.09 dex. The form of the polynomials can be written as

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}\mathrm{SFR}79 & = & a1\times x+a2\times {x}^{2}+a3\times {x}^{3}\\ & & +b1\times y+b2\times {y}^{2}+b3\times {y}^{3}\\ & & +c1\times z+d,\end{array}\end{eqnarray} \tag{ 2 }$$

where x = ${\mathrm{log}}_{10}$ EW(Hα), y = EW(Hδ_A), z = D_n(4000), and a1, a2, a3, b1, b2, b3, c1, and d are parameters determined from the fittings. The fitting parameters for different metallicities are listed in Table 1. During the fittings, we exclude mock spectra with extremely low EW(Hα) and EW(Hδ_A), since these are not encountered in the star-forming galaxies that are of interest in this work. The exclusion thresholds of these two parameters at different metallicities are listed in the last two column of Table 1. As shown in Table 1, the threshold of EW(Hα) is 1 Å for all metallicities, while the threshold of EW(Hδ_A) increases with decreasing metallicity, varying from −1 to 1.5 Å. The exclusion of quenched galaxies with very small EW(Hα) and EW(Hδ_A) is immaterial for the present purposes. We also note that very small values of these two parameters would be associated with relatively large observational uncertainties anyways.

Table 1.  Fitting Parameters of the Calibrator at Different Metallicities in Figure 4

log Z/Z⊙	a1	a2	a3	b1	b2	b3	c1	d	Scatter	EW(Hα)	EW(HδA)
0.0	1.082	−0.06909	−0.01134	−0.2922	0.01167	−5.999e−05	−1.268	1.126	0.063	>1.0 Å	>−1.0 Å
−0.2	1.014	0.006220	−0.03168	−0.2958	0.01855	−7.570e−04	−0.8483	0.5110	0.065	>1.0 Å	>−1.0 Å
−0.4	1.006	0.04190	−0.03518	−0.2576	0.01338	−5.561e−04	−0.3489	−0.3272	0.062	>1.0 Å	>0.0 Å
−0.6	0.9310	0.1042	−0.04182	−0.3572	0.03650	−0.002370	0.1615	−0.8405	0.071	>1.0 Å	>0.5 Å
−0.8	0.9319	0.1067	−0.03414	−0.4048	0.04930	−0.003331	0.7171	−1.543	0.080	>1.0 Å	>1.0 Å
−1.0	0.8920	0.1293	−0.03305	−0.4694	0.06127	−0.004017	1.244	−2.091	0.090	>1.0 Å	>1.5 Å

Download table as:  ASCIITypeset image

Figure 4 compares the SFR79 as directly measured from the mock SFHs with the results obtained from Equation (2) using the three diagnostic parameters. As shown, even though we included a huge range of possibilities in the mock SFH construction, the SFR79 can be very well determined by the combination of these three diagnostic parameters with a scatter of 0.06–0.09 dex. The scatter shows very little dependence on SFR79 for almost all of the metallicities examined. Another interesting feature is that the scatter becomes larger with decreasing metallicity. This is due to the fact that we use a fixed timescale of 800 Myr to define the change parameter for all of the metallicities. In principle, we could have increased the averaging timescale to ∼1 Gyr to reduce the scatter in the calibrator at low metallicities. However, this variable timescale would introduce more complexity in analyzing the results in Sections 4 and 5. We therefore decided to keep the timescale constant for the different metallicities. The chosen timescale of 800 Myr was selected to minimize the scatter of calibrators at the three highest metallicities (log₁₀ Z/Z_⊙ = 0.0, −0.2, and −0.4), in which most of our sample galaxies are in fact located (see Section 2.4.2).

Zoom In Zoom Out Reset image size

2.4.2. The Performance of the Calibrator

We showed in the above subsection that SFR79 could be calibrated with an overall uncertainty of 0.06–0.09 dex. In this subsection, we examine the performance of the calibrator in the ${\mathrm{log}}_{10}$ EW(Hα) versus EW(Hδ_A) diagram, to examine the performance of the SFR79 estimator in different regions of the diagram.

We show the ${\mathrm{log}}_{10}$ EW(Hα) versus EW(Hδ_A) diagram for the mock SFHs with the color-coding of the ΔSFR79 in the top group of panels in Figure 5. The ΔSFR79 is defined as the difference between the SFR79 from Equation (2) and the true SFR79 in the mock SFHs. At each metallicity, we also present the distribution of the spaxels from the MaNGA SF galaxies (of the corresponding metallicity, with bins of width of 0.2 dex) on the ${\mathrm{log}}_{10}$ EW(Hα) versus EW(Hδ_A) diagram, shown in black contours. The stellar metallicity of the sample galaxies is taken from the empirical mass–metallicity relation from Zahid et al. (2017) (see Equation (4) in Section 3.1). Most of the sample galaxies are in the three highest metallicity bins, and there are no contours in the two lowest metallicity bins. This is because the number of spaxels in the two lowest metallicity bins are quite limited (only one galaxy in the sample has a metallicity in each of these two lowest bins).

Zoom In Zoom Out Reset image size

As shown, the calibration formula can indeed give an excellent estimation of SFR79. However, we note that the estimator does not work well when EW(Hδ_A) and EW(Hα) lie beyond the threshold criteria of the calibrator (see the shaded regions of Figure 5 and also Table 1). In addition, although our mock SFHs cover a very wide range on the diagram of ${\mathrm{log}}_{10}$ EW(Hα) versus EW(Hδ_A), some regions of parameter space are still not covered (the white regions in Figure 5). The calibration polynomial is clearly not valid for the data points that are beyond the colored regions. This limitation does not affect our application to MaNGA galaxies, because almost all of the spaxels of the SF galaxies in MaNGA are within the regions of validity for the estimator (see the black contours).

At the three lowest metallicities, the estimator appears to systematically underestimate the SFR79 by up to ∼0.1 dex at the high EW(Hδ_A) end (see yellow and red colors at high EW(Hδ_A) in Figure 5). This may be due to the fact that at the edge of the parameter space, the number density of mock SFH is relative low, which contribute a small weight during the fitting. However, this systematic deviation in the calibrator is not a problem, since it can be easily corrected based on the position of the ${\mathrm{log}}_{10}$ EW(Hα)–EW(Hδ_A) diagram in the application. We perform this correction by applying the calibrator to MaNGA galaxies. We note that this correction is very minor, since in the observation, almost all of the data points from MaNGA located in the regions that the calibrator operates rather well.

In principle, we could of course simply establish a look-up table for SFR79 over the whole range of the three diagnostic parameters, so as to avoid the above correction. However, in this work, we prefer to present a simple empirical formula of SFR79 that can be used by readers.

We show the scatter in ΔSFR79 on the ${\mathrm{log}}_{10}$ EW(Hα) versus EW(Hδ_A) diagram in the bottom group of panels in Figure 5. As a whole, for all of the metallicities explored, the scatter of the calibrator is small at the high end of both ${\mathrm{log}}_{10}$ EW(Hα) and EW(Hδ_A), and increases toward the lower end of ${\mathrm{log}}_{10}$ EW(Hα) and EW(Hδ_A). This may be due to the fact that for high ${\mathrm{log}}_{10}$ EW(Hα) and EW(Hδ_A) (corresponding to high recent SFRs), the contribution from the older stellar populations in the measurement of these two parameters is correspondingly small. With increasing metallicity, the scatter decreases as a whole, which is likely due to the timescale of 800 Myr in definition of the change parameter for all of the metallicities discussed above.

Based on the bottom panels of Figure 5, we assign an uncertainty in the SFR79 according to the location on the ${\mathrm{log}}_{10}$ EW(Hα) versus EW(Hδ_A) diagram in the application of the calibrator. While the scatter of SFR79 in the mock SFHs may over-estimate the uncertainty of SFR79 for an individual galaxy because the mock SFHs may contain not found ones in the real universe, we think that this is a more reasonable approach than simply assigning a constant uncertainty of SFR79 at a given metallicity based on Figure 4. Since this uncertainty is due to the estimator itself, we refer to this uncertainty as "model uncertainty." This is distinct from the uncertainty due to measurement uncertainties of the three diagnostic parameters from the observations (see Section 4.1).

In deriving the calibration for SFR79, we assumed for each mock galaxy the same, non-evolving, IMF of Chabrier (2003). This assumption may not be the case in the real universe. In principle, a time-varying IMF and a time-varying SFH are deeply degenerate, since we have no way of knowing when stars below the stellar Main Sequence turn-off stars produced. We do not consider this issue further in the current work, but this should not be a problem unless the cosmic evolution of the IMF was significant in the last ∼1 Gyr (the timescale used in definition of the change parameter), which we consider unlikely.

However, the IMF may also be different from galaxy to galaxy, or even in different parts of the same galaxy. Based on a very sensitive index of the IMF, ¹³CO/C¹⁸O, Zhang et al. (2018) found evidence of a top-heavy stellar IMF (with respect to Chabrier IMF) in the dusty starburst galaxies at redshift ∼2–3.

In this subsection, we examine the stability of our calibrator with different IMFs, choosing to do this, for simplicity, only at a single (solar) metallicity. Based on the same approach in Sections 2.3 and 2.4, we construct a new calibration of SFR79 using a Salpeter IMF (Salpeter 1955) without changing the other settings. Then we compare the two calibrators by applying them to a sample of MaNGA galaxies with stellar mass greater than 10¹⁰M_⊙. The three diagnostic parameters of MaNGA galaxies are calculated based on the spectra binned within the effective radius. Details of the binning scheme are in Section 3.4. Note that the three diagnostic parameters are corrected for the dust attenuation, according to the approach in Section 2.6. As shown in the top panel of Figure 6, the change of IMF gives the same result but with a systematic offset of about 0.05 dex. It is to be expected that a change in the IMF produces a systematic offset in SFR79 because both IMF and SFR79 will change the relative number of stars of different masses.

Zoom In Zoom Out Reset image size

In our work, the absolute value of SFR79 is of less interest than the dispersion in SFR79. Therefore, we argue that the choice of IMF is, at least within the range of plausible possibilities, not important.

We next examine the stability of the SFR79 calibration with respect to the use of different stellar evolution model, i.e., the isochrones. In a similar way, we generate a new calibrator by adopting the MIST isochrone (e.g., Paxton et al. 2011; Choi et al. 2016; Dotter 2016; Paxton et al. 2018) with all other settings unchanged. The result is shown in the bottom panel of Figure 6. Again, we see that the use of MIST isochrones introduces a small systematic offset of −0.09 dex.

While these effects introduce systematic offsets to SFR79, this will not affect the investigation of the temporal variation of SFR in galaxy populations, because it is the scatter of SFR79 that characterizes this variation rather than the average value (Broussard et al. 2019). A significant problem would occur only if the IMF or the appropriate isochrones varied significantly from galaxy to galaxy. While the former is possible, the latter is presumably not. We note at this point that the variation with IMF (systematic offsets of 0.05 dex) is rather small compared with the observed dispersion in SFR79 within the population, which is 0.23 dex (see Section 4 below). So, we can assume that any variations in IMF are probably a negligible contributor to this scatter.

As described above, we did not include the effects of dust attenuation in the calibration of SFR79. This means that any correction for dust attenuation correction must be made to the three observational parameters before feeding them into the SFR79 calibration. As already noted, the use of equivalent widths makes them relatively insensitive to extinction. However, differential extinction between stars of different ages, or between line and continuum emission is more of a problem. Newly formed stars (<10 Myr) may have different extinction than older stars, since the extinctions of stellar continuum and of nebular emission are usually different (e.g., Calzetti et al. 2000; Moustakas & Kennicutt 2006; Wild et al. 2011; Hemmati et al. 2015). Charlot & Fall (2000) proposed a two-component dust model, where the optical depth for stellar populations older than 10 Myr is around one-third of the optical depth of stellar populations younger than 10 Myr. In this model, the regions of young stellar populations are more dusty than the regions of older stellar populations, because the newly formed stars are embedded in molecular clouds. The demarcation timescale of 10 Myr comes from the timescale of disruption of molecular clouds (Blitz & Shu 1980; Conroy et al. 2009; Murray et al. 2010). In this work, we make the same assumption of this stellar age-dependent dust model, and we adopt the Cardelli–Clayton–Mathis (CCM) dust attenuation curve (Cardelli et al. 1989).

With this dust model, we correct the dust attenuation of the three diagnostic parameters in the following way. First, we generate as before mock SFHs based on the SFHs of Illustris galaxies adding stochastic processes, but here adopt a broken power-law PSD of SFHs with α = 1.5 and τ_break = 20 Gyr (The large τ_break make the PSD close to a single power-law PSD with α = 1.5). And again, the scatter of the variations are normalized to 0.4 dex. In contrast to the situation in Section 2.3, we here try to generate SFHs that resemble the observations, rather than a huge range of all possible SFH. A power-law PSD with α = 1.5 is likely a good description of the stochasticity of SFHs without considering the intrinsic scatter of the SFMS, according to the analysis of our second paper of this series. Actually, in the second paper, we will find that, if we assume a single power-law form for the PSD of the specific SFH, a slope of 1.5 is the best to reproduce the distribution of galaxies on the ΔSFR7–ΔSFR9 plane.⁵ The slope of PSD adopted here is slightly shallower than the one assumed in Caplar & Tacchella (2019), who assumed a PSD index of 2, corresponding to a random walk process. We refer the reader to the second paper of this series for details.

We then calculate two sets of diagnostic observational parameters based on the mock SFHs with and without the dust attenuation. For the mock spectra, we know the contribution of the stellar populations at different ages, from the mock SFHs and evolving spectra of the SSP models. This enables us to apply a stellar age-dependent extinction model to obtain the reddened spectra and the values of the three diagnostic parameters. The diagnostic parameters with dust reddening are denoted as D_n(4000)_dust, EW(Hδ_A)_dust and EW(Hδ_A)_dust. In this process, we assign for each galaxy an E(B − V)_young,⁶ i.e., the E(B − V) for the stellar population younger than 10 Myr, that is randomly distributed between 0.0 and 1.0. At last, we compare the two set of observed parameters and define recipes for the correction.

After exploring many forms for the dust correction of the three observational parameters, we found good expressions with rather small uncertainties. Figure 7 presents the formulae to correct the dust attenuation (solid red lines) by comparing the two set of diagnostic parameters at solar metallicity. The correction of EW(Hα) and D_n(4000) are a function of E(B − V)_young, while the correction of EW(Hδ_A) is a function of both E(B − V)_young and EW(Hα). In the observation, the E(B − V)_young can be determined from the flux ratio of Hα and Hβ using the Balmer decrement. The adopted expressions for the correction for the three parameters are as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}\mathrm{EW}({{\rm{H}}}_{\alpha })/\mathrm{EW}{\left({{\rm{H}}}_{\alpha }\right)}_{\mathrm{dust}} & = & p1\times E{\left(B-V\right)}_{\mathrm{young}}\\ \mathrm{EW}({\rm{H}}{\delta }_{{\rm{A}}})-\mathrm{EW}{\left({\rm{H}}{\delta }_{{\rm{A}}}\right)}_{\mathrm{dust}} & = & p2\times E{\left(B-V\right)}_{\mathrm{young}}^{p3}\\ & & \times (\mathrm{EW}({{\rm{H}}}_{\alpha })+p4)\\ {\mathrm{log}}_{10}{D}_{n}{\left(4000\right)}_{\mathrm{dust}}/{D}_{n}(4000) & = & p5\times E{\left(B-V\right)}_{\mathrm{young}},\end{array}\end{eqnarray} \tag{ 3 }$$

where p1, p2, p3, p4, and p5 are parameters determined by fitting. Table 2 shows these fitting parameters to be used in Equation (3), as well as the resulting uncertainties of the observational diagnostic parameters (listed in the last three columns in Table 2) produced by the dust attenuation correction, at different metallicities. The uncertainties are determined by the scatters in the three panels of Figure 7. As can be seen in Table 2, the uncertainties due to this correction are small.

Zoom In Zoom Out Reset image size

Table 2.  Fitting Parameters of the Dust Attenuation for Three Diagnostic Parameters in Figure 7

${\mathrm{log}}_{10}Z/{Z}_{\odot }$	p1	p2	p3	p4	p5	${\mathrm{log}}_{10}$ EW(Hα)ERR	EW(HδA)ERR	Dn(4000)ERR
0.0	0.685	−0.00582	0.580	−4.29	0.0382	0.01 dex	0.098 Å	0.019
−0.2	0.687	−0.00691	0.569	8.23	0.0345	0.01 dex	0.084 Å	0.016
−0.4	0.689	−0.00595	0.576	12.5	0.0318	0.009 dex	0.072 Å	0.014
−0.6	0.691	−0.00496	0.564	14.4	0.0306	0.008 dex	0.063 Å	0.014
−0.8	0.693	−0.00407	0.565	15.1	0.0287	0.007 dex	0.052 Å	0.012
−1.0	0.692	−0.00388	0.563	16.4	0.0275	0.007 dex	0.050 Å	0.011

Download table as:  ASCIITypeset image

However, there may be additional uncertainties in the dust attenuation correction, due to possible variations in the f-factor (e.g., Kashino et al. 2013; Faisst et al. 2019; Lin & Kong 2020), defined as the E(B − V)_star/E(B − V)_nebular⁷ , from galaxy to galaxy or from regions to regions. This effect is not included in our dust model in the present work. Specifically, by using the MaNGA galaxies, Lin & Kong (2020) found that the f-factor decreases with increasing stellar mass, and slightly increases with increasing sSFR. For most of SF spaxels in galaxies above 10¹⁰M_⊙, the f-factor is in the range of 0.3–0.7 with a scatter of ∼0.1–0.15. We have examined the dependence of our SFR79 estimator on the value of E(B − V)_old/E(B − V)_young at solar metallicity. Increasing E(B − V)_old/E(B − V)_young by 0.1, the resulting SFR79 show an overall 0.03 dex offset with respect to the old ones, which is much smaller than the scatter of SFR79 (0.23 dex) we measured in Section 4. We conclude that the dust attenuation corrections of the three parameters are only a secondary effect due to the fact that they are relative values measured at fixed wavelength. The uncertainty of the dust attenuation is even much smaller than the value of the applied correction for the three parameters, and is therefore not likely to be a big concern.

MaNGA is one of the largest integral field spectroscopic surveys, aiming at obtaining the two-dimensional spectra for ∼10,000 galaxies in the redshift range of 0.01 < z < 0.15 (Bundy et al. 2015). The wavelength covered by MaNGA is 3600–10300 Å at a spectral resolution R ∼ 2000, which is sufficient to accurately measure the three diagnostic parameters (Li et al. 2015) used in this paper.

In this work, we utilize the well-defined sample of star-forming galaxy from W19. Here we therefore only briefly describe the sample definition, and refer the reader to W19 for further details.

The galaxy sample is originally selected from SDSS Data Release 14 (Abolfathi et al. 2018), excluding the mergers, irregulars, heavily disturbed galaxies, as well as galaxies for which the median S/N of the 5500 Å continuum is less than 3.0 at their effective radii. The quenched galaxies are excluded based on the stellar mass and SFR diagram. The stellar mass and SFR are measured within the effective radius, i.e., M_*(≤R_e) and SFR(<R_e),⁸ based on the MaNGA spectra. Our final sample consists of 976 SF galaxies, and is a representative sample of SF main-sequence galaxies in the low-redshift universe.

The stellar mass maps of MaNGA galaxies are derived from the public fitting code STARLIGHT (Cid Fernandes et al. 2004), using SSPs with Padova isochrones from Bruzual & Charlot (2003) and the Chabrier (2003) IMF. The SFR maps are determined by the extinction-corrected Hα luminosity adopting the conversion formula to SFR from Kennicutt (1998), again using a Chabrier (2003) IMF. The uncertainty in determining the SFR via this approach is 15% (or ∼0.06 dex), due to the variations in the electron temperature in the range T_e = 5000–20,000 K (Osterbrock & Ferland 2006). As above, we also refer to this uncertainty as model uncertainty.

Since, in this work, our purpose is primarily to investigate the SFR of galaxies on different timescales, we denote the SFR directly determined from the Hα luminosity as SFR7, i.e., SFR_{5 Myr} (see Figure 1). The intrinsic extinction for nebular emission is measured based on the Balmer decrement, assuming the CCM dust attenuation curve (Cardelli et al. 1989) and Case B recombination with ab intrinsic flux ratio of Hα/Hβ = 2.86. The E(B − V) for nebular emission is then a good estimation for the E(B − V)_young, i.e., the color excess for the stellar population younger than 10 Myr. We note that the IMF, isochrones, and dust attenuation curve that are used to obtain the stellar masses and SFR_{5 Myr} are consistent with those used in Section 2.

The strengths of the emission lines are measured based on the stellar continuum-subtracted spectrum by fitting a Gaussian profile to the lines. The D_n(4000) and EW(Hδ_A) are directly measured based on the observed spectra after subtracting emission lines, rather than from the best-fit continuum spectra. This avoids the uncertainties of the measurements due to the possible systematic offset (especially at 3800–4200 Å) between the model spectra and observed ones. For many SF galaxies, the bottom of the Hδ absorption is usually accompanied by weak Hδ emission, which make the EW(Hδ_A) difficult to measure. In this work, we take advantage of the χ² minimization spectral fitting code developed by Li et al. (2005), which can effectively mask the emission-line regions iteratively during the fitting. This is critical to accurately model the absorption troughs and also characterize the emission lines (see examples in Li et al. 2015).

Before applying our SFR79 estimator on the MaNGA galaxies, we first correct the three diagnostic parameters for dust attenuation based on Equation (3). Since both the estimator and the dust attenuation correction depend on the stellar metallicity of galaxies, we adopt the stellar mass–metallicity relation from Zahid et al. (2017) to estimate the stellar metallicity of individual galaxies. The relation can be written as:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\displaystyle \frac{Z}{{Z}_{\odot }}={Z}_{0}+{\mathrm{log}}_{10}\left[1-\exp \left(-{\left[\displaystyle \frac{{M}_{* }}{{M}_{0}}\right]}^{\gamma }\right)\right],\end{eqnarray} \tag{ 4 }$$

where Z₀ = 0.075, M₀ = 10^9.79M_⊙, and γ = 0.56. This relation is determined by modeling the galaxy spectra with a linear combination of sequential single burst model spectra. For a given set of three observational parameters at given metallicity, we first calculate the correction based on Equation (3) and Table 2 at the two closest metallicities. Then we use linear interpolation in log₁₀ Z/Z_⊙ to obtain the corrections of the observational three diagnostic parameters at the required metallicity. In the similar way, we then obtain the SFR79 using the (dust-corrected) values of the three diagnostic parameters by calculating the SFR79 at the two closest metallicities based on Equation (2) and Table 1 then obtain the SFR79 at the required metallicity via linear interpolation in log₁₀ Z/Z_⊙.

Having calculated SFR79 for all spaxels in our MaNGA sample, we can now carry out an important consistency check by calculating the total SFR of all spaxels in the sample, averaged over the last 5 Myr, and the total SFR averaged over the last 800 Myr. These should be roughly equal. To be precise, the ratio of these, i.e., $\langle \mathrm{SFR}7\rangle /\langle \mathrm{SFR}9\rangle$, should reflect the overall cosmic evolution of the SFR of the SF galaxy population over the last Gyr, i.e., the change in overall SFR of SF galaxies that is implied by integrating the change in the sSFR of the Main Sequence. The cosmic evolution of ${\mathrm{log}}_{10}\langle \mathrm{SFR}7\rangle /\langle \mathrm{SFR}9\rangle$ that is expected for the ensemble of main-sequence galaxies is calculated to be −0.025 dex based on the evolution of the sSFR of the SFMS from Lilly & Carollo (2016), and the stellar mass function for SF galaxies from Peng et al. (2010).

The ratio of these two total SFR in the MaNGA data (using our estimator of SFR79 to calculate individual SFR9) is actually −0.066 dex, which is reassuringly close (within 0.04 dex, or 10%) to the expected value. This very satisfactory agreement should be taken as a first confirmation that our calibration of SFR79 is quite accurate and certainly usable. In fact, given the assumptions that we made, however reasonably, about the effect of reddening, about the form of the IMF and about the choice of stellar evolution isochrones, the very close agreement to within 0.04 dex should probably be seen as fortuitous. This is explored further in the next subsection.

The top panel of Figure 8 shows the SFR79 for all of the spaxels in our sample galaxies as a function of deprojected stellar mass surface density (Σ_*). The contours show the overall number density distribution of spaxels in this diagram. These contours enclose 30%, 50%, 70%, and 90% of all spaxels, from the inside out. The blue dashed line shows the median SFR79 at given Σ_*. Based on SFR79, we can obtain log₁₀ SFR9, i.e., ${\mathrm{log}}_{10}{\mathrm{SFR}}_{800\mathrm{Myr}}$, for each individual spaxel of the sample galaxies.

Zoom In Zoom Out Reset image size

It is clear in this figure that the mean SFR79 appears to slightly increase with increasing Σ_*, suggesting that SF galaxies have on average a slightly negative SFR79 radial gradient. We suspect that this small effect is unlikely to be real, since it would suggest that the SFR was declining more at large radii, leading to a negative gradient in the rate of change of the sSFR of galaxies. Since galaxies generally have a positive gradient in sSFR at the current epoch, this would imply that this gradient was weakening. If anything, we would expect the opposite trend in any "inside-out" scenario of galaxy evolution. There are other reasons to question this small gradient in SFR79. Not least, although the dust correction is in principle computed locally, we did not consider any radial variation of the form of the dust correction nor of the stellar metallicity (e.g., Goddard et al. 2017; Zheng et al. 2017) nor any (possible) radial variation of IMF (Gunawardhana et al. 2011) in the estimation of SFR79. We suspect that some combination of these may be the cause of the trend in the upper panel of Figure 8.

Accordingly, we therefore perform a small ad hoc correction to the values of SFR79 as a function of Σ_* (only). This is constructed so as to eliminate the dependence of SFR79 on Σ_*, while not significantly perturbing the total rates of star formation, as discussed in the previous subsection. We first fit the median SFR79–Σ_* relation with a piecewise linear function, shown as the red line in Figure 8. The red line is in the form of

$$\begin{eqnarray}y=\left\{\begin{array}{ll}0.04581x-0.6089, & x\lt 7.5,\\ 0.1737x-1.568, & x\geqslant 7.5,\end{array}\right.\end{eqnarray} \tag{ 5 }$$

where y = log₁₀ SFR79, and x = log₁₀ Σ_* (M_⊙ kpc⁻²). For each spaxel, we then subtract this median SFR79 at the corresponding Σ_*. In order to preserve the total star formation rates, as discussed in the previous subsection, we then add a uniform −0.136 dex to the SFR79 computing this value to exactly match the ${\mathrm{log}}_{10}\langle \mathrm{SFR}7\rangle /\langle \mathrm{SFR}9\rangle$ with the value of the cosmic evolution for SF main-sequence galaxies (i.e., including the 0.04 dex offset discussed in the previous subsection). The bottom panel of Figure 8 shows the SFR79 as a function of Σ_* after this correction.

This correction to SFR79 is of course quite arbitrary. However, we stress that this correction only makes the overall SFR79 profile flat and does not change the scatter of SFR79 at given Σ_*. Since the variability of SFHs is indicated by the scatter of SFR79 across the population, rather than by its average or median value, this correction will not significantly affect any of our conclusions in the following analysis. In the remainder of this work, the SFR79 for both individual spaxels and individual galaxies will be corrected according to the above approach.

Based on the approach in Section 3, we can now obtain the SFR79 maps of each sample galaxy. Figure 9 shows an example of the measurements for one typical galaxy (MaNGA-ID: 8249-6102). The top panels of Figure 9 shows the SDSS color image, D_n(4000), EW(Hδ_A), and ${\mathrm{log}}_{10}$ EW(Hα) maps from left to right, respectively. It should be noted that the three diagnostic parameters are shown prior to the correction for dust attenuation. The bottom panels show the maps of sSFR7, sSFR9, and SFR79, as well as the profile of SFR79 for this galaxy, from left to right, respectively.

Zoom In Zoom Out Reset image size

As shown, the EW(Hα) map shows clumpy features (blue clumps), corresponding to the regions with high recent star formation within 5 Myr. However, these regions do not show high EW(Hδ_A), suggesting that the star formation in these regions were not unusually active in forming stars during the last ∼1 Gyr. Thus, the current star formation in these blue clumps of ${\mathrm{log}}_{10}$ EW(Hα) map is triggered recently (much less than 800 Myr). Consistent with this, these regions are seen as positive (blue) on the SFR79 map that blue clumps are seen in the same regions. On the other hand, some regions with high EW(Hδ_A) but relative low EW(Hα), are visible on the SFR79 map with yellow clumps. The star formation of these regions is reduced in a timescale much less than 800 Myr. This simple example indicates that the SFR79 we measured is indeed meaningful and provides the quantitative description for the above effect.

In the bottom rightmost panel of Figure 9, the small gray circles indicate the SFR79 derived spaxel-by-spaxel in this galaxy, while the red line shows the SFR79 profile obtained from the binned spectra in annular bins. It can be clearly seen that, between 0.5 and 1.0R_e the distribution of SFR79 is asymmetric: most of the spaxels have a relative low SFR79, while a small fraction of spaxels have increased SFR79. This is due no doubt to the duty cycle of star formation. At any given time, star formation in a given region of galaxy is found in a limited number of active regions, which themselves are active for only a short fraction of the time.

In the current work, we do not wish to study the variation of SFR (or SFR79) that is caused by this small scale local effect but rather focus on the variations in SFR on larger scales. Therefore, we use a binning scheme to average out these small scale spatial effects.

The diagnostic parameters for a given radial bin are calculated in the following way. For instance, for D_n(4000), we first calculate the flux density of the blue and red bandpass near 4000 Å as described in Section 2.4 for all of the spaxels (with S/Ns greater than 3 at 5500 Å) in a given radial bin. The D_n(4000) of this radial bin is then obtained from the ratio of the sum of the flux densities in the red bandpass to the sum of them in the blue bandpass, summing over all of the spaxels in this radial annular bin. In this process, we also estimate the observational error in the flux density for the blue and red bandpasses. We first assume the flux errors of all of the binned spaxels have no covariance in obtaining the binned error (σ_no,cov). This error is obviously less than the real error due to the covariance of binned spaxels that arises because the spatial-resolution (∼25) of the MaNGA survey is much larger than the size of each individual spaxel (05). To correct this, we adopt an empirical function from Law et al. (2016):

$$\begin{eqnarray}\displaystyle \frac{{\sigma }_{\mathrm{cov}}}{{\sigma }_{\mathrm{no},\mathrm{cov}}}\,=\,\left\{\begin{array}{ll}1.0+1.6{\mathrm{log}}_{10}{N}_{\mathrm{pixel}} & {N}_{\mathrm{pixel}}\leqslant 100\\ 4.2 & {N}_{\mathrm{pixel}}\gt 100\end{array}\right.\end{eqnarray} \tag{ 6 }$$

where σ_cov is the error with considering the covariance of binned spaxels, and N_pixel is the number of binned spaxels. We then obtain the observational error of D_n(4000) by error propagation. In a similar way, we obtain the binned EW(Hδ_A) and EW(Hα), as well as their error for a given bin. The E(B − V)_young in a given radial bin is also calculated based on the ratio of total Hα flux to total Hβ flux within the bin. Finally, the SFR79 in the bin is calculated with the binned diagnostic observational parameters using Equation (2), and the observational error of SFR79 is calculated based on the observational error of three diagnostic parameters via error propagation.⁹

The advantage of the binning scheme is (1) to improve the accuracy of the measurements of SFR79, (2) to reduce the variation of SFR79 caused by the small scale effect discussed above (including the duty cycle of star formation). We set the radial bin width to be 0.2R_e, large enough to eliminate the scale effect, and calculate the SFR79 and its observational error for a given radial bin via the above approach. In the current paper, we focus on the radial profiles of SFR79 for the sample galaxies, and will not consider the local variation of SFR79 within these radial bins.

We first examine the global SFR79 and SFR9 for the sample galaxies and use these to examine the SFMS when defined using the measures of star formation rate on the two timescales of 5 and 800 Myr. Consistently with the measurements of global SFR7 (i.e., SFR_{5 Myr}) and stellar mass, the global SFR9 is calculated for each individual galaxy by summing the flux in all of the spaxels within the effective radius to obtain the integrated SFR7, SFR79, and thus SFR9 for the sample galaxies.

We exclude 16 galaxies that are located in the Seyfert regions on the Baldwin–Phillips–Telervich (BPT) diagram (Baldwin et al. 1981; Kewley et al. 2006), based on the emission-line flux ratios within the effective radius, since in these galaxies the Hα emission is largely contaminated by the contribution of the AGN. We also exclude four galaxies that have the three diagnostic parameters beyond the valid range of the calibrator (see details in Section 2.4). These leave 956 galaxies.

In the top two panels of Figure 10, we present the SFMS based on SFR7 (left-hand panel) and SFR9 (right-hand panel). In both panels, the blue circles show the median SFR in different stellar mass bins. The black solid line is the best-fit straight line to the median SFR7–M_*(≤R_e) relation for galaxies with M_*(≤R_e) less than 10¹⁰M_⊙. Following the work of W19, we define the solid line as the "nominal" SFMS for SFR7. Interestingly, the line also matches the SFMS with SFR9 very well. This is consistent with the fact that the evolution of the SFMS is small within the last 800 Myr (e.g., Brinchmann et al. 2004; Pannella et al. 2009; Stark et al. 2013; Schreiber et al. 2015) and is a reflection of the fact that our average SFR79 is very close to unity, as discussed above. We therefore adopt the solid line as the "nominal" SFMS for both SFR7 and SFR9.

The typical errors of SFR7 and SFR9, including the uncertainties from the observations and calibrator, are shown in the corner box of each panel in Figure 10. The uncertainty of SFR7 is 0.06 dex due to the conversion formula from Kennicutt (1998). The measurement uncertainty of the Hα luminosity is negligible with respect to the uncertainty in the conversion formula to SFR, and therefore, we do not show it in the top left panel. The typical uncertainty of SFR79 is 0.076 dex, obtained by combining the uncertainty from the calibrator is 0.063 dex, and the measurement error is 0.042 dex. Comparing with the intrinsic scatter within the galaxy population of SFR79 (∼0.23 dex), the measurement and calibration uncertainty of SFR79 only broadens the distribution of log₁₀SFR79 by less than 10%. This indicates that the dispersion of log₁₀SFR79 in the bottom right panel of Figure 10 is real, which provides the basic condition to study the variability of the SFHs.

Comparing the top two panels in Figure 10, it is noticeable that the scatter of the SFMS is much smaller when using SFR9 than when using SFR7. This is to be expected. Averaging the SFR over 800 Myr eliminates the variation of SFH on shorter timescales. It should be noted that, as an extreme example of this, the SFMS would have zero scatter, if the SFR was computed as the average over the age of universe. By examining galaxies from the EAGLE simulation, Matthee & Schaye (2019) also found that the scatter of SFMS becomes smaller when using the SFR averaged over longer timescales.

Based on the "nominal" SFMS, we now define two parameters to quantify the deviation in logarithmic SFR space of each galaxy from the main sequence at its mass, ΔSFR7 and ΔSFR9 (in dex) of the galaxy, i.e., the vertical distance from the "nominal" SFMS (see Figure 10). The bottom left panel of Figure 10 shows the correlation between ΔSFR7 and ΔSFR9.

This plot contains additional important information about how galaxies move above and below the "nominal" SFMS. An extreme scenario in which all SF galaxies evolved parallel to the "nominal" SFMS would produce ΔSFR9 always equal to ΔSFR7, and therefore, galaxies would exactly follow the one-to-one line with zero scatter on the ΔSFR7–ΔSFR9 diagram. We could imagine the opposite extreme case, in which the scatter of the SFMS is purely due to the variation of SFR on very short timescales (≪800 Myr). In this case, we would expect that the ΔSFR9 for all SF galaxies would be close to zero, and therefore, galaxies would lie on a flat sequence with almost zero scatter on the ΔSFR7–ΔSFR9 diagram. For cases in between, galaxies would be located on a sequence with the slope between zero and one. The slope and the dispersion of the sequence on the ΔSFR7–ΔSFR9 diagram therefore indicates the relative contributions to the dispersion of the SFMS on long and short timescales. In the second paper of this series, we will constrain the PSD of the specific SFHs of galaxies based on the location of galaxies on the ΔSFR7–ΔSFR9 diagram. We do not discuss this plot further here and refer the reader to that second paper.

The bottom right panel of Figure 10 shows the relation between ΔSFR9 and SFR79. The dashed and solid red lines show the 16%, 50%, and 84% percentiles of SFR79 for galaxies at different ΔSFR9. The lack of galaxies in the bottom left corner is due to the fact that the definition of our SF galaxies was based on the SFR7-based SFMS (see the dashed line in the top left panel of Figure 10). Galaxies (if any) below the dashed black line in the bottom right panel of Figure 10 would not be included in the sample selection.

As discussed above, the SFR79 parameter determines the position of a galaxy on the SFMS (defined using the SFR over the last 5 Myr), i.e., ΔSFR7, relative to its position on the SFMS that is defined using the star formation averaged over the last 800 Myr, ΔSFR9. This latter quantity will to first order be the average ΔSFR7 over the last 800 Myr. Therefore, apart from the small offset due to the overall evolution of the SFR in star-forming galaxies with cosmic time, the sign of SFR79 indicates whether the galaxy is generally moving up or down in sSFR (i.e., its present position relative to its average position over the last 800 Myr). Galaxies to the right of the vertical dashed line in the lower right panel of Figure 10 are in this sense increasing their sSFR, or "going up," while those to the left are "going down" relative to the SFMS.

We have already commented that the value of ${\mathrm{log}}_{10}\langle \mathrm{SFR}7\rangle /\langle \mathrm{SFR}9\rangle$ is closely matched to the overall evolution of the sSFR of the SFMS. We furthermore here see in the lower right panel of Figure 10 that there is no apparent correlation between SFR79 and ΔSFR9. The absence of a correlation between SFR79 and ΔSFR9 is required if the scatter of the SFMS is to remain more or less constant over cosmic time. A strong positive correlation between ΔSFR9 and SFR79 means that galaxies in the upper (lower) part of the SFMS would tend to move up (down) with respect to the "nominal" SFMS, leading over time to an increased dispersion in the SFMS. Similarly, a strong anticorrelation between ΔSFR9 and SFR79 would lead to a reduced dispersion over time.

A roughly constant scatter of the SFMS is indeed seen over a wide range of cosmic epochs (e.g., Speagle et al. 2014; Whitaker et al. 2014; Schreiber et al. 2015; Barro et al. 2017), requiring that there should be no correlation between SFR79 and ΔSFR9. The fact that this is indeed seen in our SFR79 estimates is an important external consistency check that provides further confirmation that our estimator for SFR79 works well.

Furthermore, it is noticeable that the distribution of SFR79, at given ΔSFR9, is quite symmetrical about the median value. This symmetry is in contrast to the pronounced asymmetry in SFR79 that is visible in the bottom rightmost panel of Figure 9. In broad terms, this symmetry implies that, for an individual galaxy as a whole, the timescales of "above-average star formation" and "below average star formation" are broadly similar. As an example, short periods of highly elevated SFR superposed on longer periods of constant SFR would produce an asymmetric distribution in SFR79 with a peak at slightly negative SFR79 and a tail to high positive SFR79. However, we do not see this in the data (the bottom right panel of Figure 10).

Zoom In Zoom Out Reset image size

We return here to a point touched on earlier, namely that a galaxy with a constant sSFR will have slightly positive SFR79. This is because the SFR increases as the stellar mass increases. We show this in the left-hand panel of Figure 11, which reproduces the data from the bottom right panel of Figure 10. If all galaxies had a fixed position ΔSFR relative to the (slowly evolving) SFMS, and experienced no other variability in their (s)SFR, then they would lie precisely along the solid black line in the left panel of Figure 11, displaced to left and right only by observational scatter in determining SFR79. Because the sSFR of the SFMS is low at the current epoch, i.e., sSFR⁻¹ ≫ 1 Gyr, the SFR79 produced by constant ΔSFR is negligible for SFMS galaxies, essentially because the mass change of the galaxy during one Gyr is negligible.

Zoom In Zoom Out Reset image size

This emphasizes that the observed scatter in SFR79 within the population is completely dominated by time-variability in the SFR (and sSFR) of galaxies on timescales of 1 Gyr or less, and not by a range of (unvarying) sSFR within the population. This is further illustrated in the right-hand panel of Figure 11 in which we plot the current-epoch SFR79 of the 26,485 log-normal SFH of star-forming (sSFR7 > 10⁻² Gyr⁻¹) Illustris galaxies from Diemer et al. (2017) that were discussed above. These log-normal smooth SFHs will by construction not have short-term variability. The SFR79 are computed directly from the SFH, but we add Gaussian observational scatter to the points to simulate the real data. The dispersion in SFR79 of this simulated population of log-normal galaxies is only 0.08 dex (produced almost entirely by the addition of observational uncertainties), much less than the observed dispersion of 0.23 dex. This comparison emphasizes that the much broader scatter in SFR79 in the real MaNGA data is caused by real short-term temporal variations in the SFR of galaxies that are not present in the log-normal SFH fits given by Diemer et al. (2017).

The profiles of Σ_SFR9 for each individual galaxies, shown in Figure 13. For each galaxy, the Σ_SFR9 profile is generated by combining its SFR79 profile (shown in Figure 12) and the Σ_SFR7 profile in Figure 3 of W19. When presenting the Σ_SFR7 profiles in W19, we normalized the radius for each galaxy by dividing by the effective radius of the galaxy. This common procedure removes the effect due to the variation of galaxy size (e.g., González Delgado et al. 2016; Ellison et al. 2018; Medling et al. 2018; Guo et al. 2019). However, more originally, we also normalized the Σ_SFR with (0.2R_e)², i.e., computed the surface density of the SFR per area of (0.2R_e)². This ensures that the integration of a profile on a given surface density–radius diagram reflects the actual integrated quantity in physical terms. In the current work, we therefore present the Σ_SFR9 profiles in the same way.

Zoom In Zoom Out Reset image size

As shown in Figure 13, the Σ_SFR9 profiles vary from galaxy to galaxy. Analogous to the result of Σ_SFR7 profiles in W19, the scatter in Σ_SFR9 across the galaxy population shows a larger scatter in galactic center than the outer regions (<Re), at least for the three highest stellar mass bins. The median profile of Σ_SFR9 (red dots) can be well fitted by a straight line (the black dashed line), indicating a typical exponential star formation disk with no suppression of star formation (or quenching) in the central regions of galaxies. The slopes of the median Σ_SFR9 profile are −0.53, −0.58, −0.65, −0.65 and −0.67 from low to high stellar mass bins, which are almost the same to the slopes of median Σ_SFR7 profile computed in W19 at corresponding mass bins. This agrees well with the notion that star formation within main-sequence galaxies varies in a quasi-steady state within an exponential disk, ignoring other structural properties such as the presence or absence of a bulge.

We split the galaxies in each stellar mass bin into four quartiles of the global ΔSFR9, i.e., the position relative to the SFMS as defined by the SFR over the last 800 Myr. The thresholds of the quartiles for ΔSFR9 are listed in Table 3 for each stellar mass bin. The top set of panels in Figure 14 show the median profiles of Σ_SFR9 for the four quartiles in global ΔSFR9 for each of the five stellar mass bins. For comparison, we also evenly divide galaxies into four subsamples according to the global ΔSFR7, and present their median Σ_SFR7 profiles in the bottom panels of Figure 14. The thickness of the median profile indicates its uncertainty computed by the bootstrap method.

Zoom In Zoom Out Reset image size

It can be seen from Figure 14 that galaxies with higher (or lower) global SFR9 show enhanced (or suppressed) star formation surface density (within the last 800 Myr) at all galactic radii. This elevation (or suppression) in galactic center is more pronounced in more massive galaxies. This result is entirely consistent with and analogous to the result shown in W19 for Σ_SFR7, that galaxies with higher integrated SFR7 show higher SFR7 at all galactic radii.

Given the similar results in SFR7 and SFR9 in Figure 14, one might worry that this is somehow due to SFR9 being a derived quantity based on SFR7 and SFR79, the ratio of SFR7/SFR9. However, it is easy to see that this in fact makes the measurement of SFR9 largely independent of that of SFR7. The former is ultimately linked to the strength of Hδ absorption and the latter to the strength of Hα. Clearly, if the Hα emission due to star formation is over-estimated, for instance, because of the contribution of an AGN component, then both the SFR7 and SFR79 would be more or less equally perturbed. Ultimately, the value of SFR9 is based primarily on the EW(Hδ_A) value. The result in Figure 14 is therefore largely independent of the analogous result in W19.

W19 defined a parameter, ΔΣ_SFR7, to quantify the deviation from the median profile for a given stellar mass. In the same way, we now define the deviation from the median profile of Σ_SFR9 to be ΔΣ_SFR9. Since the median profiles of Σ_SFR7 have almost the same slopes of Σ_SFR9 (and also almost the same intercepts) at all five stellar mass bins, we use the same set of median profiles as in W19 to calculate these deviations. We note that this does not affect the measurement of the scatter of ΔΣ_SFR9 at a given galactic radius.

We show in Figure 15 the scatter of ΔΣ_SFR9 and ΔΣ_SFR7 (the second is taken from the right panel of Figure 9 in W19 for comparison) as a function of gas depletion time. As in W19, the gas depletion time (i.e., the inverse of the star formation efficiency SFE) is calculated based on an empirical formula from Shi et al. (2011), the so-called extended-Schmidt law:

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{SFE}}{{\mathrm{yr}}^{-1}}={10}^{-10.28}{\left[\displaystyle \frac{{{\rm{\Sigma }}}_{* }}{{M}_{\odot }{\mathrm{pc}}^{-2}}\right]}^{0.48}.\end{eqnarray} \tag{ 7 }$$

This formula is derived from integrated observations of individual galaxies over five orders of magnitude in stellar mass density but was found to also be valid for spiral galaxies at sub-kiloparsec resolutions and low-surface-brightness regions. In W19, we also presented another empirical formula for the gas depletion time by using the orbital timescale (Krumholz et al. 2012). Since the two give similar result, in this work we only present the one computed with the extended-Schmidt law. Each data point represents one galactic radial bin within a given galactic stellar mass bin, as in W19. The colored data points are for the radii less than R_e, while gray radii are for the radii greater than R_e. Here, we will only focus on those regions within the effective radius as done also in W19, because the outer regions are more likely affected by environmental effects, such as ram pressure stripping, and etc.

Zoom In Zoom Out Reset image size

Consistent with the result of W19, the scatter in ΔΣ_SFR9 is also evidently a decreasing function of the gas depletion time. The results on the SFR surface density on 800 Myr timescales, Σ_SFR9, presented in this subsection extend and confirm the results on the SFR surface density on 5 Myr timescales, Σ_SFR7, discussed in W19. The basic issue of whether variations in Σ_SFR reflect temporal variations or intrinsic (time-independent) differences from galaxy to galaxy however remains. To resolve this issue, we now turn to examine the SFR79 profiles more directly.

The SFR79 profiles of the galaxies in the sample were shown in Figure 12. We now divide the galaxies in each stellar mass bin into four quartiles according to their global SFR79 as measured within their effective radii in Section 4. The top group of panels of Figure 16 shows the median SFR79 profiles of the four quartiles in overall SFR79 for each stellar mass bin. The thickness of each line reflects the uncertainty of the median SFR79 profile as computed via a bootstrap approach.

Zoom In Zoom Out Reset image size

It can be seen that galaxies with high global SFR79 appear to have larger SFR79 at all galactic radii, and galaxies with low global SFR79 appear to have lower SFR79 at all galactic radii. This means that a galaxy with a certain global change in SFR, i.e., an elevation (or suppression) of the SFR with respect to the average over the last 800 Myr, experiences a corresponding elevation (or suppression) at all galactic radii. Furthermore, it can be seen that, at least for the three more massive bins of galaxy mass, this change (either elevation or suppression) is more pronounced at small galactic radii.

In the bottom panels of Figure 16, we present the analogous SFR79 profiles for the same five stellar mass bins but now dividing galaxies into four quartiles of ΔSFR9, i.e., the position of the galaxies relative to the SFMS defined in terms of the SFR averaged on 800 Myr timescales. For all of the mass bins, the SFR79 profiles appear to be overlapped together, independent of the ΔSFR9. This means that the SFR79 profile does not depend on the location of galaxies relative to the SFR9-based SFMS. This is completely consistent with the result for the integrated SFR79 that was presented above in Section 4.2.

These results in Figures 12 and 16 are consistent but quite distinct from the result in W19. W19 found that galaxies with higher SFR with respect to the SFMS show higher SFR at all galactic radii. However, there was no direct information on the temporal changes in the SFR in W19, only an inference of spatially coherent temporal variations in the SFR. The new result in the top panels of Figure 16 allows us to look directly at these temporal changes. Galaxies with a larger (temporal) change in their overall SFR, show larger (temporal) changes in their SFR at all galactic radii. However, as shown in the bottom panels of Figure 16, we do not find larger (temporal) changes in the SFR in galaxies that have different overall levels of star formation. This emphasizes that the SFR79 parameter, characterizing the change of SFR, gives a different perspective on galaxies than the SFR itself.

We now return to examine the dispersion in SFR79 from galaxy to galaxy at a given relative radius within a given stellar mass bin. This dispersion was shown in Figure 12 as the error bars on each point and was computed with a sigma-clipping algorithm by iteratively rejecting points beyond 3σ. For the two lowest stellar mass bins, it is evident that the scatter of SFR79 across the galaxy population does not change significantly with radius. However, for the more massive bins in stellar mass, the dispersion in SFR79 clearly increases toward the centers of galaxies relative to the scatter in the outer regions.

Since it is improbable that widely separated galaxies are varying coherently in their SFR, apart from the slow overall "cosmic" evolution of the SFMS discussed previously, the scatter in SFR79 across a population of galaxies (at a given epoch) is, as discussed earlier, unambiguously a measure of the variability of the SFR within individual objects. This is quite different from the scatter in Σ_SFR7, which could reflect temporal changes in individual galaxies, but also, conceivably, could reflect intrinsic (time-independent) differences from galaxy to galaxy. Whereas temporal variability could only be postulated in W19 as an explanation for the correlation of σ(ΔΣ_SFR7) with the depletion time, in this paper, we are able to measure temporal variability directly.

A clear predication of the gas-regulator model of Lilly et al. (2013) is that the variability of the SFR in response to variations in the inflow should strongly depend on the gas depletion time. Specifically, W19 drove the gas-regulator system with a sinusoidal inflow rate, and showed that the amplitude of the variation in the resulting SFR (σ_SFR) is the amplitude of the variation of the inflow (σ_Φ) multiplied by a frequency-dependent response curve. The response curve can be written as

$$\begin{eqnarray}&&f=\displaystyle \frac{1}{{\left[1+{(2\pi {\tau }_{\mathrm{eff},\mathrm{dep}}/{\tau }_{{\rm{P}}})}^{2}\right]}^{1/2}},\end{eqnarray} \tag{ 8 }$$

where τ_P is the period of the cold gas inflow, and τ_eff,dep is the "effective" gas depletion timescale, defined as the gas depletion timescale M_gas/SFR divided by the mass-loading factor of any outflow. As shown in Equation (8), for a given inflow, a shorter τ_eff,dep leads larger variations of SFR, and vice versa. Although Equation (8) was derived for an idealized sinusoidal input of inflow rate, similar shapes of the response curve are seen for a range of more complicated forms for the inflow rate (see details in W19).

Figure 17 shows the scatter of SFR79 at different galactic radii as a function of the inferred gas depletion time for different radii in the five stellar mass bins, computed as described in the previous subsection. In Figure 17, each data point represents one galactic radial bin at a given galactic stellar mass bin. According to Equation (7), the inner regions of galaxies correspond to the shorter gas depletion timescale.

Zoom In Zoom Out Reset image size

As shown, we find that there is a tight correlation between the observed dispersion of SFR79 and the inferred gas depletion timescale, with the shorter τ_dep associated with the larger range of SFR79. The scatter in SFR79 from galaxy to galaxy may be taken as a direct measure of the variability of the SFR in galaxies on timescales between 10⁷ and 10⁹ yr.

This new result is completely independent from that presented in W19 but is entirely consistent with the heuristic model presented therein. Regions with higher SFE, i.e., with shorter depletion times, do indeed appear to show a larger response in their SFR to changes in the inflow (see Figure 12 in W19). This supports the idea that the dynamical response of gas-regulator model to a time-varying inflow is the origin of the variation of SFR within and across galaxies. This further supports the idea that the narrow SFMS is indeed the result of the quasi-steady state between the inflow, outflow, and star formation (e.g., Bouché et al. 2010; Schaye et al. 2010; Lilly et al. 2013; Tacchella et al. 2016; Wang et al. 2019).

Throughout this work, we have excluded the regions that are located on the Seyfert regions on the BPT diagram in the analysis, because the contribution of Hα emission by AGN would likely lead to an over-estimation of SFR79 (see Section 4.1). However, some inner regions of galaxies may still be contaminated by a low-luminosity AGN; although, this effect must be weak for SF regions.

The connection between AGNs and instantaneous star formation has been investigated by many authors (e.g., Stanley et al. 2015, 2017; Harrison 2017; Bernhard et al. 2018; Dai et al. 2018; Scholtz et al. 2018; Schulze et al. 2019). Although both positive and negative AGN feedback have been proposed in the literature (e.g., Silk 2013; Bieri et al. 2016; Kalfountzou et al. 2017; Shin et al. 2019), convincing observational evidence for the impact of star formation by AGNs is still lacking (Stanley et al. 2015; Bernhard et al. 2018; Ramasawmy et al. 2019; Scholtz et al. 2020). Therefore, it is not clear to us what the SFR79 profile should be for an AGN-host galaxy.

In addition to AGNs, other physical processes have been proposed to play roles in changing the instantaneous SFR via different mechanisms, such as the existence of bar (e.g., Wang et al. 2012; Lin et al. 2017), and tidal/ram pressure stripping (Gunn et al. 1972; Moore et al. 1996; Abadi et al. 1999; Poggianti et al. 2017). For instance, the existence of the bar can effectively transfer the cold gas to the galactic center and lead to central enhanced star formation (Lin et al. 2017).

We emphasize that our star formation change parameter, SFR79, allows examination for the change of SFR in response to such processes. In contrast with previous studies, the SFR79 characterizes the change of star formation with respect to the level in the past, rather than to an assumed "control" population.

Many studies suggest that massive galaxies are assembled and quenched from the inside outwards (e.g., Pérez et al. 2013; Tacchella et al. 2015; Abdurro'uf 2018). Under the "inside-out" scenario, the profile of SFR79 would show a drop in the center of "quenching" or "newly quenched" galaxies (but see Lilly & Carollo 2016 for how sSFR gradients can arise without differential radial quenching). This, at first sight, seems consistent with the result in Figure 16, that galaxies with an overall lower SFR79 show a significant drop of SFR79 at the galactic center, at least for two highest stellar mass bins. However, we note that our galaxies are all SF galaxies as defined with the SFR7-based SFMS, which means that we may not expect to see the SFR79 profiles for newly quenched galaxies in our analysis. Instead, the SFR79 characterizes the movement of SF galaxies on the SFMS. Furthermore, analogously to our analysis in W19, galaxies with overall higher SFR79 are also seen to have elevated SFR79 in their centers.

Under the dynamical gas-regulator model (Lilly et al. 2013; Wang et al. 2019), the amplitude of the change of SFR is larger for regions of short gas depletion time, and the change of SFR can be both suppression and elevation. This is also in good agreement with our result. Therefore, the drop of SFR79 in galactic center for massive galaxies with overall lower SFR79 may have little to do with quenching per se, and may be purely due to the drop of overall inflow rate.

Finally, as discussed in W19, there is no need for any specific physical processes that operate in the center to quench a galaxy "inside-out," such as AGN feedback. Instead, the cut-off of the global inflow would naturally lead to an inside-out suppression of star formation, due to the short gas depletion time in galaxy centers.