Issue |
A&A
Volume 662, June 2022
|
|
---|---|---|
Article Number | A1 | |
Number of page(s) | 17 | |
Section | Extragalactic astronomy | |
DOI | https://doi.org/10.1051/0004-6361/202141338 | |
Published online | 25 May 2022 |
Polynomial expansion of the star formation history in galaxies
1
Departamento de Física Teórica, Universidad Autónoma de Madrid, Campus de Cantoblanco, Madrid 28049, Spain
e-mail: daniel.jimenezl@estudiante.uam.es
2
Australian Astronomical Optics, Macquarie University, 105 Delhi Rd, North Ryde, New South Wales 2113, Australia
Received:
18
May
2021
Accepted:
17
March
2022
Context. There are typically two different approaches to inferring the mass formation history (MFH) of a given galaxy from its luminosity in different bands. Non-parametric methods are known for their flexibility and accuracy, while parametric models are more computationally efficient.
Aims. In this work we propose an alternative, based on a polynomial expansion around the present time, that combines the advantages of both techniques.
Methods. In our approach, the MFH is decomposed through an orthonormal basis of N polynomials in lookback time. To test the proposed framework, synthetic observations are generated from models based on common analytical approximations (exponential, delayed-τ, and Gaussian star formation histories), as well as cosmological simulations for the Illustris-TNG suite. A normalized distance is used to measure the quality of the fit, and the input MFH is compared with the polynomial reconstructions both at the present time and through cosmic evolution. Our polynomial expansion is also compared with widely used parametric and non-parametric methods such as CIGALE and PROSPECTOR.
Results. The observed luminosities are reproduced with an accuracy of around 10 per cent for a constant star formation rate (N = 1) and better for higher-order polynomials. Our method provides good results on the reconstruction of the total stellar mass, the star formation rate, and even its first derivative for smooth star formation histories, but it has difficulties in reproducing variations on short timescales and/or star formation histories that peak at the earliest times of the Universe.
Conclusions. The polynomial expansion appears to be a promising alternative to other analytical functions used in parametric methods, combining both speed and flexibility.
Key words: methods: statistical / galaxies: star formation / galaxies: fundamental parameters
© ESO 2022
1. Introduction
One of the fundamental challenges in extragalactic astronomy is the reconstruction of the star formation history (SFH) of a given galaxy from the observed spectral energy distribution (SED). This is an inverse problem, where the emission at different wavelengths is considered to be the sum over simple stellar populations (SSPs) of different ages and metallicities, with some contribution from the nebular gas.
There are, broadly speaking, two different approaches to carrying out such a reconstruction. So-called non-parametric methods tackle the problem via a ‘brute-force’ approach, decomposing the SFH as a discrete sum over many individual SSPs with uncorrelated coefficients, with the only constraint that they all have to be positive (e.g. Heavens et al. 2000; Cid Fernandes et al. 2005; Ocvirk et al. 2006; Tojeiro et al. 2007; Koleva et al. 2009; MacArthur et al. 2009; Sánchez et al. 2016; Gomes & Papaderos 2017; Cardoso et al. 2019).
Alternatively, parametric methods describe the SFH in terms of a simple analytical function with only a few free parameters. A declining exponential is arguably the simplest one, and it has been widely applied to observations of early-type galaxies (e.g. Tinsley 1972). However, it was soon realized that many objects in the local Universe may display more complex behaviours (e.g. Talbot & Arnett 1971), and many are indeed better described by a monotonically increasing exponential (e.g. Maraston et al. 2010) or a delayed-τ model of the form (e.g. Lee et al. 2010). Some authors (e.g. Gladders et al. 2013) argue that additional degrees of freedom are necessary in order to capture the diversity of real SFHs and advocate the use of a Gaussian function with two free parameters, whereas other studies are based on a library of exponential SFHs with superimposed star formation bursts (e.g. Kauffmann et al. 2003; Walcher et al. 2008), variable galaxy age (e.g. Noll et al. 2009), and/or quenching of the star formation activity at a certain time (e.g. Smethurst et al. 2017).
Non-parametric methods tend to be more flexible and accurate, but they are computationally demanding and subject to numerous degeneracies due to the large number of free variables involved (see e.g. Walcher et al. 2011; Conroy 2013; Lower et al. 2020). Parametric methods, on the other hand, would be advantageous in this respect, but their simplicity makes it difficult to explore the whole range of possible SFHs.
In the present work we propose a new kind of parametric model that attempts to combine the advantages of both approaches. We model the stellar mass buildup of a galaxy in terms of the mass formation history (MFH),
where Ψ(t) represents the usual SFH in terms of the instantaneous star formation rate (SFR). We note that M(t) is the total mass of stars formed up to time t rather than the actual stellar mass, M*(t), which takes stellar evolution into account. Our main idea is to expand the MFH around the present time as a polynomial series based on a normalized lookback time,
that varies from at the present time, t = t0, to
at the Big Bang, t = 0. This provides a simple model where the free parameters have a straightforward physical interpretation (the total stellar mass formed, the current SFR, and its time derivatives) and whose number can be arbitrarily set by the user according to the amount and quality of the available data.
The details of our polynomial expansion are thoroughly presented in Sect. 2, and the method is calibrated against a set of synthetic observations based on analytical and numerical SFHs described in Sect. 3. Section 4 evaluates the accuracy of the reconstruction, and the main caveats and potential improvements are discussed in Sect. 5. Our main conclusions are briefly summarized in Sect. 6.
2. Mathematical formalism
In terms of , the MFH can be written as
We made the simplifying assumption that galaxies do not lose mass: neglecting stellar lifetimes as well as tidal stripping, must be a monotonically decreasing function of
. We expanded this function around the present time as a polynomial series in the normalized lookback time,
. Imposing that M(1) = 0 at the beginning of the Universe, we considered N polynomials
that are a linear combination of primordial terms of the form
where i varies from 0 to N − 1, and the coefficients βn, i are set by ensuring that the elements Bi form an orthonormal basis that satisfies
according to a suitable definition of the scalar product.
In the end, it is our aim to find an MFH that matches a given set of observables. More precisely, we considered Nν independent luminosity measurements, Lobs(ν), in different frequency bands, ν. Therefore, we needed to compute the luminosities, , of the N primordial polynomials,
, as well as those of the basis elements,
We then parameterized the MFH of a given galaxy as a linear combination,
of our basis elements and then minimized the distance,
between the predicted luminosities,
and the observed values, Lobs(ν). Minimizing with respect to the coefficients δi,
implies
if we define the scalar product as
We note that the number, N, of basis vectors must be equal to or smaller than the number, Nν, of observables for the system to be determined. If N > Nν, infinitely many degenerate solutions are possible, and thus the proposed method cannot be applied. When N equals Nν, the basis will uniquely cover the whole observable space, and any data will be exactly reproduced (d = 0) by one and only one model. If N < Nν, only a subspace of all possible observables will be accessible as a linear combination of the basis elements, and the normalized distance to the observational measurements will be, in general, larger than zero. As long as N ≤ Nν, the scalar product will yield the coefficients of the polynomial MFH whose luminosities are as close as possible to the observed values. However, there is no guarantee that such an optimal fit will be physically meaningful, especially if the actual MFH cannot be well approximated by a polynomial.
In particular, it is perfectly possible that the optimal reconstruction (in terms of best reproducing the observed luminosities) yields negative star formation for arbitrary periods of time. For that reason, we looked for the zeroes of and used them to divide the MFH into separate intervals with a constant sign. Then, we computed their individual luminosities (using the same polynomial coefficients, but only considering star formation within the selected time interval) and rescaled their MFHs by a constant factor to provide the best possible fit. This, of course, ensures positive star formation, but the quality of the fit necessarily worsens compared to the original result obtained for
. We repeated this process for all polynomial degrees N ≤ Nν and identified the values of
,
, and N that provide the best match to the observed data.
3. Synthetic observations
To test the ability of our formalism to reconstruct a realistic non-polynomial MFH, we simulated the observed luminosities for a set of synthetic galaxies. We first considered three different analytical SFHs that are fairly representative of real galaxies and have often been used to fit observational data: (i) a simple exponential, Ψ(t)∝e−t/τ, with characteristic time τ; (ii) a delayed-τ model of the form ; and (iii) a Gaussian,
, where α denotes the peak age of the SFH (i.e. Ψ peaks at a lookback time t0 − α) and
sets its full-width half maximum (FWHM).
We varied all our timescales from 0.1 to 13.5 Gyr. Short values of τ represent very old galaxies that formed the vast majority of their stars in the early Universe, while a Gaussian with a FWHM much shorter than t0 would describe a star formation burst of age α. On the other hand, large values of τ or FWHM yield smooth SFHs that can be well approximated by a low-order polynomial.
In addition, we also applied our algorithm to the SFHs of 22 155 galaxies selected from the IllustrisTNG-100-1 simulation, part of the IllustrisTNG magnetohydrodynamic cosmological simulations suite1 (Nelson et al. 2017; Springel et al. 2017; Pillepich et al. 2017; Naiman et al. 2018; Marinacci et al. 2018). They were run with the moving-mesh AREPO code (Springel 2010) and include a vast range of sub-grid physical processes that may play a critical role in galaxy evolution, such as gas cooling, star formation, supermassive black hole or supernova growth, and active galactic nucleus feedback (see Weinberger et al. 2018; Pillepich et al. 2018, for details). The TNG100-1 run consists of a cubic volume with box length ∼110 Mpc at z = 0 and dark mater and baryonic mass resolutions of 7.5 × 106 M⊙ and 1.4 × 106 M⊙, respectively, and we selected all galaxies with stellar masses in the range 109 < M*/M⊙ < 1012.
For every synthetic galaxy, we computed the luminosities, Lsyn(ν), in the ugriz photometric system (Nν = 5) of the Sloan Digital Sky Survey (SDSS; York et al. 2000) as an integral,
over SSPs whose specific luminosities per unit stellar mass, ℒSSP, have been obtained from the POPSTAR evolutionary synthesis models (Mollá et al. 2009; Martín-Manjón et al. 2010; García-Vargas et al. 2013; Millán-Irigoyen et al. 2021). These models provide a grid of 106 SSP ages, ranging from 105 yr to ∼16 Gyr, and we used the results appropriate for solar metallicity (Z = 0.02) and a Kroupa (2001) initial mass function (IMF).
The ugriz luminosities of our primordial MFHs are quoted in Table 1. We note that the values for n > 0 are all negative because our primordial terms are not monotonically decreasing functions of , and thus the instantaneous SFRs,
, may become negative. As mentioned above, this is an important feature of our models: in order to take advantage of the simple and computationally efficient scalar product, one cannot enforce positivity of the SFR at all times, which is a major source of unphysical results. The precise coefficients that yield orthonormal bases
of degree N up to five in this particular photometric system are quoted in Table 2.
Specific luminosities on the ugriz photometric system (in erg s−1 ) associated with the primordial polynomial MFH
up to N = 5 for the POPSTAR SSPs.
Coefficients of basis polynomials up to N = 5.
4. Results
We first evaluated the ability of our polynomial description to reproduce the luminosities and colours of realistic galaxies by comparing the synthetic observations with the best possible fits that can be achieved with polynomials of different order. Our ultimate goal is, however, to recover the true underlying MFH, which will be addressed separately in Sect. 4.2 and compared with the results obtained for the Illustris simulations by the CIGALE and PROSPECTOR inference codes in Sect. 4.3.
4.1. Fit quality
To quantify the agreement between the synthetic measurement and the best polynomial fit, we computed the normalized distance, , as
By definition, if the actual MFH were indeed a polynomial of degree Ntrue, a perfect fit with d = 0 would be obtained by our procedure as long as N ≥ Ntrue. Otherwise, since we only have five independent observables, N = 5 will always be able to obtain a perfect fit, although this does not necessarily imply that the polynomial that best fits the observation provides a good description of the actual MFH. The polynomial basis with N < 5 will only cover a subset (hyperplane) of the whole five-dimensional space of all possible observations. If the true MFH is close to a polynomial, its luminosities will be close to this hyperplane. On the other hand, if the MFH is very different from a polynomial, such as a short star formation burst, the observed luminosities may or may not be far away from the hyperplane covered by the linear combinations of the basis polynomials.
Figure 1 shows the normalized distance (Eq. (15)) as a function of the characteristic timescales of exponential and delayed-τ MFHs (the distances associated with the fifth degree polynomial are compatible with zero up to truncation errors). In general terms, higher-order polynomials always achieve an equal or better fit (lower distance) compared to previous orders. In addition, we also observe that the fit is usually better for longer timescales (i.e. a smoother MFH) than for those with narrower SFR periods. For cases where the highest polynomial degrees give rise to negative values of the SFR (typically narrow star formation bursts with small τ), the best positive-SFR fit will have distances comparable to N = 3. For smoother histories, with long values of τ, there are no periods of negative SFR, and the best fit is assigned to the N = 5 expansion.
![]() |
Fig. 1. Best fitting normalized distances associated with each polynomial degree, N (blue lines) and the best positive-SFR fit (red crosses) for different characteristic times of the exponential (upper panel) and delayed-τ (bottom panel) analytical models. |
A similar behaviour is shown for the Gaussian SFHs in Fig. 2, where we now represent the distances as a function of the two model parameters: the burst age and the FWHM. Smaller distances are found for higher-order polynomials and MFHs with wider FWHMs (smoother histories). White sections on the map represent values under (zero up to truncation errors). Once again, the best positive-SFR fit is similar to N ∼ 2 − 3 for SFHs that vary on short timescales and are equal to N = 5 for smooth functions that can be accurately described by a polynomial.
![]() |
Fig. 2. Best fitting normalized distances associated with each polynomial degree, N, and the best positive-SFR fit for different ages and FWHMs of the Gaussian synthetic SFHs. |
We note that, for some particular combinations of the input parameters, the fits obtained with N and N + 1 polynomials are basically identical for a certain N, which yields the sharp troughs observed for the delayed-τ model in Fig. 1 and the darker stripes on the Gaussian colour maps in Fig. 2. From the definition in Eq. (8), one can readily see that the MFH reconstructed by the N-th and (N + 1)-th degree polynomials are identical when the sum of the coefficients associated with vanishes, something that never occurs for the exponential law.
In order to illustrate the distances in luminosity space assigned to the results from the Illustris simulations, we plot the cumulative fraction of galaxies under a certain distance, , in Fig. 3. We observe that only 3.5% of the sample has purely positive SFRs for N = 5. The vast majority of galaxies have the best positive-SFR fit with an assigned distance close to the values expected for N = 2 − 3.
![]() |
Fig. 3. Fraction of galaxies with distances under a certain value, |
In general, we always find that for N ≥ 2. In other words, a linear SFH (a quadratic MFH) is able to reproduce the observed luminosities with an accuracy of the order of one per cent or even better. While it is obvious that this would indeed be, in most cases, an overly simplistic reconstruction of the underlying analytical model, the accuracy of the measured luminosity should be higher if we want to discriminate the difference with respect to a higher-order reconstruction.
4.2. Quality of the reconstruction
As stated above, our polynomial MFHs do not cover all possible histories, but just a subset that would be most appropriate for smooth functions. Then, one can imagine that SFHs that vary on short timescales may be harder to reproduce than galaxies that build up their mass more uniformly over cosmic time.
The reconstruction of the stellar mass, the SFR, and its first and second derivatives at the present time, (i.e. the physical quantities upon which the polynomial expansion is based), are shown in Fig. 4 for different values of the characteristic timescale, τ, for the exponential and delayed-tau models. The total mass is well reproduced (always better than 20 per cent) for N > 2, regardless of the value of τ. For timescales above a few gigayears, the recovery is also good for N = 2; only the first-order polynomial (i.e. a constant SFR) differs from the true mass by more than a few per cent. We note, however, that, according to Fig. 1, one may rule out N = 1 in favour of N = 2 as long as the relative uncertainty of the observational measurements remains below ∼10 per cent. The current instantaneous SFR may also be accurately estimated for large τ even with the simplest polynomials. For τ < 1 Gyr, N = 1 and N = 2 fail to yield an appropriate model, but N ≥ 3 is still able to qualitatively reproduce the mass and SFR of the input model at the present time, although the low SFRs are impossible to recover accurately as they approach zero. For short timescales, the results of our polynomial expansion should not be taken at face value. On the other hand, when τ becomes comparable to the age of the Universe, it is notable that N = 4 and N = 5 can even provide a reasonable estimate of the time derivative of the SFR, although their description of the second derivative is, at best, only a very crude approximation. The best positive-SFR fit sets the SFR (and therefore its derivatives) to zero for
, and thus it actually helps recover a better estimation of the derivatives at present time for low values of the characteristic timescale, τ, than any polynomial degree.
![]() |
Fig. 4. Estimates, from top to bottom, of the total stellar mass formed, the SFR, and its first and second time derivatives, reconstructed at the present time ( |
These trends are also observed for the Gaussian models in Fig. 5, where the results are shown for a range of values for the age and FWHM of the stellar burst. The polynomial description is able to capture the mass and the SFR, and it even crudely estimates its derivatives at the present time for large FWHM values. Unfortunately, when this parameter drops below a few gigayears, the MFH becomes increasingly similar to a short star formation burst, not amenable to a description in terms of smooth polynomials, and the estimates provided by our unconstrained polynomials cannot be trusted, even if the observed luminosities are accurately reproduced. The best positive-SFR fit automatically selects low-order polynomials (N ∼ 3) for short FWHMs, and it always provides meaningful results, although it struggles to reproduce the SFHs with FWHM < 1 Gyr. The mass and SFR returned by our method are fairly accurate, but the first derivative of the SFR is always difficult to recover and must be interpreted, at best, in a qualitative sense. The second time derivative is never accurately recovered.
![]() |
Fig. 5. Estimates, from left to right, of the total stellar mass formed, the SFR, and its first and second time derivatives, reconstructed at the present time ( |
In contrast to the analytical models, the SFHs of the simulated galaxies are inherently stochastic due to the finite mass of the stellar particles, and therefore the instantaneous SFR (let alone its derivatives) are not trivially defined and are severely affected by Poisson noise. Thus, we assessed the quality of our reconstruction in terms of the relative difference
between the mass formed since a certain
in our model and the simulation. We note that −1 ≤ δM ≤ 1, where δM = 0 indicates a perfect reconstruction, while δM = −1 for ΔMpoly ≤ 0 and δM = 1 for ΔMsim = 0 and ΔMsim > 0. When our best positive-SFR fit correctly identifies that no mass was formed in the simulation during the requested interval, we set δM = 0.
The probability distribution of δM at is represented in Fig. 6. These time intervals trace stellar populations of very different ages, and we think they capture most of the relevant information about the SFH of the galaxy. We find that the mass formed in each interval is recovered with an accuracy of the order of 0.2 dex, although our unconstrained polynomials yield negative masses for a sizeable fraction of the objects. By construction, the best positive-SFR fit completely eliminates the peak at δM = −1, which is associated with these null or negative masses, and it helps detect quenched galaxies, bringing them from δM = 1 to δM = 0.
![]() |
Fig. 6. Probability densities for different values of |
One can directly extrapolate our reconstruction towards the past in an attempt to infer the full time evolution of the MFH from the observed luminosities. From Figs. A.1 and A.2, we find that the polynomial expansion provides a good description for exponential and delayed-τ models with high values of the characteristic timescale, τ (above several gigayears), whereas the fit is not as good for values much shorter than the age of the Universe, even for the highest N degrees.
We can also appreciate a similar behaviour for Gaussian SFHs. Figures A.3, A.4, and A.5 show the MFH, the SFR, and its first derivative as a function of for different values of the FWHM and the peak age of the Gaussian. In general, the accuracy of the polynomial reconstructions is very sensitive to these parameters. Very narrow SFHs can never be well reproduced, not even by the highest degree polynomials. At best, one can recover qualitative information, such as an approximate estimate of the height and width of the star formation peak. The polynomial fit becomes a better description as the FWHM increases above a few gigayears, and the SFH is fairly well reproduced, especially for the higher-order polynomials (perhaps excluding the extrapolation towards the earliest moments of the Universe, although this problem is significantly alleviated for the best positive-SFR fit).
For a small sample of randomly selected Illustris galaxies covering a broad mass range, Fig. A.6 exemplifies how a constant SFR (N = 1) tends to underestimate the mass-to-light ratio and thus the reconstructed mass. The overall buildup of stellar mass in recent times is nicely recovered for N > 1; however, extrapolating back to the first gigayear of cosmic evolution is always problematic, and our simple polynomials often feature unphysical oscillations. The best positive-SFR prescription fully corrects for this undesirable behaviour and yields MFHs that faithfully reproduce the simulation results at all times. The instantaneous SFR (Fig. A.7) is also accurately recovered, even for early cosmic times. At recent epochs, the smooth polynomial solution is representative of the average discrete SFR of the simulated galaxies, and our best positive-SFR fit is not only able to successfully identify quenched galaxies, but also to provide an order-of-magnitude estimate of their death time (Corcho-Caballero et al. 2021).
4.3. Comparison with other methods
In order to put our results into context, we compared them with two other algorithms representative of state-of-the-art parametric and non-parametric algorithms. More precisely, we fitted our Illustris galaxies with the publicly available codes CIGALE and PROSPECTOR, which are widely used in the literature.
CIGALE (Code Investigating GALaxy Emission; Boquien et al. 2019) is an energy-balance code developed for modelling the observed galaxy SEDs from X-rays to radio wavelengths. This package applies a Bayesian strategy through the analysis of the likelihood distribution that takes into account the uncertainties on the observations and the effect of intrinsic degeneracies between physical parameters. It is characterized by a series of modules that (i) build stellar population models, (ii) consider the dust absorbing and re-emitting in the infrared wavelengths, (iii) add non-thermal sources of dust emission, (iv) include the interstellar lines, (v) take redshift effects into account, (vi) calculate best-fit models, and (vii) estimate weighted galaxy properties such as the SFR, attenuation, dust luminosity, stellar mass, and many other physical quantities.
We first generated synthetic photometry in the ugriz bands from the simulated SFHs using the Bruzual & Charlot (2003) stellar basis with a Chabrier (2003) IMF and metallicity Z = 0.02. Then, the synthetic luminosities were fit to a grid of models with delayed-τ SFHs, typically used to reconstruct the SFHs of large observational galaxy samples (e.g. Noll et al. 2009; Ciesla et al. 2015). The characteristic timescale was varied from 0.1 to 7 Gyr in linear steps of 0.1 Gyr (a total of 70 models), and the age of the stellar population (i.e. the onset of star formation) was 13.7 Gyr.
The code spent about 0.1 − 0.6 s fitting each galaxy on a standard computer using four cores. Most of that time was actually invested in input and output operations. The fit itself is negligible in comparison, and in this sense there in no practical difference in terms of computational time with respect to our method, both methods being virtually instantaneous for the simple setup considered in the present work. In terms of accuracy, one can clearly see in Fig. 6 that results obtained with the delayed-τ model are comparable to our polynomial fits for the total stellar mass, although the dispersion around δM = 0 appears to increase for the mass formed on shorter timescales.
While parametric models may indeed provide a reasonable description of the SFHs, their limited number of degrees of freedom prevents them from fully capturing the diversity found in individual galaxies (see e.g. Lower et al. 2020, and references therein). On the other hand, non-parametric models enable a more thorough exploration of the possible SFHs at the expense of a significantly higher computation time. Here, we calibrated the performance of our algorithm against the Bayesian inference code PROSPECTOR (Johnson & Leja 2017; Johnson et al. 2019, 2021), which performs a Monte Carlo sampling to compute the posterior parameter distribution of relatively complex models. We used the Flexible Stellar Population Synthesis (FSPS) stellar population synthesis models (Conroy et al. 2009; Conroy & Gunn 2010) to compute the synthetic luminosities of the Illustris galaxies, setting a default signal-to-noise ratio for each photometric band of S/N = 100. These data were then fit with the default continuity non-parametric model, which estimates the difference Δ log(SFR) between adjacent time bins, weighted by a Student’s-t distribution to prevent sharp changes (see Leja et al. 2019, for details). We chose N = 7 adaptive time bins and the same default parameters as discussed in Leja et al. (2019). Dust extinction and metallicity were fixed to AV = 0 and Z = Z⊙, respectively.
PROSPECTOR spends of the order of ∼20 ± 10 minutes fitting an individual galaxy. Running the code on the full sample would have been computationally unfeasible, and therefore we randomly selected 373 objects with the intent of sampling the whole mass (9 ≤ log10(M*/M⊙)≤11) and specific star formation rate (sSFR) range (− 12 ≤ log10(sSFR/yr−1) < − 9) as uniformly as possible, verifying by visual inspection that continuous SFHs, as well as suddenly quenched galaxies with different death times, are duly represented. One can see in Figs. A.6 and A.7 that, in general, PROSPECTOR provides a good fit to the MFH, and it has sufficient flexibility to accurately recover variations in the instantaneous SFR in recent epochs. As shown in Fig. 6, the differences δM with respect to the mass formed in different time intervals is comparable to the results obtained by our polynomial reconstruction. If anything, the best positive-SFR fit seems to be slightly more successful in correctly identifying recently quenched galaxies.
5. Discussion
Our results show that the proposed polynomial expansion is able to successfully reproduce representative smooth SFHs from photometric observables, as long as most of the stellar mass is formed over long characteristic timescales (τ or FWHM comparable to the age of the Universe). On the other hand, the method has difficulties when the MFH features abruptly changes or peaks in the early Universe. In any case, it offers a compromise between the stability, computational efficiency, and straightforward physical interpretation of parametric methods and the flexibility of the non-parametric approach. Technically, the fit only involves a scalar product, providing a significant advantage with respect to non-parametric methods in terms of computing time. The coefficients of the expansion are trivially related to the current mass, the SFR, and its derivatives with respect to time at the moment the object was observed.
Extrapolating backwards, the fit provides a fairly accurate reconstruction of the SFH when it is indeed smooth, although the unconstrained polynomials may provide unphysical negative values, especially at higher orders. For that reason, we also computed the best positive-SFR fit, where a finite time interval is identified where star formation is positive and definite. No stars are formed for t < tmin nor t < tmax. Our results show that this model is capable of reproducing not only simple analytical SFHs, but also the complex MFH found in cosmological numerical simulations. The accuracy of the best positive-SFR fit rivals, or even surpasses, that of state-of-the-art parametric and non-parametric algorithms such as PROSPECTOR and CIGALE, involving a negligible computational burden. A particularly relevant aspect where this prescription has proven to excel is the identification of quenched galaxies from the synthetic photometry. As discussed by Corcho-Caballero et al. (2021), many of the simulated galaxies have not formed any stars over the last few gigayears, while it is at present unclear whether such ‘quenching’ of the star formation activity is observed in the real Universe (cf. Ciesla et al. 2018; Aufort et al. 2020; Corcho-Caballero et al. 2020). Based on the distribution of the relative difference δM shown in Fig. 6, we argue that the best positive-SFR fit provides a promising alternative for tackling this particular problem.
Nevertheless, there are also important caveats and potential improvements to the scheme proposed here. First and foremost, observational uncertainties must be taken into account in order to provide error estimates on the recovered quantities. Typical errors in the SDSS fluxes are of the order of 2% in the u band and 1% in g, r, i, and z. As noted above, this is larger than the values of obtained for N ≥ 2. Therefore, it is unlikely that meaningful results will be obtained by increasing the polynomial degree, N, further than this limit. This, in turn, suggests that constraining even the first derivative of the SFR will be difficult in practice. A more quantitative assessment will be carried out in future work.
In addition, our tests do not consider the effects of metallicity and dust extinction. Treating these non-linear parameters represents a major improvement to the proposed framework, but it would be necessary to apply it to real observational data. A Bayesian approach is arguably the most appropriate way of quantifying the posterior probability distribution of the model parameters, at the expense of computing time. A significant advantage would be that the priors would make it possible to supplement our best positive-SFR fit with additional constraints, such as assigning less (or zero) weight to the individual polynomial terms that provide a significant amount of negative mass formation. More elaborate priors based on chemical evolution constraints could also be explored, and the initial and final times, tmin and tmax, could be left to vary as two additional free parameters. These would again behave non-linearly, but it is straightforward to include them within a Bayesian framework.
Another trivial, albeit important, extension would be to consider multi-wavelength photometric observations in other bands (e.g. infrared and/or ultraviolet) that contain valuable information, as well as spectroscopic data. Including additional observables, appropriately weighted to yield a combined likelihood, would make it possible to increase and discriminate between models that are very similar in the optical colours considered in the present work. This is arguably the most promising avenue for increasing the maximum degree, N, of the polynomial basis in a practical application of the algorithm.
6. Conclusions
Here we propose an analytical description of the MFH of a given galaxy in terms of a polynomial expansion around the present time. Our tests against a set of synthetic observations in the SDSS ugriz photometric system of exponential, delayed-τ, and Gaussian SFHs, as well as a set of galaxies from the Illustris-TNG suite of cosmological simulations, show that:
-
The polynomial reconstruction of order N = 1 (constant SFR) is able to reproduce the synthetic luminosities with an accuracy of the order of 10 per cent or better. The accuracy improves with N, with higher-order polynomials reaching values that are always below 1 per cent.
-
For smooth SFHs that vary on timescales comparable to the age of the Universe, the polynomial expansion does indeed provide a faithful reconstruction of the actual SFH, especially at recent times. On the other hand, our approach is not well suited to describe an SFH that varies on short timescales or peaks in the very early Universe, yielding negative stellar masses and/or SFRs.
-
Such unphysical results can be avoided by selecting the best positive-SFR fit, where star formation is always positive within the time interval
and vanishes outside this range. This prescription helps polynomial models recover a more accurate estimation of the mass formed at the most recent times and identify quenched galaxies in the numerical simulations.
-
The results of the best positive-SFR fit are highly competitive with state-of-the-art parametric and non-parametric methods, both in terms of accuracy and computational speed.
We conclude that the proposed polynomial expansion represents a promising alternative to the analytical functions that have been traditionally used in parametric methods. The polynomial description is more flexible than the exponential or delayed-τ models, and the number of free parameters can be trivially adapted to the amount and quality of the available data. Its physical interpretation in terms of the current total mass, SFR, and higher-order time derivatives is straightforward, and its linear nature is amenable to an extremely fast computation, which could in principle be applied to large photometric surveys on a pixel-by-pixel basis.
Nonetheless, further work is still necessary before our approach can be applied to real measurements. Most importantly, the effects of metallicity, dust extinction, and observational uncertainties must be taken into account, and we will explore a Bayesian extension of the present formalism in the near future. This powerful statistical framework will enable us to test the ability of the polynomial expansion to tackle the problem under realistic conditions.
Acknowledgments
We thank the anonymous referee for a helpful, constructive report, that has encouraged us to develop the best positive-SFR fit, add numerical simulations to our test suite, and compare our results with other codes in the literature. This work has been funded by the Spanish Research Agency (grant PID2019-107408GB-C42/AEI/10.13039/501100011033). S.Z. also acknowledges support from contract BES-2017-080509 associated to grant AYA2016-79724-C4-1-P from the Spanish Ministry of Economy, Industry and Competitiveness through the MINECO-FEDER research.
References
- Aufort, G., Ciesla, L., Pudlo, P., & Buat, V. 2020, A&A, 635, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
- Cardoso, L. S. M., Gomes, J. M., & Papaderos, P. 2019, A&A, 622, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
- Cid Fernandes, R., Mateus, A., Sodré, L., Stasińska, G., & Gomes, J. M. 2005, MNRAS, 358, 363 [Google Scholar]
- Ciesla, L., Charmandaris, V., Georgakakis, A., et al. 2015, A&A, 576, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ciesla, L., Elbaz, D., Schreiber, C., Daddi, E., & Wang, T. 2018, A&A, 615, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Conroy, C. 2013, ARA&A, 51, 393 [NASA ADS] [CrossRef] [Google Scholar]
- Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833 [Google Scholar]
- Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486 [Google Scholar]
- Corcho-Caballero, P., Ascasibar, Y., & López-Sánchez, Á. R. 2020, MNRAS, 499, 573 [NASA ADS] [CrossRef] [Google Scholar]
- Corcho-Caballero, P., Ascasibar, Y., & Scannapieco, C. 2021, MNRAS, 506, 5108 [NASA ADS] [CrossRef] [Google Scholar]
- García-Vargas, M. L., Mollá, M., & Martín-Manjón, M. L. 2013, MNRAS, 432, 2746 [CrossRef] [Google Scholar]
- Gladders, M. D., Oemler, A., Dressler, A., et al. 2013, ApJ, 770, 64 [NASA ADS] [CrossRef] [Google Scholar]
- Gomes, J. M., & Papaderos, P. 2017, A&A, 603, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Heavens, A. F., Jimenez, R., & Lahav, O. 2000, MNRAS, 317, 965 [NASA ADS] [CrossRef] [Google Scholar]
- Johnson, B., & Leja, J. 2017, https://doi.org/10.5281/zenodo.1116491 [Google Scholar]
- Johnson, B. D., Leja, J. L., Conroy, C., & Speagle, J. S. 2019, Prospector: Stellar population inference from spectra and SEDs Astrophysics Source Code Library, [record ascl:1905.025] [Google Scholar]
- Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22 [NASA ADS] [CrossRef] [Google Scholar]
- Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 33 [Google Scholar]
- Koleva, M., De Rijcke, S., Prugniel, P., Zeilinger, W. W., & Michielsen, D. 2009, MNRAS, 396, 2133 [CrossRef] [Google Scholar]
- Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
- Lee, S.-K., Ferguson, H. C., Somerville, R. S., Wiklind, T., & Giavalisco, M. 2010, ApJ, 725, 1644 [NASA ADS] [CrossRef] [Google Scholar]
- Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3 [Google Scholar]
- Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33 [NASA ADS] [CrossRef] [Google Scholar]
- MacArthur, L. A., González, J. J., & Courteau, S. 2009, MNRAS, 395, 28 [NASA ADS] [CrossRef] [Google Scholar]
- Maraston, C., Pforr, J., Renzini, A., et al. 2010, MNRAS, 407, 830 [NASA ADS] [CrossRef] [Google Scholar]
- Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113 [NASA ADS] [Google Scholar]
- Martín-Manjón, M. L., García-Vargas, M. L., Mollá, M., & Díaz, A. I. 2010, MNRAS, 403, 2012 [CrossRef] [Google Scholar]
- Millán-Irigoyen, I., Mollá, M., Cerviño, M., et al. 2021, MNRAS, 506, 4781 [CrossRef] [Google Scholar]
- Mollá, M., García-Vargas, M. L., & Bressan, A. 2009, MNRAS, 398, 451 [CrossRef] [Google Scholar]
- Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206 [Google Scholar]
- Nelson, D., Pillepich, A., Springel, V., et al. 2017, MNRAS, 475, 624 [Google Scholar]
- Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ocvirk, P., Pichon, C., Lançon, A., & Thiébaut, E. 2006, MNRAS, 365, 74 [Google Scholar]
- Pillepich, A., Nelson, D., Hernquist, L., et al. 2017, MNRAS, 475, 648 [Google Scholar]
- Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077 [Google Scholar]
- Sánchez, S. F., Pérez, E., Sánchez-Blázquez, P., et al. 2016, Rev. Mex. Astron. Astrofis., 52, 21 [NASA ADS] [Google Scholar]
- Smethurst, R. J., Lintott, C. J., Bamford, S. P., et al. 2017, MNRAS, 469, 3670 [NASA ADS] [CrossRef] [Google Scholar]
- Springel, V. 2010, MNRAS, 401, 791 [Google Scholar]
- Springel, V., Pakmor, R., Pillepich, A., et al. 2017, MNRAS, 475, 676 [Google Scholar]
- Talbot, R. J., & Arnett, W. D. 1971, ApJ, 170, 409 [NASA ADS] [CrossRef] [Google Scholar]
- Tinsley, B. M. 1972, ApJ, 178, 319 [NASA ADS] [CrossRef] [Google Scholar]
- Tojeiro, R., Heavens, A. F., Jimenez, R., & Panter, B. 2007, MNRAS, 381, 1252 [NASA ADS] [CrossRef] [Google Scholar]
- Walcher, C. J., Lamareille, F., Vergani, D., et al. 2008, A&A, 491, 713 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, Ap&SS, 331, 1 [NASA ADS] [CrossRef] [Google Scholar]
- Weinberger, R., Springel, V., Pakmor, R., et al. 2018, MNRAS, 479, 4056 [Google Scholar]
- York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Reconstructed histories
The reconstructed MFHs, M(t), SFHs, Ψ(t) ≡ Ṁ(t), and their time derivatives, and
, are illustrated in Figs. A.1 and A.2 for the exponential and delayed-τ models, respectively, with characteristic timescales τ = {0.5, 3.5, 7.5, 13.5} Gyr. Figures A.3, A.4, and A.5 show M, Ψ, and
, respectively, for the Gaussian models with representative values of their peak age and FWHM. Finally, Figs. A.6 and A.7 show the reconstructed histories from the Illustris-TNG sample, as well as the comparison with the results obtained with CIGALE and PROSPECTOR.
![]() |
Fig. A.1. Sample of a MFHs reconstructed with our parametric model from the luminosities of a synthetic exponential SFR with different timescales, τ. With blue lines we plot the stellar mass formed (first row), the SFR (second row), the SFR first derivative (third row), and the SFR second derivative (fourth row), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N, respectively). Red crosses correspond to the best positive-SFR fit. |
![]() |
Fig. A.2. Sample of MFHs reconstructed with our parametric model from the luminosities of a synthetic exponential-delayed SFR with different timescales, τ. Using blue lines, we plot the stellar mass formed (first row), the SFR (second row), the SFR first derivative (third row), and the SFR second derivative (fourth row), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N, respectively). Red crosses correspond to the best positive-SFR fit. |
![]() |
Fig. A.3. Mass formation history reconstructed with our parametric model from the luminosities of a synthetic Gaussian SFR with different peak ages (increasing from top to bottom) and FWHMs (increasing from left to right), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N). Red crosses correspond to the best positive-SFR fit. |
![]() |
Fig. A.4. Star formation rate reconstructed with our parametric model from the luminosities of a synthetic Gaussian SFR with different peak ages (increasing from top to bottom) and FWHMs (increasing from left to right), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N). Red crosses correspond to the best positive-SFR fit. |
![]() |
Fig. A.5. Time derivative of the SFR reconstructed with our parametric model from the luminosities of a synthetic Gaussian SFR with different peak ages (increasing from top to bottom) and FWHMs (increasing from left to right), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N). Red crosses correspond to the best positive-SFR fit. |
![]() |
Fig. A.6. Mass formation history for different galaxies of the Illustris sample (black line). The polynomial expansion reconstruction is shown as the different solid (odd degree) and dashed (even degree) blue lines. Results from PROSPECTOR are plotted as green circles, and the results from CIGALE are plotted as yellow triangles. Red crosses correspond to the best positive-SFR fit. |
![]() |
Fig. A.7. Star formation rate for different galaxies of the Illustris sample (black line). The polynomial expansion reconstruction is shown as the different solid (odd degree) and dashed (even degree) blue lines. Results from PROSPECTOR are plotted as green circles, and the results from CIGALE are plotted as yellow triangles. Red crosses correspond to the best positive-SFR fit. |
All Tables
Specific luminosities on the ugriz photometric system (in erg s−1 ) associated with the primordial polynomial MFH
up to N = 5 for the POPSTAR SSPs.
All Figures
![]() |
Fig. 1. Best fitting normalized distances associated with each polynomial degree, N (blue lines) and the best positive-SFR fit (red crosses) for different characteristic times of the exponential (upper panel) and delayed-τ (bottom panel) analytical models. |
In the text |
![]() |
Fig. 2. Best fitting normalized distances associated with each polynomial degree, N, and the best positive-SFR fit for different ages and FWHMs of the Gaussian synthetic SFHs. |
In the text |
![]() |
Fig. 3. Fraction of galaxies with distances under a certain value, |
In the text |
![]() |
Fig. 4. Estimates, from top to bottom, of the total stellar mass formed, the SFR, and its first and second time derivatives, reconstructed at the present time ( |
In the text |
![]() |
Fig. 5. Estimates, from left to right, of the total stellar mass formed, the SFR, and its first and second time derivatives, reconstructed at the present time ( |
In the text |
![]() |
Fig. 6. Probability densities for different values of |
In the text |
![]() |
Fig. A.1. Sample of a MFHs reconstructed with our parametric model from the luminosities of a synthetic exponential SFR with different timescales, τ. With blue lines we plot the stellar mass formed (first row), the SFR (second row), the SFR first derivative (third row), and the SFR second derivative (fourth row), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N, respectively). Red crosses correspond to the best positive-SFR fit. |
In the text |
![]() |
Fig. A.2. Sample of MFHs reconstructed with our parametric model from the luminosities of a synthetic exponential-delayed SFR with different timescales, τ. Using blue lines, we plot the stellar mass formed (first row), the SFR (second row), the SFR first derivative (third row), and the SFR second derivative (fourth row), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N, respectively). Red crosses correspond to the best positive-SFR fit. |
In the text |
![]() |
Fig. A.3. Mass formation history reconstructed with our parametric model from the luminosities of a synthetic Gaussian SFR with different peak ages (increasing from top to bottom) and FWHMs (increasing from left to right), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N). Red crosses correspond to the best positive-SFR fit. |
In the text |
![]() |
Fig. A.4. Star formation rate reconstructed with our parametric model from the luminosities of a synthetic Gaussian SFR with different peak ages (increasing from top to bottom) and FWHMs (increasing from left to right), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N). Red crosses correspond to the best positive-SFR fit. |
In the text |
![]() |
Fig. A.5. Time derivative of the SFR reconstructed with our parametric model from the luminosities of a synthetic Gaussian SFR with different peak ages (increasing from top to bottom) and FWHMs (increasing from left to right), compared to the input model (solid black line). Colour gradients and line styles indicate the degree of the polynomial reconstruction from N = 5 (solid) to lower degrees (more diffuse lines; solid and dashed for odd and even N). Red crosses correspond to the best positive-SFR fit. |
In the text |
![]() |
Fig. A.6. Mass formation history for different galaxies of the Illustris sample (black line). The polynomial expansion reconstruction is shown as the different solid (odd degree) and dashed (even degree) blue lines. Results from PROSPECTOR are plotted as green circles, and the results from CIGALE are plotted as yellow triangles. Red crosses correspond to the best positive-SFR fit. |
In the text |
![]() |
Fig. A.7. Star formation rate for different galaxies of the Illustris sample (black line). The polynomial expansion reconstruction is shown as the different solid (odd degree) and dashed (even degree) blue lines. Results from PROSPECTOR are plotted as green circles, and the results from CIGALE are plotted as yellow triangles. Red crosses correspond to the best positive-SFR fit. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.