Resolving cosmic star formation histories of present-day bulges, disks, and spheroids with ProFuse

The tool that we use to conduct the simultaneous spectral and photometric decomposition of our sources is ProFuse (Robotham et al., 2022). This combines the SED-fitting capabilities of ProSpect Robotham et al. (2020), the structural modelling capabilities of ProFit (Robotham et al., 2017), and the source-finding capabilities of ProFound (Robotham et al., 2018). While ProFuse as a tool in its entirety is relatively new, it has been developed over many years, through the gradual creation and application of its individual building blocks.

The first of these tools was ProFit (Robotham et al., 2017), which was developed to facilitate light profile fitting of galaxies in a generative, Bayesian manner. Previous codes had been prone to model error, making 2D light modelling of large samples of galaxies unwieldy, and the Bayesian nature of ProFit improved this behaviour. ProFit was wielded by studies such as Cook et al. (2019) and Cook et al. (2020) to assess the impact of galaxy bulges in xGASS (Catinella et al., 2010) on the scatter of HI scaling relations and the star forming main sequence, and later by Hashemizadeh et al. (2022) to study the evolution of bulges and disks through cosmic time in the DEVILS survey (Davies et al., 2018).

Separately to ProFit, there was a need to develop a generative and Bayesian tool that was capable of conducting SED fitting in a flexible and physically motivated manner. This inspired the development of ProSpect (Robotham et al., 2020). With a particular desire to flexibly extract unbiased star formation histories (which in turn would ensure more accurate properties such as total stellar mass), ProSpect was applied to a volume-limited sample by Bellstedt et al. (2020b), where the capacity for this code to forensically reconstruct the cosmic star formation history was demonstrated. This was only possible due to the physically motivated implementation of metallicity evolution within ProSpect, shown to realistically evolve metals in galaxies as demonstrated by the correct inferred evolution of the mass–metallicity relation (Bellstedt et al., 2021). This SED fitting tool was soon shown to work just as effectively at higher redshifts using the DEVILS survey by Thorne et al. (2021), with the link between metallicity and star formation histories explored by Thorne et al. (2022b). These studies demonstrated that ProSpect could self-consistently model the local Universe and infer the correct evolution of stellar mass, star formation rate and metallicity. This is as yet undemonstrated by other SED fitting codes, of which numerous exist with various degrees of similarity, e.g Prospector (Leja et al., 2017; Johnson et al., 2021), BAGPIPES (Carnall et al., 2018), MAGPHYS (da Cunha et al., 2008), CIGALE (Noll et al., 2009; Boquien et al., 2019), and BEAGLE (Chevallard & Charlot, 2016) (see Thorne et al. 2021 for a direct comparison of their utility and capabilities). This physically motivated implementation lends itself incredibly well to the incorporation of additional phenomena, as shown through the inclusion of an AGN component by Thorne et al. (2022a), which was even extended to include radio data in the fitting process (Thorne et al., 2023). Through application of this self-consistent stellar and AGN modelling, the link between star formation and AGN activity over 12.5 Gyr of cosmic time could be studied using GAMA and DEVILS by D’Silva et al. (2023).

The final implement in the ProFuse toolbag is a source detection routine, capable of identifying and isolating the flux of the source object, but also nearby stars to facilitate a local extraction of the image point spread function (PSF). This source extraction software is ProFound (Robotham et al., 2018), which has been readily used since its development to extract photometry of large, multi-wavelength imaging datasets for surveys (Bellstedt et al., 2020a; Davies et al., 2021).

While ProFuse is described in great detail by Robotham et al. (2022), we describe our implementation of ProFuse in the following sections.

It is essential to be explicit about the nomenclature adopted in our work, as the usage of structural terms varies across the literature. The three terms that we adpot are disks, bulges, and spheroids.

Disk terminology is consistent with the common use, describing any flattened, circular structure. Disks can either be one part of a two-component system, or galaxies can be purely disks. Bulge is used to describe the ellipsoidal structure at the centre of a two-component system. Spheroid, in our work, is used to describe the single-component structures that are ellipsoidal. Note that with this usage, bulges are not deemed to be a subset of spheroids (which is how the term is often used in the literature). In this sense, disks, bulges, and spheroids are not treated as overlapping categories, and are instead viewed as “eigenstructures” of galaxies. We note that our spheroids are a much broader class than the “elliptical” galaxy class.

Regardless of the manner in which a galaxy decomposition is conducted, an essential part of the process is determining which model configuration best describes a galaxy’s two-dimensional distribution of light. This can be the most challenging part of structural decomposition, as in many cases a single-component model can describe the light distribution of a galaxy just as well as a bulge plus disk model. What should then determine the choice? This choice has frequently been made using visual classifications or inspections (such as in component modelling by Gadotti, 2009). While a completely numerical quantifier is desired (to remove subjective visual decisions), this approach has remained elusive (although it has been described at length in works like Hashemizadeh et al., 2021).

Because reliable and purely numerical discriminators remain controversial, visual classifications have remained important in much morphology/structure-based work. As datasets are increasing in size dramatically though, requiring visual inspections of galaxies is becoming increasingly unfeasible. An alternate approach increasingly being used is machine learning, where effectively a computer is trained to quickly and efficiently replicate a visual classification on a very large scale. Many galaxy morphological classifications have been made in this manner (see for example Walmsley et al., 2022; Li et al., 2023; Fang et al., 2023; Tian et al., 2023). At this stage these classifications have not yet readily been used for much in the way of follow-on science, although an exception to this is Cavanagh et al. (2023), who used machine learning morphologies to study the evolution of lenticular galaxies with cosmic time. The benefit from a machine learning approach is that larger sample sizes can be processed. Machine learning morphologies will suffer the same potential biases that any visual classification scheme will, due do its inherently qualitative nature, and because nearly all machine learning classifictions have been trained on visually classified data.

For our work, we apply a numeric quantifier from our ProFuse outputs to quantitatively categorise our sample into bulges, disks, and spheroids. Acknowledging however that any classification scheme (whether numerical or visual) still bears some ambiguity, a key aim of this work is to demonstrate the inherent uncertainties that still accompany any attempt to classify galaxies into structural classes. We therefore also incorporate a number of different visual morphological classifications that have been conducted for this sample of galaxies. These different classification schemes are outlined in the following subsections.

A visual re-classification of the 6,664 galaxies of this work has been conducted (by SB) into structural classes that are equivalent to those defined by the automatic ProFuse classifications described in the previous section. The classes used in this classification are much broader and simpler than previously-used schemes, and are limited to three categories: disk, spheroid, and two-component galaxies. These have been called “$\rm{Visual}_{\rm{SB}}$ classes” throughout the paper.

These classes complement visual classifications that have been previously conducted on this sample by Moffett et al. (2016a) (“$\rm{Visual}_{\rm{Hubble}}$ Classes” hereafter), sorting the galaxies into their Hubble types. The Hubble classes used are ellipticals (E), S0-Sa, SB0-SBa, Sab-Scd, SBab-SBcd, Sd-Irr, and Little Blue Spheroid (LBS) classes.

As part of GAMA DR4, Driver et al. (2022b) conducted an additional visual classification that was based on the visual structure of these galaxies, rather than the visual morphology. We refer to these classifications as the “$\rm{Visual}_{\rm{DR4}}$ Classes” hereafter. The DR4 classes separate galaxies into pure disks (D), ellipticals (E), compact sources (C), and then two-component systems, either with a compact bulge (cBD), or with a diffuse bulge (dBD). Galaxies that are too difficult to classify due to irregular structure such as from mergers are labelled “hard” (H), or “hard elliptical” (HE). These account for a tiny fraction of the sample, and as such are not significant in this analysis.

When presenting results from visual classification methods throughout this paper, we use properties derived from our ProFuse modelling. When visual classes have described a galaxy as having a single-component fit, we assign the FS properties to the galaxy. When a two-component system has instead been selected, then we can either assign the BD or PD properties to the galaxy. Wherever relevant for the duration of this paper, we present results with both options, to indicate the plausible uncertainty originating from model selection. We highlight that this is the most pessimistic uncertainty, as a likelihood criterion (as used in the automatic ProFuse classes) can be used to isolate the better-fitting description to shrink this uncertainty. This same approach is taken when presenting the bulge and disk properties using the $\rm{Visual}_{\rm{DR4}}$ classes.

We show a comparison between our ProFuse model classifications and the $\rm{Visual}_{\rm{SB}}$ Classes, $\rm{Visual}_{\rm{DR4}}$ Classes, and $\rm{Visual}_{\rm{Hubble}}$ Classes in Fig. 3. The bijective agreement (which combines the fractional occupation of every class along both the row and the column) between the ProFuse classes and each of the visual classifications is shown in Fig. 3. The darker the matrix element is coloured (as indicated by the colourbar), the better the agreement. For each matrix class combination, we plot the $\log(M_{*}/\rm{M}_{\odot})$ – $\log(n)$ distribution of the subset, coloured by the $g-i$ colour (the relevant subset of the left-hand panel). For each matrix entry, a dashed box indicates regions that we would reasonably expect to be populated if the classification schemes were working perfectly as intended. The fact that the majority of these regions are highly populated generally suggests that the classification schemes are working fairly consistently.

The lower panel of Fig. 3 shows the equivalent matrix comparison, where the individual agreement percentages are shown in each category (showing greater detail than simply the coloured bijective agreement).

These figures are designed to provide insight into the potential differences in the classification schemes, and include an enormous amount of information. Therefore, an exhaustive analysis of each intersection in these comparison matrices is not the aim of this section. Rather, the below subsections summarise some of the major conclusions apparent from this complex comparison.

An example of full ProFuse outputs can be seen in figure 8 of Robotham et al. (2022), which demonstrates how the SED and 2D fitting elements of this technique come together to produce an accurate model with a wealth of data for an individual galaxy and its structural components.

As a brief demonstration of the galaxy fits achieved our ProFuse implementation in this work, we show three-band ($g$/$r$/$Z$) thumbnails of the best-fitting model and corresponding residual for a handful of sources in Figs. 7-8 (with a further selection presented in Appendix A). Residuals shown are the absolute residual, highlighting both the negative and positive features. For each galaxy, the image, model and residual thumbnails are identically scaled, to make all visible features directly comparable.

Fig. 7 shows a handful of examples for which the ProFuse model results in little residual flux, largely due to the lack of visual substructures within the galaxies (seen in the right-most column of the figure). Because these galaxies are smaller on the sky, the lower effective resolution likely smooths over any substructure in the galaxies, making them easier to model.

Structure types that have not been explicitly modelled in this work include components such as bars and spiral arms. Galaxies that have very strong bar features are shown in Fig. 8. Here, it is clear in the right column that there is a notable portion of the galaxy light that has not been well modelled by ProFuse. In particular, the bar is clearly visible in the residual images, and the red colour is dramatically different from the blue spiral arms, hinting at very distinct stellar populations within these structures. While we find that the total SED of the galaxy is still well approximated on average by ProFuse (and hence the global properties are not likely to be dramatically affected), the 2D model is non-optimal. It has been shown by Erwin et al. (2021) that a classical bulge plus disk decomposition of a barred galaxy can overestimate the spheroid components by factor of 4-100. Using the SDSS (Sloan Digital Sky Survey, Stoughton et al., 2002), Nair & Abraham (2010) estimate that 26% of galaxies with $g<16$ have bars, which is higher than the results that come from a visual inspection of our colour residuals, in which we identify bars in 9% of our sample in the same magnitude range. The resulting impact of neglecting bars in our fitting therefore has the potential to affect a significant fraction of our sample, although the light fraction of a bar in a galaxy can vary dramatically, so the true impact is difficult to estimate. Interestingly, while the residual structure for the galaxies in Fig. 8 seems to be similar, a mix of models (BD, FS and PD) were selected for the galaxies themselves. This highlights that the presence of a bar does not tend to influence the automatic classifications from ProFuse, as bars are badly modelled by all model configurations in general.

The residual features have a very strong colour signature. This indicates that an SED analysis of the residual features with ProSpect will already be possible with existing ProFuse outputs, allowing the stellar populations of features such as bars and spiral arms to be studied. For regular structures such as bars, it may be possible to include this as extra modelled structures in ProFuse itself, although care would need to be taken that there is sufficient quality in the imaging to constrain these extra free parameters. While it is beyond the scope of this work to conduct such an analysis, we feel it is a natural future extension. Additional samples of galaxies are presented in the same manner in Appendix A, to demonstrate some of the interesting galaxy phenomena that are identifiable from this modelling method.

The stellar mass distributions of the full sample, as well as the disk, bulge, and spheroid populations are shown in Fig. 9, for each of the three classification schemes used in this work. While the sample itself is identical in each of the schemes, the overall stellar mass distribution displays minor variation, thanks to the different combinations of ProFuse models that were selected to best describe each galaxy in each scheme. As mentioned in Sec. 4.1, either PD or BD configurations could be assigned to galaxies that were visually deemed to be two-component systems. Rather than attempt assign a better model numerically to each individual galaxy, we conduct the analysis with both, which gives the most pessimistic estimate of the error. The uncertainty ranges in the bulge and two-component disk components in the middle and bottom panels of Fig. 9 are the extremes bounded by the scenarios where either all two-component systems are described by PD, or all two-component systems are described by BD. The use of PD generally results in a lower-mass bulge (because of the limited size that a bulge can have), and conversely the use of BD can produce higher-mass bulges. This is what is responsible for the substantial uncertainty in the bulge mass distribution at the low- and high-mass extremes. From the automatic ProFuse classes, it seems that lower-mass bulges are generally best described by the PD configuration, whereas high-mass bulges are best described by BD.

There are significant differences that demonstrate the fundamentally different approaches used to classify galaxies. The double-peaked distribution of disks separates relatively cleanly between the pure disk and the two-component disk systems. This is true for all three classification schemes, although to different degrees. The $\rm{Visual}_{\rm{DR4}}$ classes produce the highest low-mass peak in two-component disks, whereas the $\rm{Visual}_{\rm{SB}}$ classes have the lowest number of low-mass two-component disks (with the ProFuse classes somewhere in the middle). This behaviour is mirrored in the bulge population, with the $\rm{Visual}_{\rm{DR4}}$ classes having the most low-mass bulges, the $\rm{Visual}_{\rm{SB}}$ classes the fewest, with the ProFuse classes somewhere in the middle. In general, it seems that the ProFuse classes produce a result that is the compromise between the different visual schemes at play in the $\rm{Visual}_{\rm{DR4}}$ and $\rm{Visual}_{\rm{SB}}$ classes.

To compare the overall mass distribution of our sample with previous work, we present the total mass function in Fig. 10. This is compared against the recent mass function from Driver et al. (2022a), and the previous GAMA mass function from Wright et al. (2017). We see that our mass function appears to be shifted slightly towards higher stellar masses than that of both Wright et al. (2017) and Driver et al. (2022a). This is to be expected, given that the ProFuse-derived stellar masses are on average 0.1 dex larger than that of ProSpect (which is where the stellar masses originated from for the Driver et al. 2022a analysis). This offset is shown in fig. 11 of Robotham et al. (2022), with more detailed discussion than we repeat here.

The mass function we derive when dividing by structure types is presented in Fig. 11, where the indicated uncertainty ranges convey the error produced through different classification approaches (compressing the information conveyed in the three panels of Fig. 9). Here, a completeness correction was conducted on a per-galaxy basis using the 1/$V_{\rm{max}}$ values for each galaxy computed using ProSpect. For ease of interpretation, the top panel of Fig. 11 focuses on disk components, whereas the bottom panel focuses on the bulge and spheroid components. We include the Schechter fits to equivalent structure mass functions from Moffett et al. (2016b) in coloured lines. The total disk mass functions are generally consistent, as is the two-component disk distribution. Equivalently, the spheroid distribution is also entirely consistent. Some differences exist when comparing the pure disk systems, where we note that our ProFuse pure disk mass function extends to much higher stellar masses than that of Moffett et al. (2016b). Similarly, our mass function for bulges is much less peaked, and contains a higher density of low-mass bulges than noted by Moffett et al. (2016b).

We present the contribution of bulges, disks, and spheroids to the overall mass budget as a function of stellar mass in Fig. 12, alongside the corresponding curves from Moffett et al. (2016b). It is remarkable that our mass fractions agree well with those measured in the past using Hubble-type classifications alone. At high stellar masses ($>10^{10}\rm{M}_{\odot}$), our values are indeed entirely consistent. It is just at lower stellar masses where we tend to assign more mass to both bulges and spheroids, and less mass to disks, than Moffett et al. (2016b). Note that LBS objects are not explicitly considered in the curves of Moffett et al. (2016b) (potentially accounting for the difference in the spheroidal populations).

Using our volume-limited sample, we derive the star formation rate density (SFRD) with the $z<0.06$ volume to representing the Universe. The SFRD is derived by summing the star formation histories of all the galaxies in the sample, normalised by the volume. To account for incompleteness at the lowest stellar masses (below $10^{9}\,\rm{M}_{\odot}$), we scale the contribution of each galaxy to the SFRD by $V/V_{\rm{max}}$,³³3Often referred to simply as a $1/V_{\rm{max}}$ correction to account for the portion of the studied volume over which the galaxy is undetectable within the selection limits. Here, $V$ is the volume of the obervational sample, and $V_{\rm{max}}$ is the lower between the maximum volume within which the object is observable, or the full volume of the sample.

The resulting total SFRD of our sample derived using our ProFuse method is shown in Fig. 13 as the thick black line. The standard compilation of data points⁴⁴4Measurements included in this compilation come from the following studies: Sanders (2003); Takeuchi et al. (2003); Wyder et al. (2005); Schiminovich et al. (2005); Dahlen et al. (2007); Reddy & Steidel (2009); Robotham & Driver (2011); Magnelli et al. (2011); Magnelli et al. (2013); Cucciati et al. (2012); Bouwens et al. (2012a, b); Schenker et al. (2013); Gruppioni et al. (2013) by Madau & Dickinson (2014) is shown as faint grey data points, as well as the more recent observational measurements of the CSFH by D’Silva et al. (2023).

The forensically derived CSFH using only broadband SED fitting via ProSpect by Bellstedt et al. (2020b)⁵⁵5We present here the SFRD presented in appendix B of Bellstedt et al. (2020b) rather than the main body of the paper, to compare results using the equivalent metallicity evolution prescription. is shown in Fig. 13 as a dashed magenta line. While both this derivation and the in-situ measurement by D’Silva et al. (2023) uses the modelling power of ProSpect, the implemented methods of deriving the SFRD are entirely different. Although the same parametrisation for the SFH and metallicity is used within both ProSpect and ProFuse, the fact that the SFHs can effectively be multi-modal using ProFuse means that these CSFH derivations were by no means destined to look the same. Despite this, the results do agree for the most part, except for slight differences in the early Universe and small discrepancies at low redshift.

This lowest redshift range has been shown in the Fig. 13 inset, highlighted in orange, so that discrepancies are clearer to see. At a lookback time of 1 Gyr, the ProSpect SFRD is larger than that derived by ProFuse by 13% (although the two are still entirely consistent within the uncertainty range introduced by different classification schemes). A possible explanation for the recent dip in the ProFuse CSFH is that star-forming clumps in galaxies (such as in their spiral arms) are not explicitly fitted. This is evident in the residuals shown in Fig. 8. If recent star formation is systematically under-fitted, because the structure is clumpy and not well-described by our model configurations, then we may be systematically underestimating the recent SFR in some cases. An additional contribution to the potential difference at recent times is the manner in which a completeness correction is conducted. A more simplistic correction was conducted by Bellstedt et al. (2020b), which may account for some differences in the lowest lookback time bins of the SFRD.

The other difference between the total CSFH in this work and that derived using ProSpect by Bellstedt et al. (2020b) can be seen at lookback times of greater than 6 Gyr, where the added flexibility of multiple components in this work has resulted in slightly more star formation being recovered (although considering the uncertainties, the discrepancy is only small). Through a comparison of SED-fitting outputs from GAMA, galaxies from the semi-analytic model Shark (Lagos et al., 2018a), and SED-fitted outputs of the Shark galaxies, Bravo et al. (2022) conducted a very thorough analysis to determine the reliability of forensic properties derived from ProSpect. They determined that the evolution of galaxy properties could be accurately reproduced to lookback times of $\sim$6 Gyr, but that beyond this epoch, the results start to be dominated by the modelling assumptions. Given that the light from the oldest stellar populations is much fainter than from younger ones, it is unsurprising that the constraint from SED-fitting is much lower here. For this reason, the slight deviation we see in the CSFH between the ProSpect-only and ProFuse derivations may well simply be the consequence of low constraints from these stellar populations. While the absolute SFHs at these older epochs likely have greater uncertainty from SED fitting that limit our derivation of absolute ages, the relative differences between histories for individual populations is still meaningful.

The measurements by D’Silva et al. (2023) are also derived using ProSpect with a linear metallicity evolution model, consistent with the SED modelling approach in this paper. D’Silva et al. (2023) also include AGN in their modelling as described by Thorne et al. (2022a), resulting in an improved estimation of galaxy properties, especially at high-$z$. Our forensic derivation is broadly consistent, however we note that our forensic SFRD is perhaps underestimated in the 8-10 Gyr range, and then overestimated beyond 10 Gyr. The agreement between each of the ProSpect-derived approaches in measuring the SFRD is remarkable within the most recent 8 Gyr of lookback time.

The SFRD with contributions by different components is presented in Fig. 14, with the total SFRD shown in black, and in colour the contributions toward the total SFRD by different galaxy structures, including disks (blue), bulges (red), and spheroids (orange). Values plotted in Fig. 14 are presented in Table 5, with full data available in the supplementary material.

The corresponding evolution in the stellar mass density inferred by our star formation history is shown in Fig. 15. The impact of the slight difference in total SFRD between ProSpect (from Bellstedt et al. 2020b) and ProFuse is a 16% factor increase in the total stellar mass density. Note that despite this increase, the present-day stellar mass density is still entirely consistent with the compilation data⁶⁶6This compilation includes SMD measurements by Arnouts et al. (2007); Gallazzi et al. (2008); Pérez-González et al. (2008); Kajisawa et al. (2009); Li & White (2009); Marchesini et al. (2009); Yabe et al. (2009); Pozzetti et al. (2010); Caputi et al. (2011); González et al. (2011); Bielby et al. (2012); Lee et al. (2012); Reddy et al. (2012); Ilbert et al. (2013); Labbé et al. (2013); Moustakas et al. (2013); Muzzin et al. (2013) of Madau & Dickinson (2014). Values plotted in Fig. 15 are presented in Table 6, with full data available in the supplementary material.

We discuss the results for each component individually in the following subsections, reiterating that all trends refer to the stars that are within each component at the present day, as opposed to the component that the stars may have been in at previous epochs.

Despite the substantial effort that has been invested in recent years to better constrain the CSFH prior to cosmic noon, there is still significant uncertainty. This is driven by the inherent difficulty in measuring SFRs of galaxies in this epoch, due to uncertainties in the substantial dust corrections required at this epoch (examples of the numerous studies that have attempted to make these measurements using both rest-frame UV and FIR include Kistler et al., 2009; Gruppioni et al., 2013; Bouwens et al., 2015; Rowan-Robinson et al., 2016; Novak et al., 2017; Khusanova et al., 2020, 2021).

The analysis presented in Fig. 13 suggests that the capacity for codes like ProSpect to forensically measure the SFRD is excellent over the last $\sim$8 Gyrs. We therefore speculate that in future, one of the most promising methods used to constrain the $z>3$ SFRD will be to conduct a forensic analysis of galaxies at $z\sim 0.8$. This is realistically the greatest distance at which forensic-style CSFH extractions are possible, as the stellar mass of the completeness limits becomes too high at greater lookback times, and morphologies become too irregular to model the galaxies well. This is also the epoch at which the uncertainties of the derived stellar mass function increase dramatically in the analysis by Wright et al. (2018). With this forensic approach at higher redshifts, there should then be enough constraining power within high-quality data to access cosmic dawn. The dataset that will make this possible is the complete spectroscopic redshift survey WAVES-deep (Driver et al., 2019) to be conducted on 4MOST (de Jong et al., 2019), and other redshift surveys conducted with facilities like MOONS (Cirasuolo et al., 2020), in combination with high-quality imaging from facilities such as Euclid (Laureijs et al., 2012).

The overall trends presented for bulges, spheroids and disks are qualitatively similar to the parametric model presented by Driver et al. (2013), where star formation in present-day spheroids (referring in that work to all ellipsoidals) was deemed responsible for the bulk of star formation at and before cosmic noon, followed by an increase in star formation in disks. The subtle differences we identify are more prolonged star formation in spheroids, with a slight increase in the fraction of earlier formed stars that end up in present-day disks.

Comparing directly with bulge/disk star formation histories from IFU-based work is difficult, because of the substantially smaller (and incomplete) samples currently studied in that level of detail. The main takeaway from an analysis of lenticular galaxies by Johnston et al. (2022) is that bulges in general seem to be pretty old, but that the spread in ages for disks is much broader (and on average younger than bulges). This is consistent with the picture we see from our complete sample of $\sim$7,000 galaxies. Even for individual subsets of galaxies however, a quantitative comparison is hampered by the fundamental differences derived in the star formation histories. For example the difficulty in distinguishing very old populations means that the star formation histories by Johnston et al. (2022) begin at epochs earlier than the start of the Universe. Different stellar populations analysis tools, whether SED fitting tools like ProSpect or spectral fitting tools like ppxf (Cappellari & Emsellem, 2004), will all suffer from similar difficulties in constraining older populations. While there is more information in spectra than broadband photometry, there is also therefore more modelling complexity, and sensitivity to modelling degeneracies and assumptions (for example, ppxf is very sensitive to implementation treatments such as the use of regularization). Given the different approaches used to overcome these challenges, the quantitative histories are not yet directly comparable.

Various observational studies have identified “young” bulges, the existence of which seems to be at odds with our results at first glance. Such young bulge-like structures have been seen at the low-mass end (such as the star-forming nuclear star clusters identified by Johnston et al., 2020), which would not impact the overall bulge CSFH due to the small amount of stellar mass involved. Young star-forming bulges were identified by Barsanti et al. (2021a) in an analysis of cluster galaxies. Such galaxies with bulges younger than disks do exist in our sample (as shown in fig. 21 of Robotham et al., 2022, and discussed in Sec. 5.3.2), however as this only accounts for $\sim$12% of the two-component sample, the impact is insufficient to have an overall effect on the mass-weighted bulge CSFH. Nonetheless in Fig. 14 we separate the contribution of these younger bulges, and show that it is distinct to that of the rest of the bulge population.

Recent work by Costantin et al. (2021) applied SED fitting to multi-wavelength photometric decompositions of 156 high-$z$ galaxies to study the stellar populations of bulges, disks and spheroids. In their work, a bimodality was identified in the formation epoch of bulges, which contributed to the interpretation by Costantin et al. (2022) of two distinct phases of bulge formation in cosmic time. The bulge SFRD we present in Fig. 14 demonstrates no hint of two distinct phases. Instead, it suggests that bulges are simply a continuously old population, with only a small tail of star formation at recent times (contributed to by the subset of bulges that are younger than their disks). Even when analysing the age distribution of bulges identified in our work (as presented in detail by Robotham et al., 2022), we do not see any clear bimodality. The different redshift epochs studied hamper a direct comparison though, as our forensic approach may not be sensitive to two separate epochs that are both quite old. If so, this suggests that any potential high-$z$ bimodality is not reflected in the properties of present-day bulges. Furthermore, the two studies use different SFH parametrisations (Costantin et al. 2021 use a declining delayed exponential parametrisation), which can also contribute to differences in recovered ages.

As demonstrated by the need to clarify our nomenclature in Sec. 3.2, the physical interpretation of the “ellipsoidal” structures in the Universe (bulges and spheroids) varies across the literature. Are they in fact distinguishable structures with distinct evolutionary pathways (true “eigenstructures”), or are their structural and stellar population properties similar enough to warrant interpretation as a single entity?

Our analysis with ProFuse consistently suggests that spheroids and bulges are distinct populations, whether that be due to distinct size–mass relations as shown by Robotham et al. (2022), the stellar mass distributions, or the overall age of the populations as shown by their distinct contributions to the CSFH⁸⁸8Despite taking care to indicate these CSFH resulting from a range of classification schemes, any contamination of disk-like sources into the spheroidal population could still be influencing the distinction between the bulge and spheroid CSFH. in Fig. 14. This conclusion is similar to that of Gadotti (2009), who finds that high-mass bulges and ellipticals are offset in the size–mass relation, also concluding that they are in fact separate populations⁹⁹9Despite the similar conclusion, the size–mass relation for bulges measured by Robotham et al. (2022) is different to Gadotti (2009).. It is difficult to make clear comparisons, as the spheroidal population that we analyse is a broader population than ellipticals studies by Gadotti. This impact of this is seen when comparing (for example) the spheroid CSFH in Fig. 14 with that of the elliptical CSFH presented in fig. 8 of Bellstedt et al. (2020b), where the elliptical contribution to the CSFH dominates much earlier than the spheroids studies in this work generally. Our spheroids contain a larger portion of younger galaxies with more active star formation, typically at lower stellar masses. In the past, such galaxies have been separated out into distinct classes, such as Kormendy & Bender (2012) who treat elliptical and dwarf elliptical (dE) galaxies as separate classes with distinct luminosity functions (see their fig. 3). This was also presented in the luminosity functions of Binggeli et al. (1988), where elliptical galaxies at the high-luminosity end are held distinct from dE and irregular (Irr) galaxies at the low-luminosity end.

Contrastingly, studies of SMBH-bulge correlations treat ellipticals and bulges not as separate populations, but as a continuous population described collectively as “spheroids” (for example Savorgnan et al., 2016; Savorgnan & Graham, 2016). Ellipsoidals were also treated as a common-origin structure in the simple model by Driver et al. (2013) (also referred to collectively as simply spheroids in that work).

This complexity of interpretation is not clearer in simulations. In analysing the IllustrisTNG simulation, Tacchella et al. (2019) study spheroids and bulges as one continuous class. This is contrasted by Ortega-Martinez et al. (2022), who separate classical bulges and spheroids in their EAGLE-based analysis. Despite separating these populations into categories that seem more consistent with our observational approach, however, the spheroids in their analysis tend to be much more extended, and with a higher Sérsic index, suggesting that this divide is still not quite the same. In semi-analytic models like Galform (Cole et al., 2000) and Shark (Lagos et al., 2018b), pure spheroids aren’t explicitly formed. Instead, bulge-like structures are formed via two main mechanisms: mergers or disk instabilities. For this reason, ellipsoidals in these models are usually divided into disk-instability and merger-origin bulges, which are again different to any other definition of this structure. Yet the labels assigned to these structures are often “pseudobulges” and “classical spheroids” respectively (as applied by Huško et al., 2023), which again makes direct comparisons between studies difficult.

The analysis from this paper, and how it relates to the literature, suggests that finding comparable definitions for disks is relatively straightforward, but that bulges/spheroids are much less clearly defined at this stage. This highlights that our definitions are far from universal, and it is essential that more care is taken in future to characterise the specific elements of the “bulge complex”, to facilitate more direct comparisons.

Our forensic CSFH suggests that around 60% of all stars formed at $z\sim 3$ end up in present-day spheroidal systems, with the other 40% having ended up pretty evenly in both disks and bulges. Because our analysis is a forensic one, it simply reflects which structures stars end up in at the present day, not what type of structures they were formed in. Early results from JWST indicate that substantial morphological evolution has occurred since cosmic noon, with morphological classification of 850 galaxies by Kartaltepe et al. (2023) suggesting that $>$56% of galaxies at $z>3$ have a visually identifiable disk, with only 38% of galaxies visually identified as spheroids. To really quantify this evolution however, complete studies at high redshift are needed to estimate the mass fraction of each structure, not just the fraction of galaxies that include these structures.

With the quality of imaging now available from JWST, it is possible to conduct bulge/disk decomposition (as Chen et al., 2022, did for a sample of 7 galaxies, demonstrating that a bulge component was already in place for each of them). The contrast between high-$z$ measurements and our forensic view of the early Universe clearly indicates that some morphological change must be occurring. The exact mechanism by which this tranformation is occurring is hard to tell from our forensic analysis.

We use the outputs from the semi-analytic model Shark (Lagos et al., 2018a) to further explore the differences between “true” mass fractions of bulges, disks and spheroids with cosmic time, as compared with forensically-derived mass fractions like the ones measured in this work. Seeing how these measurements differ provides some intuition for how better to interpret forensic analyses. This analysis is presented in Fig. 16, where our results from the bottom panel of Fig. 15 are renormalised to show the change in mass fraction for each structural component, relative to the $z=0$ epoch. The absolute fractions of mass in disks and bulges is quite different in Shark to that seen in GAMA, which is quite possibly caused by the tuning of the semi-analytic model. For this reason we do not focus on the absolute fractions, but on the relative build-up of disk and bulge mass, which is more fundamentally linked to the modelled physics. The ProFuse results and the equivalent forensic view from Shark shown as solid lines. The “true” view from Shark is shown as dashed lines.

What the comparison of the true and forensic mass fractions from Shark shows, is that the historical mass in ellipsoidals is likely overestimated with a forensic approach, whereas the historical disk mass is likely underestimated. In the middle panel of Fig. 16, the dark solid line shows the mass fraction in bulges assumed if Shark galaxies were to be analysed forensically, where the age distributions of stars in the present-day structures are analysed (as we have done observationally using ProFuse). This implies that the fraction of mass in bulges was higher at higher lookback times than in the present day (consistent with our forensic measurement). When the relative mass fraction is extracted from Shark by actually measuring this value in-situ at different epochs (shown by the dashed line), it can be seen that the fraction of stars held in bulges is actually relatively constant. This indicates that some of the old stars in present-day bulges must not always have been in a bulge in the past. Given that the reverse trend is shown in the disk component in the bottom panel of Fig. 16, this suggests that it is possible that the old stars forensically attributed to bulges, may have actually come from disks. This points to a morphological evolution pathway, in which stars that are formed in disks likely end up in ellipsoidal structures. This is consistent with the picture of morphological transformation built up based on both observations (for example Dressler, 1980; Hashemizadeh et al., 2022; quilley2022; Marasco et al., 2023) and simulations (for example Bournaud et al., 2007; Clauwens et al., 2018; Martin et al., 2018; Jackson et al., 2019).

We caution that in different semi-analytic models, the “true” picture can differ dramatically to that shown in Shark. For example, Huško et al. (2023) uses Galform to present the fraction of mass in bulges with cosmic time, and these curves differ substantially. In particular, a much lower fraction of mass exists within bulges at the present day, but rises substantially beyond $z=2$.