© The Authors 2023

1. Introduction

Constraining the processes that drive galaxy evolution in different cosmic epochs is still a work in progress. At higher redshifts, external processes, such as mergers, galaxy interactions, and gas inflows, have an important impact on how galaxies evolve (e.g., Schreiber et al. 2006; Genzel et al. 2008; Law et al. 2009; Dekel et al. 2009; Oser et al. 2010). As the Universe expands and interactions are less frequent, internal processes begin to play an important role in the evolution of galaxies (Kormendy & Kennicutt 2004).

Among the most important internal drivers of the evolution of disc galaxies are stellar bars, which efficiently redistribute angular momentum, as well as both stars and gas (e.g., Lynden-Bell & Kalnajs 1972; Combes & Gerin 1985; Athanassoula 2003; Munoz-Tunón et al. 2004; Sheth et al. 2005; Romero-Gómez et al. 2007; Di Matteo et al. 2013; Halle et al. 2015; Fragkoudi et al. 2016, 2017). Bars have been linked to a global quenching of star formation in galaxies (e.g., Masters et al. 2012; Schawinski et al. 2014; Haywood et al. 2016; Géron et al. 2021) and are also responsible for inducing bursts of star formation in central regions (e.g., Ishizuki et al. 1990; Ellison et al. 2011; Coelho & Gadotti 2011). Bars are common structures in the local Universe: many studies have found that between ∼30% and 70% of disc galaxies host bars (e.g., Eskridge et al. 2000; Menéndez-Delmestre et al. 2007; Barazza et al. 2008; Sheth et al. 2008; Aguerri et al. 2009; Nair & Abraham 2010; Buta et al. 2015; Erwin 2018). Therefore, understanding the time of their formation and how they affect their host galaxies is key to understanding the late-stage evolution of galaxies themselves.

By investigating how galaxies evolve dynamically at different redshifts, studies find that a large fraction of galaxies at high redshift are rotationally supported (e.g., Shapiro et al. 2008; Schreiber et al. 2009; Epinat et al. 2012; Wisnioski et al. 2015; Rizzo et al. 2020; Lelli et al. 2021), with fractions varying from 70 − 90% at z ∼ 1 and 47 − 74% at z ∼ 2. Using the James Webb Spatial Telescope. Early Release Observations of the galaxy cluster SMACS 0723, Ferreira et al. (2022) show that the fraction of disc galaxies in the early Universe is still an open question. They find that disc galaxies dominate the morphology at z ∼ 1.5 by more than a factor of ∼10 higher compared with previous results from the Hubble Space Telescope. Nevertheless, these rotationally supported discs at higher redshifts often present higher velocity dispersions, suggesting they are turbulent, unsettled, and thick (e.g., Elmegreen & Elmegreen 2006; Cresci et al. 2009; Newman et al. 2013; although see also Rizzo et al. 2020). Since analytical and numerical work indicates that the bar can only form once the disc is dynamically settled, the moment of bar formation marks this transition epoch, during which the discs are at least partially dynamically cold (e.g., Kraljic et al. 2012 and references therein). Exactly how and when the switch from external to internal processes as the driver of galaxy evolution happens is still not clear, but it is associated with the settling of galactic discs. Therefore, timing the epoch of bar formation is a major step forwards in piecing together the different phases of galaxy evolution.

Different studies have attempted to time when the bar formed for a handful of galaxies using different approaches. Gadotti et al. (2015) analysed stellar populations in the nuclear disc (ND) of NGC 4371 using high quality Multi Unit Spectroscopic Explorer (MUSE) data, and found a bar age of about 10 Gyrs. Pérez et al. (2017) used the formation of the boxy/peanut-shaped bulge, which originated from instabilities of the bar, to estimate the age of the bar in NGC 6032 as 10 Gyrs. de Lorenzo-Cáceres et al. (2019) analysed the star formation histories (SFHs) of nuclear structures, including inner bars, and find that the inner bars in NGC 1291 and NGC 5850 must have formed at least 6.5 Gyrs and 4.5 Gyrs ago, respectively. These studies provide observational evidence that bars can be long-lived, which is in agreement with cosmological simulations that find bars that formed between redshifts 1 and 2 and survive down to z = 0 (e.g., Kraljic et al. 2012; Rosas-Guevara et al. 2020; Fragkoudi et al. 2020, 2021).

One of the immediate effects that follow bar formation – and the ensuing onset of tangential forces in the disc – is the gas inflow to the central parts of the galaxy. This gas inflow builds the rotation-supported structures known as NDs, which are commonly referred to as pseudo-bulges or disc-like bulges. Simulations and theoretical studies find that, after the bar forms, it only takes ∼10⁸ yr to form the ND (e.g., Athanassoula 1992a,b; Emsellem et al. 2015; Seo et al. 2019; Baba & Kawata 2020). Thus, one can use the time of the formation of the ND – obtained through galactic archaeology – to time the age of the bar (Gadotti et al. 2015). The study by Gadotti et al. (2015) was offered as a proof of concept, using NGC 4371, to demonstrate the feasibility of timing bar formation with archaeological evidence from the ND itself. However, the ND in NGC 4371 has no star formation and is mainly dominated by old stars, which makes it a special case for this type of analysis. More generally, with more complex SFHs, estimating stellar population properties in the central region of disc galaxies is not a trivial task, since the observed light carries tangled information of the ND with the underlying main disc (MD) – and possibly other structural components that were already present when the ND formed. Therefore, in order to reliably detect the oldest stars in the ND, we need a way to disentangle the light of the ND from that of the underlying MD.

This is precisely the goal of this work, in which we present a new methodology for disentangling the light of the ND from other structures in the central region of the galaxy. This in turn enables us to derive the bar formation epoch. With the disentangled light, we derive independent SFHs for different stellar structures in order to time the moment the ND formed and, therefore, when the bar formed. We selected NGC 1433, one of the galaxies in the Time Inference with MUSE in Extragalactic Rings (TIMER) survey (Gadotti et al. 2019), as a pilot galaxy to present this methodology, which will in the future be applied to all TIMER galaxies.

This paper is organised as follows. In Sect. 2 we describe the data and explain the reason we selected NGC 1433 for this study. In Sect. 3 we describe the methodology that we developed to disentangle the SFH of the ND from that of the underlying MD, and how this can be used to obtain the bar age, as well as the tests of the methodology carried out using hydrodynamic simulations. In Sect. 4 we present our results for the age of the bar in NGC 1433 and the detailed buildup of its ND. In Sect. 5 we discuss the implications of our results in the context of galaxy evolution and the formation of NDs. We summarise and conclude in Sect. 6. Additional tests of the methodology, including on two control galaxies, are presented in the appendix.

2. Sample and data description

In order to present the methodology we have developed – which separates the light of NDs from the underlying population – we select NGC 1433 for this pilot study, since it has a ND with similar properties to most NDs in the TIMER sample (e.g., Gadotti et al. 2020; Bittner et al. 2020 – which will be future targets of study). Also, we ensure that our methodology does not produce artificial bar ages for unbarred galaxies by applying it to two control bar-less galaxies: NGC 1380 and NGC 1084 (see Appendix A).

Buta et al. (2015) classified NGC 1433 morphologically as strongly barred with a nuclear ring/lens and a nuclear bar. In addition, the galaxy is at a distance of 10 Mpc, has an inclination of 34°, and stellar mass of 2 × 10¹⁰ M_⊙ (see references in Gadotti et al. 2019). Following Gadotti et al. (2020), the radius of the ND is defined as the peak in the v/σ radial profile, which is 440 pc (see Fig. 3). We used data from the TIMER project: the observations were carried out using the MUSE instrument at the European Southern Observatory Very Large Telescope (ESO-VLT), in Period 97, from 2016 March to October, in Wide Field Mode. Considering the point spread function full width at half maximum, the spatial resolution of these observations is about 50 pc. Further details on the observations and data reduction can be found in Gadotti et al. (2019).

In Fig. 1 we display colour composites of our target along with maps of stellar kinematics and population properties: stellar velocity and velocity dispersion, the higher-order moments of the line of sight velocity distribution h₃ and h₄, stellar age, metallicity, and α enhancement, which were produced applying the GIST pipeline (Bittner et al. 2019) to the TIMER MUSE data (more details on the following section). One can see that NGC 1433 has clear signatures of a ND, that is, a rapidly rotating structure detached from the MD and coincident with a drop in velocity dispersion and mean ages. In addition, it also shows an anti-correlation between velocity and h₃, and an increase in h₄, as expected for NDs. These results are all in agreement with the findings in Gadotti et al. (2020) and Bittner et al. (2020), respectively.

Fig. 1. NGC1433 data and derived maps: The top two images are the colour composites of NGC 1433 from the Carnegie-Irvine Galaxy Survey (top left; Ho et al. 2011) and TIMER (top right; Gadotti et al. 2019). Below are maps that show (from left to right and top to bottom) the stellar kinematics and population properties derived from the TIMER data using GIST: stellar velocity, velocity dispersion, h3, h4, age, metallicity, and α enhancement. The ND radius is shown with a solid black contour and displays a faster rotation, a drop in velocity dispersion, an anti-correlation between velocity and h3, and an increase in h4. In addition, the same region corresponds to a drop in mean age and α enhancement and an increase in metallicity. It is clear that NGC 1433 hosts a younger ND with more rotational support in the central region than the underlying population, in agreement with Gadotti et al. (2020) and Bittner et al. (2020).

3. Methodology

In this section we describe our methodology step-by-step, as illustrated in Fig. 2. We describe how we create the underlying stellar population contribution (Sect. 3.1) and how we isolate the SFH of the ND (Sect. 3.2). In Sect. 3.3 we describe the details in our data analysis and our criterion for obtaining the time of bar formation. In Sect. 3.4 we test our methodology by applying it to a simulated barred galaxy.

Fig. 2. Illustration of the methodology described in Sect. 3 for a galaxy that hosts a ND. In the top left we highlight the field of view from MUSE in the centre of the galaxy and the position of the representative region just outside the ND. The position of the representative ring is chosen taking the v/σ radial profile into consideration. From it we derive the representative spectrum used to build the MD. In the bottom left we display the output from the light disentangling: the original data cube (red), the ND data cube (blue), and the representative MD data cube (green). All three data cubes have spaxels corrected to v = 0 km s−1 and σ = σmax km s−1. Lastly, we collapse each data cube into a mean spectrum (as illustrated) and derive SFHs for each one. The steps in our methodology are described to the right.

3.1. Building the underlying population contribution and disentangling the nuclear disc light

Recent studies deriving radial profiles of stellar ages found that galaxies hosting NDs display a drop in average ages towards the centre. This drop in average ages is coincident with the radius of the peak in stellar v/σ (e.g., Falcón-Barroso et al. 2002; Bittner et al. 2020), which implies that these galaxies have a central stellar structure that is younger than the main underlying population. This agrees with the scenario whereby NDs are structures formed by a relatively late gas inflow induced by stellar bars (e.g., Gadotti et al. 2015, 2019; Bittner et al. 2020). Hereafter, we consider the ‘underlying main disc’ to be every stellar population that was present before the formation of the ND. Since the observed light carries combined information of the younger ND and other underlying central structures, the ages derived for the central region of the galaxy can only be considered as an upper limit to the mean stellar age of the ND.

In this section we describe a new strategy to build the underlying MD and disentangle the light of it from that of the ND. In summary, this methodology consists of deriving a spectrum that represents the main underlying disc (hereafter referred to as the ‘representative spectrum’), using it to build a representative MD data cube, and later subtract it from the observed data. We make two hypothesis: (1) the MD extend all the way to the centre, and its surface brightness profile follows an exponential law and (2) the stellar populations of the MD do not change significantly in the inner regions and, therefore, we can assume is similar to the spectrum extracted from an aperture surrounding the nuclear ring. The result from this subtraction is expected to be the isolated ND. A step-by-step outline and an illustration of our methodology are shown in Fig. 2.

Before building the MD, we treated the original data cube, taking the kinematic properties into consideration. We first used the data analysis pipeline (DAP) from the PHANGS-MUSE survey (Emsellem et al. 2022) for the original data cube, deriving kinematic and line emission properties. DAP is a module-based pipeline based on the Galaxy IFU Spectroscopy Tool (GIST; Bittner et al. 2019) able to extract properties from data cubes such as kinematic information, emission lines fluxes and more. For the emission lines fitting, DAP uses pPXF (Cappellari 2012), considering emission lines as extra gaussian elements added to the stellar continuum. We measured the kinematic maps for Voronoi-binned (Cappellari & Copin 2003) data with signal to noise of 100 and emission lines spaxel-by-spaxel, between the wavelengths 4800 and 7000 Å, employing the E-MILES simple stellar population (SSP) model library (Vazdekis et al. 2016) to remove the stellar continuum. From the results, we have kinematic properties such as stellar velocity and velocity dispersion and emission line fluxes. Applying the kinematic information, we shift all spectra to have a final velocity of v = 0 km s⁻¹, accounting for both the galaxy recession velocity and internal kinematics. In addition, we also convolved the original data to ensure that every spectra has σ = σ_max, where the latter is the highest velocity dispersion in the central region of the original data cube. For NGC 1433, σ_max is 122 km s⁻¹. We point out that we verified that the convolved data and the original data result in the same stellar population, so such a procedure does not affect our results. Our aim in doing so is to ensure we do not create artificial wings or artificial emission/absorption lines when subtracting the spectra (see below). Once the original data were shifted to v = 0 km s⁻¹ and σ = σ_max km s⁻¹, we built our underlying MD.

The first step is to define the representative region from which we will extract the representative spectrum. This region is chosen as a ring surrounding the ND – from now on referred to as the ‘representative ring’. The representative ring is expected to be as close as possible to the ND without being contaminated by its light. For galaxies hosting a ND, such as NGC 1433, we take into consideration the stellar v/σ radial profile (Fig. 3). We define the inner radius of the representative ring as the first minimum point outside the ND. That is the radius outside of which we expect the light from the MD to start dominating and, as expected, it is also the radius with the oldest mean age (see Fig. 3). Nevertheless, we tested different positions for the representative ring to assess the systematic error linked to this decision (see Appendix B.1). In addition, the representative ring has a width of 2″, so one does not expect strong age and flux gradients inside of it. Nonetheless, we performed tests with widths between 1″ and 4″ with little variation in the outcome.

Fig. 3. v/σ radial profile, displayed in dashed red contours (values in the left-axis), and light-weighted average age, in solid blue contours (values in the right axis) for NGC 1433. We display the median values together with the first and third quartiles. The vertical dotted black line marks the ND radius and the hatched area the representative ring. Note that the representative ring is placed in the first v/σ minimum outside the ND, which matches the oldest mean age.

Secondly, we aim to mask spaxels dominated by active galactic nucleus (AGN) emission, based on the Baldwin, Phillips & Tervelich (BPT) diagram classification (Baldwin et al. 1981), with AON ≥ 20. We used the emission lines Hα, Hβ, [OIII], and [NII] extracted by DAP, as described above, to build the BPT diagram in Fig. 4.

Fig. 4. NGC1433 BPT classification: Left: BPT diagram (Baldwin et al. 1981) classification of each spaxel for the ND in NGC 1433. The BPT classification is continued in the middle and right panels, which display the physical position of the different spaxels and the light radial profile, respectively. We also display the radius of the inner mask dominated by AGN contributions (dotted contours), the radius of the ND (solid contours), and the representative ring (dashed contours).

In the third step, we normalise the flux in each spaxel of the representative ring to account for the radial flux increase in the MD towards the centre. We normalised each spectra to r = 0 pc, assuming that the disc light profile follows the exponential function

$$\begin{array}{c}\hfill I(r)=I_0\times e^{-r/h},\end{array}$$(1)

where h is the disc scale-length, r is the distance of a given spaxel from the centre (corrected for inclination effects), and I₀ is the flux at r = 0 pc. The value of h taken for NGC 1433 comes from Salo et al. (2015) and is 3100 pc. Considering we have I(r) for each spaxel in the representative ring, we divided it by e^−r/h, extrapolating the observed flux to the centre of the galaxy. We derived the representative spectrum as the mean flux per wavelength of all spectra from the representative ring with v = 0 km s⁻¹, σ = σ_max km s⁻¹ and the corresponding flux at r = 0 pc. Finally, to re-construct a data cube of the MD, we extrapolated the representative spectrum back to a range of radii, again taking an exponential light profile into account for the MD. Although the MD is often described with an exponential light profile, recent studies show that discs may not follow an exponential light profile all the way to the centre (e.g., Zhu et al. 2018; Breda et al. 2020; Papaderos et al. 2022). In light of this, we tested the implications of assuming an exponential profile in Appendix B.2, by applying the same methodology to a flat light profile MD. We show that the choice of a profile for the main underlying population does not affect the results noticeably.

Next, we subtracted the reconstructed MD from the original data. As mentioned, in order to prevent the creation of artificial wings from the subtraction, we also used the original data convolved and shifted to v = 0 km s⁻¹ and σ = σ_max km s⁻¹. With that, we can disentangle our original data into the MD and the ND data cubes, as exemplified in Fig. 5. Lastly, we collapse each data cube into a single mean spectra, deriving SFHs for each one. To assure that the collapsed spectra is the correct description of the central region for the three data cubes, we masked contributions from the AGN in the centre. Using the BPT classification (with AON ≥ 20) together with the light radial profile (Fig. 4, right), we delimited a central region to be masked of around 140 pc (15 spaxels), since it is dominated by AGN emission and can contaminate the total flux of the ND.

Fig. 5. Illustration of different data products derived from the methodology described in Sect. 3. We display the sum of the original MUSE fluxes between 4800 and 5800 Å (top row), and the derived representative MD (bottom row, left) and the ND data (bottom row, right), which is the result of subtracting the representative MD from the original data cube. For the representative MD, it is possible to notice the exponential increase in flux towards the centre. All data cubes were masked for AGN contributions with AON ≥ 20.

Due to the high quality of the TIMER data (Gadotti et al. 2019), the decision to collapse the data cube is not motivated by the increase of signal-to-noise ratio. Instead, we collapse the data cube to save computational time and guarantee low statistical errors. Nevertheless, we applied the same methodology for the non-collapsed data cube and the bar age is not strongly affected (see Appendix B.3).

3.2. Deriving star formation histories and mass assembly

For each collapsed spectrum – MUSE original, MD, and clean ND data – we ran the GIST pipeline (Bittner et al. 2019) to derive stellar population properties and SFHs. To guarantee consistency in our analysis with previous TIMER work, we used the same GIST configuration as in Bittner et al. (2020). Firstly, GIST employs an un-regularised run of pPXF (Cappellari & Emsellem 2004; Cappellari 2017) to derive stellar kinematic properties. We also include a low-order multiplicative Legendre polynomial in the fit to account for small differences between the shape of the continuum templates and the observed spectra. Next, GIST employs pyGandALF (see Bittner et al. 2019) to model emission lines as Gaussians, simultaneously fitting the stellar continuum while the stellar kinematics remains fixed from the previous step. With this, we obtain the emission-subtracted spectra. Lastly, GIST performs a regularised pPXF run in the emission-subtracted spectra, in order to fit a combination of stellar populations and derive mean properties. Since metallicity and stellar velocity dispersion can both be responsible for absorption line broadening, causing possible degeneracies (e.g., Sanchez-Blazquez et al. 2011), we keep the stellar kinematics fixed from previous steps. In addition, to account for extinction and continuum mismatch effects, we apply an eighth-order multiplicative Legendre polynomial in the fit. For the last step, we employ the MILES library (Vazdekis et al. 2015), light-weighted, with [M/Fe] between −1 and +0.4, ages in the range 0.03–14 Gyr, and [α/Fe] enhancements of +0.0 and +0.4. Lastly, we use the regularisation error value of 0.15 derived for TIMER data by Bittner et al. (2020). As described in Cappellari (2017), the regularisation of the SFHs allows one to derive the smoothest and still physically meaningful result. To assess how much our final bar age relies on the regularisation error, we tested different values in Appendix B.4. For further details on the data analyses, we refer the reader to previous TIMER papers (e.g., Bittner et al. 2020).

pPXF also estimates different weights for SSPs with different ages, allowing us to build a non-parametric SFH. Each weight represents the fraction of the light formed in the different SSPs. In order to convert the light-weighted SFHs to mass-weighted SFHs, we consider the distance to the galaxy to derive the intrinsic luminosity. We then use the mass-to-light ratios¹ predicted from the BaSTI isochrones (Pietrinferni et al. 2004, 2006, 2009, 2013) to convert luminosity into stellar mass. The mass-to-light ratios assume a Kroupa revised initial mass function (Kroupa 2001), and the MILES template library (Vazdekis et al. 2015) with [α/Fe] enhancements of +0.0 and +0.4. These mass-to-light ratios depend on the age, metallicity and [α/Fe] enhancement that best describe the observed data, and account for both stellar and remnants masses.

3.3. Age-dating bar formation

Considering the scenario in which the ND is formed from a late gas inflow due to the bar, in an ideal case, one could expect the original SFH characterised by an older star formation event followed by a later burst. Once we are able to disentangle the MD and the ND, we can expect the former’s SFH to map the oldest star formation event and the latter’s to map the youngest burst.

If we could perfectly disentangle the light of the ND from the underlying MD, one could simply use the first peak in the ND SFH as the formation time of the ND and, consequently, the bar. However, due to gradients in the stellar population properties, the region around the ND that we use to obtain the SFH of the MD might not be identical to the real SFH of the main underlying disc within the ND region. Thus, it might not fully remove the contamination of the MD from the ND light. With that in mind, and by testing our methodology on N-body+hydrodynamic simulations (see Sect. 3.4), we employed a criterion to time the bar formation epoch as the moment when the ND dominates the star formation, as a signature of the bar bringing gas towards the centre. That corresponds to the first time in which the ratio between the star formation in the ND and that in the MD, ND/MD, rises above 1, with a positive slope towards younger ages. In order to verify whether this is a reliable criterion, we test our methodology using an N-body+hydrodynamic simulation of a barred galaxy below.

3.4. Testing the method using hydrodynamic simulations

In order to test the robustness of the methodology developed here and some of the assumptions employed, we used an N-body and hydrodynamic simulation of an isolated Milky-Way-like disc galaxy, which forms a bar and a ND self-consistently (the simulation is part of a suite of models developed to study the evolution of barred galaxies; Fragkoudi & Bieri, in prep.). The simulation has two collisionless components (a stellar disc and a dark matter halo) and a collisional component (gaseous disc), which is able to form stars that subsequently return mass, energy, and metals to the interstellar medium via supernova feedback. We refer the reader to Appendix C for technical details of the simulation.

The simulation is evolved for a total of 3.3 Gyr. The axisymmetric stellar disc that is in place at the start of the simulation (which we refer to as the ‘old stellar component’) rapidly forms a bar after ∼0.3 Gyr (we define the bar as being fully formed when the m = 2 Fourier mode of the surface density, A₂ > 0.3). During and after the formation of the bar, gas piles up at the leading edges of the bar, where it shocks, loses angular momentum, and is funnelled to the centre, where it forms a dense gaseous ND (see panel (a) of Fig. 6), which proceeds to form stars. These new stars – formed out of gas pushed to the centre by the bar – form a highly rotating stellar ND (see panel (b) of Fig. 6), whose size is set by the bar orbits in the inner regions (see Athanassoula 1992b), similar to those observed in the local barred galaxies (e.g., Gadotti et al. 2020).

Fig. 6. Testing the methodology using hydrodynamic simulations. a: face-on projection of the gas temperature. The dust lanes on the leading edge of the bar and the gaseous ND are clearly visible as low temperature regions. b: kinematic map showing the velocity of stars in the y direction with stellar isodensity contours over plotted, which outline the shape of the bar (note the different scale compared to panel a). A highly rotating stellar component in the central kiloparsec (i.e. the ND) is clearly evident in the kinematic map. c: age map of all the stars in the simulated galaxy. The grey box shows the inset that is represented in panel d. d: zoomed-in view around the ND region. The inner and outer dashed lines denote the region used to obtain the SFH of the ‘representative disc region’. The solid black line denotes the radius within which the SFH of the ND is estimated. The scatter points indicate the locations of the pixels used for deriving the SFHs in panel e. The points are colour-coded by the mean age in the pixel. e: the top panel shows the SFHs of the original ND region (solid red), of the representative SFH of the MD (dot-dashed green), and of the ND with the representative SFH subtracted (dashed blue). The bottom panel shows the ratio of the subtracted SFH to the representative SFH. The vertical orange line indicates the time of bar formation in the simulation (3 Gyr), and the vertical dashed line indicates the time at which ND/MD is above 1 (3.3 Gyr).

We highlight that, while the new stars formed out of the gas have self-consistent ages – according to when they are formed in the simulation, the old stellar component can have any SFH we assign to it. In order to model the age gradient often found in galaxies, we assigned ages to the old stellar component at the initial snapshot such that a negative age gradient was obtained² (i.e. with older stars in the centre and younger stars at the edge of the disc). This old underlying population, together with the new stars formed out of the gas in the simulation, give rise to the age map of the galaxy at the end of the simulation shown in panel (c) of Fig. 6.

We can now extract the SFH for a given ‘pixel’ in the ND region – as we do in the observations – which will contain stars born from gas pushed to the centre by the bar, as well as old stars that were present before the bar formed. We also extract the SFH for pixels in the region just outside the ND, which give us the representative SFH of the underlying main stellar disc (MD). Therefore, as in the methodology used for the observations, we can extract the SFH of both the ND region and the MD. We then subtract the SFH of the MD from the SFH in the ND region, in order to obtain the SFH of the ‘clean’ ND itself. These SFHs are shown in panel (e) of Fig. 6, with solid red for the original total SFH within the ND region, dot-dashed green for the SFH of the main underlying disc, and dashed blue lines for the clean ND SFH. This can be compared to the ‘true’ SFHs of the ND and MD, which are shown in the left panel of Fig. C.1.

We find that due to the age gradient in the underlying disc as well as the gradient in the light profile, there can be a contamination of the ‘old component’ in the subtracted ND SFH, which cannot be fully removed by subtracting the underlying MD (see Fig. C.1, where we show how assuming a negative age gradient versus a flat age gradient affects the methodology). The extent to which there is ‘contamination’ by the oldest stars in the ND will depend on both how steep the age gradient is, and on the location of the ring used to obtain the SFH of the MD (see also Fig. C.2). This indicates that we cannot simply use the oldest peak in the subtracted SFH (ND) in order to obtain the time at which the ND formed, but we rather should use the comparison between the subtracted ND and the representative SFH of the MD (i.e. ND/MD). In practice, this is the time when the SFH of the ND increases above that of the MD (i.e. when ND/MD rises above 1), with a positive slope towards younger ages.

As can be seen from panel (e) of Fig. 6, this method allows us to recover the time at which the bar formed (which in the simulation occurs at t_lookback ∼ 3 Gyr), which is marked with the vertical solid orange line, while the ratio of the ND over the representative MD SFH (ND/MD) gives a bar age of t_lookback = 3.3 Gyr. Therefore, we find that, even with the contamination of older ages in the subtracted spectrum, we can recover the time of bar formation, as the first time at which ND/MD rises above one, with an accuracy that will depend on the width of the age bins in the SFHs. The main limitation in obtaining the bar age using this methodology therefore stems from uncertainties in deriving stellar ages – and therefore the SFH – which are typically of the order of ∼1 Gyr (see e.g., Bittner et al. 2020).

4. Results

After following the methodology to disentangle SFHs described in Sect. 3, we present our results for our pilot study galaxy, NGC 1433. We remind the reader that the chosen criterion to time bar formation epoch is the first moment when the ratio ND/MD increases above 1 with positive slope towards younger ages.

Figure 7 shows our main results: the mass of stars formed with different age bins of the stellar templates – analogous to SFHs – for the MUSE original (solid-red line and contour), the MD (dot dashed green line) and the ND (dashed blue line and contour) for NGC 1433, together with the ratio ND/MD in the bottom panel. The age is related to the SSP template combination that best fits the observed spectra and can be understood as a ‘look-back time’ (i.e. t = 0 is the present time).

Fig. 7. Mass-weighted SFHs for NGC 1433 from collapsed spectra of the MUSE original data, ND data, and representative MD data, as illustrated in Fig. 2. The y-axis corresponds to the mass formed, in solar masses, for each age bin and the x-axis to the age in Gyr, that is, the look-back time. The original data are shown with a solid red line, the subtracted data with a dashed blue line, and the representative spectrum with the dot-dashed green line. The lines are the SFHs derived from the data, while the shaded regions are results from 100 Monte Carlo runs. In the bottom panel we display ND/MD in black, and highlighted with an orange circle is the age where ND/MD rises above one for the first time, together with the statistical uncertainty: ${}_{-0.5}^{+0.2}(\mathrm{stat})$ Gyr.

We measure the bar formation epoch of NGC 1433 to occur $7.5_{-1.1}^{+1.6}{(\mathrm{sys})}_{-0.5}^{+0.2}(\mathrm{stat})$ Gyrs ago, corresponding to a redshift of z ≈ 1. In order to quantify the statistical error of the methodology, we perform 100 Monte Carlo runs for each of the collapsed spectra (original data, ND data and representative spectra) to derive variations on the bar age. We use the noise information to sample a distribution of fluxes for each wavelength, creating 100 artificial spectra. We then run pPXF on each of these spectra, to obtain the different SFHs and the subsequent bar ages. This is shown in Fig. 7 as the shaded area of each SFH. From it, we derive a statistical uncertainty corresponding to ${}_{-0.5}^{+0.2}(\mathrm{stat})$ Gyr. The statistical error is subdominant, since, by collapsing the data cubes into single spectra, we achieve signal-to-noise values over 2000. To further quantify uncertainties in the derived bar age, which can be introduced due to various aspects of the methodology, we perform multiple tests with different configurations (see Appendix B): different locations for the representative ring, different light profiles to describe the increase in density of the main underlying population towards the centre, and different regularisation errors for the pPXF run that results in the derivation of the SFHs. From these tests, we find that there is a systematic uncertainty in our measurements of the bar age of the order of ${}_{-1.1}^{+1.6}(\mathrm{sys})$ Gyr, which we quote in addition to our statistical errors.

From the SFHs, we can derive estimates of the total stellar mass in the isolated ND and the underlying MD within the ND radius, by summing the mass formed through time following each curve. However, if NGC 1433 has significant age gradients within the central kiloparsec, the very oldest population in the underlying disc may still be partly present in the isolated ND, as discussed above. Therefore, these would be, respectively, an upper limit to the mass of the ND, and a lower limit to the mass of the underlying population within the ND region. For NGC 1433, we measure the underlying MD mass within the ND radius as 2.95 × 10⁸ M_⊙ and the ‘cleaned’ ND mass as 4.05 × 10⁸ M_⊙. To explore whether these values are in agreement with the literature, we extrapolate the underlying MD mass to obtain the total mass of the galaxy assuming an exponential function, following Eq. (2):

$$\begin{array}{c}\hfill M_h=2\pi {\displaystyle {\int }_0^{\infty }\mathrm{\Sigma }(r)r\mathrm{d}r=2\pi {\mathrm{\Sigma }}_0{h}^2,}\end{array}$$(2)

where Σ₀ is the mass density at the centre and h is the disc scale-length. Using the mass of the underlying MD within the ND radius, we measure the mass density as 512 M_⊙ pc⁻². This value gives the extrapolated mass for the entire galaxy of 2.74 × 10¹⁰ M_⊙, which is consistent with the total stellar mass of 2 × 10¹⁰ M_⊙ derived by Muñoz-Mateos et al. (2015), considering the uncertainties involved. This indicates that our measurements for the total stellar mass of the ND and the underlying main population are reliable. In addition, this also shows that the methodology described above to disentangle the light of the nuclear and underlying discs is trustworthy. Interestingly, our mass estimates indicate that the ND dominates the stellar mass budget in the central region, with the ND being ∼40% more massive than the underlying MD in the same region.

In addition, we analyse the SFH in different radial bins inside the ND region for the original data (left panel) and the ND cleaned data (right panel; Fig. 8). At the top of each panel we display the mean age for each radius, colour-coded according to distance from the centre. As one can see for the cleaned ND results, the SFH and the mean ages gradually get younger at larger radii, in agreement with the inside-out growth scenario (Bittner et al. 2020). We discuss the implications of these findings further in Sect. 5.2.

Fig. 8. SFHs for NGC 1433 in different galactocentric radii for the original data (left) and the ND data (right). In the y-axis we display the fraction of mass formed. In the upper part of each panel we use arrows to display the mean age for each SFH, following the same colour coding with respect to the distance to the centre. For both data cubes, it is clear that the farther from the centre, the younger the SFH is (with the exception of the last curve), meaning that the fraction of younger stars increases. This is shown even more strikingly with the subtracted data cube: it is clear that the outskirts of the ND are in fact younger, in accordance with the inside-out growth scenario of Bittner et al. (2019).

5. Discussion

In this section we discuss the results from applying our methodology to the galaxy NGC 1433 and what they reveal about the formation and evolution of bars and NDs.

5.1. The old bar in NGC 1433

Applying the methodology described in Sect. 3 and considering our chosen criterion, we find a bar age of $7.5_{-1.1}^{+1.6}{(\mathrm{sys})}_{-0.5}^{+0.2}(\mathrm{stat})$ Gyr for NGC 1433, which hosts a ND with a radius of 440 pc (e.g., Gadotti et al. 2019; Bittner et al. 2020). This corresponds to a redshift z ≈ 1 and is consistent with other observational studies, such as Simmons et al. (2014), who also find massive galaxies with bars at z ≈ 2. In addition, Gadotti et al. (2015) and Pérez et al. (2017) also inferred similar bar ages for other local galaxies, using independent methodologies. Furthermore, this is comparable to predictions from cosmological simulations, which also find bars that form between redshifts 1 − 2 and that survive down to z = 0 (e.g., Kraljic et al. 2012; Fragkoudi et al. 2020). Our results indicate that NGC 1433 hosts an old and long-lived bar, which is in accordance with several observed characteristics of this galaxy and our current understanding of bar formation and bar-driven evolutionary processes.

As the bar evolves, the shape of its light profile may change. Kim et al. (2015) suggested that the bar Sérsic index (n_bar) evolves from an exponential (n_bar ≈ 1 − 2) to a flat (n_bar ≈ 0.2) profile (but see also Anderson et al. 2022). They also measured n_bar for 144 local galaxies from the S⁴G sample and found that the oldest bars have a n_bar index close to 0.2. In addition, simulations demonstrate that, as the bar grows older, it becomes more elongated and stronger (e.g., Athanassoula et al. 2013), which can be measured by the bar length, R_bar, and the A₂ index. The latter is the second component of the Fourier light/mass distribution in the galaxy and is associated with the bar non-axisymmetry. Lastly, the bar-to-total (Bar/T) luminosity ratio also increases as the bar grows longer and more massive, and the bar morphology evolves to a more rectangular or boxy shape (e.g., Kim et al. 2015). The bar morphology is measured by the index c, where c = 2.0 represents a perfect ellipse and c ≥ 2.7 a strongly boxy bar (see Gadotti 2011). In summary, an old bar is expected to have a flat light profile (n_bar ≤ 0.7), high strength (A₂ ≥ 0.4), a strongly boxy shape (c ≥ 2.7), relatively large bar-to-total ratio (Bar/T ≥ 0.09) and bar length normalised by the disc scale-length R_bar/h ≥ 1.5. Apart from A₂, these values are the median values of the corresponding distributions as found in Gadotti (2011) for a sample of about 300 massive barred galaxies.

NGC 1433 has n_bar = 0.3, A₂ = 0.56, a deprojected, normalised bar length of 1.2, a boxy shape with c = 2.9, and lastly, a Bar/T = 0.08 (see Kim et al. 2014; Díaz-García et al. 2016). Most of these characteristics are consistent with the expectation of an old bar, also in accordance with our bar age estimate of $7.5_{-1.1}^{+1.6}{(\mathrm{sys})}_{-0.5}^{+0.2}(\mathrm{stat})$ Gyr. Lastly, NGC 1433 has similar characteristics as NGC 4371, which was studied by Gadotti et al. (2015) also using the archaeological evidence present on its ND. The authors investigated the oldest stars in the ND to estimate a bar age of ≈10 Gyr. In addition, NGC 4371 has n_bar = 0.2, A₂ = 0.62, deprojected normalised bar length of 1.3, c = 2.7, and Bar/T = 0.08 (see Gadotti et al. 2015 and Díaz-García et al. 2016). Both the age and the characteristics of NGC 1433 and NGC 4371 are very consistent with each other and with the scenario whereby they are amongst the first galaxies to form bars.

The conditions to form a bar are connected with the MD properties, and depend on the galaxy assembling enough mass and settling in more ordered motion. With this in mind, one can expect that the more massive galaxies will reach the necessary conditions first, following a downsizing picture (e.g., Sheth et al. 2012). Despite that, another plausible possibility is the one where early interactions and/or minor mergers, which happened after the disc has settled, could have triggered the bar formation (e.g., Noguchi 1987; Gerin et al. 1989; Miwa & Noguchi 1998; Peschken & Łokas 2019; Łokas 2021). Gadotti et al. (2015) argued how this could be possible for NGC 4371, which is member of the Virgo cluster. Similarly, NGC 1433 is part of the Dorado group (Maia et al. 1989). Following the pre-processing picture, galaxies enter clusters with already some level of processing due to earlier interactions (Haines et al. 2015) while in groups, which could also be responsible for early bar formation ∼7.5 Gyrs ago. Méndez-Abreu et al. (2010, 2012) found that galaxies in the Coma and Virgo clusters that host bars are mainly massive (10⁹ ≤ M/M_⊙ ≤ 10¹¹).

5.2. The inside-out scenario for the growth of nuclear discs

Bittner et al. (2020) showed that in the TIMER sample, the derived age profiles of NDs follow a negative gradient, with the outer parts of the ND being younger. This is in accordance with an inside-out growth scenario, in which one can have a growing gaseous ND that forms stars, and/or star formation concentrated in rings of increasing size, on the external borders of the ND, as shown by Hα maps (Bittner et al. 2020). In this scenario, the gas funnelled by the bar towards the centre piles up near the inner Lindblad resonance of the bar. As the bar grows longer, the radius of the inner Lindblad resonance increases, which leads to the buildup of a ND in an inside-out fashion.

In order to test this scenario, we analyse the SFHs of the subtracted ND in different radial bins, as shown in Fig. 8, and derive the mean ages at each radius (shown as arrows on the top part of the panels). For the original data, the inside-out evidence is subtle, with little change in the mean ages in each radial bin. Nonetheless, with the subtracted data, the inside-out growth of the disc becomes strikingly evident, since in the outer part of the ND the older stellar populations are almost completely absent from the clean ND. This is reflected in the mean age at different radii for the subtracted ND, in which there is a clear gradient towards younger ages at larger radii. This result is also testament that the methodology we develop to subtract the underlying disc component is robust. It shows that the ND of NGC 1433 is consistent with the inside-out growth picture for NDs. This highlights how inner structures such as NDs might be assembled in a self-similar way to the corresponding larger-scale structure of the MD.

Finally, we note that Fig. 8 shows a prominent peak in the SFHs at all radii at young ages (∼2.5 Gyr), which is also evident in Fig. 7: there is a late burst of star formation, which occurs at younger ages than the first burst associated with the formation of the ND at ∼7.5 Gyr. This implies an event that leads to a renewed inflow of gas at late times, which gives rise to such a burst of star formation. Various mechanisms could give rise to such a late gas inflow event through the bar, such as mechanisms that remove angular momentum from gas, for example an interaction or flyby, or – interestingly for this galaxy – this late inflow could be related to the buckling of the bar and the formation of the boxy/peanut bulge (e.g., Pérez et al. 2017).

6. Summary and concluding remarks

In this study we have developed and presented a new method for dating the bar formation epoch of observed disc galaxies. We summarise this work as follows:

• We present a new methodology that allows us to disentangle the light from NDs, which are formed by the bar, from the underlying MDs of galaxies, using high-resolution integral field spectroscopic data from MUSE on the VLT. This allows us to find the time at which the ND formed, by isolating the moment when its star formation starts to dominate over the star formation in the underlying population. As NDs are formed due to bar-driven inflow – which is concurrent with the formation time of the bar – this allows us to determine the age of the bar.

• We performed a number of tests of our methodology, both on observed data and on a hydrodynamic simulation of a barred galaxy (which self-consistently forms a ND), in order to validate the robustness of the methodology.

• As a pilot study, we applied our methodology to the barred galaxy NGC 1433 from the TIMER survey (Gadotti et al. 2019) and find a bar age of $7.5_{-1.1}^{+1.6}{(\mathrm{sys})}_{-0.5}^{+0.2}(\mathrm{stat})$ Gyr. This implies that NGC 1433 has an old bar that formed around z ∼ 1. This aligns with a number of observational characteristics of the galaxy, such as its mass, bar strength, and light profile (see Sect. 5.1 for a more detailed discussion on this). Our results are consistent with other studies in the literature that find old bars (e.g., Simmons et al. 2014; Gadotti et al. 2015), as well as with results from cosmological simulations (e.g., Kraljic et al. 2012; Fragkoudi et al. 2020) that imply that bars can be old, long-lived structures.

• By examining the SFH of the ND of NGC 1433 at different radii, we find that the disc grows inside out, with younger stars forming at progressively larger radii and the youngest stars forming at the edge of the ND (i.e. at the location of the nuclear ring). This is in agreement with an inside-out growth scenario for NDs (see Bittner et al. 2020, and Sect. 5.2).

We will apply the methodology presented here to the full TIMER sample of barred galaxies, which will provide, for the first time, robust age determinations for bars in a sizeable sample of disc galaxies. This will enable us to compare the age of bars with various galaxy properties, such as the total mass, the bar length, and the pattern speed, enabling us to place important constraints on the evolution of dynamical properties of disc galaxies with time, the epoch of disc settling, and the effects of bar-driven evolution.

Acknowledgments

We thank the anonymous referee for the constructive report. This work was supported by STFC grant number ST/T000244/1. Based on observations collected at the European Southern Observatory under ESO programmes 296.B-5054, 097.B-0640 and 099.B-0242. Raw and reduced data are available at the ESO Science Archive Facility. A.d.L.C. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through the Spanish State Research Agency, under Severo Ochoa Centres of Excellence Programme 2020-2023 (CEX2019-000920-S). K.F. acknowledges support through the ESA research fellowship program. P.C. acknowledges support from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP, 2021/08813-7) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, 310555/2021-3). G.G. acknowledges support from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). T.K. was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) grants (No. 2019R1I1A3A02062242) funded by the Ministry of Education, grants (No. 2022R1A4A3031306) funded by the Korean government (MSIT), and grants (WISET 2021-541) funded by the Korea Foundation for Women In Science, Engineering and Technology. J.M.A. acknowledges the support of the Viera y Clavijo Senior program funded by ACIISI and ULL. P.S.B. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MCIUN) under grant number PID2019-107427GB-C31. G.v.d.V. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 724857 (Consolidator Grant ArcheoDyn). J.N. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343) and the Science and Technology Facilities for support through the Consolidated Grant Cosmology and Astrophysics at Portsmouth, ST/S000550/1.

References

Appendix A: Control galaxies: NGC 1380 and NGC 1084

In order to assess whether our methodology creates spurious results for galaxies that do not host a ND, we applied the same methodology for two control galaxies: NGC 1380 and NGC 1084.

Considering NGC 1380, Gao et al. (2019) describe the galaxy as an inclined system with a classical bulge and no clear presence of a bar. The galaxy is at a distance of 21.2 Mpc with inclination of 47° (see references in Sarzi et al. 2018 and Gao et al. 2019). We consider the effective radius of the bulge to be 1080 pc (12.1 arcsec) from the photometric decomposition performed in Gao et al. (2019) in lieu of the ND radius. The data for NGC 1380 comes from the ESO archive, PI: Sarzi, M., programme ID 296.B-5054, using MUSE in Wide Field Mode. Further details about the galaxy and the observations can be found in Sarzi et al. (2018).

As a second control galaxy, NGC 1084 is a bar-less galaxy (Gao et al. 2019) with inclination³ of 49.9° that also does not host a ND. Gao et al. (2019) also describe its morphology with a broken inner disc and a bulge with effective radius of ≈ 150 pc (4.3 arcsec). The data for NGC 1084 comes from the ESO archive, PI: Carollo, C. M., programme ID 099.B-0242, using MUSE in Wide Field Mode.

Figure A.1 shows the outcome of these tests. It is clear that the subtracted spectra show SFHs similar to those in the original spectra. Therefore, our methodology does not artificially produce differences in the SFHs of the regions where the nuclear and main underlying discs dominate, which would wrongly be attributed to the formation of the bar.

Fig. A.1. Same as Fig. 7 but for NGC 1380 (top) and NGC 1084 (bottom). It is clear that for both galaxies neither the representative MD nor the ‘cleaned’ data show important deviations in their SFHs.

Appendix B: Further tests of the methodology

In order to test our methodology and assess the effects of systematic errors, we tested different configurations to assess the possible range of bar ages thus deduced. Our tests include varying the position of the representative ring, using different light profiles to build the underlying MD, not collapsing the data cubes, but rather deriving a mean SFH from the SFHs corresponding to each individual spaxel in the data cubes, and lastly, changing the pPXF regularisation error parameter.

Considering all the possible different configurations, we estimate the final systematic errors as $7.5_{-1.1}^{+1.6}(\mathrm{sys})$ Gyr. In the following we describe each of the tests individually.

B.1. The location of the representative spectrum

Ideally, the representative ring should be located immediately after the end of the ND so that it is not contaminated by light from it and represents, as closely as possible, the underlying population in the region where the ND dominates. However, the end of the ND is not trivial to pinpoint, as this structure may gradually fade into the MD. If the region is too close to the ND, it may be affected by its young star formation, but if it is too far, it may not map the underlying MD old star formation. Our methodology used the physical justification based on kinematic maps to select the representative region. Nevertheless, to constrain how much the position of the representative ring could affect our final result, we produce a number of tests with the representative ring at different positions closer to the centre than the position we employed. The results are summarised in Fig. B.1. We note that the mask used that the mask used for the representative ring of Fig. B.1 is slightly different than the one used in the main result. For this case, we masked all spaxels with AGN contribution, with no threshold for AON (as opposite of AON ≥ 20 used for the main result). With this, the bar age is slightly older, 8 Gyrs (Fig. B.1, last row, center). With this test we aim to assess the relative systematic error due to the location of the representative ring, therefore this does not affect the conclusions. We further note that the change in bar age due to the different mask is well within the systematic error uncertainties.

Fig. B.1. Results for different configurations in the methodology presented in Sect. 3, varying the position of the representative ring (row by row) and the assumed light profile of the main underlying disc. To be continued on the next page.

Fig. B.2. Results for different configurations in the methodology presented in Sect. 3, varying the position of the representative ring (row by row) and the assumed light profile of the main underlying disc. In the first column, we display the stellar velocity dispersion map, indicating the representative ring position with the dashed lines (the solid line corresponds to the peak in v/σ at the outskirts of the ND). In the middle and left columns we display results for exponential and flat light profiles, respectively. One can see that the different configurations result in bar ages in the range 7.0–9.5 Gyr, with the exception of the representative ring closest to the centre. The middle panel of the bottom row corresponds to the standard configuration of our methodology. We note that the mask used that the mask used for the representative ring is slightly different than the one used for the main result. With this, the bar age is slightly older, 8 Gyrs.

With the exception of the most central representative ring, for underlying MDs following an exponential light profile, the derived bar ages are constrained to the range 7–9 Gyr.

B.2. Different main disc profiles: Flat versus exponential

In order to build the underlying MD, one has to assume a light profile that describes it. Although galaxy discs are usually assumed as having exponential light profiles, recent studies demonstrate that some cases may follow a flat light profile (e.g. Zhu et al. 2018; Breda et al. 2020; Papaderos et al. 2022). In order to assess how much such a decision may affect our final result, we compared exponential and flat profiles as extreme possible cases. The results are summarised in Fig. B.1.

For each representative ring position, one sees that the final bar age derived using a flat profile is 0.5-1.0 Gyr older than that derived with an exponential profile (with the exception of the most central representative ring).

Considering the light profile choice together with the position of the representative ring, the derived bar ages vary between 7.0-9.5 Gyr.

B.3. Spaxel-by-spaxel analysis

Considering the elevated signal-to-noise ratio in the TIMER data (Gadotti et al. 2019), it would be possible to run the same methodology spaxel-by-spaxel, deriving spatially resolved SFHs, and then obtaining a mean SFH for each data cube. This would be in contrast to our regular methodology, in which the individual spectra in the data cubes are collapsed into a single spectrum, to then proceed to the derivation of the SFH. The disadvantages of doing so spaxel by spaxel are that the statistical uncertainties for each SFH would be larger and the computational time would increase significantly. Since our goal is to apply the same methodology for the whole TIMER sample, it is too detrimental to require substantial computational time for the analysis of an individual galaxy. Nevertheless, to verify that the results do not differ significantly with the different approaches, we applied the alternative methodology using our data for NGC 1433. For the spatially resolved products, we considered the mean SFHs to measure the bar age. The result is shown in Fig. B.3 and, as one can see, the derived bar age is not substantially different; in fact, the difference is within typical uncertainties in the derivation of stellar ages. With that in mind, we decided to keep our regular methodology employing collapsed spectra to ensure smaller statistical errors and lower computational requirements.

Fig. B.3. Same as Fig. 7, but for the test with spaxel-by-spaxel SFHs instead of a single collapsed spectrum produced for each data cube. The presented SFHs are mean SFHs calculated over all spaxels.

B.4. pPXF regularisation

As described in Cappellari (2017), regularisation allows one to have the smoothest SFH solution, without affecting the physical reality of it. In our methodology, we used the regul_err value of 0.15 derived by Bittner et al. (2020) for the TIMER sample following the procedure developed by McDermid et al. (2015). Nevertheless, in order to assess how much our final result depends on it, we produce tests with values of regul_err in the range 0.1-5.0. The results are summarised in Fig. B.4. Since regul = 1/regul_err, the larger the regularisation error, the less regularised the solution is.

Fig. B.4. Same as Fig. 7 but for different values of the regularisation error parameter applied by pPXF. The original configuration uses a regularisation error of 0.15. The systematic effect of employing different regularisations is constrained to 0.5 Gyr.

As one can see, by changing the regularisation error parameter our bar age has variations of ±0.5 Gyr. Considering the variations due to the representative ring position and the MD light profile, this leaves us with systematic errors of ${}_{-1.1}^{+1.6}$ Gyr, after summing in quadrature the maximum variations we find with the four tests on systematic effects.

Finally, the regularisation error of 1.00 leads to a SFH that looks similar to that of the simulation presented in Fig. 6 (panel e), which could lead to the incorrect interpretation that the first peak (∼8 Gyrs) is a contamination of the MD and the second peak (∼5 Gyrs) is the real bar formation. We point out however that, even though they might look similar at a first glance, one should expect that contamination from the MD would be from the oldest stars, as showed from the different tests for the simulated galaxies (Appendix C.1). Since the oldest peak of the ND is younger than the one from the MD, we do not expect this to be a sign contamination. In addition, this feature is not a well converged feature (as it is not evident for regularisation errors of 0.45 and 5) and is particular to a singular regularisation error. As noted in the main body of the paper, the peak at younger ages could arise through a number of mechanisms that induce gas inflow, such as the formation of the boxy/peanut bulge (e.g. Pérez et al. 2017) as well as through an interaction or flyby that removes angular momentum from gas in the disc, causing a burst in gas inflow to the central regions through the bar.

Appendix C: Hydrodynamic simulations: Description and tests

The simulation used in Sect. 3.3 has both a collisionless (stellar disc and dark matter halo) and collisional (gaseous disc) component and is part of a suite of isolated disc simulations developed to study the evolution of barred galaxies (these will be presented in detail in Fragkoudi & Bieri, in prep.). Here we describe the main properties of the simulation that are relevant for this study. The simulation is run using the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). The AMR grid is refined using a quasi-Lagrangian strategy, and has seven refinement levels, with the maximum resolution reached in the simulation being 10 pc. The initial conditions of the model are created with the Markov chain Monte Carlo code DICE (Perret et al. 2014; Perret 2016). The total mass of the system is Ṁ_tot = 2 × 10¹² M_⊙, with 98.5%, 1.425% and 0.075% of this distributed in the dark matter, stellar and gaseous components, respectively. The dark matter and stellar components have particles with masses 3.7 × 10⁴ M_⊙ and 4.3 × 10⁵ M_⊙, respectively. The dark matter halo has a Navarro-Frenk-White profile (Navarro et al. 1997) with a scale-length of 3 kpc, while the stars (gas) are modelled as an exponential disc, with a scale-length of 3 kpc (4 kpc) and a scale-height of 150 pc (50 pc). Gas in the simulation cools via atomic and metal-dependant cooling processes.

Star formation is modelled as a Schmidt law,

$$\begin{array}{c}\hfill {\dot{\rho }}_{\star }={ϵ}_{\star }{\rho }_{\mathrm{gas}}/t_{\mathrm{ff}},\end{array}$$(C.1)

where ${\dot{\rho }}_{\star }$ is the local star formation rate, and $t_{\mathrm{ff}}=\sqrt{3\pi /(32G{\rho }_{\mathrm{gas}})}$ is the free-fall time. Star formation is triggered when the gaseous density ρ_gas is larger than 1 cm⁻³ and the temperature is less than 100 K, with an efficiency of ϵ_⋆ = 1%. Core-collapse supernova feedback is implemented by assuming that a fraction of the stellar population, η_SN = 0.2, will explode as supernovae. The explosion itself is modelled using the mechanical feedback implementation presented in Kimm & Cen (2014) and Kimm et al. (2015). In this implementation, the supernova explosion is injected into the surrounding interstellar medium according to the phase of the explosion (energy conserving or momentum conserving); in other words, it is injected as momentum or thermal energy, depending on whether the cooling radius is resolved. More details on the feedback implementation can be found in (Kimm & Cen 2014; Kimm et al. 2015).

C.1. Testing the methodology with hydrodynamic simulations

In this section we describe some of the tests carried out on the simulations that allowed us to assess the effects of some of the assumptions in our methodology. In particular, we tested how an assumed underlying age gradient of the old disc, as well as how the location of the ring used to obtain the representative spectrum, affected the results.

In the left panel of Fig. C.1, we show what are the true SFHs of the ND and the underlying MD in the simulations. Then, we show that we can accurately reproduce the time the bar formed using these two values and our chosen criterion for the time of bar formation (i.e. the first moment when the ratio of ND/MD increases above 1 with positive slope towards younger ages; see Sect. 3). We define the true SFH of the ND as the SFH of all stars formed from the gas pushed to the centre after the bar forms within the radius corresponding to the ND (r_ND). The true SFH of the MD is defined as the SFH of all the ‘old’ stars that were present before the bar and ND formed in the MD, within the same radius r_ND. As we can see, the bar age is well constrained by the moment when ND/MD is above 1 with positive slope towards younger ages for the first time. We can compare true SFHs to the recovered SFHs from our methodology in Fig. 6, where we see that the SFH of the ND is well recovered, apart from at the oldest ages.

Fig. C.1. Testing the bar age criterion and how the age gradient in the underlying MD affects our methodology for a simulated galaxy. In the left panel we show the ‘true’ SFH of the ND and MD. The true SFH of the ND is defined as the SFH of all stars formed from the gas pushed to the centre after the bar forms within the ND radius, rND. The true SFH of the MD is defined as the SFH of all the ‘old’ stars in the disc that were present before the bar and ND formed in the disc, within the same radius (i.e. rND). As one can see, in the bottom panel, the bar age criterion – the first moment at which ND/MD is above 1 with a positive slope towards younger ages (the vertical grey dashed line) – successfully times the moment the bar is formed (given by the vertical orange line). By comparing this and Fig. 6 we see that the methodology we present in this paper – which involves using a representative region around the ND to model the MD – is able to recover the SFH of the ND and, therefore, the bar age. In the middle and right panels we explore how the age gradient of the old underlying MD affects our methodology: in the middle panels we show a case where there is a negative age gradient applied to the underlying MD, while in the right panels we show the effects of a flat age gradient on the obtained results. The top panels show the average age gradient as a function of radius, and the bottom panels show how this affects the SFHs of the original ND region (red), of the representative SFH of the underlying MD (green), and of the cleaned ND (blue), and the lower panels show the ND/MD ratio.

In addition, we test how the underlying age gradient of the MD in the galaxy can affect the determination of the bar age. In the right panel of Fig. C.1, we show a case where there is no age gradient in the underlying MD of the galaxy. By subtracting the main underlying disc (as obtained from the representative ring around the ND) from the original data, we obtain a peak in ND/MD that corresponds to the time of bar formation, while before this peak, ND/MD is consistently below < 1. When a negative age gradient is added to the old disc in the simulation (middle panel), we find similar results, with the difference that, at old ages, ND/MD rises again above one. This happens because there are old stars that are present in the ND region, which are not present in the ‘representative’ MD region which is at larger radii. Therefore, due to the age gradient, there is a lack of old stars in the ring, which leads to a rise in ND/MD for the oldest populations. This effect will be larger, the larger the age gradient (which will vary from galaxy to galaxy). Due to this, for our criterion of determining the bar age, we explicitly require that ND/MD be above 1 with a positive slope towards younger ages (rather than simply ND/MD > 1, as this could occur at old ages due to a strong age gradient in the inner galaxy).

In Fig. C.2 we test how the location of the ring affects our methodology, in a case where there is a negative age gradient (if there is no age gradient, the location of the ring does not have a significant effect on the methodology, as long as the SFH of the MD is appropriately rescaled to take into account the higher surface density in the centre of the galaxy, due to its exponential profile). When the ring is placed at larger radii, the representative SFH will be lacking more of the old stars that are present at smaller radii, which can be seen by the increase in ND/MD > 1 for old ages, as we move to larger radii (i.e. the middle and right panels of Fig. C.2). Therefore, the representative ring should not be placed at too large radii from the ND, as the farther away it is placed, the more different the stellar populations will be from those of the underlying MD within the ND region. As described in Sect. 3, we therefore select the radius just outside the ND, where the average age is oldest (see Fig. 3). This ensures we select a region outside the influence of the ND itself, but without going to larger radii where the stellar population properties of the representative region will be significantly different from those of the underlying MD in the ND region.

Fig. C.2. Testing the effects of ring location on the methodology for a simulated galaxy. From left to right we show cases where the ring used to reconstruct the representative SFH of the underlying MD is placed at larger radii, in a case where the model has a negative age gradient in the underlying MD. The top rows show the locations of the rings with respect to the ND region, and the bottom panels show the SFHs. We see that when the ring is placed at larger radii, the representative SFH of the MD does not fully match the ND region at the oldest ages.

We also note that we applied our methodology to other simulations with the same initial conditions, but with modifications in the star formation and stellar feedback prescriptions, and we were able to recover the bar age reliably also with these changes (i.e. the methodology does not depend on the details of the star formation and stellar feedback prescriptions in the simulation).

All Figures

Fig. 1. NGC1433 data and derived maps: The top two images are the colour composites of NGC 1433 from the Carnegie-Irvine Galaxy Survey (top left; Ho et al. 2011) and TIMER (top right; Gadotti et al. 2019). Below are maps that show (from left to right and top to bottom) the stellar kinematics and population properties derived from the TIMER data using GIST: stellar velocity, velocity dispersion, h3, h4, age, metallicity, and α enhancement. The ND radius is shown with a solid black contour and displays a faster rotation, a drop in velocity dispersion, an anti-correlation between velocity and h3, and an increase in h4. In addition, the same region corresponds to a drop in mean age and α enhancement and an increase in metallicity. It is clear that NGC 1433 hosts a younger ND with more rotational support in the central region than the underlying population, in agreement with Gadotti et al. (2020) and Bittner et al. (2020).

In the text

Fig. 2. Illustration of the methodology described in Sect. 3 for a galaxy that hosts a ND. In the top left we highlight the field of view from MUSE in the centre of the galaxy and the position of the representative region just outside the ND. The position of the representative ring is chosen taking the v/σ radial profile into consideration. From it we derive the representative spectrum used to build the MD. In the bottom left we display the output from the light disentangling: the original data cube (red), the ND data cube (blue), and the representative MD data cube (green). All three data cubes have spaxels corrected to v = 0 km s−1 and σ = σmax km s−1. Lastly, we collapse each data cube into a mean spectrum (as illustrated) and derive SFHs for each one. The steps in our methodology are described to the right.

In the text

Fig. 3. v/σ radial profile, displayed in dashed red contours (values in the left-axis), and light-weighted average age, in solid blue contours (values in the right axis) for NGC 1433. We display the median values together with the first and third quartiles. The vertical dotted black line marks the ND radius and the hatched area the representative ring. Note that the representative ring is placed in the first v/σ minimum outside the ND, which matches the oldest mean age.

In the text

Fig. 4. NGC1433 BPT classification: Left: BPT diagram (Baldwin et al. 1981) classification of each spaxel for the ND in NGC 1433. The BPT classification is continued in the middle and right panels, which display the physical position of the different spaxels and the light radial profile, respectively. We also display the radius of the inner mask dominated by AGN contributions (dotted contours), the radius of the ND (solid contours), and the representative ring (dashed contours).

In the text

Fig. 5. Illustration of different data products derived from the methodology described in Sect. 3. We display the sum of the original MUSE fluxes between 4800 and 5800 Å (top row), and the derived representative MD (bottom row, left) and the ND data (bottom row, right), which is the result of subtracting the representative MD from the original data cube. For the representative MD, it is possible to notice the exponential increase in flux towards the centre. All data cubes were masked for AGN contributions with AON ≥ 20.

In the text

Fig. 6. Testing the methodology using hydrodynamic simulations. a: face-on projection of the gas temperature. The dust lanes on the leading edge of the bar and the gaseous ND are clearly visible as low temperature regions. b: kinematic map showing the velocity of stars in the y direction with stellar isodensity contours over plotted, which outline the shape of the bar (note the different scale compared to panel a). A highly rotating stellar component in the central kiloparsec (i.e. the ND) is clearly evident in the kinematic map. c: age map of all the stars in the simulated galaxy. The grey box shows the inset that is represented in panel d. d: zoomed-in view around the ND region. The inner and outer dashed lines denote the region used to obtain the SFH of the ‘representative disc region’. The solid black line denotes the radius within which the SFH of the ND is estimated. The scatter points indicate the locations of the pixels used for deriving the SFHs in panel e. The points are colour-coded by the mean age in the pixel. e: the top panel shows the SFHs of the original ND region (solid red), of the representative SFH of the MD (dot-dashed green), and of the ND with the representative SFH subtracted (dashed blue). The bottom panel shows the ratio of the subtracted SFH to the representative SFH. The vertical orange line indicates the time of bar formation in the simulation (3 Gyr), and the vertical dashed line indicates the time at which ND/MD is above 1 (3.3 Gyr).

In the text

Fig. 7. Mass-weighted SFHs for NGC 1433 from collapsed spectra of the MUSE original data, ND data, and representative MD data, as illustrated in Fig. 2. The y-axis corresponds to the mass formed, in solar masses, for each age bin and the x-axis to the age in Gyr, that is, the look-back time. The original data are shown with a solid red line, the subtracted data with a dashed blue line, and the representative spectrum with the dot-dashed green line. The lines are the SFHs derived from the data, while the shaded regions are results from 100 Monte Carlo runs. In the bottom panel we display ND/MD in black, and highlighted with an orange circle is the age where ND/MD rises above one for the first time, together with the statistical uncertainty: ${}_{-0.5}^{+0.2}(\mathrm{stat})$ Gyr.

In the text

Fig. 8. SFHs for NGC 1433 in different galactocentric radii for the original data (left) and the ND data (right). In the y-axis we display the fraction of mass formed. In the upper part of each panel we use arrows to display the mean age for each SFH, following the same colour coding with respect to the distance to the centre. For both data cubes, it is clear that the farther from the centre, the younger the SFH is (with the exception of the last curve), meaning that the fraction of younger stars increases. This is shown even more strikingly with the subtracted data cube: it is clear that the outskirts of the ND are in fact younger, in accordance with the inside-out growth scenario of Bittner et al. (2019).

In the text

	Fig. A.1. Same as Fig. 7 but for NGC 1380 (top) and NGC 1084 (bottom). It is clear that for both galaxies neither the representative MD nor the ‘cleaned’ data show important deviations in their SFHs.
In the text

	Fig. B.1. Results for different configurations in the methodology presented in Sect. 3, varying the position of the representative ring (row by row) and the assumed light profile of the main underlying disc. To be continued on the next page.
In the text

Fig. B.2. Results for different configurations in the methodology presented in Sect. 3, varying the position of the representative ring (row by row) and the assumed light profile of the main underlying disc. In the first column, we display the stellar velocity dispersion map, indicating the representative ring position with the dashed lines (the solid line corresponds to the peak in v/σ at the outskirts of the ND). In the middle and left columns we display results for exponential and flat light profiles, respectively. One can see that the different configurations result in bar ages in the range 7.0–9.5 Gyr, with the exception of the representative ring closest to the centre. The middle panel of the bottom row corresponds to the standard configuration of our methodology. We note that the mask used that the mask used for the representative ring is slightly different than the one used for the main result. With this, the bar age is slightly older, 8 Gyrs.

In the text

	Fig. B.3. Same as Fig. 7, but for the test with spaxel-by-spaxel SFHs instead of a single collapsed spectrum produced for each data cube. The presented SFHs are mean SFHs calculated over all spaxels.
In the text

	Fig. B.4. Same as Fig. 7 but for different values of the regularisation error parameter applied by pPXF. The original configuration uses a regularisation error of 0.15. The systematic effect of employing different regularisations is constrained to 0.5 Gyr.
In the text

Fig. C.1. Testing the bar age criterion and how the age gradient in the underlying MD affects our methodology for a simulated galaxy. In the left panel we show the ‘true’ SFH of the ND and MD. The true SFH of the ND is defined as the SFH of all stars formed from the gas pushed to the centre after the bar forms within the ND radius, rND. The true SFH of the MD is defined as the SFH of all the ‘old’ stars in the disc that were present before the bar and ND formed in the disc, within the same radius (i.e. rND). As one can see, in the bottom panel, the bar age criterion – the first moment at which ND/MD is above 1 with a positive slope towards younger ages (the vertical grey dashed line) – successfully times the moment the bar is formed (given by the vertical orange line). By comparing this and Fig. 6 we see that the methodology we present in this paper – which involves using a representative region around the ND to model the MD – is able to recover the SFH of the ND and, therefore, the bar age. In the middle and right panels we explore how the age gradient of the old underlying MD affects our methodology: in the middle panels we show a case where there is a negative age gradient applied to the underlying MD, while in the right panels we show the effects of a flat age gradient on the obtained results. The top panels show the average age gradient as a function of radius, and the bottom panels show how this affects the SFHs of the original ND region (red), of the representative SFH of the underlying MD (green), and of the cleaned ND (blue), and the lower panels show the ND/MD ratio.

In the text

Fig. C.2. Testing the effects of ring location on the methodology for a simulated galaxy. From left to right we show cases where the ring used to reconstruct the representative SFH of the underlying MD is placed at larger radii, in a case where the model has a negative age gradient in the underlying MD. The top rows show the locations of the rings with respect to the ND region, and the bottom panels show the SFHs. We see that when the ring is placed at larger radii, the representative SFH of the MD does not fully match the ND region at the oldest ages.

In the text