<figure><img src="image_1.png" /></figure>

# A Duality in the Origin of Bulges and Spheroidal Galaxies

Luca Costantin¹, Pablo G. Pérez-González¹, Jairo Méndez-Abreu²,³,⁴,⁵, Marc Huertas-Company²,³,⁶, Paola Dimauro⁷, Belén Alcalde-Pampliega⁸, Fernando Buitrago⁹, Daniel Ceverino¹⁰,¹¹, Emanuele Daddi¹², Helena Domínguez-Sánchez¹³,¹⁴, Néstor Espino-Briones¹⁵, Antonio Hernán-Caballero¹⁶, Anton M. Koekemoer¹⁷, and Giulia Rodighiero¹⁸

¹Centro de Astrobiología (CSIC-INTA), Ctra de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain; lcostantin@cab.inta-csic.es

²Instituto de Astrofísica de Canarias, E-38200, La Laguna, Tenerife, Spain

³Departamento de Astrofísica, Universidad de La Laguna, E-38205, La Laguna, Tenerife, Spain

⁴Dpto. de Física y del Cosmos, Campus de Fuentenueva, Edificio Mecenas, Universidad de Granada, E-18071, Granada, Spain

⁵Instituto Carlos I de Física Teórica y Computacional, Facultad de Ciencias, E-18071, Granada, Spain

⁶LERMA, Observatoire de Paris, CNRS, PSL, Université de Paris, France

⁷Observatório Nacional, Rua General José Cristino, 77, São Cristóvão, 20921-400, Rio de Janeiro, Brazil

⁸European Southern Observatory (ESO), Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile

⁹Instituto de Astrofísica e Ciências do Espaço, Universidade de Lisboa, Tapada da Ajuda, 1349-018, Lisbon, Portugal

¹⁰Universidad Autonoma de Madrid, Ciudad Universitaria de Cantoblanco, E-28049, Madrid, Spain

¹¹CIAFF, Facultad de Ciencias, Universidad Autonoma de Madrid, E-28049 Madrid, Spain

¹²AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191, Gif-sur-Yvette, France

¹³Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Magrans, E-08193, Barcelona, Spain

¹⁴Institut d'Estudis Espacials de Catalunya (IEEC), E-08034, Barcelona, Spain

¹⁵Departamento de Física de la Tierra y Astrofísica, Faultad de CC Físicas, Universidad Complutense de Madrid, E-28040, Madrid, Spain

¹⁶Centro de Estudios de Física del Cosmos de Aragón, Plaza San Juan 1, planta 2, E-44001, Teruel, Spain

¹⁷Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁸Department of Physics and Astronomy, University of Padova, Vicolo Osservatorio 3, I-35122, Padova, Italy

Received 2020 December 17; revised 2021 February 26; accepted 2021 March 15; published 2021 June 2

## Abstract

Studying the resolved stellar populations of the different structural components that build massive galaxies directly unveils their assembly history. We aim at characterizing the stellar population properties of a representative sample of bulges and pure spheroids in massive galaxies ($M_{*} > 10^{10} M_{\odot}$) in the GOODS-N field. We take advantage of the spectral and spatial information provided by SHARDS and Hubble Space Telescope data to perform the multi-image spectrophotometric decoupling of the galaxy light. We derive the spectral energy distribution separately for bulges and disks in the redshift range $0.14 < z \le 1$ with spectral resolution $R \sim 50$. Analyzing these spectral energy distributions, we find evidence of a bimodal distribution of bulge formation redshifts. We find that 33% of them present old mass-weighted ages, implying a median formation redshift $z_{form} = 6.2^{+1.5}_{-1.7}$. They are relics of the early universe embedded in disk galaxies. A second wave, dominant in number, accounts for bulges formed at median redshift $z_{form} = 1.3^{+0.6}_{-0.6}$. The oldest (first-wave) bulges are more compact than the youngest. Virtually all pure spheroids (i.e., those without any disk) are coetaneous with the second-wave bulges, presenting a median redshift of formation $z_{form} = 1.1^{+0.3}_{-0.3}$. The two waves of bulge formation are distinguishable not only in terms of stellar ages but also in star formation mode. All first-wave bulges formed fast at $z \sim 6$, with typical timescales around 200 Myr. A significant fraction of the second-wave bulges assembled more slowly, with star formation timescales as long as 1 Gyr. The results of this work suggest that the centers of massive disk-like galaxies actually harbor the oldest spheroids formed in the universe.

*Unified Astronomy Thesaurus concepts: Galaxy bulges (578); Galaxy structure (622); Galaxy disks (589); Galaxy spheroids (2032); Galaxy formation (595); Galaxy evolution (594); Galaxy quenching (2040); Galaxy stellar content (621); Galaxy properties (615); Galaxy ages (576)*

*Supporting material: machine-readable table*

## 1. Introduction

Integrated or spatially resolved observations are typically used to infer the stellar content of galaxies at low redshift, reconstructing their star formation history (SFH) with an “archaeological” approach (Thomas et al. 2005; Rogers et al. 2010; González Delgado et al. 2015). Complementarily, to trace back the formation of stars in the universe, it is possible to compare the stellar content of similar samples of galaxies at different redshifts (the so-called “look-back” approach; Schiavon et al. 2006; Sánchez-Blázquez et al. 2009; Gallazzi et al. 2014). In these studies, stellar population models with different SFHs are usually compared to the best model either obtained from the fit of the spectral energy distribution (SED) or derived

using key spectral features that are sensitive to fundamental physical parameters such as age, metallicity, or $\alpha$ enhancement (Kriek et al. 2011; Domínguez Sánchez et al. 2016). Both methods converge to a common picture, which shows that the most massive galaxies form at earlier epochs and in short timescales in the so-called “downsizing” scenario (Cowie et al. 1996; Pérez-González et al. 2008; Thomas et al. 2010; McDermid et al. 2015).

One fundamental piece of evidence is the size evolution of galaxies. Massive quenched galaxies at high redshift ($z \sim 2$) are up to 6 times smaller in size than their counterpart of the same mass today (Daddi et al. 2005; Trujillo et al. 2006b, 2006a; Buitrago et al. 2008; van der Wel et al. 2014). Similarly, at all
fixed galaxy masses, the more extended galaxies are younger, more metal-poor, and less α-enhanced compared to the more compact galaxies (Scott et al. 2017). Morphologically, at $z \sim 2$ there is a bimodal distribution of galaxies: the universe is composed of quiescent compact pure spheroids with no clear evidence of a disk component and star-forming clumpy disks, which gradually transform from disturbed to normal disk galaxies by $z \sim 1$. By contrast, at lower redshift massive galaxies ($M > 10^{10.8} M_\odot$) present almost no morphological evolution. The abundance of bulgeless galaxies decreases with redshift until almost the entire population presents a significant bulge component (Huertas-Company et al. 2016).

Indeed, this global view misses a crucial piece of information: galaxies are complex systems, which generally consist of multiple morphological components (i.e., bulges, disks, bars, etc.). In the simplest scenario, disk galaxies are composed by a central bulge and an outer disk. In particular, according to the photometric definition, the bulge is identified as the light excess in the central part of a galaxy, over and above the light associated with the exponential profile of the external stellar disk. In this work we not only focus on pure spheroids but also study the stellar populations of bulges up to redshift $z=1$ in the context of their assembly history.

Spheroids are traditionally thought to arise from either a violent and dissipative collapse of protogalaxies (Eggen et al. 1962; Larson 1976) or accumulation and rearrangement of stars in merger events (Cole et al. 2000; Hopkins et al. 2009). Recently, an alternative scenario emerged, since cosmological simulations point toward a rapid spheroid formation while high-redshift galaxies go through a gas-compaction phase (Dekel & Burkert 2014; Zolotov et al. 2015; Tacchella et al. 2016). Considering the high density and high gas fraction of the universe at high redshifts, multiple linked mechanisms (i.e., minor mergers, violent disk instabilities, clump migration, counterrotating streams, etc.) seem to conspire to quickly fuel the gas into the central region of the disky galaxy (Mihos & Hernquist 1996; Noguchi 1999; Immeli et al. 2004a, 2004b; Hopkins et al. 2006; Dekel et al. 2009; Ceverino et al. 2012; Forbes et al. 2014; Jog 2014; Renaud et al. 2014; Danovich et al. 2015; Dyda et al. 2015; Wellons et al. 2015; Bournaud 2016). Torques in the inner regions facilitate gas inflow and inward clump migration, compacting the stellar mass in the form of a “blue nugget” (Barro et al. 2013, 2014), growing a massive rotating stellar spheroid (Genzel et al. 2008; Ceverino et al. 2010) or a classical bulge (Ceverino et al. 2015), and feeding the central black hole (Bournaud et al. 2011).

In this first paper of a series, we provide robust estimations of some key properties of bulges and pure spheroids (i.e., spheroids without an extended stellar disk), such as the stellar mass and the mass-weighted age, discussing the implications for their assembly history across cosmic time. Key questions still need to be addressed: (1) When does the spheroidal population assemble? (2) Why are some dark matter halos evolving more rapidly and forming massive spheroids at earlier cosmic time? (3) Do spheroids present different properties at different cosmic epochs? To date, answering these questions for the separate galaxy components remains very elusive, since spatially resolved studies at high redshift are very scarce because very limited data (mainly broadband only) can be gathered with spectral information.

We investigate the mechanisms that drive the evolution of spheroids (both bulges and pure spheroids) taking advantage of

the Survey for High-z Absorption Red and Dead Sources (SHARDS; Pérez-González et al. 2013), an ESO/GTC Large Program that provided an ultradeep ($m < 26.5$ AB mag) imaging survey in 25 filters covering the wavelength range of 500–950 nm. The state-of-the-art photometric multifilter surveys, like SHARDS, allow us to easily study large samples of galaxies at different redshifts with observations typically deeper than spectroscopic ones. Indeed, the “snapshots” across cosmic time of spheroid properties offer a more continuous view of their evolution. Thanks to its 25 medium-band filters, the SHARDS data set permits us to accurately determine the main properties of the stellar populations of galaxies, providing a smooth SED with a resolution $R \sim 50$. The observed photometry provides a pseudospectrum at each pixel on the sky, which allows for a 2D spatial analysis treating each filter independently. Moreover, SHARDS photometry offers accurate estimates of absorption and emission features (i.e., $M_{g,UV}$, D4000, [O II], and [O III]) that are visible within the SHARDS wavelength range at high redshift. SHARDS data are significantly deeper than spectroscopic surveys and grant a consistent improvement in the spectral resolution with respect to broadband studies, mitigating the typical degeneracies that affect the inference of individual SFHs of galaxies (Hernández-Caballero et al. 2013; Domínguez Sánchez et al. 2016).

We propose a novel approach based on the multiband bulge/disk spectrophotometric decoupling of the galaxy light to derive the SED separately for both the bulge and disk of each galaxy. This methodology is based on recently developed techniques of spectrophotometric decomposition (Johnston et al. 2017; Tabor et al. 2017; Méndez-Abreu et al. 2019a) and takes advantage of the high spatial resolution images from the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) and Wide Field Camera 3 (WFC3), which are used as priors for the decoupling of the bulge and disk light in SHARDS images. Moreover, we complement the information provided by SHARDS and HST with that at longer wavelengths by means of the Canada-France-Hawaii Telescope WIRCam data. This allows us to constrain the stellar mass and the stellar populations of each individual component in the sample with unprecedented accuracy. Thus, the spectrophotometric analysis of SHARDS data represents a significant step forward in reaching a more exhaustive picture of galaxy formation and understanding the interplay between baryons and their dark matter hosts.

The paper is organized as follows. The sample of galaxies is characterized in Section 2. The spectrophotometric decoupling of the galactic components, the error analysis, and the retrieval of the stellar population of bulges and pure spheroids are presented in Section 3. We present and discuss our results in Sections 4 and 5, respectively. Our conclusions are summarized in Section 6. Throughout the paper we assume a flat cosmology with $\Omega_m = 0.3$, $\Omega_\lambda = 0.7$, and a Hubble constant $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2. Sample

We analyze galaxies in the North field of the Great Observatory Origins Deep Survey (GOODS-N), which presents photometric data provided both by HST and SHARDS. In particular, taking advantage of the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011), we use the seven filters for HST images (i.e., ACS F475W, F606W, F775W, and
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Properties of the sample galaxies; from left to right: galaxy redshift and stellar mass (Barro et al. 2019), and galaxy Sérsic index (van der Wel et al. 2012). Top panels: the gray histogram defines the mother sample of 478 galaxies in the SHARDS field of view, while the green histogram stands for the final sample of 156 spheroids, namely, 65 pure spheroids and 91 bulges (red histogram). Bottom panels: correspondent cumulative distributions.</figcaption></figure>

F850LP and WFC3 F105W, F125W, and F160W) and combine them with the 25 filters for SHARDS data in the optical wavelength range 500–941 nm (see Pérez-González et al. 2013; Barro et al. 2019, for all details).

Our mother sample consists of 478 massive ($M_{*} > 10^{10} M_{\odot}$) and luminous ($m_{F160W} < 21.5$ mag) galaxies at redshift $z \le 1$ in the 141 arcmin$^2$ area surveyed by SHARDS (Barro et al. 2019). In Figure 1 (gray histogram) we present their mass, redshift, and Sérsic index distributions (van der Wel et al. 2012; Barro et al. 2019). We cross-match our sample galaxies with the Dimauro et al. (2018) catalog, which provides photometric bulge-disk decompositions in this field of view, and the van der Wel et al. (2012) catalog, which provides photometric parameters for the galaxies. In particular, by means of an unsupervised feature learning (deep learning) technique, Dimauro et al. (2018) measured the best analytic model describing the surface brightness distribution of a given galaxy in the WFC3 F160W band: 62 are modeled with a single Sérsic component (i.e., pure spheroids with bulge-over-total luminosity ratio $B/T > 0.8$), 41 galaxies are classified as pure disks (i.e., modeled with a single exponential profile), while 228 are modeled with two components (i.e., Sérsic for inner bulge and exponential for outer disk). We discard 147 galaxies, which were poorly fitted in both Dimauro et al. (2018) and van der Wel et al. (2012). Galaxies with very low and thus unphysical values of their Sérsic index were removed from the sample.

This first paper of a series focuses on the spheroidal component, i.e., either bulges in disk galaxies or pure spheroids (not accompanied by any other structural component), limiting the sample selection to 290 galaxies. Moreover, since it is difficult to discriminate the nature of the central light prominence in nearly edge-on galaxies, and in order to avoid the contamination from boxy-peanut structures, we also restrict the sample of galaxies with two components to the 192 with

inclination $i < 70^{\circ}$. We further model the 1D surface brightness profile of galaxies with two components (see Section 3.1.1), classifying 42 galaxies as pure disks and 15 as pure spheroids. This criterion prevents unphysical or poor solutions derived from a blind photometric decomposition procedure due to parameter degeneracy. At this level, the sample is composed of 135 two-component galaxies and 77 pure spheroids. Finally, we visually classify these galaxies looking for signs of interaction and/or foreground/background contaminating objects, removing 46 of these galaxies from our sample.

The final sample is composed of 156 galaxies, 65 pure spheroids, and 91 galaxies with a bulge and disk component, as presented in Figure 1 (green and red histograms). Ten more galaxies were discarded from the 166 galaxies analyzed because the 2D modeling of their individual SEDs provides no constraints for their stellar populations (see Sections 3.1.1 and 3.2). The representativeness of the final sample is checked by means of a Kolmogorov–Smirnov test, assuring that each subsample does not introduce any substantial bias in the ($z, M_{*}$) parameter space ($p$-value$_{z} > 10\%$, $p$-value$_{M_{*}} > 30\%$). Moreover, it is worth nothing that our selection clearly discards galaxies modeled with low values of the Sérsic index ($n \lesssim 1.5$). This is mainly because we focus on the spheroidal component, partially due to our inclination threshold, but also because we checked the goodness of our selection for each individual galaxy, discarding unphysical or poorly constrained solutions.

### 3. Analysis

In this section we present the preprocessing of the SHARDS images, the spectrophotometric decoupling of the bulge and disk components combining HST and SHARDS information, and the analysis of statistical errors and degeneracies.
### 3.1. Spectrophotometric Decoupling: Bulge and Disk

For this work, we started with the v1.14.5 SHARDS images from Barro et al. (2019). The reduction in all the bands redder than 800 nm was repeated (v1.15.0) to include larger masks for extended objects and improve the fringing correction and sky subtraction for galaxy outskirts. As a result of this, all v1.15.0 images probing wavelengths redder than 800 nm improved their reliability at surface brightness levels fainter than 24.5 mag arcsec⁻² and down to approximately 26 mag arcsec⁻². This turns out to be, in fact, a very important feature for our analysis of galaxy structural components.

We perform a spectrophotometric decoupling of the different galactic structures (i.e., bulge and disk) in the sample of massive galaxies presented in Section 2 as a function of wavelength. Our method consists of a robust and accurate 2D photometric decomposition of the 25 medium-band SHARDS images, jointly with HST WFC3 and K-band data, covering the wavelength range between 400 and 2000 nm. This analysis mimics the strategy used by the C2D code (Méndez-Abreu et al. 2019a, 2019b), where high spatial resolution broadband images are used to robustly constrain spectroscopic data with lower spatial resolution. At the core of C2D, the GASP2D algorithm (Méndez-Abreu et al. 2008, 2014) performs the actual photometric decompositions.

The high flexibility of GASP2D allows us to represent the different galactic structures with a variety of analytical functions. We model the light of the bulge component with a Sérsic function (Sérsic 1968):

$$I_{\text{bulge}}(r_{\text{b}}) = I_{\text{e}} 10^{-b_{\text{n}} \left[ \left( \frac{r_{\text{b}}}{R_{\text{c}}} \right)^{1/n} - 1 \right]}, \quad (1)$$

where $r_b$ is the radius measured in the reference system of the bulge (see Costantin et al. 2017), $R_c$ is the effective radius, $I_c$ is the intensity at the effective radius, $n$ is the Sérsic index describing the curvature of the profile, and $b_n \simeq 0.868n - 0.142$ (Caon et al. 1993). The light of the disk component is parameterized by a single exponential (Freeman 1970):

$$I_{\text{disk}}(r_{\text{d}}) = I_0 e^{-r_{\text{d}}/h}, \quad (2)$$

where $r_d$ is the radius measured in the reference system of the disk (see Costantin et al. 2017), while $I_0$ and $h$ represent the central intensity and scale length of the disk, respectively. Both the bulge and the disk are assumed to have elliptical isophotes centered on the galaxy center ($x_0, y_0$), with constant position angle PA (counterclockwise measured from the north) and constant axial ratio $q$, defined as the minor axis divided by major axis.

Ideally, we would like to apply a standard 2D structural decomposition to each SHARDS image. This would imply fitting a total of 11 free parameters for each image: five for the bulge ($I_c, R_c, n, q_b, PA_b$), four for the disk ($I_0, h, q_d, PA_d$), and the galaxy center ($x_0, y_0$). Given that the spatial resolution of the seeing-limited SHARDS images is around 0.9, and considering that we are targeting galaxies at high redshift, the degeneracies involved in determining those 11 parameters independently for each SHARDS bands are quite significant. To account for this, we take advantage of the high spatial information of multiple optical and near-infrared HST images to derive the structural parameters of the bulge (i.e., $R_c, n, q_b, PA_b$) and disk (i.e., $h, q_d, PA_d$) in our sample galaxies. Then, we allow only the corresponding intensity to vary in the SHARDS

filters (see Méndez-Abreu et al. 2019a, 2019b, for a similar application using a combination of CALIFA and SDSS data).

### 3.1.1. Synergy between HST and SHARDS Data

The key strength of this work consists in combining HST and SHARDS data. Indeed, we take full advantage of the spectral information provided by SHARDS medium-band images and the spatial resolution given by HST. This allows us to retrieve individual SEDs of bulges and disks with spectral resolution $R \sim 50$ up to redshift $z=1$ and down to $\log(M_{\text{gal}}/M_{\odot}) = 10$. In Figure 2 we show an example of the different spatial information provided by SHARDS and HST for the galaxy GDN 18522, where it is possible to appreciate the differences of the bulge and disk component as a function of wavelength.

As a first step in our method, we retrieve initial guesses for the bulge and disk structural parameters from the analysis of HST data presented in Dimauro et al. (2018). We used these initial guesses as priors to characterize the bulge and disk parameter space by means of Bayesian inference (see Section 3.1.2, for a full description). The most probable values of the bulge and disk structural parameters (i.e., $R_c, n, h$) are linearly interpolated over wavelengths, transferring the HST information to the SHARDS images. Then, the galaxy image in each of the SHARDS and HST WFC3 filters is fitted using GASP2D. The intensity of the two components varies independently at each wavelength, while their structure is kept frozen (to the values obtained from the interpolation in the HST results). This is possible thanks to the exquisite astrometric calibration and distortion correction of the SHARDS images. In order to better constrain the stellar mass of each individual component, we extend the analysis to the K band using WIRCam data (Hsu et al. 2019). For the K-band images, we use as an initial guess the structural parameters determined for the WFC3 F160W band. Even if this is an extrapolation of the bulge and disk properties, we expect mild wavelength variation of the structural parameters for our intermediate-redshift galaxy sample.

As a result, we compute a pseudodata cube with all the spatial and spectral information for each component included in the fit, as well as a data cube with the galaxy model and the residuals at each wavelength. As an example, in Figure 3 we show the radial surface brightness profile, the 2D model, and the residual map of GDN 18522 in the WFC3 F160W band. For comparison purposes, in Figure 4 we present the radial surface brightness profiles obtained in three different SHARDS bands (see Appendix A for more examples). Dimauro et al. (2018) classified the galaxy as a two-component system (probability of 0.88) fitted with an exponential disk and a Sérsic bulge profile with light-weighted bulge-over-total ratio $B/T_{\text{F160W}} = 0.49$. Accordingly, Huertas-Company et al. (2015) classified the galaxy as spheroidal-like, with a probability of having a spheroid equal to 0.99 and a probability of having a disk equal to 0.18. The galaxy was flagged as a single-component “bad fit” in van der Wel et al. (2012), reinforcing the need of modeling the galaxy with more components to properly reproduce its structure.

The data cubes of the separated components allow us to calculate the relative contribution on the light coming from the bulge compared to the disk in each galaxy. Using our bulge and disk model data cube, we integrate the light of each component up to a radius where the galaxy surface brightness reaches
<figure><img src="image_3.png" /><figcaption>**Figure 2.** From left to right: postage RGB images for a typical galaxy in our sample (GDN 18522) and corresponding images in the three different filters: HST ACS F435W, F600W, and F850LP (top panels) and SHARDS f500w17, f602w17, and f840w17 (bottom panels). The galaxy is oriented north up, east to the left, and the field of view is 18 × 18 arcsec².</figcaption></figure>

26.5 mag arcsec⁻² in the WFC3 F160W filter, in order to maximize signal-to-noise ratio and avoid problems linked to the extrapolation of the model light. In Figure 5 we show an example of *B/T* for GDN 18552. First, in Figure 5 the spectral information provided by SHARDS medium-band filters could be appreciated in the distribution of *B/T* through wavelength. Second, it is remarkable how our analysis is sensitive to spectral features, displayed as an abrupt change in the *B/T* value in correspondence of the D4000 break at 400 nm or the [O III] emission line at 500.7 nm. For instance, the bulge is dimmer than expected based on a smooth wavelength interpolation exactly where an emission line can be present (most probably, revealing significant star formation in the disk). Finally, we find an overall good agreement between our values obtained from SHARDS data and the values estimated using broadband HST images presented in Dimauro et al. (2018). A quantitative test of the robustness of our decomposition of SHARDS images can be carried out by calculating *B/T* ratios for SHARDS bands and comparing them with those obtained for HST images at similar wavelengths (see Figure 14 in Appendix A). Considering all the galaxies in the sample, the average relative difference in the mass-weighted *B/T* measured from the two data sets is <16% and totally consistent with the statistical errors (see Tables 3 and 4 in Appendix A). The color dependence of the *B/T* ratio observed in the HST data is also highly consistent with that obtained with the SHARDS bands.

### 3.1.2. Error Analysis

It is well known that the minimization algorithms implemented in available routines for photometric decomposition do not usually provide a comprehensive representation of the real errors (Häussler et al. 2007; Costantin et al. 2017; Méndez-Abreu et al. 2017). Thus, for each galaxy in our sample and

each SHARDS band, we build a set of mock galaxies that take care of mimicking two sources of errors: (a) each image is perturbed pixel by pixel according to the background noise, and (b) mock galaxies are simulated varying their structural properties starting from the best-fitted value in the WFC3 F160W band provided in Dimauro et al. (2018). It is worth noting that the definition of the *n*-dimensional parameter space for each galaxy is critical, since the structural parameters (Sérsic index, size, and shape) are kept frozen throughout our spectrophotometric analysis, allowing only for variations of the relative intensity of the two components.

We use the No-U-Turn Sampler (NUTS; Hoffman & Gelman 2014) to find the posterior distribution of the model best fitting a set of data.<fnref n="19" /> Given an analytic model, this Markov Chain Monte Carlo (MCMC) algorithm, which closely resembles a Hamiltonian Monte Carlo (HMC) method, allows us to estimate the unknown posterior probability distribution. For each galaxy, our model

$$I_{\text{model}} = I_{\text{bulge}} + I_{\text{disk}} \quad (3)$$

is built to represent the 1D surface brightness radial profile. In this case, the intensity of the bulge component depends on the effective radius, the intensity measured at the effective radius, and the Sérsic index parameter, i.e., $I_{\text{bulge}}(I_e, R_e, n)$; on the other hand, the intensity of the disk component depends on the central intensity and the disk scale length, i.e., $I_{\text{disk}}(I_0, h)$. In detail, a normal prior probability distribution is assumed for each variable in the model, with standard deviation derived in Dimauro et al. (2018), and a normal likelihood is assumed for the total surface brightness $I_{\text{model}}$. We sample two independent

<footnote n="19">¹⁹ This analysis was done using PyMC3, a Python open-source probabilistic programming framework (Salvatier et al. 2016).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** 2D photometric decomposition of the galaxy GDN 18522 in the WFC3 F160W band as obtained from GASP2D. The left panel shows the ellipse-averaged radial profile of the surface brightness measured in the observed (black circles with gray error bars) and PSF-convolved modeled image (green solid line) and their corresponding difference. The surface brightness radial profiles of the best-fitting bulge (red dashed line) and disk (blue dashed-dotted line) are also shown in both linear and logarithmic scale for the distance to the center of the galaxy. The right panels (from top to bottom) show the map of the observed, modeled, and residual (observed-modeled) surface brightness distributions. The field of view is oriented with north up and east to the left. The vertical black dotted lines in the left panels and the black dotted ellipses in the right panels mark the radius where the galaxy surface brightness reaches 26.5 mag arcsec⁻².</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Same as Figure 3, but for SHARDS f500w17 (left), f602w17 (middle), and f840w17 (right) bands.</figcaption></figure>

chains, drawing 1000 samples and tuning with 3000 iterations. We evaluate the MCMC convergence by means of the Gelman–Rubin diagnostic $\hat{R}$ (Gelman & Rubin 1992), imposing a strict criterion for convergence $\hat{R} < 1.01$.

As a consequence of Bayes's theorem, the posterior probability distribution is evaluated, and the probability distributions of each parameter in the model are retrieved as the end result of MCMC. For each galaxy, we sample from this probability distribution to generate mock galaxies with physical properties that took into account correlations between parameters and intrinsic degeneracies. Each mock galaxy is

perturbed pixel by pixel according to the background noise and fitted using GASP2D as described in Section 3.1.1, using these new guesses for the structural parameters of the bulge and disk component. We perform 100 Monte Carlo (MC) realizations for each galaxy and each SHARDS band. This allows us to compute statistical errors for the flux of the two independent components, taking into account degeneracies and correlations between parameters, as well as possible biases in the photometric decomposition procedure.

The final values for the flux (and the structural parameters) of the bulge and the disk are retrieved as the median value of
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Bulge-over-total luminosity ratio as a function of wavelength for the galaxy GDN 18522. Diamonds, circles, and squares represent values for SHARDS, WFC3, and WIRCam bands, respectively. Errors are reported as the 16th–84th percentile interval. Blue crosses stand for values derived in Dimauro et al. (2018). From left to right, vertical dashed lines represent the location of [O II], D4000, Hβ, and [O III] features, respectively.</figcaption></figure>

the 100 MC realizations, while the statistical errors are computed as half of the 16th–84th percentile range, which corresponds to the standard deviation for normally distributed errors. In the final SED, for each component and each SHARDS band, we discard all the values of the flux that are compatible with zero according to their errors at 1σ level. This is mostly the case of shorter wavelengths, where the signal-to-noise ratio of the images does not allow us to consistently measure the flux, even for the total galaxy in the HST ACS F435W filter.

## 3.2. Stellar Populations

The reliable estimation of the stellar mass for a large and representative sample of bulges and disks, as well as a proper characterization of their SFH, is critical to quantify the evolutionary process. The typical degeneracies that affect the study of stellar populations in nearby and distant galaxies could be mitigated by the use of photometry with higher spectral resolution than broadband data (Pacifici et al. 2013). Using SHARDS data, the statistical significance of the best solution improves up to 10%–20% compared to the fit of broadband data alone (Pérez-González et al. 2013), granting estimations of stellar masses with typical uncertainties around 0.2 dex (Pérez-González et al. 2008; Barro et al. 2011).

Shortly, photometric data at different wavelengths are considered to be the end product of the galaxy SFH, assumed to be a declining delayed exponential:

$$SFR(t) \propto t/\tau^2 e^{-t/\tau}, \quad (4)$$

where τ runs from 200 Myr to a roughly constant SFH (τ = 100 Gyr). Increasing observational evidence and recent simulations justify the choice of an SFH described by a rising followed by a declining phase, irrespective of the specific parameterization (Behroozi et al. 2013; Pacifici et al. 2016; López Fernández et al. 2018; Costantin et al. 2019). Using the

synthesizer fitting code (see Pérez-González et al. 2003, 2008, for all details), we compare the measured SEDs with the Bruzual & Charlot (2003) stellar population library, assuming a Chabrier (2003) initial mass function integrated in the range 0.1 < M/M☉ < 100.

We allow the metallicity of the models to take discrete values Z/Z☉ = [0.4, 1, 2.5] (i.e., subsolar, solar, and super-solar). This allows us to explore the trend between metallicity and stellar mass, with more massive galaxies being more metal-rich than less massive ones (Gallazzi et al. 2005, 2014). The extinction is parameterized with a V-band attenuation assuming the extinction law of Calzetti et al. (2000), with values ranging from 0 to 3 mag. A minimization of the reduced χ² maximum likelihood estimator is used to search for the best-fitting model, allowing the stellar age to range from 1 Myr to the age corresponding to the age of the universe at the galaxy redshift.

We run 500 MC simulations for each galaxy’s component to estimate the uncertainties in the stellar population parameters, also accounting for possible degeneracies in the solutions, as comprehensively described in Domínguez Sánchez et al. (2016). Briefly, we model each photometric data point in the SED with a Gaussian distribution of width equal to the photometric uncertainty and randomly vary it, repeating the fit again and sampling from all possible models. The final solution is computed from the analysis of the (possible) different clusters in the τ-age parameter space, assigning a statistical significance to each of them according to the fraction of solutions belonging to a particular cluster. By means of the multidimensional age-τ-AV-Z space we compute the best solutions and the corresponding uncertainties as the median and the 68% confidence interval of the most significant cluster’s values, respectively. It is worth noting that the estimations of stellar masses from different clusters of solutions are totally compatible within their uncertainties (with average scatter of 0.1 dex), meaning that this parameter is robustly retrieved and not strongly affected by the degeneracies.

In this work, we individually fit the UV-to-near-IR SEDs of bulges and disks. Moreover, since we are assuming that the SFH of the galaxy is the sum of the bulge and disk one, we also fitted the total galaxy SED. The new physical parameters of the galaxy are mostly consistent with the ones provided in Barro et al. (2019) but could be considered as an improvement in terms of the complexity of the SFH of the galaxy. While the observed and model values for the galaxy’s SED are totally consistent, the little discrepancies in the galaxy physical parameters could arise from the different stellar models and parameter space. As an example, in Figure 6 we present the best model for the bulge, disk, and galaxy SED for the galaxy GDN 18522 (see also Appendix A).

## 4. Results

For the first time, we have applied a structure decomposition method to obtain the SEDs with spectral resolution R ~ 50 for bulges and disks in a representative sample of massive galaxies at redshift 0.14 < z ≤ 1. Furthermore, the spectral resolution and depth of the SHARDS data allow us to measure absorption indices (such as D4000) that are closely correlated to stellar ages (see Figure 16). In this paper, we focus our analysis of these SEDs on the characterization of the spheroids (either bulges surrounded by a disk or pure naked spheroids) at the mentioned redshifts. Indeed, we were able to properly reconstruct the SFHs
<figure><img src="image_7.png" /><figcaption>**Figure 6.** SED of the bulge (red), disk (blue), and galaxy (green) GDN 18522. Diamonds represent the individual photometric results of our decoupling analysis, while black circles represent the measured integrated photometry of the galaxy in Barro et al. (2019). Errors are reported as the 16th–84th percentile interval. The best models for the bulge, disk, and galaxy are shown as red, blue, and green lines, respectively. From left to right, vertical gray dashed lines represent the location of [O II], D4000, Hβ, and [O III] features, respectively.</figcaption></figure>

<table><caption>**Table 1**

Best Parameters for the Sample of Bulges and Pure Spheroids</caption><thead><tr><th>ID<br>(1)</th><th>log(M<sub>*</sub>)<br>(M<sub>☉</sub>)<br>(2)</th><th>τ&#773;<sub>M</sub><br>(Gyr)<br>(3)</th><th>z<sub>form</sub><br>(4)</th><th>τ<br>(Myr)<br>(5)</th><th>R<sub>c</sub><br>(kpc)<br>(6)</th><th>log(Σ<sub>1.5</sub>)<br>(M<sub>☉</sub> kpc<sup>-1.5</sup>)<br>(7)</th><th>Type<br>(8)</th></tr></thead><tbody><tr><td>750</td><td>10.6<sup>+0.1</sup><sub>-0.1</sub></td><td>2.5<sup>+0.3</sup><sub>-0.4</sub></td><td>1.8<sup>+0.2</sup><sub>-0.2</sub></td><td>193<sup>+31</sup><sub>-17</sub></td><td>1.0 ± 0.1</td><td>10.65 ± 0.07</td><td>B</td></tr><tr><td>775</td><td>10.7<sup>+0.1</sup><sub>-0.1</sub></td><td>1.4<sup>+0.1</sup><sub>-0.2</sub></td><td>1.4<sup>+0.1</sup><sub>-0.1</sub></td><td>194<sup>+27</sup><sub>-21</sub></td><td>0.67 ± 0.08</td><td>10.98 ± 0.08</td><td>B</td></tr><tr><td>912</td><td>11.0<sup>+0.1</sup><sub>-0.1</sub></td><td>3.7<sup>+0.4</sup><sub>-0.6</sub></td><td>1.5<sup>+0.2</sup><sub>-0.2</sub></td><td>2490<sup>+320</sup><sub>-270</sub></td><td>0.79 ± 0.04</td><td>11.20 ± 0.03</td><td>B</td></tr><tr><td>2104</td><td>10.1<sup>+0.1</sup><sub>-0.1</sub></td><td>1.1<sup>+0.2</sup><sub>-0.1</sub></td><td>1.0<sup>+0.1</sup><sub>-0.0</sub></td><td>200<sup>+27</sup><sub>-21</sub></td><td>0.40 ± 0.08</td><td>10.7 ± 0.1</td><td>B</td></tr><tr><td>5131</td><td>10.6<sup>+0.1</sup><sub>-0.1</sub></td><td>5.7<sup>+0.3</sup><sub>-0.7</sub></td><td>7.3<sup>+4.2</sup><sub>-3.1</sub></td><td>203<sup>+27</sup><sub>-24</sub></td><td>0.9 ± 0.2</td><td>10.7 ± 0.1</td><td>B</td></tr><tr><td>6098</td><td>10.7<sup>+0.1</sup><sub>-0.1</sub></td><td>2.2<sup>+0.3</sup><sub>-0.3</sub></td><td>1.7<sup>+0.2</sup><sub>-0.1</sub></td><td>198<sup>+27</sup><sub>-22</sub></td><td>1.3 ± 0.1</td><td>10.55 ± 0.06</td><td>B</td></tr><tr><td>6379</td><td>10.4<sup>+0.1</sup><sub>-0.1</sub></td><td>1.9<sup>+0.3</sup><sub>-0.3</sub></td><td>1.2<sup>+0.1</sup><sub>-0.1</sub></td><td>203<sup>+33</sup><sub>-23</sub></td><td>1.30 ± 0.06</td><td>10.22 ± 0.03</td><td>B</td></tr></tbody></table>

**Note.** Column (1): CANDELS ID of the galaxy. Column (2): stellar mass. Column (3): mass-weighted age. Column (4): redshift of formation. Column (5): timescale of exponentially declined SFH. Column (6): effective radius. Column (7): mass surface density. Column (8): type (B = bulge; PS = pure spheroid).
(This table is available in its entirety in machine-readable form.)

of our sample of spheroids, deriving fundamental physical quantities that constrain their stellar populations: the stellar mass ($M_*$), the age $t_0$, the star formation timescale $\tau$, the metallicity (Z), and the dust attenuation ($A_V$). In Table 1 we report the main properties of the sample spheroids.

### 4.1. Mass-weighted Ages and Formation Redshift

We consistently characterize the SFH of both bulges and pure spheroids, pushing the stellar population analysis of individual structural components of galaxies to an unexplored redshift range. In particular, we compute their mass-weighted age

$$ \bar{t}_M = t_0 - \frac{\int_0^{t_0} \mathrm{SFR}(t) \times \mathrm{tdt}}{\int_0^{t_0} \mathrm{SFR}(t) \mathrm{dt}}, \qquad (5) $$

where time runs from the start of the star formation in the galaxy onward, i.e., $t_0$ Gyr before the universe age corresponding to the galaxy redshift. This mass-weighted age represents a better approximation to the average age of the stellar population, which takes into account the extent of the star formation and allows us to mitigate the age-$\tau$ degeneracy (see Table 1).
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Mass-weighted stellar ages of bulges (red circles) and pure spheroids (black triangles) as a function of their stellar mass. Errors are reported as the 16th–84th percentile interval. The red and black shaded regions show the distribution of the 500 MC realizations for each bulge and pure spheroid, respectively. The gray dashed–dotted horizontal line marks $\bar{t}_M = 3$ Gyr. Red arrows mark upper limits for four systems with mass-weighted ages $\bar{t}_M < 100$ Myr. The red and black histograms represent the frequency of the mass-weighted ages of bulges and pure spheroids, respectively. The dotted black histogram shows the distribution for the total spheroidal population.</figcaption></figure>

In Figure 7 we present the mass-weighted ages of spheroids as a function of their stellar mass. The first result is that we find a very old population of bulges ($\bar{t}_M > 6–7$ Gyr). In particular, bulges at redshift $0.14 < z \le 1$ probed by our survey clearly display a bimodal distribution: 52% of them have $\bar{t}_M < 3$ Gyr with median value $\bar{t}_M = 1.8_{-0.8}^{+0.6}$ Gyr, while 48% have $\bar{t}_M > 3$ Gyr with median value $\bar{t}_M = 6.3_{-1.1}^{+2.0}$ Gyr. This behavior could either reflect a sharp difference in the formation mode of the two types of spheroids or be due to some kind of selection bias that we cannot identify (see Appendix B). There are not massive bulges ($M_* > 3 \times 10^{10} M_\odot$) with young ages ($\bar{t}_M < 1$ Gyr). Second, we find that pure spheroids are on average younger than bulges: they have median mass-weighted ages $\bar{t}_M = 1.3_{-0.2}^{+1.2}$ Gyr and $\bar{t}_M = 2.7_{-1.6}^{+3.9}$ Gyr, respectively. Considering the global spheroidal population (bulges and pure spheroids), they have a median mass-weighted age $\bar{t}_M = 1.7_{-0.7}^{+3.5}$ Gyr. When combining all spheroids, the age bimodality holds, and it is dominated by systems less than 3 Gyr old (68% of the entire sample): young spheroids ($\bar{t}_M < 3$ Gyr) have a median mass-weighted age $\bar{t}_M = 1.3_{-0.6}^{+0.7}$ Gyr, while older spheroids ($\bar{t}_M > 3$ Gyr) present $\bar{t}_M = 5.3_{-1.3}^{+2.1}$ Gyr. At $M_* \gtrsim 7 \times 10^{10} M_\odot$ the population of old spheroids dominates, while less massive systems show similar masses when divided into old and young ones. Four spheroids are caught in the middle of their formation process ($\bar{t}_M < 100$ Myr) at the redshift of observation.

It is worth remembering that our sample of galaxies spans along a range in redshift $0.14 < z \le 1$ (see Figure 1), which translates to a wide range in cosmic time ($t \sim 6$ Gyr). Thus, since we are dealing with galaxies observed at different epochs,

we derive the redshift corresponding to mass-weighted ages (i.e., formation redshift $z_{form}$) to better understand their evolutionary pathways (see Table 1). In Figure 8 we show the correlation between formation redshift and stellar mass of both bulges and pure spheroids. In this way, we confirm the observation of a very old population of bulges. The bimodality in formation redshift reflects the one in mass-weighted ages: 67% of bulges have $z_{form} < 3$ with median value $z_{form} = 1.3_{-0.6}^{+0.6}$, while 33% have $z_{form} > 3$ with median value $z_{form} = 6.2_{-1.7}^{+1.5}$. At $M_* > 7 \times 10^{10} M_\odot$, 67% of bulges are formed at redshift $z_{form} > 3.8$. We find that pure spheroids are formed (on average) later than bulges. The totality of pure spheroids builds up at redshift $z_{form} \lesssim 2$, with a median formation redshift $z_{form} = 1.1_{-0.3}^{+0.3}$. Considering the global spheroidal population (bulges and pure spheroids), 19% of the systems formed in a first wave at redshift $z_{form} > 3$ (with median $z_{form} = 6.2_{-1.8}^{+1.5}$), while the majority of spheroids (81%) assemble in a second wave at redshift $z_{form} < 3$ (with median $z_{form} = 1.2_{-0.4}^{+0.4}$). Finally, there is a positive trend between formation redshift and stellar mass. For $M_* > 7 \times 10^{10} M_\odot$ the spheroidal population is dominated by systems formed at $z_{form} > 3$ (77%), while at lower mass 88% of the population is formed at $z_{form} < 3$.

We tested that the bimodality in the mass-weighted age and formation redshift of our bulges is not biased by the redshift of the observed galaxies. The two waves of bulge formation are also observed when dividing the sample into two redshift ranges, $z \le 0.75$ and $z > 0.75$ (see Appendix B). On the other hand, pure spheroids do not display a bimodal distribution and have similar mass-weighted ages and formation redshifts to the second wave of bulges.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Formation redshift of bulges (red circles) and pure spheroids (black triangles) as a function of their stellar mass. Errors are reported as the 16th–84th percentile interval. The red and black shaded regions show the distribution of the 500 MC realizations for each bulge and pure spheroid, respectively. The gray dashed–dotted horizontal line marks the formation redshift $z_{\text{form}} = 3$. Red arrows mark upper limits for two systems with masses $M_{*} < 2 \times 10^{10} M_{\odot}$. The red and black histograms represent the frequency of the formation redshift of bulges and pure spheroids, respectively. The dotted black histogram shows the distribution for the total spheroidal population.</figcaption></figure>

In summary, we find an old population of bulges, with 33% of them being already in place by $z_{\text{form}} > 3$. Spheroids in the redshift range $0.14 < z \le 1$ display a bimodal distribution: systems with median mass-weighted age $\bar{t}_{\text{M}} = 1.3$ Gyr, formed at median redshift $z_{\text{form}} = 1.2$, and considerably older bulges (“embedded” spheroids), with a median mass-weighted age of 5.3 Gyr, formed at median redshift $z_{\text{form}} = 6.2$. Pure spheroids (which do not present a disk) belong to the first type, having median mass-weighted age of $\bar{t}_{\text{M}} = 1.3_{-0.2}^{+1.2}$ Gyr and all of them being formed at $z \lesssim 2$.

Hints for an early formation of bulges in disk galaxies were already provided by Morelli et al. (2016) by direct age measurements in a sample of 12 local galaxies. Despite the low statistics, they identify a double population of bulges: seven young bulges (2–7 Gyr) with solar metallicity and five old bulges ($\gtrsim$13 Gyr) with a large spread in metallicity. Nonetheless, regardless of the actual bimodality, one of the main results of this work is that a fraction of bulges formed at very high redshift. Despite studying a small sample of 10 star-forming galaxies at redshift $0.45 < z < 1$, Mancini et al. (2019) already saw hints that quiescent bulges present ages approaching the age of the universe at the time of observation.

Our results agree and provide a deeper insight about the trend recently found by independent spectroscopic observation of galaxies from redshift $z = 0.7$ to redshift $z = 2.5$ (Gallazzi et al. 2014; Belli et al. 2019; Carnall et al. 2019). In particular, Gallazzi et al. (2014) studied a sample of ~70 galaxies at redshift $z \sim 0.7$, finding average r-band light-weighted ages from ~2.3 Gyr ($z_{\text{form}} \sim 1.3$) at $\log(M_{\star}/M_{\odot}) = 10.5$ to ~2.9 Gyr ($z_{\text{form}} \sim 2.2$) at $\log(M_{\star}/M_{\odot}) = 11.4$. Considering the average bulge-over-total mass ratio ($B/T_{\text{mass}} = 0.51$ of our sample

galaxies, the mass range proven by our second-wave bulges overlaps with their sample galaxies. On the other hand, our first wave of spheroid formation could be compared with higher-redshift observations. Carnall et al. (2019) derived the formation redshift of 75 massive ($\log(M_{\star}/M_{\odot}) > 10.3$) UVJ-selected galaxies at redshifts of $1.0 < z < 1.3$, finding a population of old systems ($z_{\text{form}} > 3$) that have higher stellar masses compared to a younger population that dominates their sample. This trend is also consistent with the one provided by Belli et al. (2019) studying 24 quiescent galaxies at $1.5 < z < 2.5$. It is worth noting that the scenario holds besides the differences in the details of the SFHs and the different spectrophotometric data sets, providing really strong constraints for cosmological simulations aiming at reproducing the observed universe.

## 4.2. Morphological Properties

In order to understand what drives the formation redshift of the substructures of massive galaxies and the possible differences in the morphological properties of older and younger spheroids, we plot in Figure 9 their mass–size relations. For pure spheroids we use sizes according to their half-light radius in the WFC3 F160W band from Dimauro et al. (2018), while we characterize the size of the bulge using the best value of its half-light radius obtained in the spectrophotometric decoupling using the WFC3 F160W filter (see Section 3.1.2). The stellar masses of spheroids are derived from the stellar population analysis described in Section 3.2.

We compare the position of bulges and pure spheroids in the mass–size plane with the trend expected for early- and late-type
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Mass–size relations for bulges (circles) and pure spheroids (triangles), color-coded according to their formation redshift $z_{\text{form}}$. Errors are reported as the 16th–84th percentile interval. The dark-red and blue dashed lines correspond to the best-fitting trend for early-type and late-type galaxies at redshift $z = 0.75$ in van der Wel et al. (2014), respectively. The light-red dashed–dotted line stands for the best-fitting trend of bulges at redshift $z \sim 0$ in Méndez-Abreu et al. (2021). Black arrows mark upper limits for two systems with masses $M_{\star} < 2 \times 10^{10} M_{\odot}$. The red and black histograms represent the frequency of the mass-weighted ages of bulges and pure spheroids, respectively. The dotted black histogram shows the distribution for the total spheroidal population.</figcaption></figure>

galaxies described in van der Wel et al. (2014). The two populations actually correspond to star-forming and quiescent galaxies classified by means of a rest-frame color selection (see Wuyts et al. 2007; Williams et al. 2009). In particular, we use the best-fitting relation at redshift $z = 0.75$, but similar trends are observed when considering the evolving mass–size relationships described in other works (Buitrago et al. 2008; Damjanov et al. 2009; van der Wel et al. 2014). Both bulges and pure spheroids lie in the region of early-type galaxies, having masses and sizes that are not compatible with late-type systems. The bulge population is more compact than its counterpart at redshift $z \sim 0$, consistent with an increased star formation activity in the galaxy’s most internal regions at earlier times. We note in Figure 9 that they are offset from the best-fitting relation in Méndez-Abreu et al. (2021) toward lower sizes and higher masses.

Interestingly, at fixed stellar mass, bulges are on average smaller than pure spheroids in the same redshift range. The first ones have median sizes $R_e = 1.0_{-0.4}^{+0.9}$ kpc, while the latter show $R_e = 1.9_{-0.6}^{+1.4}$ kpc. Indeed, they are two different populations in terms of size, as proven by means of the Kolmogorov–Smirnov statistics ($p$-value < 1%). Bulges extend to lower sizes, that is, pure spheroids are not found among the smallest spheroids: 92% of them have $R_e > 1$ kpc, while 48% of bulges have $R_e \le 1$ kpc. Although it seems that there is a continuity between these two classes, a multivariate Wald–Wolfowitz test reinforces the result that they are not the same population. At low redshift, various studies confirm that elliptical galaxies and

bulges are two different populations in terms of their size (Gadotti 2009; Méndez-Abreu et al. 2021). Local bulges and massive elliptical galaxies seem to follow offset mass–size relations and occupy different loci of the fundamental plane (Djorgovski & Davis 1987), providing clues to their different evolutionary processes (Gadotti 2009). Nevertheless, it has to be noted that this result is less strict for low-luminosity elliptical galaxies, which fits in the picture that elliptical galaxies actually harbor a more compact (and probably older) component in their core (de la Rosa et al. 2016).

Regarding their evolution, young spheroids ($z_{\text{form}} < 3$) have median sizes $R_e = 1.4_{-0.6}^{+1.5}$ kpc and older ones ($z_{\text{form}} > 3$) have $R_e = 1.3_{-0.6}^{+0.8}$ kpc. In particular, we find a lack of extended ($R_e > 2.5$ kpc) spheroids among the population that formed at the highest redshift.

The most important result we can infer from the mass–size relation in Figure 9 is that first-wave bulges are more compact than second-wave ones: they are not only slightly smaller but also more massive. Since we find no clear trend either with the Sérsic index of the galaxy or with the Sérsic index of the bulge, to properly quantify the evolution of spheroids, we plot in Figure 10 their formation redshift as a function of their mass surface density $\Sigma_{\alpha} = MR_e^{-\alpha}$ (see Table 1). According to the definition in Barro et al. (2013), we define $\Sigma_{1.5}$ ($\alpha = 1.5$), such as it lies between the mass surface density $M/R_e^2$ and $M/R_e$, both of which strongly correlated with color and star formation rate (SFR) up to high redshifts (Franx et al. 2008). Following Barro et al. (2013), we could identify our spheroids with their
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Mass surface density $MR_e^{-1.5}$ of bulges (circles) and pure spheroids (triangles) as a function of their mass, color-coded according to their formation redshift $z_{form}$. Errors are reported as the 16th–84th percentile interval. Systems are separated between compact and extended ones at $log(\Sigma_{1.5}) = 10.3 M_{\odot} \text{ kpc}^{-1.5}$ (dashed gray horizontal line; Barro et al. 2013).</figcaption></figure>

compact population of quiescent galaxies, defined by $log(\Sigma_{1.5}) > 10.3 M_{\odot} \text{ kpc}^{-1.5}$. Young bulges ($z_{form} < 3$) have median mass surface density $log(\Sigma_{1.5}) = 10.2_{-0.4}^{+0.5} M_{\odot} \text{ kpc}^{-1.5}$, while older ones have median mass surface density $log(\Sigma_{1.5}) = 10.6_{-0.4}^{+0.4} M_{\odot} \text{ kpc}^{-1.5}$. While bulges span a wider range of mass surface densities (median values $log(\Sigma_{1.5}) = 10.4_{-0.6}^{+0.5} M_{\odot} \text{ kpc}^{-1.5}$), 86% of pure spheroids have $log(\Sigma_{1.5}) < 10.3 M_{\odot} \text{ kpc}^{-1.5}$.

Summarizing, the first-wave bulges are more compact than the second-wave bulges, and bulges are more compact than pure spheroids.

# 5. Discussion

In this first paper of a series we focus on studying independently the separate structures in disk galaxies up to redshift $z=1$, which allows us to compare the properties of bulges with pure spheroids. Although a plethora of high-resolution data are collected at redshift $z \sim 0$, it still remains difficult to resolve SFHs when galaxies are older than 5 Gyr owing to the similarity of their stellar spectra. On the contrary, the window $0.14 < z \le 1$ gives us the advantage of discerning different stellar populations of galaxies up to higher formation redshift, simply because the stellar ages are bounded by the age of the universe at the epoch of observation. Moreover, the SHARDS spectral resolution and wavelength coverage range permit us to probe young stellar populations owing to the presence of optical spectral features in their rest-frame SEDs.

## 5.1. Bulges Form in Two Waves

In this work we prove that the population of more massive spheroids ($M_{\star} > 7 \times 10^{10} M_{\odot}$) is mainly composed by systems already in place at very early cosmic time ($z_{form} \gtrsim 4$), while the low-mass end is dominated by spheroids formed in a prolonged time span of $\sim 8$ Gyr. Our result suggests that in the early phase of galaxy evolution the conditions were very favorable for rapidly assembling and quenching the building blocks of

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Distribution of timescales as a function of formation redshift for our bulges (first wave: purple circles; second wave: orange circles) and pure spheroids (black triangles). The red shaded region shows the density distribution of the 500 MC realizations for each bulge. The purple and orange histograms represent the frequency of the timescales for all the MC realizations of first- and second-wave bulges, respectively.</figcaption></figure>

today's galaxies. Indeed, the population of bulges formed at $z_{form} > 4$ can be considered as the relic of the early universe.

Star formation can be very rapidly quenched at earlier cosmic time, as suggested by the presence of quiescent galaxies at redshifts $z \gtrsim 4$ (Straatman et al. 2014). They are typically highly compact elliptical-like galaxies (Daddi et al. 2005; Damjanov et al. 2009). We further extend this picture, finding a bimodal distribution of bulges in terms of their formation redshift. Indeed, these spheroids present an extended stellar disk component at the time of observation. Until today, great effort was put into studying the quiescent population at high redshift, selecting it by means of the galaxy color. Although the UVJ selection has proven efficient in selecting genuinely quiescent galaxies, it could miss a significant fraction of mass hidden within spheroids in luminous star-forming disk galaxies (Schreiber et al. 2018).

## 5.2. Spheroids Evolve in Two Modes

The question that arises is whether the bimodal distribution in the ages of the spheroidal population is due to a different channel of evolution. Looking at Figure 8, our analysis seems to suggest that in the redshift range $2 \lesssim z \lesssim 4$ the processes responsible for assembling the more massive spheroids were less efficient, even though it remains difficult to state whether it is a reflection of different pathways and transitional epochs in galaxies' evolution. In order to quantify the formation timescale of the spheroids in our sample, we show in Figure 11 their $\tau$ distributions. The buildup of the spheroidal component not only comes in two waves but also seems to be characterized by two different modes of formation: a fast and a slow one. The majority of our systems form in short timescales, having a median $\tau = 203_{-4}^{+160}$ Myr. However, first-wave bulges have average timescale $\tau = 232 \pm 16$ Myr, while second-wave bulges present $\tau = 716 \pm 203$ Myr. The scatter of the $\tau$ distribution is $\sim 90$ Myr for first-wave bulges and $\sim 1.6$ Gyr for second-wave ones. Among second-wave bulges, 15 out of 61 (25%) have $\tau > 500$ Myr and 31% have $\tau > 300$ Myr. On the other hand, 97% (29 out of 30) of first-wave bulges have $\tau < 500$ Myr and 90% have $\tau < 300$ Myr. Furthermore, we
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Averaged SFHs of the first-wave bulges (purple solid line), second-wave bulges (orange dashed line), and pure spheroids (black solid line). The red and black shaded curves represent the corresponding 16th–84th percentile interval. The blue shaded area indicates the redshift studied in this work.</figcaption></figure>

present in Figure 12 the averaged SFHs separately for bulges and pure spheroids. The fast mode is mostly the only channel of evolution at earlier cosmic epochs ($z \gtrsim 5$). On the other hand, the slow mode starts to be relevant only in the last $\sim 10$ Gyr ($z \lesssim 2–3$). This picture is consistent with an increasing number of results, where the fast mode dominates at redshift $z \gtrsim 2.5$ but the two channels have to be in place simultaneously and with similar quenching rates at $z \sim 2$ (Schawinski et al. 2014; Wild et al. 2016; Belli et al. 2019).

### 5.3. The Old Population of Compact Bulges

Another piece of information is provided by the mass-size evolution of our spheroidal population. Comparing systems at different formation redshift (see Figure 9), we find that older bulges are more compact than younger spheroids (both bulges and pure spheroids; Estrada-Carpenter et al. 2020). This finding agrees with the fact that not only do two waves of spheroid formation seem to be in place but also fast-mode spheroids appear more compact than slow-mode ones (see also Belli et al. 2019). On the other hand, first- and second-wave bulges have similar sizes of $R_e = 1.3_{-0.6}^{+0.8}$ kpc and $R_e = 1.0_{-0.4}^{+0.8}$ kpc, respectively. Despite being a result already observed in Whitaker et al. (2012), dividing the red sequence of their quiescent galaxies into blue and red halves, it is worth noting that this lack of evolution is seen here for the first time in the bulge component.

As stated before, our results provide the evidence of a population of bulges that formed by compaction in the early universe. They are fast-track spheroids, which present similar sizes but higher masses compared to their younger counterparts. In this context, we favor a scenario in which older spheroids in our sample were formed by violent disk instabilities and clump migration, rather than mergers. The higher abundance of cold gas available in the early universe could provide the ideal conditions to rapidly generate very efficient starbursts, which result in very compact and massive systems. Cosmological simulations presented by Zolotov et al. (2015) show that their most massive galaxies start their

compaction phase at redshift $z \gtrsim 4$ and successfully quench by redshift $z \sim 2$. Furthermore, recent results from cosmological simulations presented in Ceverino et al. (2018) show that high SFRs are driven by high gas accretion rates and the successive compaction of the stellar systems. Indeed, in the case in which the compaction event induces an intense star formation burst at $z = 10$ with maximum specific SFR $\sim 20$ Gyr$^{-1}$, this translates to SFR $\sim 200$ $M_{\odot}$ yr$^{-1}$ for $10^{10}$ $M_{\odot}$ at $z = 10$, mimicking the picture we are providing in Figure 12.

The compaction phase of galactic disks into spheroids naturally leads to quenching the systems stabilizing the violent disk instabilities. The direct result of building this mass concentration is rapidly increasing the angular velocity of the gas (morphological quenching; Martig et al. 2009), diminishing the surface density of the cold gas component (i.e., star formation and outflows; Forbes et al. 2014), and increasing the gas radial velocity dispersion (i.e., stellar or AGN feedback; Krumholz & Thompson 2013). The consequent suppression of the star formation is followed by a gradual growth and expansion into a larger spheroidal galaxy that could develop an extended stellar disk surrounding the spheroid owing to late and slower gas accretion. Indeed, at the time of observation, the bulge population presents an extended stellar disk, which fits the more efficient redistribution of the gas driven by violent disk instabilities rather than mergers, which could lead to a disruption of the disky kinematics, diminishing more abruptly the angular momentum of the forming system.

### 5.4. The Second Dominant Wave of Spheroids

The totality of pure spheroids in our sample is formed at the epoch of the star formation activity peak (Lilly et al. 1996; Madau et al. 1996; Madau & Dickinson 2014). Likewise, more than half of bulges were also formed in the last 10 Gyr of the universe life. This second (dominant) wave, which comprises a younger population of spheroids, is more challenging to interpret. While the majority of spheroids are still formed in a fast mode, a slower mode starts to take place. Again, there is a subpopulation of bulges that is more compact and formed in shorter timescales, pointing toward a formation by compaction at delayed times with respect to the older population. On the other hand, there are pure spheroids and bulges that are on average less compact and formed with a variety of timescales. It is difficult to discriminate any scenario, and probably different ones start to contribute at some level in shaping the properties of those spheroids.

Finally, we find that bulges are on average more compact than pure spheroids. This made us question the true nature of these supposedly pure spheroids: are they hosting an inner (older) component that belongs to the same family of older bulges? Theoretical arguments favor minor mergers being responsible for their size growth, since it would be the physical processes allowing them to efficiently increase the size of the spheroid compared to the growth of its stellar mass (Naab et al. 2009; Hopkins et al. 2010). After the compaction phase, the high-redshift spheroid is likely to experience further accretion, becoming more star dominated as it evolves (Oser et al. 2010; Porter et al. 2014). This wet and/or dry growth could be substantial, leading the spheroid to acquire an extended stellar envelope and making it appear like a typically observed elliptical galaxy in the local universe (van Dokkum et al. 2014; Buitrago et al. 2017).
<figure><img src="image_14.png" /><figcaption>Figure 13. Illustration of the proposed formation scenario for the spheroidal component of massive galaxies at redshift $0.14 < z \le 1$.</figcaption></figure>

<table><caption>Table 2
Median Physical Properties of the Spheroidal Population at Redshift $0.14 < z \le 1$</caption><thead><tr><th>Type</th><th>log(M<sub>*</sub>)<br>(M<sub>☉</sub>)</th><th>I&#773;<sub>M</sub><br>(Gyr)</th><th>z<sub>form</sub></th><th>τ<br>(Myr)</th><th>R<sub>e</sub><br>(kpc)</th><th>log(Σ<sub>1.5</sub>)<br>(M<sub>☉</sub> kpc<sup>-1.5</sup>)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th></tr></thead><tbody><tr><td>First-wave bulges</td><td>10.8<sup>+0.2</sup><sub>-0.3</sub></td><td>6.5<sup>+1.5</sup><sub>-1.4</sub></td><td>6.2<sup>+1.5</sup><sub>-1.7</sub></td><td>203<sup>+15</sup><sub>-3</sub></td><td>1.3<sup>+0.8</sup><sub>-0.6</sub></td><td>10.6<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>Second-wave bulges</td><td>10.3<sup>+0.4</sup><sub>-0.5</sub></td><td>1.7<sup>+2.0</sup><sub>-0.8</sub></td><td>1.3<sup>+0.6</sup><sub>-0.6</sub></td><td>213<sup>+790</sup><sub>-13</sub></td><td>1.0<sup>+0.8</sup><sub>-0.4</sub></td><td>10.2<sup>+0.5</sup><sub>-0.4</sub></td></tr><tr><td>Pure spheroids</td><td>10.5<sup>+0.3</sup><sub>-0.2</sub></td><td>1.3<sup>+1.2</sup><sub>-0.2</sub></td><td>1.1<sup>+0.3</sup><sub>-0.3</sub></td><td>201<sup>+114</sup><sub>-3</sub></td><td>1.9<sup>+1.4</sup><sub>-0.6</sub></td><td>10.1<sup>+0.2</sup><sub>-0.3</sub></td></tr></tbody></table>

Note. Column (1): spheroidal type. Column (2): stellar mass of each component. Column (3): mass-weighted age. Column (4): redshift of formation. Column (5): timescale of exponentially declined SFH. Column (6): effective radius. Column (7): mass surface density.

We qualitatively illustrate in Figure 13 the proposed scenario for the formation and evolution of bulges and pure spheroids in massive galaxies at $0.14 < z \le 1$. Moreover, we summarize the main properties of the different spheroids in Table 2. The bulges in $z \le 1$ massive disk galaxies form in two waves. A first wave of bulges builds up fast in the early universe ($z > 3$) through very dissipative processes such as violent disk instabilities. They undertake a compaction phase, probably evolving to the well-known $z \sim 2$ red nuggets, and then acquire a disk. A second wave of spheroids forms at $z = 1-2$; some of them accrete new material and acquire a disk by $z \sim 0.5$, becoming a bulge, and some of them reach our sample as a pure spheroid. The first wave is characterized by short formation timescales and large peak SFR. A relatively slower mode of formation starts to be relevant for the second-wave bulges and pure spheroids, presenting longer timescales and smaller peak SFRs. As mentioned earlier, first-wave bulges are characterized by their compactness: they present similar sizes but are more massive than second-wave bulges. Pure spheroids are larger than bulges and present similar masses, i.e., they are not as compact.

## 6. Conclusions

In this work we investigate the assembly history and evolutionary pathways of the spheroidal structures within

massive galaxies at redshift $z \le 1$. For that purpose, we use a sample of massive galaxies selected by the SHARDS spectrophotometric survey in GOODS-N. We decomposed spectrophotometrically the light of our galaxies into two main stellar components (i.e., the central bulge and the outer extended disk), retrieving their separate SEDs with a resolution $R \sim 50$. We fit the photometry for the spheroids to stellar population synthesis models, characterizing the SFH of the 91 bulges and comparing their properties with the 65 pure spheroids in our sample, all at redshift $0.14 < z \le 1$.

By deriving the mass-weighted ages and the formation redshifts of spheroids, we find that they form in two waves. We distinguish a fraction of very old bulges, which are formed at redshift $z_{form} > 3$, and a dominant population of spheroids (bulges and pure spheroids) that formed at median redshift $z_{form} = 1.2^{+0.4}_{-0.4}$. Furthermore, spheroids not only form in two waves but also form in two modes. At higher redshift, a fast mode (timescale around 200 Myr) is driving the rapid evolution of those systems, while at lower redshift a slower mode starts to become relevant (with timescales ranging from 200 Myr to 1 Gyr and longer). Finally, the old population of bulges is more compact than the young population of spheroids. We propose that the first wave of formation is characterized by a violent compaction phase, which builds up distinctively dense spheroids in short timescales.
Considering the rarity of high-z red nuggets surviving in the local universe (Trujillo et al. 2009; Tortora et al. 2016; Charbonnier et al. 2017; Buitrago et al. 2018), our results have important implications for the evolution of this population. The population of red nuggets, formed at high redshift through rapid compaction (Damjanov et al. 2009; Dekel et al. 2009), could possibly evolve in today's early-type galaxies or settle in the center of a disk galaxy. Recently, Costantin et al. (2020) suggested that disk galaxies hide the remnants of the compact and quiescent population observed at high redshift. The results of this work reinforce the increasing evidence that the central regions of early-type galaxies actually harbor the population of spheroids formed at high redshift (MacArthur et al. 2003; Graham 2014; Graham et al. 2015), broadening this comprehensive and compelling picture to later Hubble types.

We would like to thank the anonymous referee for improving the content of the manuscript. L.C. wishes to thank Cristina Cabello for the support provided while this project was devised and Michele Perna for the useful discussion.

L.C. acknowledges financial support from Comunidad de Madrid under Atracción de Talento grant 2018-T2/TIC-11612. L.C. and P.G.P.-G. acknowledge support from Spanish Ministerio de Ciencia, Innovación y Universidades through grant PGC2018-093499-B-I00. J.M.-A. acknowledges support from the MCIU by the grant AYA2017-83204-P and the Programa Operativo FEDER Andalucía 2014-2020 in collaboration with the Andalucian Office for Economy and Knowledge. D.C. is a Ramon-Cajal Researcher and is supported by the Ministerio de Ciencia, Innovación y Universidades (MICIU/FEDER) under research grant PGC2018-094975-C21.

This work has made use of the Rainbow Cosmological Surveys Database, which is operated by the Centro de

Astrobiología (CAB/INTA), partnered with the University of California Observatories at Santa Cruz (UCO/Lick, UCSC).

# Appendix A

## Spectrophotometric Decoupling: Robustness of the Method

In this appendix we describe the robustness of our decomposition process and provide some more examples of the spectrophotometric decoupling performances.

As described in Section 3.1, we combined the HST and SHARDS data to get SEDs for bulges and pure spheroids with a spectral resolution $R \sim 50$. In this sense, we did not fit two components blindly and independently to the SHARDS data, but we use the decomposition of the HST images as priors. We fit the characteristic surface brightness levels of bulges and disks in SHARDS bands, keeping their structural parameters (i.e., Sérsic index, size, and shape) frozen from the 2D decomposition of the galaxy light in HST filters. Since only the relative intensity of the bulge and disk component is allowed to vary in the decomposition of SHARDS images, we compare in Figure 14 the trend of $(B/T)_{HST}$ and $(B/T)_{SHARDS}$, dividing the galaxies according to the size of their bulges. The sensibility of our decomposition process is presented in Tables 3 and 4. We find that the average relative error between $(B/T)_{HST}$ and $(B/T)_{SHARDS}$ is less than 16% for galaxies hosting small bulges ($R_{e,b} < 0.2''$), and it is less than 11% for galaxies hosting large bulges ($R_{e,b} > 0.2''$).

We randomly pick four galaxies out of our sample (i.e., GDN 3360, GDN 9386, GDN 17242, and GDN 20441) in order to demonstrate how the 2D photometric decomposition performed on HST images transfers to SHARDS data, as detailed in Section 3.1. We show in Figure 15 how the synergy between the two data sets allows us to fit the individual SEDs of bulges and disks (right panels), characterizing their SFHs, as detailed in Section 3.2.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Luminosity-weighted *B/T* as a function of wavelength (circles) for galaxies with *R*<sub>e,b</sub> < 0.2 (left panel) and *R*<sub>e,b</sub> > 0.2 (right panel). Circles are color-coded according to the bulge effective radius. Diamonds represent the mean values measured from our SHARDS photometry, while squares stand for mean values retrieved from Dimauro et al. (2018). Error bars stand for the standard error of the mean.</figcaption></figure>

<table><caption>**Table 3**
Comparison of *B/T* from SHARDS and HST Photometry for Galaxies with *R*<sub>e,b</sub> < 0.2</caption><thead><tr><th>λ<sub>rest</sub><br>(nm)<br>(1)</th><th>(B/T)<sub>HST</sub><br>(2)</th><th>σ<sub>HST</sub><br>(3)</th><th>(B/T)<sub>SHARDS</sub><br>(4)</th><th>σ<sub>SHARDS</sub><br>(5)</th><th>Relative<br>Difference<br>(6)</th></tr></thead><tbody><tr><td>0.26–0.32</td><td>0.28 ± 0.02</td><td>0.16</td><td>0.24 ± 0.01</td><td>0.14</td><td>14%</td></tr><tr><td>0.32–0.38</td><td>0.29 ± 0.03</td><td>0.12</td><td>0.28 ± 0.01</td><td>0.15</td><td>4%</td></tr><tr><td>0.38–0.44</td><td>0.37 ± 0.02</td><td>0.13</td><td>0.31 ± 0.01</td><td>0.17</td><td>16%</td></tr><tr><td>0.44–0.50</td><td>0.40 ± 0.02</td><td>0.13</td><td>0.36 ± 0.01</td><td>0.19</td><td>11%</td></tr><tr><td>0.50–0.56</td><td>0.41 ± 0.02</td><td>0.10</td><td>0.34 ± 0.02</td><td>0.19</td><td>15%</td></tr><tr><td>0.56–0.62</td><td>0.34 ± 0.03</td><td>0.12</td><td>0.34 ± 0.02</td><td>0.16</td><td>&lt;1%</td></tr></tbody></table>

**Note.** Column (1): rest-frame wavelength range. Column (2): mean bulge-over-total luminosity ratio and standard error (HST; Dimauro et al. 2018). Column (3): standard deviation. Column (4): mean bulge-over-total luminosity ratio and standard error (SHARDS). Column (5): standard deviation. Column (6): relative difference between $(B/T)_{HST}$ and $(B/T)_{SHARDS}$.

<table><caption>**Table 4**
Same as Table 3, but for Galaxies with *R*<sub>e,b</sub> > 0.2</caption><thead><tr><th>λ<sub>rest</sub><br>(nm)<br>(1)</th><th>(B/T)<sub>HST</sub><br>(2)</th><th>σ<sub>HST</sub><br>(3)</th><th>(B/T)<sub>SHARDS</sub><br>(4)</th><th>σ<sub>SHARDS</sub><br>(5)</th><th>Relative<br>Difference<br>(6)</th></tr></thead><tbody><tr><td>0.26–0.32</td><td>0.25 ± 0.05</td><td>0.24</td><td>0.25 ± 0.02</td><td>0.16</td><td>4%</td></tr><tr><td>0.32–0.38</td><td>0.33 ± 0.05</td><td>0.19</td><td>0.33 ± 0.02</td><td>0.19</td><td>&lt;1%</td></tr><tr><td>0.38–0.44</td><td>0.41 ± 0.05</td><td>0.24</td><td>0.39 ± 0.01</td><td>0.20</td><td>6%</td></tr><tr><td>0.44–0.50</td><td>0.48 ± 0.05</td><td>0.22</td><td>0.43 ± 0.02</td><td>0.22</td><td>11%</td></tr><tr><td>0.50–0.56</td><td>0.41 ± 0.05</td><td>0.21</td><td>0.42 ± 0.02</td><td>0.22</td><td>4%</td></tr><tr><td>0.56–0.62</td><td>0.51 ± 0.05</td><td>0.22</td><td>0.47 ± 0.02</td><td>0.19</td><td>7%</td></tr></tbody></table>

**Note.** Column (1): rest-frame wavelength range. Column (2): mean bulge-over-total luminosity ratio and standard error (HST; Dimauro et al. 2018). Column (3): standard deviation. Column (4): mean bulge-over-total luminosity ratio and standard error (SHARDS). Column (5): standard deviation. Column (6): relative difference between $(B/T)_{HST}$ and $(B/T)_{SHARDS}$.
<figure><img src="image_17.png" /><figcaption>Figure 15. 2D photometric decomposition on HST (left panels) and SHARDS (middle panels) images, similarly to those proposed in Figures 3 and 4. The SEDs and the best-fit stellar models of the galaxy, bulge, and disk are also presented (right panels), similarly to those proposed in Figure 6. From top to bottom, the galaxies GDN 3360, GDN 9386, GDN 17242, and GDN 20441 are shown as an example.</figcaption></figure>

## Appendix B
### Mass-weighted Ages: Degeneracies

It is well known that stellar age, dust content, and metallicity are strongly degenerate, introducing uncertainties in physical parameter estimates based on SED fitting. However, key

spectral features in the optical rest-frame spectra of galaxies (i.e., the 4000 Å break and the Balmer break) better constrain the typical degeneracies that affect the SFH of each galaxy (Whitaker et al. 2011; Ferreras et al. 2012). In particular, since the break in the stellar continuum at 4000 Å is a good age
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Measured D4000 values as a function of mass-weighted ages for bulges (red circles) and pure spheroids (black triangles). The solid blue line shows the evolution of D4000 Bruzual & Charlot (2003) models at solar metallicity, while the dashed and dashed-dotted orange lines correspond to supersolar and subsolar metallicities, respectively.</figcaption></figure>

indicator (Bruzual 1983; Kauffmann et al. 2003), it has been successfully used to infer redshifts and ages of stellar populations in massive galaxies at high z (Saracco et al. 2005) and can also be successfully measured using medium-band photometry (Kriek et al. 2011). It is worth noting that, since we are dealing with estimations from photometric data alone, we decide to use the definition in Bruzual (1983) instead of the narrow index $D_n$4000 (Balogh et al. 1999), in order to reduce the uncertainties in the measurement. The SHARDS data provide a combination of depth and spectral resolution, which allows us to detect such features not only in the integrated SED of the galaxy but even in the bulge one (see Figure 6), providing an independent validation of its age. Thus, to assess the robustness of the mass-weighted ages retrieved from the SED fitting proposed in Sections 3.2 and 4.1, we show in Figure 16 our measurements of D4000 for the sample bulges and pure spheroids, comparing them with the trend expected from Bruzual & Charlot (2003) stellar population models.

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Distribution of mass-weighted age considering the 500 MC simulations for each bulge in our sample. The yellow histogram represents bulges in galaxies at redshift $z \le 0.75$, while the green histogram stands for bulges in galaxies at redshift $z > 0.75$.</figcaption></figure>

In Figure 17 we show the distribution of mass-weighted ages of the bulges in our sample, dividing them into those observed at redshift $z \le 0.75$ (46 out of 91) and those at redshift $z > 0.75$ (45 out of 91). We see that the trends presented in Figure 7 (and Figure 8) hold even when we separate bulges into lower- and higher-redshift ones. This allows us to rule out that the bimodality of bulge ages is biased because of their redshift distribution and/or large-scale structure (see Figure 1). Additionally, to further quantify that the age bimodality reported in Section 4.1 is not driven by degeneracies in the cluster of solutions resulting from the SED fitting, we show in Figure 18 the distribution of the 500 MC simulations performed for every bulge in our sample. In particular, we separate first-wave bulges ($z_{form} > 3$; left panel) from second-wave ones ($z_{form} < 3$; right panel). We show how the age bimodality holds independently of choosing the best cluster solution, implying a physical separation between the two waves of bulges.
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Density distribution of the mass-weighted ages of our bulges retrieved from the MC simulations described in Section 3.2. Bulges are separated between first-wave (purple distribution in the left panel) and second-wave ones (orange distribution in the right panel). Darker colors stand for higher density.</figcaption></figure>

## ORCID iDs

Luca Costantin https://orcid.org/0000-0001-6820-0015

Pablo G. Pérez-González https://orcid.org/0000-0003-4528-5639

Jairo Méndez-Abreu https://orcid.org/0000-0002-8766-2597

Marc Huertas-Company https://orcid.org/0000-0002-1416-8483

Paola Dimauro https://orcid.org/0000-0001-7399-2854

Belén Alcalde-Pampliega https://orcid.org/0000-0002-4140-0428

Fernando Buitrago https://orcid.org/0000-0002-2861-9812

Daniel Ceverino https://orcid.org/0000-0002-8680-248X

Emanuele Daddi https://orcid.org/0000-0002-3331-9590

Helena Domínguez-Sánchez https://orcid.org/0000-0002-9013-1316

Néstor Espino-Briones https://orcid.org/0000-0001-6426-3844

Antonio Hernán-Caballero https://orcid.org/0000-0002-4237-5500

Anton M. Koekemoer https://orcid.org/0000-0002-6610-2048

Giulia Rodighiero https://orcid.org/0000-0002-9415-2296

## References

Balogh, M. L., Morris, S. L., Yee, H. K. C., Carlberg, R. G., & Ellingson, E. 1999, ApJ, 527, 54

Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2013, ApJ, 765, 104

Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2014, ApJ, 791, 52

Barro, G., Pérez-González, P. G., Cava, A., et al. 2019, ApJS, 243, 22

Barro, G., Pérez-González, P. G., Gallego, J., et al. 2011, ApJS, 193, 30

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17

Bournaud, F. 2016, in Galactic Bulges, ed. E. Laurikainen, R. Peletier, & D. Gadotti (Berlin: Springer), 355

Bournaud, F., Dekel, A., Teyssier, R., et al. 2011, ApJL, 741, L33

Bruzual, A. G. 1983, ApJ, 273, 105

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Buitrago, F., Ferreras, I., Kelvin, L. S., et al. 2018, A&A, 619, A137

Buitrago, F., Trujillo, I., Conselice, C. J., et al. 2008, ApJL, 687, L61

Buitrago, F., Trujillo, I., Curtis-Lake, E., et al. 2017, MNRAS, 466, 4888

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Caon, N., Capaccioli, M., & D'Onofrio, M. 1993, MNRAS, 265, 1013

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2019, MNRAS, 490, 417

Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151

Ceverino, D., Dekel, A., Mandelker, N., et al. 2012, MNRAS, 420, 3490

Ceverino, D., Dekel, A., Tweed, D., & Primack, J. 2015, MNRAS, 447, 3291

Ceverino, D., Klessen, R. S., & Glover, S. C. O. 2018, MNRAS, 480, 4842

Chabrier, G. 2003, PASP, 115, 763

Charbonnier, A., Huertas-Company, M., Gonçalves, T. S., et al. 2017, MNRAS, 469, 4523

Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000, MNRAS, 319, 168

Costantin, L., Iovino, A., Zibetti, S., et al. 2019, A&A, 632, A9

Costantin, L., Méndez-Abreu, J., Corsini, E. M., et al. 2017, A&A, 601, A84

Costantin, L., Méndez-Abreu, J., Corsini, E. M., et al. 2020, ApJL, 889, L3

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680

Damjanov, I., McCarthy, P. J., Abraham, R. G., et al. 2009, ApJ, 695, 101

Danovich, M., Dekel, A., Hahn, O., Ceverino, D., & Primack, J. 2015, MNRAS, 449, 2087

de la Rosa, I. G., La Barbera, F., Ferreras, I., et al. 2016, MNRAS, 457, 1916

Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870

Dekel, A., Sari, R., & Ceverino, D. 2009, ApJ, 703, 785

Dimauro, P., Huertas-Company, M., Daddi, E., et al. 2018, MNRAS, 478, 5410

Djorgovski, S., & Davis, M. 1987, ApJ, 313, 59

Domínguez Sánchez, H., Pérez-González, P. G., Esquej, P., et al. 2016, MNRAS, 457, 3743

Dyda, S., Lovelace, R. V. E., Ustyugova, G. V., Romanova, M. M., & Koldoba, A. V. 2015, MNRAS, 446, 613

Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748

Estrada-Carpenter, V., Papovich, C., Momcheva, I., et al. 2020, ApJ, 898, 171

Ferreras, I., Pasquali, A., Khochfar, S., et al. 2012, AJ, 144, 47

Forbes, J. C., Krumholz, M. R., Burkert, A., & Dekel, A. 2014, MNRAS, 438, 1552

Franx, M., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2008, ApJ, 688, 770

Freeman, K. C. 1970, ApJ, 160, 811

Gadotti, D. A. 2009, MNRAS, 393, 1531

Gallazzi, A., Bell, E. F., Zibetti, S., Brinchmann, J., & Kelson, D. D. 2014, ApJ, 788, 72

Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41

Gelman, A., & Rubin, D. B. 1992, StaSc, 7, 457

Genzel, R., Burkert, A., Bouché, N., et al. 2008, ApJ, 687, 59

González Delgado, R. M., García-Benito, R., Pérez, E., et al. 2015, A&A, 581, A103

Graham, A. W. 2014, in ASP Conf. Ser. 480, Structure and Dynamics of Disk Galaxies, ed. M. S. Seigar & P. Treuthardt (San Francisco, CA: ASP), 185

Graham, A. W., Dullo, B. T., & Savorgnan, G. A. D. 2015, ApJ, 804, 32

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Häussler, B., McIntosh, D. H., Barden, M., et al. 2007, ApJS, 172, 615
Hernán-Caballero, A., Alonso-Herrero, A., Pérez-González, P. G., et al. 2013, MNRAS, 434, 2136

Hoffman, M. D., & Gelman, A. 2014, J. Mach. Learn. Res., 1593, 15

Hopkins, P. F., Bundy, K., Hernquist, L., Wuyts, S., & Cox, T. J. 2010, MNRAS, 401, 1099

Hopkins, P. F., Hernquist, L., Cox, T. J., Robertson, B., & Springel, V. 2006, ApJS, 163, 50

Hopkins, P. F., Somerville, R. S., Cox, T. J., et al. 2009, MNRAS, 397, 802

Hsu, L.-T., Lin, L., Dickinson, M., et al. 2019, ApJ, 871, 233

Huertas-Company, M., Bernardi, M., Pérez-González, P. G., et al. 2016, MNRAS, 462, 4495

Huertas-Company, M., Pérez-González, P. G., Mei, S., et al. 2015, ApJ, 809, 95

Immeli, A., Samland, M., Gerhard, O., & Westera, P. 2004a, A&A, 413, 547

Immeli, A., Samland, M., Westera, P., & Gerhard, O. 2004b, ApJ, 611, 20

Jog, C. J. 2014, AJ, 147, 132

Johnston, E. J., Häußler, B., Aragón-Salamanca, A., et al. 2017, MNRAS, 465, 2317

Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kriek, M., van Dokkum, P. G., Whitaker, K. E., et al. 2011, ApJ, 743, 168

Krumholz, M. R., & Thompson, T. A. 2013, MNRAS, 434, 2329

Larson, R. B. 1976, MNRAS, 176, 31

Lilly, S. J., Le Fevre, O., Hammer, F., & Crampton, D. 1996, ApJL, 460, L1

López Fernández, R., González Delgado, R. M., Pérez, E., et al. 2018, A&A, 615, A27

MacArthur, L. A., Courteau, S., & Holtzman, J. A. 2003, ApJ, 582, 689

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Madau, P., Ferguson, H. C., Dickinson, M. E., et al. 1996, MNRAS, 283, 1388

Mancini, C., Daddi, E., Juneau, S., et al. 2019, MNRAS, 489, 1265

Martig, M., Bournaud, F., Teyssier, R., & Dekel, A. 2009, ApJ, 707, 250

McDermid, R. M., Alatalo, K., Blitz, L., et al. 2015, MNRAS, 448, 3484

Méndez-Abreu, J., Aguerri, J. A. L., Corsini, E. M., & Simonneau, E. 2008, A&A, 478, 353

Méndez-Abreu, J., de Lorenzo-Cáceres, A., & Sánchez, S. F. 2021, MNRAS, 504, 3058

Méndez-Abreu, J., Debattista, V. P., Corsini, E. M., & Aguerri, J. A. L. 2014, A&A, 572, A25

Méndez-Abreu, J., Ruiz-Lara, T., Sánchez-Menguiano, L., et al. 2017, A&A, 598, A32

Méndez-Abreu, J., Sánchez, S. F., & de Lorenzo-Cáceres, A. 2019a, MNRAS, 484, 4298

Méndez-Abreu, J., Sánchez, S. F., & de Lorenzo-Cáceres, A. 2019b, MNRAS, 488, L80

Mihos, J. C., & Hernquist, L. 1996, ApJ, 464, 641

Morelli, L., Parmiggiani, M., Corsini, E. M., et al. 2016, MNRAS, 463, 4396

Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJL, 699, L178

Noguchi, M. 1999, ApJ, 514, 77

Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312

Pacifici, C., Kassin, S. A., Weiner, B., Charlot, S., & Gardner, J. P. 2013, ApJL, 762, L15

Pacifici, C., Kassin, S. A., Weiner, B. J., et al. 2016, ApJ, 832, 79

Pérez-González, P. G., Cava, A., Barro, G., et al. 2013, ApJ, 762, 46

Pérez-González, P. G., Gil de Paz, A., Zamorano, J., et al. 2003, MNRAS, 338, 525

Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234

Porter, L. A., Somerville, R. S., Primack, J. R., et al. 2014, MNRAS, 445, 3092

Renaud, F., Bournaud, F., Kraljic, K., & Duc, P. A. 2014, MNRAS, 442, L33

Rogers, B., Ferreras, I., Peletier, R., & Silk, J. 2010, MNRAS, 402, 447

Salvatier, J., Wieckiå, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, e55

Sánchez-Blázquez, P., Jablonka, P., Noll, S., et al. 2009, A&A, 499, 47

Saracco, P., Longhetti, M., Severgnini, P., et al. 2005, MNRAS, 357, L40

Schawinski, K., Urry, C. M., Simmons, B. D., et al. 2014, MNRAS, 440, 889

Schiavon, R. P., Faber, S. M., Konidaris, N., et al. 2006, ApJL, 651, L93

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, A&A, 618, A85

Scott, N., Brough, S., Croom, S. M., et al. 2017, MNRAS, 472, 2833

Sérsic, J. L. 1968, Atlas de Galaxias Australes (Cordoba: Observatorio Astronomico de Cordoba)

Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2014, ApJL, 783, L14

Tabor, M., Merrifield, M., Aragón-Salamanca, A., et al. 2017, MNRAS, 466, 2024

Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 458, 242

Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673

Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775

Tortora, C., La Barbera, F., Napolitano, N. R., et al. 2016, MNRAS, 457, 2845

Trujillo, I., Cenarro, A. J., de Lorenzo-Cáceres, A., et al. 2009, ApJL, 692, L118

Trujillo, I., Feulner, G., Goranova, Y., et al. 2006a, MNRAS, 373, L36

Trujillo, I., Förster Schreiber, N. M., Rudnick, G., et al. 2006b, ApJ, 650, 18

van der Wel, A., Bell, E. F., Häussler, B., et al. 2012, ApJS, 203, 24

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28

van Dokkum, P. G., Bezanson, R., van der Wel, A., et al. 2014, ApJ, 791, 45

Wellons, S., Torrey, P., Ma, C.-P., et al. 2015, MNRAS, 449, 361

Whitaker, K. E., Kriek, M., van Dokkum, P. G., et al. 2012, ApJ, 745, 179

Whitaker, K. E., Labbé, I., van Dokkum, P. G., et al. 2011, ApJ, 735, 86

Wild, V., Almaini, O., Dunlop, J., et al. 2016, MNRAS, 463, 832

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Wuys, S., Labbé, I., Franx, M., et al. 2007, ApJ, 655, 51

Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327