© The Authors 2026

1. Introduction

The star formation main sequence (SFMS) and the mass-metallicity relation (MZR) are two of the most widely studied relations in galaxy evolution. Their very existence implies a link between the current properties of the galaxies and their cumulative history as both the stellar mass (M_★) and metallicity are the direct consequence of the star formation history (SFH).

Many authors have interpreted the low scatter among the SFMS (∼0.3 dex), observed at different redshifts, as evidence that star formation in disk galaxies is regulated by feedback processes and gas accretion (e.g., Whitaker et al. 2012; Ilbert et al. 2015; Salmon et al. 2015; Popesso et al. 2019). Alternatively, Kelson (2014) argues that the relation and its scatter can be reproduced using stochastic prescriptions for the SFH as a random walk of the star formation rate (SFR), and that it arises from the central limit theorem.

Several works have studied secondary correlations of the SFMS to assess which mechanisms preferentially regulate the growth of galaxies. González Delgado et al. (2016) and Cano-Díaz et al. (2019) propose that morphology is partly responsible for the scatter of the SFMS in both its global and resolved versions. Saintonge et al. (2016) find similar results regarding the gas fraction (which strongly correlates with morphology), but Hall et al. (2018) find no effect of morphology or gas fraction in the resolved SFMS, but rather a critical mass above and below which different mechanisms contribute to the scatter. Barrera-Ballesteros et al. (2021) and Ellison et al. (2024) favor instead gas pressure as a factor contributing to the scatter and, therefore, regulating the efficiency of the star formation. On global scales, Berti et al. (2021) tie the scatter to the large-scale environment, and Gladders et al. (2013), Rodríguez-Puebla et al. (2016), Matthee & Schaye (2019), and Blank et al. (2021) propose that most of the scatter is due to differences in the formation time of the parent halos. This means that it is not intrinsically a difference in SFR values that creates the scatter. Instead, some galaxies have lower or higher M_★ values due to having had different lifetimes to accumulate M_★ for a given halo mass (which drives gas accretion rates and therefore the SFR).

These mechanisms are expected to act on a variety of temporal scales, meaning that one of the ways to discriminate between them is to check whether the scatter results from short-term SFR variations or whether the current position of a galaxy relative to the SFMS is maintained over longer timescales in its lifetime. If galaxies are dominated by short-scale fluctuations (< 1 Gyr), the SFMS is intrinsically tighter than what is observed and the scatter originates in short stochastic deviations, such as random bursts or decreases in star formation (e.g., Peng et al. 2010; Behroozi et al. 2013; Speagle et al. 2014).

Longer-timescale fluctuations of the SFMS (several gigayears) are typically associated with larger physical scales in the Universe, either the aforementioned delays in halo formation time in galaxies or the environment via differences in the clustering of galaxies (Berti et al. 2021). In practice, they produce offsets or differences in shape between SFHs such that some galaxies spend most of their lifetime in the main sequence (MS) above the SFMS and others preferentially stay below it. Alternatively, Morselli et al. (2017) show that variations perpendicular to the SFMS correlate with the morphology of the galaxies.

Measuring the motions of galaxies on the parameter space of the SFMS at different timescales is very difficult. The variability in the SFMS has been mainly studied in cosmological simulations where the SFH of the galaxies can be readily obtained without worsening reliability at older ages, as happens with spectral energy distribution and full spectrum fitting techniques. Tacchella et al. (2020), Iyer et al. (2020), and Wang & Lilly (2020b), among others, have studied the variability of the SFHs on different timescales and have found results similar to those of Matthee & Schaye (2019), showing that the temporal variation in the SFR in galaxies is a combination of longer timescales with changes on gigayear scales and shorter fluctuations with periods of 10–100 Myr, consistent with the ages of star-forming regions and starbursts. Observationally, SFR variability, or “burstiness”, has only been studied using indirect methods, for example by comparing SFRs measured using different tracers such as Hα and UV fluxes, which are associated with different timescales and therefore the difference between their values can be considered to arise from short-term variability (e.g., Weisz et al. 2012; Guo et al. 2016; Emami et al. 2019; Caplar & Tacchella 2019; Wang & Lilly 2020a; Patel et al. 2023; Clarke et al. 2024).

In this study we approached this topic from a purely observational viewpoint by measuring the number of times each galaxy crossed the SFMS in its recent history. We used SFHs obtained via full spectral fitting and focused on the better-sampled and more reliable last gigayear to determine how many times they crossed the SFMS, how much they deviate from it after crossing them, and how much time passes before they cross it again. The crossings are associated with short time variability (10–100 Myr) in the galaxies, and we used their properties and those of the host galaxy to discern which physical processes they could be associated with. We also assessed whether these oscillations can account for the scatter in the relation or if galaxies maintain their position relative to the SFMS over longer timescales. In Sect. 2 we present the data employed and in Sect. 3 how they were processed to obtain the number of crossings. In Sect. 4 we present the results, which we discuss in Sect. 5. We adopt a standard Λ cold dark matter cosmology (H_O = 73 km s⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7) and a Salpeter initial mass function (Salpeter 1955) throughout this work.

2. Data

The Mapping Nearby Galaxies at APO (MaNGA) survey, as described by Bundy et al. (2015) and Abdurro’uf et al. (2022), comprises integral field unit spectroscopic observations conducted on a group of roughly 10 000 local galaxies selected based on their luminosity in the redshift range 0.01 < z < 0.15 (⟨z⟩∼0.03). Data were acquired using the Baryon Oscillation Spectroscopic Survey (BOSS) spectrographs (Smee et al. 2013) on the Sloan 2.5 m telescope at Apache Point Observatory (Gunn et al. 2006). In addition to the fiber bundles used for the observation of each galaxy, varying in size depending on the observed object (Drory et al. 2015), additional fibers and fiber bundles were used for the flux calibration and sky subtraction tasks (Yan et al. 2016). The data cubes, spanning a spectral range from 3600 Å to 10 300 Å with a spectral resolution of approximately R ∼ 2000 and a spatial resolution of roughly 2.5″/FWHM, were derived through the reduction and calibration of the observations using the Data Reduction Pipeline (Law et al. 2016).

We refined the full sample (DR17, 10 245) following Camps-Fariña et al. (2022), i.e., by selecting galaxies that passed the quality control described in Sect. 4.5 of Sánchez et al. (2022). Our sample selection and subsequent analysis make use of the pyPipe3D data produces, which are described in detail in Sect. 3.1. As a brief summary, first an automatic check is done to make sure that all data-products have been created and all the relevant scripts have been successfully executed, followed by comparing the values for the redshift and M_★ that we recover with those of the NASA-Sloan Atlas (NSA) catalog¹. These are expected to be stable values regardless of the method employed to measure them, so galaxies exceeding a 30% discrepancy are discarded. This is followed by a multiple-round human inspection of the central spectrum and its fit, mass assembly and chemical enrichment histories, mock photometry versus Sloan Digital Sky Survey (SDSS) images, and maps of emission lines, ages, stellar abundances and diagnostics such as the Baldwin-Phillips-Terlevich diagram (Baldwin et al. 1981). Since we measure the metallicity at the effective radius R_e (see the next section), we additionally select only galaxies below 70° in inclination and remove optically detected active galactic nuclei (types I and II; see Lacerda et al. 2022) to avoid contamination of their emission in the stellar spectra. The final sample includes 8960 galaxies spanning a range of stellar masses between log M_★/M_⊙ ∼ 8.5 − 12 as well as all morphological types.

3. Analysis

We use full spectral fitting to derive the SFH of the galaxies and trace their position in the SFMS during the last gigayear. Full spectral fitting is used to determine what fraction of light emitted today can be assigned to a set of stellar populations of given ages and metallicities. These fractions of light can be converted into the amount of M_★ that was formed at a given age and the average abundance of the interstellar medium at that time.

3.1. pyPipe3D

We derived SFHs using pyPipe3D (Sánchez 2006; Sánchez et al. 2016b,a; Lacerda et al. 2022) together with the MaStar-sLOG stellar population template library (Sánchez et al. 2022). pyPipe3D is a full spectrum fitting software capable of fitting both the absorption and emission components of galaxy spectra and analyzing them separately. The individual spaxels of the MaNGA data cubes are first spatially binned to ensure a good S/N ≳ 50. Each stellar spectrum is then analyzed using a nonparametric fitting procedure to find the linear combination of spectral models that best fits the observed spectrum. What this means is that no shape of the SFH is assumed and the pipeline only returns the light fractions for each individual population, which are used to reconstruct the SFH. Therefore, any variability in the histories can only arise due to the intrinsic physical changes in the values or be due to noise and uncertainties in the measurements or fitting.

MaStar-sLOG is a single stellar population (SSP) library constructed using spectra from the MaStar (Yan et al. 2019) stellar library and Charlot-Bruzual 2019 models² made with the GALAXEV software (Bruzual & Charlot 2011). MaStar spectra were observed using the BOSS spectrographs, which are the same as those employed in MaNGA, ensuring that the resolution and spectral coverage of the fitted and template spectra are identical. This setup is ideal for the purposes of stellar population fitting techniques and will reduce the systematic errors in the results. The library is comprised of 39 ages between 1 Myr and 13.5 Gyr using a pseudo-logarithmic spacing and 7 metallicity values (Z = 0.0001, 0.0005, 0.002, 0.008, 0.017, 0.03, and 0.04; see Sect. 4.2 in Sánchez et al. 2022).

pyPipe3D performs the fitting in several steps meant to improve the precision of the recovered results. First, the nonlinear parameters (redshift, dust attenuation and velocity dispersion of the stellar component) are determined using a reduced set of SSP and spectral range to avoid degeneracies (see Sánchez-Blázquez et al. 2011), with the emission lines masked. We employed a Cardelli et al. (1989) extinction law for the dust attenuation. Once these parameters are fixed, the resulting stellar continuum is subtracted and the emission lines fitted. The derived emission lines are then subtracted from the spectrum, leaving only the stellar component, which is then fitted with the full SSP templates. The fitting procedure includes a random exploration of the parameter space for the line-emission fitting as well as Monte Carlo iterations with perturbations for the derivation of the SSP light fractions, ensuring the results are not affected by local minima in the parameter space and that they are as robust and reliable as possible. For more details, we refer the reader to Lacerda et al. (2022) for a description of the code and Sánchez et al. (2022), where its application to MaNGA and its results are shown and discussed in detail. Our SFHs are derived from the light fractions obtained in Sánchez et al. (2022).

For the analysis presented in this work, we derived global values of the SFH for each galaxy, obtained by adding the SFHs of all the bins as well as the metallicity [Z/H] at the effective radius (R_e) to derive the MZR. The [Z/H] at the R_e was calculated by fitting the radial gradient of metallicity at each age and selecting the value at the R_e. The metallicity at the R_e has been shown to be a robust estimation of a galaxy’s global metallicity (González Delgado et al. 2014; Sánchez 2020).

Although the spatial coverage of the data is not exactly the same for all the galaxies, the MaNGA sample selection guarantees that we reach out to 1.5 R_e and 2.5 R_e for 80% of the Primary and Secondary subsamples, respectively, which does not significantly change our measurement of [Z/H] at the R_e. For the SFR and M_★, however, this will impact their total value, which can potentially affect the results in two ways, either in terms of the fluctuations of the SFH or by inducing an artificial offset to the SFMS. The fluctuations in the SFH should not be affected significantly unless the bulk of star formation occurs at the outskirts of the galaxies, which should be rare. The loss of light will decrease both the SFR and M_★, meaning that the offset to the SFMS will be canceled out to a large extent. In Appendix B of Camps-Fariña et al. (2024) it is shown that the lost flux is relatively small (0–20%) and fairly independent of M_★ and in Appendix A of this article we show that the Primary and Secondary subsamples do not show a significant offset in the SFMS and MZR.

In Fig. 1 we show the SFH for two galaxies as an example, one of them located within 1σ of the SFMS and one below this value. The galaxy currently below the SFMS used to be star-forming but eventually retired with a steep drop in SFR.

Fig. 1. Example SFHs obtained from the full spectral fitting performed on the MaNGA sample. We show one galaxy within the MS of the SFMS (MS) and one located more than 1σ below it (RET). The corresponding PLATEIFU identifiers of the data cubes are 12514-12701 for the MS galaxy and 10001-1902 for the RET galaxy.

The emission line spectra are not used to measure the SFR or [Z/H] at any point in the article and are only used for comparison purposes. We did, however, use the EW_Hα to select star-forming galaxies and define the star formation sequence in the SFMS (Cano-Díaz et al. 2016). We also used the emission lines for a consistency check on the reliability of the SFR values by comparing the current values derived from emission lines to equivalent measurements using the SFH in Fig. 2. In making this comparison it is very important to ensure that the timescales associated with each estimate of the abundance and the SFR match properly, which for the case between full spectral fitting and optical emission lines is about 30 Myr (Asari et al. 2007; González Delgado et al. 2016). Therefore, the full spectral fitting values we show in the figure are measured by selecting only the light fractions for templates with ages under 30 Myr. The emission line SFR values are calculated from the dust-corrected H_α emission using the Hα/Hβ ratio as described in Sánchez et al. (2022). The SFR values are taken as the global value for the galaxy and they follow the 1:1 relation quite well (σ ∼ 0.35 dex at log SFR EL > −2) except at very low values, which is also the case in which both the light fractions and the emission line fluxes will be very low and closer to the noise.

Fig. 2. Bottom: Comparison between the SFR measured using emission lines (x-axis) and the full spectral fitting we employed, selecting the last 30 Myr (y-axis) for all 8960 galaxies in our sample. The contours are set to encompass 95%, 65%, and 35% of the sample. Top: Distribution of data points.

3.2. Crossing the SFMS

To quantify the evolution of galaxies in the SFMS, we calculated the number of times they cross it during the last gigayear, transitioning from having an SFR higher than predicted for its mass to values below, or vice versa. The first step is to determine the SFMS as well as the MZR, which will be used to see if the crossings of the SFMS correlate with where galaxies are in the MZR. We are mostly interested in currently star-forming galaxies since we probe the recent history. We therefore employed a cut, selecting galaxies with EW_Hα > Å (Cano-Díaz et al. 2016) to define the star-forming sample, which we used to measure the SFMS and MZR. Since the SFMS evolves over time and we are measuring the position of galaxies relative to it, we recalculated it at each age for the comparison. As such, we ensured that each crossing represents a change in the position of the galaxy with respect to the average value for currently star-forming galaxies. As seen in Fig. 3, the SFMS for the most recent bin agrees well with the measurements from Speagle et al. (2014) and Popesso et al. (2023) at z ∼ 0.

Fig. 3. SFMS and MZR measured from the most recent values of the SFHs and the [Z/H] we employed (∼20 Myr) obtained from full spectral fitting. The full distribution (ALL) is shown as black points, with MS (blue) and RET (brown) shown in colored points in the left and right columns, respectively. The contours follow the star-forming subsample as defined in Sect. 3.2, which is used to measure the SFMS and MZR (red lines). In the top two panels we show the SFMS at z ∼ 0 from Speagle et al. (2014, left) and Popesso et al. (2023, right) for comparison with ours.

Due to the redshift range in MaNGA, the age values at which we measured the SFR correspond to different look-back times (LBTs) when corrected for light-travel time. Ideally, we would apply this correction so that the histories are on the same cosmic clock, but this would introduce a selection bias in the sample as the LBT bins would not be equally populated. In fact, a significant number of galaxies are located far enough that their most recent measurement is older than 1 Gyr in LBT. Because of these inconsistencies it is better to use ages so that the results are comparable between galaxies at the expense of measuring galaxies at different cosmic times for the same age bin. Due to the selection strategy of MaNGA (which observes more luminous galaxies at higher z) this will mainly correlate with M_★ in our results, meaning that massive galaxies are measured up to ≳1 Gyr earlier. We do not expect this to cause a significant difference in the results, as galaxies with similar mass will have a very similar redshift. As a result, the SFMS and MZR will be largely consistent within M_★ bins, which is good enough for us to measure the crossings for each galaxy individually.

The number of crossings is determined as the number of times the sign of the difference between a galaxy’s SFR and that of the SFMS changes within the past gigayear of each galaxy’s history, in other words, each time that a galaxy goes from being above the relation to being below or vice versa. An example of this is shown in Fig. 4. We selected the last gigayear due to the sampling of the SSP library we employed, which has 24 values within the last gigayear but only 3 between 1 and 2 Gyr. This is because the spectra of SSPs change much more slowly with time at ages greater than 1 Gyr. The first few bins have very low age intervals (∼10 − 10 Myr) on which galaxies are not expected to significantly vary their SFR, which is why we averaged the first 20 Myr of the SFH (see Table 1).

Fig. 4. Example of how crossings of the SFMS are detected. We show the position relative to the SFMS in dex of two galaxies with similar masses but very different numbers of crossings, obtained by subtracting the value of the SFMS at each age from the SFR at each age. The age bins where a crossing is detected (i.e., the previous age bin has the opposite sign for Δlog SFR) are indicated with a dot.

We selected three bins for the galaxies based on their position in the SFMS: (i) all galaxies (ALL), (ii) galaxies within 1 σ_SFMS of the measured SFMS that are considered to be MS galaxies, and (iii) a band of the same width (1 σ_SFMS) as the MS but immediately below it in the SFMS (RET). The RET bin is meant to represent galaxies located below the SFMS. We imposed the same width as MS for similarity between the bins and to avoid galaxies that have very low fractions of light in young populations, impacting the quality of the SFH and [Z/H] measured in the last 1 Gyr. There are 8960 galaxies in ALL, 4654 in MS and 2650 in RET.

In Fig. 3 we show the measured SFMS and MZR for the galaxies and how MS and RET are distributed in the relations. In the SFMS, MS galaxies are located close to the relation for star-forming galaxies and RET are located immediately below them, as would be expected since these are the definitions employed to select them. In the MZR the MS galaxies are distributed over the entire parameter space that ALL occupy, though centered around the measured relation for star-forming galaxies. RET galaxies, on the other hand, are located in a tighter sequence above the MZR, showing that quenched galaxies have higher metallicities. However, the age-metallicity degeneracy can make it hard to distinguish between whether a galaxy has old ages or high abundances in their stellar populations, in Fig. 2 the galaxies with very low SFR (which will populate mostly the RET sample) appear to be offset to higher SFR as measured from their stellar populations compared to using emission lines. It is possible then that the higher metallicity for RET galaxies comes from these galaxies being measured with younger stellar populations than they should be. Since the measurements use very young stellar populations (≲20 Myr), it is unlikely that the entire offset is due to degeneracy since the correlation between age and metallicity largely breaks down for younger ages. This is at the expense of the metallicity being harder to measure due to its associated absorption lines being less visible in spectra dominated by very young stars (e.g., Vazdekis 2008; Conroy 2013).

4. Results

Using the methodology outlined in the previous section, we could obtain the age values at which each galaxy crosses the SFMS. We measured three parameters for each galaxy, (i) the total number of times it is crossed (𝒞_SFMS), (ii) the average maximum deviation from the relation between two crossings measured in dex (ΔSFR), and (iii) the average time interval between two crossings (ΔAge). These are proxies for the frequency of the variability, the intensity of the oscillations, and the duration of the deviation or period of oscillation, respectively. We considered both increases and decreases in the SFR relative to the SFMS, taking the absolute value of ΔSFR when a galaxy goes below the relation. See Fig. 4 for an example indicating how the values are measured.

Because we measured the deviation relative to the SFMS, it is not reflective of how much a galaxy deviates from its steady state but instead how much it deviates from the average. This introduces a slight inconsistency in the case where galaxies have different long-term tracks. Consider a galaxy that is typically 0.1 dex below the SFMS and which has a burst of star formation that increases its SFR by 0.2 dex. We would detect this burst as a pair of crossings of the SFMS but measure its ΔSFR as 0.1 dex instead of 0.2 dex. Our deviations reflect how far from the average for its M_★ a galaxy goes when it oscillates, rather than the intrinsic amplitude of its oscillations. This can also affect the measurement of the intervals, if the time it takes for a galaxy to reach the relation is longer than the age sampling.

To properly interpret the results shown here, it is very important to understand how reliable the results obtained are. Given that we are measuring fluctuations in the values of the SFR and the abundance, we are uniquely affected by noise and random uncertainties that produce fluctuations in value as well. Because of this, we tested how the results would change if random noise was the origin of the scatter in SFR. The results, shown in Appendix B prove that the distribution of values and the resulting correlations with M_★ due to stochastic noise are distinct from those we measured. Our results are still very likely affected by random uncertainties but not to the point that we cannot find underlying correlations.

In Fig. 5 we show the distribution of values for the number of crossings and deviation from the SFMS for ALL galaxies, calculated as the average per galaxy and separated into three M_★ bins. We do not include the interval ΔAge due to how the uneven age sampling affects this measurement; we study this parameter in the following subsection. The distribution is shown using histograms as well as smoothed distributions calculated using kernel density estimation, which adaptively selects the smoothing parameters based on the data. We employed the gaussian_kde method from scipy with the bandwidth estimation method from Scott (2015).

Fig. 5. Distribution of the parameters for ALL galaxies for three M★ bins. Top: Average number of crossings per galaxy 𝒞SFMS. Bottom: Average deviation from the SFMS ΔSFR. Each distribution is shown as both a histogram and the kernel density estimator-smoothed curve for each M★ bin. The kernel bandwidth was estimated using the default method from SciPy.

𝒞_SFMS shows zero as the most common number of crossings and less massive galaxies tend to have a higher number of crossings compared to more massive ones. In addition to the peak at 𝒞_SFMS = 0, there is a wide distribution centered around 𝒞_SFMS ∼ 4 that constitutes most of the galaxies in the lowest M_★ bin. The distribution of 𝒞_SFMS is very similar in shape for the two lower M_★ bins in contrast to the most massive one, which shows a decline toward higher 𝒞_SFMS rather than a wide peak. There is also a preference for an even number of crossings with about twice the fraction compared to and odd number of crossings, meaning that galaxies tend to deviate occasionally from their steady state before going back to it. We can infer that there is a certain value of SFR for each galaxy that represents its steady state, with intense bursts or quenching episodes occurring occasionally, as opposed to a scenario where they are constantly oscillating around a mean value. In the latter case there would be no preference for an even number of crossings, just as them being noise-dominated would also produce a smooth distribution (see Appendix B). Regarding ΔSFR, the distribution becomes wider and shifts to higher deviations for the lower M_★ bins but the distributions overlap almost completely.

4.1. Measuring the time intervals

A third parameter to consider is the interval between crossings (ΔAge), which can be taken as an approximate proxy for the period of the variability or the duration of bursts and/or small quenching episodes. Similar to ΔSFR, this measurement is affected by the offset between the typical position of a galaxy in the SFMS, as a galaxy needs to reach the average value for its M_★ bin before its interval is measured. This offset can affect the results if the time it takes for an oscillating galaxy to reach the MS is longer than the sampling in age. Measuring the intervals using our dataset is particularly tricky due to the highly uneven sampling in age of our SFHs (pseudo-logarithmic). Due to this, very short intervals can only be measured if they occur at a very recent age, with age steps ranging from a few to ∼200 Myr. This makes it difficult to properly assess whether there are intrinsic timescales associated with the variability of galaxies, as the age at which an oscillation occurs will play a significant role in determining the value of its interval.

To circumvent this problem, we employed a statistical approach that will allow us to detect the presence of intrinsic timescales and approximate their value. Assuming these timescales exist, they would be overrepresented in the distribution of ΔAge values compared to random fluctuations. Thus, we generated a new set of SFHs such that the distance to the SFMS for each galaxy at each age bin is randomly sampled from a normal distribution whose amplitude matches the scatter. For more details see Appendix B. This new set of histories reproduces the SFMS at each age perfectly but erases any real, physical, information related to variability in the SFH. To estimate uncertainties and to improve the robustness of the calculation, we repeated the calculation for the noise-generated distribution 100 times and used the average at each ΔAge bin.

In Fig. 6 we show the distribution of ΔAge resulting from our measurements as well as the one resulting from imposing random noise as the origin of the scatter in the SFMS. The general shape of the distributions is similar, but the real results have about a third of the crossings, showing that the detected fluctuations are mostly stochastic in nature but also that random uncertainties do not dominate our results. In the bottom panel of the figure, we show the ratio between the measured and noise-generated distributions. Any intrinsic timescales that are more likely than would be expected from random fluctuations should show up here as peaks due to these ΔAge values being overrepresented.

Fig. 6. Distribution of ΔAge intervals for galaxies within the MS, sampled in ΔAge bins 10 Myr wide. Top: Our measured distribution compared to that resulting from imposing noise as the separation from the SFMS at each age value (see Appendix B). Bottom: Ratio between the measured and noise-generated intervals for each bin. The peaks correspond to statistically overrepresented intervals in the measurements. Black bars show the error estimated assuming that the number of intervals measured follows a Poisson distribution. Dotted lines show the age intervals between age bins within the last 1 Gyr, and red shaded areas indicate the bins for which no combination of consecutive age bins is possible and therefore cannot be produced.

Indeed, in the ratio plots the larger peak centered at lower ΔAge values disappears (comparatively) and a few very prominent peaks appear. Of particular relevance are the values obtained for the intrinsic timescales in the SFMS: 135, 305, 595, 720–740, and 900–940 Myr. These follow very closely a progression of multiples of ∼150 Myr: 150, 300, 450, 600, 750, and 900 Myr with 450 Myr skipped, though there is a small peak at precisely this value. The significance of the peaks is also very high, being between 10 and 120 times higher than the values expected for a noise-generated distribution (despite this distribution generating ∼3 times more crossings in general). It bears mention that the height of the peaks in the ratio between distributions is a measure of how significant the signal is but does not mean that the highest peak is the most likely interval to be observed. We show the number of crossings measured for each bin next to the peaks and it is clear that the 135 Myr peak is the most important timescale.

The appearance of a linear sequence of ΔAge_SFMS values is striking and strongly suggests the existence of an intrinsic timescale for the variability of the SFMS of around ∼150 Myr. In this scenario the peaks at higher timescales would result from the combination of two or more bursts or quenching episodes where we did not measure the valley between them. For an example of sorts, see Fig. 4 at the location where ΔSFR is indicated, this peak and the one to its left are detected as separate fluctuations, but if the data point in between them was (due to noise or the sampling) a bit higher they would instead be detected as one fluctuation of double the interval.

The peaks in the ratio of the observed data to the simulated data in the left hand panels of Fig. 6 appear unnaturally sharp; however, this is largely due to the use of evenly sampled histograms even though not all these bins can be populated equally or at all. The allowed values of ΔAge are either the intervals of the sampling (indicated in the figure) or the sum of several consecutive intervals. This is the reason that, despite averaging 100 instances, the noise-generated distribution exhibits many peaks and null values (indicated as red zones in the figure) instead of a smooth distribution. Some of the bins have more opportunities to be populated and others are impossible to fill as there is no combination of intervals that lies in them. The expectation of a wider distribution around each peak is also not necessarily possible for all peaks. For each of the age bins within a ΔAge interval, there is a chance to artificially cross the SFMS due to noise, which breaks up the interval into pieces of arbitrary lengths. For a regular sampling of the ages the adjacent ones would be expected to be populated more, but this is not our case. Populating the adjacent bins requires the lower age bins to be involved since their separation in age is smaller, and for the larger intervals, which span a large portion of the last 1 Gyr, the involvement of younger, smaller age bins is much more likely and that is why these do show wider distributions. That being said, it is interesting to note that none of the prominent ΔAge values coincide with the sampling, meaning that they must arise from a combination of age bins. Even the 135 Myr peak lies between sampling intervals, meaning that it is populated by oscillations that have several SFR measurements between the two crossings.

4.2. Crossings versus M_★

In Fig. 7 we show the correlation between M_★ and the parameters of the crossings for the three selected groups described in Sect. 3.2: ALL, MS, and RET. For 𝒞_SFMS and ΔSFR we averaged the values for each galaxy within 0.25 dex wide M_★ bins and the errors correspond to the error of the mean (scatter divided by $\sqrt{N}$) computed within each bin. For the ΔAge values we followed an approach similar to what was laid out in the previous section. In each bin, we calculated the distribution of ΔAge values for the galaxies included in it as well as the distribution if those same galaxies had noise-like offsets from the SFMS over time (see Appendix B). We then determined the ratio between the distributions and selected the peaks where the S/N reaches above 10. Each selected ΔAge value is included in a weighed sum whose weights are how many times each value is measured. In this manner we removed those ΔAge values that cannot statistically be distinguished from random fluctuations and then obtained a representative average μ such that the more a ΔAge value appears the more weight it has on the average. For the error we estimated the variance using the second moment with the same weights:

$$\begin{array}{cc}\hfill \mu & =\frac{{\sum }_i{N}_i·\mathrm{\Delta }{\mathrm{Age}}_i}{{\sum }_i{N}_i}\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\sigma }^2& =\frac{{\sum }_i{N}_i·{(\mathrm{\Delta }{\mathrm{Age}}_i-\mu )}^2}{{\sum }_i{N}_i}·\hfill \end{array}$$(2)

Fig. 7. Properties related to the crossings of the SFMS compared to M★, for the three groups of galaxies considered: ALL (black), MS (blue), and RET (red). Top: Distribution of the galaxies in M★. Second pane: Average number of crossings (𝒞SFMS) for each M★. Third panel: ΔSFR in dex after a crossing has happened. Fourth pane: Average interval (ΔAge) between crossings for the galaxies, in Myr. For the bottom three panels, the shaded areas correspond to the error of the mean.

The value of σ is then used to calculate the error of the mean dividing it by the square root of the total number of ΔAge values for the selected peaks with S/N above 10 ($\sqrt{\sum _i{N}_i}$).

In Fig. 7 the number of crossings of the SFMS strongly varies between MS and RET galaxies, the former of which show a fairly flat ∼4.5 average crossings below log M_★/M_⊙ ∼ 10 at which point it starts to decline. RET galaxies, on the other hand, have significantly fewer crossings on average for all masses (∼1 − 2) and the number of crossings declines steadily with mass from 2 to almost zero for the most massive galaxies. The average ΔSFR is fairly constant for MS but varies irregularly for RET, likely due to the lower number of both galaxies and the average number of crossings for this bin. Still, there is a clear offset compared to MS such that the deviations are greater.

Since RET galaxies have on average 𝒞_SFMS ≲ 2 and are selected to currently be 1σ below the SFMS, a significant portion of this sample are galaxies that have crossed the SFMS one last time within the last 1 Gyr and are currently on their way to becoming retired. In Fig. 4 the RET galaxy almost reaches the SFMS one last time before its SFR drops below –1.25 dex at the lowest point, a much larger offset than the typical deviations shown by the MS galaxy.

The averaged values for ΔAge_SFMS do not appear to change significantly with M_★, other than perhaps having slightly higher values at 10⁹ M_⊙ and below and shorter timescales above 10^10.75 M_⊙. An offset compared to RET can also be seen, but these features are not very clear compared to the uncertainties so in general we do not detect strong trends in this parameter with M_★ or between MS and RET.

In Fig. 7 the average values for Δ_SFR show an offset between RET and MS galaxies, but we do not observe two peaks in the global distribution in Fig. 5 despite the difference in value (0.4–0.6) being significant enough that it would be detectable. The contribution of RET galaxies to the distribution in Fig. 5 is therefore the extended wing toward higher ΔSFR that makes the distribution asymmetric. For MS galaxies their average values coincide with the peak in the distribution at ∼0.4 dex, suggesting that the intensity of a burst of star formation or a quenching episode for a galaxy on the MS is, on average, a fundamental value ≳0.4 dex. This means that, regardless of their mass, when galaxies have a burst of star formation or a quenching episode they typically vary their SFR by ∼2.5 times. This value is lower than the typically reported offsets of the starburst population of 0.6–1 dex (Rodighiero et al. 2011; Sargent et al. 2012; Rodighiero et al. 2014; Schreiber et al. 2015), but we need to take into consideration the fact that we required a crossing of the SFMS to occur before we measured the deviation. As such, if we use Δ_SFR as a proxy for the position of the starbursts, we are selecting galaxies that typically lie below the SFMS and ignoring those that typically lie above it (as the latter will not cross despite the increase). This biases Δ_SFR to lower values by a factor that should be of order ∼0.5 ⋅ σ_SFMS, which, given that σ_SFMS = 0.4, makes a total ΔSFR of 0.6; this matches the lower end of the expected values. Also, Δ_SFR includes the deviation from galaxies temporarily dipping below the SFMS, which might have a lower amplitude for the MS subsample. Finally, the SSP sampling can effectively “smooth” the peaks in the SFH compared to emission line measurements. Because of this, we consider our results to be roughly compatible with the reported offset for starbursts accounting for a systemic offset of 0.2 dex due to how it is measured. Therefore, the narrow distribution of Δ_SFR that does not depend on M_★ indicates that individual bursts of star formation occur in a similar, self-regulated manner, in all galaxies with a similar relative increase in SFR whose value is not strongly affected by halo-level parameters.

In general, we can characterize the main results for the crossings as such: (i) more massive galaxies have fewer crossings in general no matter what, (ii) MS galaxies have more crossings in general, (iii) ΔSFR is fairly constant with M_★ but offset to larger values for RET galaxies, (iv) we find no clear correlation between ΔAge and M_★. These results support star formation bursts or quenching episodes as self-regulated incidents where the relative intensity and the duration do not depend on the global parameters of the host galaxy, whose role is to instead regulate whether and how often the bursts occur.

4.3. Crossings versus current position in the SFMS and MZR

Another way to probe how the number of crossings can be tied to the evolution of the galaxies is to check how the averaged parameters vary depending on where the galaxies are currently located in the MZR and SFMS. In other words, we measured the average parameters for galaxies that are currently above or below either the MZR or SFMS. For this analysis we employed only the MS group of galaxies, as the RET galaxies cannot be currently located above the SFMS by definition. As can be seen in the figures, the number of galaxies above and below the relations is not exactly the same even selecting only MS galaxies. This is because the SFMS and MZR are calculated using galaxies selected with EW_Hα > 6 Å, while the MS sample requires galaxies to be within ±1σ of the SFMS. In theory, this should be roughly equivalent to selecting MS after determining the SFMS but, as with all empirical relations, there is scatter such that some galaxies are currently measured as within the MS using full spectral fitting but have lower EW_Hα so they are not selected for measuring the SFMS. The result is an excess of galaxies below the MS and above the MZR; the latter is due to how galaxies with lower SFR are located above the MZR (see Fig. 3). The difference in galaxies for the SFMS is about 12% of the MS sample.

In Figs. 8, 9, and 10 we show the results for the three parameters, 𝒞_SFMS, ΔSFR, and ΔAge, respectively. The position of the numbers in the diagrams indicates which bin of galaxies is considered. In general, the differences between bins are smaller than those observed for M_★ and some parameters appear largely unaffected by the current position of the galaxies. One thing to consider when interpreting these results is that, ideally, we would select galaxies based on their “preferred” location within the relations rather than their current one which could be due to a temporary oscillation.

Fig. 8. Average number of crossings (𝒞SFMS) for several bins of galaxies selected based on the galaxies’ currently observed position in the MZR and SFMS. In the upper panel, the diagonal lines represent the relations; the position of the numbers relative to them informs which positional bin (above or below the relation) is being considered. In the bottom panel the bins considering combinations of positions in both relations are shown, with the shaded areas indicating the position in each relation. The error is computed as the error of the mean and is indicated in parentheses: the last digit shown in the error corresponds to the last digit shown in the value. In the corners in red we indicate the number of galaxies included in each bin.

Fig. 9. Same as Fig. 8 but for the deviation from the SFMS (ΔSFR). For galaxies to be included in a particular bin, they need to have 𝒞SFMS > 0, and therefore the number of galaxies per bin is different from that in Fig. 8.

Fig. 10. Same as Fig. 8 but for the interval (in Myr) between crossings from the SFMS. For galaxies to be included in a particular bin, they need to have 𝒞SFMS > 0, and therefore the number of galaxies per bin is different from that in Fig. 8.

One of the most relevant correlations appears to be in 𝒞_SFMS, which is significantly higher for galaxies that are either currently above the SFMS (> SFMS) or below the MZR (< MZR) when considering only one relation. Considering the position in both relations at once, this difference appears to be driven by the bin located at < SFMS and > MZR, which has the lowest 𝒞_SFMS. This result fits well with pristine gas accretion sustaining star formation in galaxies, such that those that are deficient in gas accretion will not be able to fuel as many bursts of star formation. Since the gas has very low abundances, galaxies that accrete less of it will also tend to be located at higher average abundances for their M_★. We do not observe the opposite where the > SFMS and < MZR bin has a higher average 𝒞_SFMS, however, which might indicate that changes in SFR and [Z/H] are not simultaneous in general, allowing us to find galaxies with significant gas accretion in any combination of positions of the SFMS and MZR. The accretion deficient galaxies, on the other hand, do not shift around the SFMS as much and therefore correlate with their current position in both relations much more.

For ΔSFR we find higher deviations for the bins below the SFMS and MZR, but the differences are much smaller relative to the values of ΔSFR and closer to compatibility with the error. ΔAge follows similar correlations to 𝒞_SFMS, with galaxies < SFMS and > MZR having longer intervals.

Another interesting parameter to compare is the number of galaxies per bin. The distribution of galaxies is uneven: most galaxies lie on opposite sides of the two relations, that is, either the > SFMS and < MZR bin or vice versa. This suggests a coevolution of both relations, as proposed by the fundamental metallicity relation (Mannucci et al. 2010). The number of galaxies per bin changes for Figs. 9 and 10, as galaxies with no crossings do not contribute to measuring the deviation and interval. The greatest difference in the number of galaxies appears for the aforementioned bin that is expected to be deficient in gas accretion (< SFMS and > MZR) with 42% of the galaxies showing no crossings of the SFMS. Strikingly, in the bin opposite to it (> SFMS and < MZR) only one galaxy out of 1247 shows no crossings, and for the other two bins the fraction is 7% and 13%. This result agrees very well with the above hypothesis that the < SFMS and > MZR bin is the location for galaxies deficient in gas accretion. Since these galaxies still belong to the MS sample it is likely that they are just now starting to retire and will soon cross into the Green Valley. The fact that the position in the MZR is the determining factor to detect them shows that, as proposed in Camps-Fariña et al. (2023), galaxy quenching generally occurs after the interstellar medium stops being significantly diluted, likely due to gas accretion being shut off.

4.4. Fraction of time above the SFMS

So far, we have characterized the motions of the galaxies around the SFMS, showing that at least some of its scatter originates in short-term variability in the SFR (several to hundreds of megayears). The question remains whether these fluctuations can account for most if not all the scatter. If they do not, it means that galaxies of the same mass have different average values of the SFR over time. It would result in some galaxies being intrinsically more star-forming than others. These different evolutionary tracks are expected to come from differences in halo-level properties, such as the halo formation time (Matthee & Schaye 2019; Blank et al. 2021). To assess this, we measured the percentage of time that star-forming galaxies currently above and below the SFMS have spent above it.

In Fig. 11 we show the fraction of time above the SFMS versus M_★ for star-forming galaxies selected using the EW_Hα > 6 Å criterion (see Sect. 3.2). As explained in Sect. 4.3, the MS sample includes 12% more galaxies that do not fulfill the EW_Hα criterion used to determine the SFMS. These are preferentially located below the SFMS and are likely galaxies whose SFR is beginning to decline. This makes them less suitable for the purpose of this section, which is checking whether actively star-forming galaxies tend to stay on one side of the relation. As such, instead of the MS sample, we applied the EW_Hα criterion, which is independent of the full spectral fitting and therefore makes for more robust results. This sample contains 4253 galaxies.

Fig. 11. Fraction of time spent above the SFMS in the last 1 Gyr for currently star-forming galaxies. The shaded areas correspond to the error of the mean.

We show only the fraction of time above the SFMS and not that of time below it because these are interchangeable, and we have ensured that for all galaxies the two fractions add up to 1. There is a clear offset in the fraction such that, while the sample as a whole averages to galaxies spending half their time above the SFMS as would be expected, galaxies currently located above or below the relation tend to have spent, on average, most of their time on the same side (∼60%).

This result shows that a significant part of the scatter of the SFMS originates in galaxies having intrinsically different values of SFR at the same M_★. This could be due to a difference in the efficiency at which molecular gas is converted into stars (SFR/M_H2), the same but for HI to H₂ conversion, a difference in the gas fraction of the galaxies, or a difference in the efficiency at which the gas is condensed into molecular gas. As mentioned above, differences in halo-formation time are a likely explanation for this though in that case a better interpretation of the scatter of the SFMS is differences in M_★ originating in some halos having had less time than others to form stars (Matthee & Schaye 2019; Blank et al. 2021). In other words, two galaxies with the same SFH would always have similar SFR, but if we offset them in time, due to one dark matter halo forming earlier, then for a typical t-tau shape that declines after cosmic noon there will be an offset between them such that the one that formed later will have higher SFR. Environment is also a likely suspect as it can affect most of the mechanisms listed above and acts on long timescales (e.g., McGee et al. 2009; Haines et al. 2015; Paccagnella et al. 2016; Maier et al. 2019; Kim et al. 2023). The first explanation could be checked by using the halo-mass instead of M_★ for Fig. 11, which should bring the fractions closer to 0.5.

5. Discussion

Full spectral fitting with stellar populations is rife with degeneracies and systematic uncertainties (e.g., Walcher et al. 2011; Cid Fernandes et al. 2014). This makes quantitative determinations of individual galactic parameters less reliable, but they work well for statistical comparisons between galaxies such as the results we present here. The actual number of crossings and intervals shown in each figure are heavily model-dependent, primarily on the number of age values present in the SSP library employed. These represent the number of times a galaxy can be detected to have crossed the SFMS and therefore we are largely blind to crossings that occur with higher frequency. Individual histories can be rather noisy too and these surely contribute to the total number of crossings detected on each galaxy. The reliable results are therefore the statistically significant differences between bins of galaxies that we report.

Perhaps one of the most surprising results is the detection in Sect. 4.1 of what appears to be a statistically preferred timescale of ∼135–150 for SFR fluctuations such that this interval and its multiples appear much more frequently than would be expected for random fluctuations. The fraction of crossings associated with the ∼135 − 150 Myr timescale and its multiples is quite small (∼400 out of ∼22 000), but it has very high statistical significance. This possibly indicates either the existence of a regulatory mechanism that acts on this timescale or that this timescale is intrinsic to the triggering and propagation of star formation over a galaxy, such as a crossing time for self-propagated star formation. This feature must arise from either real physics or from systematic method-related uncertainties that correlate with this particular timescale of 135–150 Myr. It is never wise to assert that systematic effects are not present in any measurement obtained using full spectral fitting techniques, as they are numerous and difficult to fully disentangle from the results. That being said, our age sampling is pseudo-logarithmic in nature, and changes in the shape of the SSP templates do not generally follow linear progressions such as the one seen here. We are not aware of any systematic effect capable of artificially producing this feature, especially since the values are not perfect multiples of each other, which would be expected of a systematic, method-induced, intrinsic timescale. We consider this a promising result, but proper verification is required from other authors using different SSP templates and fitting algorithms for its detection to be certain.

From Fig. 6 it might be tempting to dismiss most of the detected crossings as noise, since the general shape of the distribution vanishes when divided by a noise-like distribution. However, this is not actually an inference that can be made from this figure. What we are measuring is how distinct from random fluctuations the distribution is, which allows us to detect the presence of the overrepresented ΔAge values we report. However, just because they follow a random-like distribution does not mean they originate from uncertainties, as they can be the result of physical stochastic processes that would produce the same distribution. In other words, the similarity between the two distributions can also be interpreted as the triggering of star formation being mostly a stochastic process. In Sect. 4.4 we show clearly that galaxies have a preferred side of the SFMS in which to reside, which we would not be able to measure if we were dominated by random uncertainties. The much lower number of crossings for the MS compared to the noise-like distribution in Fig. 6 (∼1/3) also reflects this. The best evidence that the smaller ΔAge values are not dominated by noise is in Fig. 5, which shows that the distribution of crossings is strongly biased to even numbers. In Fig. B.1 no such feature can be seen in the noise-like distribution of 𝒞_SFMS.

In conclusion, our study of the distribution of ΔAge_SFMS values points to a scenario where star formation is triggered stochastically in sub-100 Myr timescales, but these are modulated by longer bursts of star formation regulated to last about 135–150 Myr. This scenario of larger bursts of star formation on ∼100 Myr scales within which there are several shorter, ∼10 Myr individual bursts fits well with previous studies (Förster Schreiber et al. 2003; Di Matteo et al. 2008; McQuinn et al. 2009; French et al. 2018), lending credibility to our results.

One of the main questions we wanted to answer is whether short-term oscillations can explain the scatter of the SFMS. We find that the average deviation for MS galaxies is of a similar magnitude to the scatter of the relation (0.3–0.4 dex), but we also show that galaxies currently above and below the SFMS tend to spend most of their time on the same side for the past gigayear. In our results the excess is moderate, with a 60–40% split rather than the 50–50% that would arise from purely stochastic fluctuations, but as shown in Fig. B.3 uncertainties strongly affect this value bringing it closer to 50–50. The 60–40% split is reproduced when noise equal to σ_SFMS is introduced to SFHs that have no intrinsic fluctuations and their offset reproduces the scatter, so the offset between SFHs is likely a significant contributor to σ_SFMS given that we are strongly affected by noise. Consequently, our results favor the mixed scenario predicted by Matthee & Schaye (2019) and Tacchella et al. (2020) where galaxies are mostly spread around the SFMS because of long-term variation star formation and short-term variability contributes to increasing the scatter. In our results, it interesting that the number of crossings declines for high mass galaxies but the scatter of the SFMS increases, going from 0.34 at log M_★/M_⊙ ∼ 9 to 0.56 at log M_★/M_⊙ ∼ 11. Since the scatter increases even though the oscillations become less frequent it should mean that the long timescale differentiations in the SFH account for a larger portion of the scatter of the SFMS in more massive galaxies. Oddly, such a M_★ dependence is not observed in Fig. 11 for the SFMS. The dependence of the scatter of the SFMS with M_★ is still debated, several observational studies (Brinchmann et al. 2004; Guo et al. 2013; Popesso et al. 2019) find an increase with mass, especially a jump in scatter at masses approaching log M_★/M_⊙ ∼ 11 just like us, and Sparre et al. (2015) finds the same in the Illustris simulations. However, Matthee & Schaye (2019) find the opposite in EAGLE simulations, in agreement with the observations of Chang et al. (2015).

In Fig. 3 we show where MS and RET galaxies are located in both relation highlighting that, while MS galaxies occupy most of the parameter space that the full sample covers, RET galaxies are strongly shifted to higher metallicity values. This fits very well with the expectation that pristine gas accretion is required for galaxies to stay in the MS (Camps-Fariña et al. 2023). Vaughan et al. (2022), however, find that this offset in the MZR between star-forming and passive galaxies correlates with the size of the galaxies in observational data, such that due to evolutionary tracks smaller galaxies tend to be the most metallic for their M_★.

6. Summary and conclusions

In this work we used full spectral fitting to obtain the SFH and [Z/H] for 8960 MaNGA galaxies in order to study how the SFR changes over time and check whether the scatter of the SFMS can be attributed to short-term fluctuations or if instead galaxies have offsets in these parameters that are sustained over time. Because of how close to the uncertainties our measurements lie, we needed robust statistical measurements to ensure that our results are significant, which we achieved mostly by ensuring that our results cannot have arisen from random uncertainties and that they do not follow expected trends arising from fitting degeneracies.

In general, we find that this type of study is feasible with large datasets and find results that fit well with what is known in the field and predictions from simulations, making this a promising methodology especially for application to future spectroscopic surveys of galaxies with improved spectral resolution and S/N.

Some of the most important results are:

1. SFHs obtained from full spectral fitting are capable of yielding useful results on short timescales. We used this to measure the variability in the SFR by measuring how galaxies crossed the SFMS over the last 1 Gyr.

2. Variability in the SFH appears to be a combination of purely stochastic sub-100 Myr fluctuations and longer bursts or quenching episodes, with a possible slight excess on a 135–150 Myr timescale and its multiples.

3. The number of crossings of the SFMS (∼4.5) does not appear to change much with M_★ in the 10^{8.5 − 10.5} M_⊙ range. Above this range, galaxies tend to fluctuate less.

4. The degree to which galaxies deviate from the SFMS does not appear to correlate strongly with M_★ and is similar to the scatter of the relation at 0.3–0.4 dex.

5. The current locations in the SFMS and the MZR have a significant effect on how galaxies cross the former, pointing to a tie between the enrichment history and the SFH. Galaxies located below the SFMS and above the MZR show the lowest number of crossings, likely because 42% of such galaxies have not had a burst of star formation in the last 1 Gyr. These galaxies will likely soon fall into the Green Valley.

6. Galaxies appear to have preferred locations in the SFMS and have spent 60% of the last gigayear there. This shows that the scatter of the SFMS is a combination of long- (> 1 Gyr) and short-term (≲100 Myr) fluctuations, in agreement with models.

Acknowledgments

We thank the referee for their useful comments and S. Pascual for discussion on the analysis. ACF acknowledges financial support by the Spanish Ministry of Science and Innovation through the research grants PID2019-107427-GB-31 and PID2022-138855NB-C31, funded by MCIN/AEI/10.13039/501100011033/FEDER, EU. SFS thanks the PAPIIT-DGAPA AG100622 project and CONACYT grant CF19-39578. This work was supported by UNAM PASPA – DGAPA. MCC acknowledges the support of AC3, a project funded by the European Union’s Horizon Europe Research and Innovation programme under grant agreement No 101093129, PID2022-138621NB-I00, funded by MCIN/AEI/10.13039/501100011033/FEDER, EU. This work has been partially funded by “Tecnologías Avanzadas para la exploración del universo y sus componentes” (PR47/21 TAU-CM) project funded by Comunidad de Madrid, by the Recovery, Transformation and Resilience Plan from the Spanish State, and by NextGenerationEU from the European Union through the Recovery and Resilience Facility. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian (CfA), the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

References

Appendix A: Primary and secondary samples in MaNGA

The MaNGA sample is composed of three subsamples: Primary, Secondary, and Color-enhanced. The first two account for most of the galaxies (84%) and the Color-enhanced one was added to increase statistical significance in areas of the NUV-i versus M_i parameter space that are more sparsely populated. The Primary and Secondary samples are defined such that at least 80% of the galaxies contained in them are covered by the integral field unit out to 1.5 and 2.5 times the effective radius, respectively. Consequently, our SFR and M_★ determinations for the Primary sample are measured using a lower fraction of light compared to the Secondary sample. The metallicity values are less affected by this due to being measured at the effective radius.

The variability of the SFR should not be strongly affected by this except in galaxies where the bulk of star formation occurs at the outskirts, which should be relatively unusual. That being said, if there is a significant statistical difference in the total SFR measured this could manifest as an offset between the samples in the SFMS such that galaxies in the Secondary sample lie preferentially above or below it. In Fig. A.1 we show the SFMS and MZR as in Fig. 3 but with contours for the Primary and Secondary samples, with no significant offset between the samples observed.

Fig. A.1. Comparison between the Primary and Secondary samples in MaNGA regarding their distribution in the SFMS (top) and MZR (bottom) relations. In both panels the black points correspond to all galaxies in the sample, while the red lines (averaged SFMS and MZR) and contours (35%, 65%, and 95%) represent the EWHα-selected star-forming galaxies used to measure the relations. Orange contours indicate no further selection based on the MaNGA target selection, while the blue and green contours represent, respectively, the Primary and Secondary samples.

Appendix B: How distinct are our results from pure noise?

There is a variety of ways in which our results are affected by the method, notably the sampling of the stellar population ages and the intrinsic uncertainties in the determination of the parameters. Additionally, the uncertainties are not well constrained. As such, it is of interest to check how our results change if we impose random noise as a dominant factor, to see if we find significantly distinct values.

B.1. Change in our parameters if crossings are dominated by noise

If we substitute the residuals of the SFMS we used to detect crossings (see Sect. 3.2) with random values sampled from a Gaussian distribution of width equal to the measured width of the SFMS at the corresponding M_★ bin, we reproduce the overall distribution of the relation but erase any physical information in the galaxy histories. The specific age sampling will produce characteristic average values for each of the parameters when processing noise. For example, we would expect the average number of crossings to be close to half the number of age values since there is a 50% chance for the galaxy to switch sides on the SFMS if its position is determined randomly.

In Fig. B.1 we show the distribution for two of the parameters of the crossings for MS galaxies as well as what results from imposing noise as the origin of the scatter. The methodology to measure representative values for ΔAge shown in Sect. 4.1 already includes a comparison to the noise-like distribution and therefore cannot be included here. For 𝒞_SFMS and the ΔSFR the distributions for MS and noise are different, though the position of the peaks are somewhat similar in ΔSFR. The average number of crossings, which is the most important parameter since the other ones follow from the detection of the crossings, has very different distributions. Noise produces a Gaussian-like distribution centered at 𝒞_SFMS ∼ 8 while the actual data have a bi-modal distribution with a larger peak at 0 crossings and a wider component centered around 𝒞_SFMS ∼ 4 instead. The reported preference for even values of crossings is not observed for noise-like fluctuations.

Fig. B.1. Top: Distribution of values for 𝒞SFMS. Bottom: ΔSFR distributions. The actual measured values are shown in black, and the predicted values when assuming purely random fluctuations as the origin of the scatter are shown in green.

In Fig. B.2 we reproduce Fig. 7 from the main text but instead of the ALL and RET galaxy groups we only show the results for MS and the same galaxies but with noise as the scatter of the SFMS. The results are similar to those of Fig. B.1, with the differences between the actual results and noise are very apparent as the dependence on M_★ either has large offsets or differences in the shape. The number of crossings remains virtually constant except for edge cases, while the deviation is mainly modulated by the scatter at that particular M_★.

Fig. B.2. Properties related to the crossings of the SFMS compared to M★. The actual measured values for MS galaxies are shown in black, and the predicted values when assuming purely random fluctuations as the origin of the scatter are shown in green.

B.2. How noise affects the fraction of time above the SFMS

In Sect. 4.4 we show that at least part of the scatter for the SFMS is due to the galaxies preferentially spending their time on one side of the relation. Galaxies appear to spend ≳60% of the time on the same side of the SFMS but this value could be significantly affected by noise, lowering it. As should be expected, using random fluctuations the percentage of time spent on each side is almost perfectly 50% for each mass bin, with the median value of the fraction for the M_★ bins of 0.488 and a standard deviation of 0.05, which drops to 0.007 if we consider only the 10^9 − 11 M_⊙ range, removing the least populated bins. This reinforces the validity of the offset seen in Fig. 11 as a physical result, showing that it cannot originate in noise or random uncertainties from the method.

Given that noise tends to produce even amounts of time above and below the SFMS, we can also assess how much lower the fraction of time spent above it we measure will be depending on the noise. To this end, we performed a different inference regarding the residuals, assuming that the scatter of the SFMS is 100% due to an offset in SFR values over time, with no crossings. As such, we took the separation from the SFMS that the galaxies currently have and imposed the same separation for all age values in the past gigayear, thus ensuring fractions of time above the relation of 1 and 0 for galaxies currently above and below, respectively. Then we added noise to see how much this value is reduced depending on the noise level.

In Fig. B.3 we show the result of adding Gaussian noise of different amplitudes to flat residuals where the current separation from the SFMS is preserved over time. We used five values of the amplitude scaled to the scatter of the relation at each M_★, finding that higher noise amplitude makes the fraction of time spent above the relation closer to 0.5. A noise amplitude equal to the scatter reduces the fraction for galaxies currently above from 1 to ≳0.6, similar to the values we observe. Of course, it is impossible for the separation and the noise to simultaneously have the same measured scatter, as that would widen the scatter to $\sqrt{2}{\sigma }_{\mathit{SFMS}}$. Still, these results show that our measured ≳60% value for the fraction of time spent on one side of the SFMS is a lower limit. Our scatter measurements of 0.3-0.5 are about 0.1 dex wider than the range reported by the studies with the lowest measured scatters in the Local Universe (e.g., Guo et al. 2013; Popesso et al. 2019), though they are compatible with other determinations (e.g., Brinchmann et al. 2004; Ilbert et al. 2015; Chang et al. 2015). As such, we expect a possible underestimation of the fraction of time galaxies spend on one side but the offset should not be very large.

Fig. B.3. Fraction of time spent above the SFMS for currently star-forming galaxies when the current separation between each galaxy and the relation is imposed for all age values before adding noise. The colors correspond to the width of the random fluctuations that are added to this separation as a factor of the scatter of the SFMS.

All Tables

All Figures

	Fig. 1. Example SFHs obtained from the full spectral fitting performed on the MaNGA sample. We show one galaxy within the MS of the SFMS (MS) and one located more than 1σ below it (RET). The corresponding PLATEIFU identifiers of the data cubes are 12514-12701 for the MS galaxy and 10001-1902 for the RET galaxy.
In the text

	Fig. 2. Bottom: Comparison between the SFR measured using emission lines (x-axis) and the full spectral fitting we employed, selecting the last 30 Myr (y-axis) for all 8960 galaxies in our sample. The contours are set to encompass 95%, 65%, and 35% of the sample. Top: Distribution of data points.
In the text

Fig. 3. SFMS and MZR measured from the most recent values of the SFHs and the [Z/H] we employed (∼20 Myr) obtained from full spectral fitting. The full distribution (ALL) is shown as black points, with MS (blue) and RET (brown) shown in colored points in the left and right columns, respectively. The contours follow the star-forming subsample as defined in Sect. 3.2, which is used to measure the SFMS and MZR (red lines). In the top two panels we show the SFMS at z ∼ 0 from Speagle et al. (2014, left) and Popesso et al. (2023, right) for comparison with ours.

In the text

Fig. 4. Example of how crossings of the SFMS are detected. We show the position relative to the SFMS in dex of two galaxies with similar masses but very different numbers of crossings, obtained by subtracting the value of the SFMS at each age from the SFR at each age. The age bins where a crossing is detected (i.e., the previous age bin has the opposite sign for Δlog SFR) are indicated with a dot.

In the text

	Fig. 5. Distribution of the parameters for ALL galaxies for three M★ bins. Top: Average number of crossings per galaxy 𝒞SFMS. Bottom: Average deviation from the SFMS ΔSFR. Each distribution is shown as both a histogram and the kernel density estimator-smoothed curve for each M★ bin. The kernel bandwidth was estimated using the default method from SciPy.
In the text

Fig. 6. Distribution of ΔAge intervals for galaxies within the MS, sampled in ΔAge bins 10 Myr wide. Top: Our measured distribution compared to that resulting from imposing noise as the separation from the SFMS at each age value (see Appendix B). Bottom: Ratio between the measured and noise-generated intervals for each bin. The peaks correspond to statistically overrepresented intervals in the measurements. Black bars show the error estimated assuming that the number of intervals measured follows a Poisson distribution. Dotted lines show the age intervals between age bins within the last 1 Gyr, and red shaded areas indicate the bins for which no combination of consecutive age bins is possible and therefore cannot be produced.

In the text

Fig. 7. Properties related to the crossings of the SFMS compared to M★, for the three groups of galaxies considered: ALL (black), MS (blue), and RET (red). Top: Distribution of the galaxies in M★. Second pane: Average number of crossings (𝒞SFMS) for each M★. Third panel: ΔSFR in dex after a crossing has happened. Fourth pane: Average interval (ΔAge) between crossings for the galaxies, in Myr. For the bottom three panels, the shaded areas correspond to the error of the mean.

In the text

Fig. 8. Average number of crossings (𝒞SFMS) for several bins of galaxies selected based on the galaxies’ currently observed position in the MZR and SFMS. In the upper panel, the diagonal lines represent the relations; the position of the numbers relative to them informs which positional bin (above or below the relation) is being considered. In the bottom panel the bins considering combinations of positions in both relations are shown, with the shaded areas indicating the position in each relation. The error is computed as the error of the mean and is indicated in parentheses: the last digit shown in the error corresponds to the last digit shown in the value. In the corners in red we indicate the number of galaxies included in each bin.

In the text

	Fig. 9. Same as Fig. 8 but for the deviation from the SFMS (ΔSFR). For galaxies to be included in a particular bin, they need to have 𝒞SFMS > 0, and therefore the number of galaxies per bin is different from that in Fig. 8.
In the text

	Fig. 10. Same as Fig. 8 but for the interval (in Myr) between crossings from the SFMS. For galaxies to be included in a particular bin, they need to have 𝒞SFMS > 0, and therefore the number of galaxies per bin is different from that in Fig. 8.
In the text

	Fig. 11. Fraction of time spent above the SFMS in the last 1 Gyr for currently star-forming galaxies. The shaded areas correspond to the error of the mean.
In the text

Fig. A.1. Comparison between the Primary and Secondary samples in MaNGA regarding their distribution in the SFMS (top) and MZR (bottom) relations. In both panels the black points correspond to all galaxies in the sample, while the red lines (averaged SFMS and MZR) and contours (35%, 65%, and 95%) represent the EWHα-selected star-forming galaxies used to measure the relations. Orange contours indicate no further selection based on the MaNGA target selection, while the blue and green contours represent, respectively, the Primary and Secondary samples.

In the text

	Fig. B.1. Top: Distribution of values for 𝒞SFMS. Bottom: ΔSFR distributions. The actual measured values are shown in black, and the predicted values when assuming purely random fluctuations as the origin of the scatter are shown in green.
In the text

	Fig. B.2. Properties related to the crossings of the SFMS compared to M★. The actual measured values for MS galaxies are shown in black, and the predicted values when assuming purely random fluctuations as the origin of the scatter are shown in green.
In the text

	Fig. B.3. Fraction of time spent above the SFMS for currently star-forming galaxies when the current separation between each galaxy and the relation is imposed for all age values before adding noise. The colors correspond to the width of the random fluctuations that are added to this separation as a factor of the scatter of the SFMS.
In the text