© ESO 2021

1. Introduction

A two-phase formation scenario is becoming increasingly favoured to explain the formation of massive early-type galaxies (ETGs; e.g., Naab et al. 2009; Oser et al. 2010; Hilz et al. 2013; Rodriguez-Gomez et al. 2016; Zibetti et al. 2020). According to this scenario, a first phase characterised by an intense and fast dissipative series of processes occurring at early cosmic times (z >  2), forms blue, disk-dominated, and centrally concentrated galaxies (van der Wel et al. 2011) that form stars ‘in situ’ at a high rate (star formation rate, SFR ≥ 10³ M_⊙ yr⁻¹): the so-called blue nuggets (Barro et al. 2013; Dekel & Burkert 2014; Zolotov et al. 2015; Tacchella et al. 2016). Subsequently, the star formation quenches, for reasons that are still heavily debated (see e.g., Martín-Navarro et al. 2019 and references therein) and the blue nuggets quickly and passively evolve into compact ‘red nuggets’ (Daddi et al. 2005; Damjanov et al. 2009).

These high-z red nuggets can typically reach masses similar to those of local giant elliptical galaxies (M_⋆ >  10¹¹M_⊙), which indicates that a large portion of the mass is assembled during this first formation phase. However, their sizes are only about a fifth of those of local ETGs of similar masses (Daddi et al. 2005; Trujillo et al. 2007; van Dokkum 2008; Werner et al. 2018). Thus, during a second phase, at lower redshifts, red nuggets must increase in size, becoming, over billions of years, the present-day giant ETGs.

After the second and more time-extended phase, the stars formed during the first phase become the ‘core’ of today’s ETGs and occupy their innermost regions. Stars ‘accreted’ during the second phase, instead, live preferentially at larger radii (e.g., Spavone et al. 2017), although it has also been shown that major mergers and accretion of more massive satellites can deposit stars at the centres of the hosts (Amorisco 2017).

Since theoretical evidence suggests that massive (M_⋆ > 10¹¹ M_⊙) halos at z = 2 typically experience no more than one or two major mergers until z = 0 (Genel et al. 2008; Khochfar & Silk 2006), and this is not enough to account for the required mass and size growth; a series of minor dry mergers (i.e. accretions of small gas-poor satellites) seems to be the preferred channel for size evolution, building extended stellar envelopes around compact galaxy cores (Naab et al. 2009; van Dokkum et al. 2010; Oser et al. 2010, 2012; Tortora et al. 2018a; Stockmann et al. 2020). Together with a continuous addition of larger, newly quenched (blue) galaxies to the quiescent population (e.g., van der Wel et al. 2009; Saglia et al. 2010; Carollo et al. 2013; Cassata et al. 2013), this scenario may explain the observed increasing average size of quiescent galaxies with decreasing redshift (e.g., van der Wel et al. 2014). The early domination of in situ star formation and the subsequent growth by stellar mergers is also in good agreement with predictions from semi-analytical models (e.g., De Lucia et al. 2006; Guo & White 2008; Guo et al. 2011; Shankar et al. 2010) and with observational results in the local Universe (Coccato et al. 2010; Longobardi et al. 2015; Iodice et al. 2017; Pulsoni et al. 2018; Ferré-Mateu et al. 2019; Zibetti et al. 2020).

Unfortunately, in local giant galaxies, the accreted material contaminates the in situ component¹ that encodes the information about high-z baryonic processes, affecting, along the line-of-sight (LOS), its spatial and orbital distributions, and thus irreversibly limiting our resolving power (e.g., Naab et al. 2014; Ferré-Mateu et al. 2019). This makes it very difficult to observationally separate star formation history from assembly history. Fortunately, given the stochastic nature of the merging processes, a small fraction of massive quiescent and very compact galaxies must exist at all redshifts, although this fraction should decrease with time.

Relic galaxies, which are red nuggets that survived intact to the present-day Universe without experiencing any merger or interaction, provide a unique opportunity to track the formation of the pristine in situ galaxy stellar component, which is otherwise mixed with the accreted one in normal massive ETGs (Trujillo et al. 2009). A detailed study of the stellar population of a large, statistical sample of high-z red nuggets would require prohibitive integration times with the current available facilities. Thus, relics are the only systems that allow us to study the physical processes that shaped the mass assembly of massive galaxies in the high-z Universe with the amount of details currently reachable only in the nearby Universe (e.g., high spatial resolution and spectroscopic signal-to-noise ratio).

However, the exact fraction of such systems at lower redshift critically depends on the physical processes shaping the size and mass evolution of galaxies, such as the relative importance of major and minor galaxy merging (e.g., Naab et al. 2009; Oser et al. 2010, 2012; Oogi & Habe 2013; Quilis & Trujillo 2013; Wellons et al. 2016). Moreover, many questions have still to be answered: how many relics exist in the Universe at each redshift, and how many survive until the present day? How can they passively evolve through cosmic time without experiencing any interaction? What can we learn on the merging history and the evolution of the most massive ETGs? Is there a physical scenario, predicted by hydrodynamical cosmological simulations, able to explain all the current observational results at all z? These topics are currently heavily debated.

At high redshift (z >  1), the number density of compact quiescent galaxies is measured to be in the range ρ ∼ 10⁻⁵–10⁻⁴ Mpc⁻³ (e.g., Bezanson et al. 2013; Saracco et al. 2010; Cassata et al. 2010, 2013; van der Wel et al. 2014). In the local Universe, so far, only three relics have been spectroscopically confirmed: Mrk 1216, PGC 032873, and NGC 1277, (Trujillo et al. 2014; Ferré-Mateu et al. 2017, hereafter F17). Based on these three confirmed relics, F17 obtained a lower limit to the number density (ρ) of relics in the local Universe of 6 × 10⁻⁷ Mpc⁻³. This is also consistent with the number density of massive object that formed early and evolved varying their stellar mass by less than 30% predicted from either semi-analytical models (Millennium simulations – Quilis & Trujillo 2013) or hydrodynamical simulations (Illustris – Wellons et al. 2016).

The three local relics have large rotation velocities (V ∼ 200–300 km s⁻¹) and high central stellar velocity dispersions (σ_⋆ >  300 km s⁻¹); moreover, their morphology and density profiles match those found in high-z red nuggets and in nearby fossil groups remarkably well (e.g., NGC 6482, see Buote & Barth 2019). A further confirmation of the genuine relic nature of one of these three systems, NGC 1277, came from the uni-modal and uniquely red colour distribution of its globular cluster population (Beasley et al. 2018). For Mrk 1216, Buote & Barth (2019) found a highly concentrated dark matter (DM) halo, reflecting the extremely quiescent and isolated evolution of its stars, and a larger DM fraction within the effective radius.

From a stellar population point of view, the three objects found by F17 have a single stellar population with super-solar metallicities and old ages out to several effective radii. Their stars are also characterised by a large [Mg/Fe] over-abundance, which is consistent with early and short star formation episodes (within timescales < 1 Gyr, Thomas et al. 2005). Finally, the three relics also have a bottom-heavy stellar initial mass function (IMF), with the fraction of low-mass stars being at least a factor of 2 larger than that found in the Milky Way (Martín-Navarro et al. 2015; F17). This result is in good agreement with the proposed tight relation between IMF slope and galaxy velocity dispersion (e.g., Conroy & van Dokkum 2012; Cappellari et al. 2012; Tortora et al. 2013; La Barbera et al. 2013; Spiniello et al. 2014). However, while for ETGs with ‘normal sizes’ the bottom-heavy IMF is concentrated only in the very centre (Sarzi et al. 2018; Parikh et al. 2018; La Barbera et al. 2019; Barbosa et al. 2021), for the F17 relics the IMF remains virtually invariant with radius. A possible working hypothesis may be that the excess of dwarf stars originates in the first phase of the mass assembly, the only one that relics have experienced (Smith 2020).

Clearly, in order to validate these various scenarios and test the predictions for the size evolution from different models, more relic statistics are needed, extending the redshift boundaries to bridge the gap with the high-z Universe. This constitutes one of the main goals of the INvestigating Stellar Population In Relics (INSPIRE) Project, which we present in this paper.

The starting point of the INSPIRE project is an approved ESO Large Program (LP; ID: 1104.B-0370, PI: C. Spiniello) of 154 hr of observations on the X-shooter spectrograph (XSH, Vernet et al. 2011) at the ESO Very Large Telescope. The LP is spread over two years and targets 52 confirmed ultra-compact massive galaxies (UCMGs) with redshifts 0.1 <  z <  0.5 and with red colours, found in the Kilo Degree Survey (KiDS, Kuijken 2011) Data Release 3 (DR3, de Jong et al. 2017) from a dedicated project (Tortora et al. 2016, 2018b; Scognamiglio et al. 2020).

Ultra-compact massive galaxies are outliers in the mass-size relation (e.g., Shen et al. 2003) since they have very small sizes (R_e  ≤  1.5 kpc) despite their relatively large stellar masses (M_⋆ ≥ 8 × 10¹⁰ M_⊙)². They are therefore the perfect relic candidates, to test whether they are populated by a very old single stellar population, formed in situ, likely with a super-solar [α/Fe] and metallicity, and with very little or no sign of subsequent star formation events. Proving this, and thus constraining the stellar population of these objects, is, indeed, the main goal of the INSPIRE Project.

By means of a detailed stellar population analysis on the sample of 52 KiDS-confirmed UCMGs (relic candidates), once observations are completed, we will build the largest sample of relics at z <  0.5 obtaining a precise relic number density and its evolution in the low-to-intermediate redshift window. The INSPIRE catalogue will represent a fundamental benchmark for ongoing and future cosmological simulations, which will have to be able to reproduce the number density and characteristics of these systems in order to predict the formation and evolution of local massive galaxies. Here in this pilot paper we present the preliminary results on the first three INSPIRE targets, demonstrating that our observation and data analysis strategy works, and therefore that we will be able to assess the true (relic) nature of the candidates.

This paper is organised as follows. In Sect. 2, we give an overview of the INSPIRE ESO LP, including details on the sample definition and the current status of the observations and analysis. In Sect. 3, we introduce the pilot study based on the first three fully observed and analysed galaxies. We describe how we obtained the final UVB+VIS spectrum for each system and how we calculated the line-of-sight velocity distribution and the stellar population parameters from full spectral fitting. Finally, in Sect. 4 we present the results of this study, including a comparison to F17, and in Sect. 5 we draw our conclusions. Throughout the paper, we assume a standard ΛCDM cosmology with H₀ = 69.6 km s⁻¹ Mpc⁻¹, Ω_vac = 0.714 and Ω_M = 0.286 (Bennett et al. 2014).

2. The ESO Large Program INSPIRE

The INSPIRE project aims to characterise the kinematical, dynamical, and stellar population properties of 52 confirmed UCMGs at redshift 0.1 <  z <  0.5 with unprecedented detail. These objects have been selected because their multi-band photometry is consistent with having old stellar populations, from multi-band photometry (see Sect. 2.1 for more details).

INSPIRE has been designed to obtain spectra with a signal-to-noise (S/N) high enough for detailed spectroscopic stellar population analysis, which is the only way to confirm the relic nature of the stellar populations of these UCMGs. Furthermore, the large wavelength range of XSH is particularly suitable to infer the IMF slope at its low-mass end, which represents one of the main pillars of the project. In fact, a very wide wavelength interval (also covering the near-infrared where cold dwarf stars contribute up to 30% of the stellar light (Worthey 1994)) is crucial in breaking the stellar population degeneracies between variations in age, metallicity, elemental abundances, and IMF slope (Spiniello et al. 2014; La Barbera et al. 2019).

In the following sections, we describe how the relic candidates targeted by the INSPIRE project have been selected, and we give some details regarding the strategy and the current status of the observations. The full sample will be presented in the first INSPIRE Data Release, expected at the beginning of 2021, where we will also describe in more detail the main goals and scientific objectives of the project. Constraints on the IMF will be presented in a dedicated forthcoming publication, since this requires very careful treatment of the telluric atmospheric absorption lines in the infrared, and possibly a more sophisticated approach to the stellar population analysis.

2.1. Sample definition

With a dedicated collaboration project, our team has started a systematic census of UCMGs in the redshift range 0 <  z <  0.5 within the KiDS footprint. Thanks to its exquisite image quality (angular scale of 0.21″/pixel and a median r-band seeing of ∼0.65″; Kuijken 2011; de Jong et al. 2015), KiDS is optimally suited for the purpose of carrying on a complete census of these ultra-compact systems. In Tortora et al. (2018b, hereafter T18), we presented ∼1000 candidates, with R_e ≤ 1.5 kpc and M_⋆ ≥ 8 × 10¹⁰, collected from the first 333 deg² of KiDS (DR3, de Jong et al. 2017).

We also recently started a programme to spectroscopically confirm these candidates using public spectra from the literature (2DFLenS, SDSS, and GAMA) as well as new spectroscopy. We currently confirmed ∼100 of them (T18 and Scognamiglio et al. 2020, hereafter S20). Although these spectra allowed us to estimate the redshifts and thus confirm the extragalactic, ultra-compact, and massive nature of these objects, their S/Ns were too low to use them for any stellar population analysis.

We built the largest sample of spectroscopically confirmed UCMGs at redshifts 0.1 <  z <  0.5. This redshift range is an important link between the very small sample of ultra-compact galaxies in the local Universe and the large population of high-redshift red nuggets. In particular, it represents a sweet spot where it is possible to survey larger volumes with large-sky photometric surveys but also to target the candidates with spectroscopic surveys/facilities (Damjanov et al. 2013, 2014, 2015a; Buitrago et al. 2018), as they are brighter than the higher-z counterparts. Currently, this is the only way to enable statistical studies of structural, dynamical, and environmental properties of UCMGs (see e.g., Damjanov et al. 2015b; Tortora et al. 2020).

Among the ∼100 spectroscopically confirmed UCMGs, half of them have broad-band colours consistent with old stellar populations, typical of red and dead galaxies. To show this, and thus select the optimal relic candidates, we compared the UCMG’s g − i colours to those of single-burst stellar population models from Bruzual & Charlot (2003). We then only selected objects consistent with a formation epoch that occurred at least 8 Gyr ago (also allowing for super- and sub-solar metallicities), as visible from Fig. 1, where the g − i colour (in the observed-frame) is plotted against the redshift. Blue points in the figure are the ∼100 UCMGs for which we have spectroscopic redshifts available. The systems selected for the INSPIRE ESO LP are plotted in green, while the three magenta diamonds show the position of the three systems that we study in this pilot paper, as described in Sect. 3. Finally, the yellow triangles show the positions of two of the three confirmed local relics in F17 (NGC1277 and PGC 032873) that are also in the SDSS footprint. For these objects, the g − i colour is computed using the ‘modelMag’ from the SDSS DR16 (Ahumada et al. 2020). Unfortunately, the third relic in F17 (Mrk 1216) is not in the SDSS footprint and thus we do not have magnitudes available for this object.

Fig. 1. g − i colour plotted versus redshift for the spectroscopically confirmed UCMGs in KiDS (∼100, T18; S20). The relic candidates, targets of the XSH INSPIRE LP, are plotted as green points. The three objects, which are the focus of this study, are plotted as magenta diamonds. Blue points are UCMGs from T18 and S20 consistent with younger stellar ages. Two out of the three confirmed relics of F17 are plotted as yellow triangles (colours are taken from the SDSS DR16). Lines show single-burst old stellar population models (Bruzual & Charlot 2003) with solar (Z = 0.02, solid lines), sub-solar (Z = 0.008, dotted), and super-solar (Z = 0.05, dashed) metallicity and ages of 12 Gyr (black) and 8 Gyr (red).

The 52 INSPIRE targets span a very wide range in right ascension and declination (i.e. the optimum observability from Paranal spread over a full year) and are all brighter than 20mag in the r-band. Since the effective radii of the galaxies are very small and there is no possibility to spatially resolve them from the ground (without adaptive optics), we required relatively ‘poor’ conditions for the observations (seeing up to 1.2″, CLR nights, grey lunar phase up to a Moon FLI of 0.6). This increased the efficiency of the service mode under a LP.

For all 52 systems, we have multi-band optical photometry available from the KiDS (DR3) catalogue, structural parameters derived in Roy et al. (2018) and spectroscopic redshifts retrieved from the literature (from SDSS, Ahn et al. 2012, 2014; 2dFLenS, Blake et al. 2016; GAMA, Driver et al. 2011) or calculated through new spectroscopic observations presented in previous papers of our group (T18 and S20).

These UCMGs are the optimal relic candidates: they have been confirmed to be massive (M_⋆ ≥ 6 × 10¹⁰ M_⊙) and compact (R_e  ≤ 2 kpc) extragalactic objects and are consistent with having passive stellar populations. We note that the size and mass selection criteria used here to select relic candidates are slightly more relaxed than those used for UCMGs in T18 and S20. In fact, these systems were originally selected to have R_e  ≤  1.5 kpc and M_⋆ ≥ 8 × 10¹⁰ M_⊙ on the basis of their photometric redshifts. However, after recomputing these quantities using the spectroscopic redshifts instead, some sizes (stellar masses) resulted to be slightly larger (smaller) for some objects, which could still be interesting, on the basis of their stellar population properties. Because the validation of UCMGs is based merely on the measured values, without taking into account any uncertainties (which in some extreme cases, i.e. for very small radii, can be as high as 50–80%: see e.g., T18), and we need to build a statistically large sample, we relaxed both the compactness and massiveness thresholds. Thus we selected, for the XSH Proposal, galaxies with R_e <  2 kpc³ and M_⋆ ≥ 6 × 10¹⁰ M_⊙.

All these objects have similar or slightly smaller sizes and masses than the three confirmed local relics of F17⁴ and high-z red nuggets, and they lie below the mass-size relation defined from normal-size ETGs of equal masses at similar redshift (Shen et al. 2003; Roy et al. 2018), as can be seen from Fig. 2. In terms of their morphology, like the confirmed local relics, the majority of the INSPIRE targets show compact elongated disk-like shapes typical of the massive galaxy population at redshifts z >  2 (F17; van der Wel et al. 2011). Finally, the KiDS images, which usually go as deep as 26 mag arcsec⁻² in the r-band surface brightness (see Fig. 3 in T18), do not seem to reveal signatures of tidal tails, asymmetries, or any other features that could indicate any past or current interaction or stripping (see Fig. 2 in T18 and Fig. 1 in S20)⁵.

Fig. 2. Stellar mass-size relation for high-z nuggets (black points, from Szomoru et al. 2012), the F17 local relics (yellow triangles), and the relic candidates targeted by INSPIRE (green points and magenta diamonds), compared to the relation for normal SDSS ETGs (grey line, from Shen et al. 2003). By selection, all UCMGs are clear outliers of the local stellar mass-size relation, and high- and low-z systems lie on the same region of the diagram.

2.2. Observation strategy and current status

The observing strategy was set to be the same for all the systems: we chose slit widths of 1.6″, 1.5″, and 1.5″, respectively for UVB, VIS, and NIR to ensure minimal slit loss. The XSH slit length is 11″, which ensures that we cover within the slit a region far enough from the galaxy, allowing us to perform a proper sky subtraction. For this purpose, we used a dithering scheme (in the NODDING MODE) with multiple frames where the galaxy is offset by a small amount from the centre of the slit (jitter box of 0.2″ along the slit width and a shift of 2″ along the slit length). The position angle (PA) was always aligned with the major axis of the galaxies, which was inferred from the r-band KiDS images.

We scheduled from 1 to 4 OBs (of one hour each) per system, based on their surface brightness luminosities (obtained performing a Seŕsic fit) averaged within the effective radius (⟨μ_e⟩) and aperture magnitudes in the r-band (obtained from the KiDS catalogue). To determine the integration times for each target, we used the ESO XSH exposure time calculator with a template spectrum of an elliptical galaxy and testing both an extended spatial distribution and a point-like source, since our objects are very compact on the sky (r-band R_e  ∼ 0.2−0.6″). In particular, the integration times are driven by the quality of the spectra we need to reach in terms of S/N, in order to achieve the scientific goal of accurately constraining the integrated stellar population parameters, including a precise estimate of the stellar IMF slope. Based on our previous experience and on literature results (Cid Fernandes et al. 2014; Ge et al. 2018; Costantin et al. 2019), an integrated S/N ≥ 15 per Å allows the recovery of stellar age, metallicity, and [α/Fe] with uncertainties smaller than 0.1 dex (using both full-spectral fitting with single stellar population modelling and line index measurements), and thus the confirmation of the relic’s nature. However, a much higher S/N (> 40 per Å, based on our experience) is required, at least in the NIR, to also properly infer the low-mass end slope of the IMF (see e.g., La Barbera et al. 2019). For this purpose, XSH is the ideal instrument since its very large wavelength baseline allows one to break the degeneracies between the stellar population parameters and at the same time obtain a secure and precise constraint (with 10% accuracy) on the low-mass end slope of the IMF. We will constrain the IMF, fitting the full UVB-to-NIR wavelength range (at least for the objects with the highest S/N spectra) in a forthcoming paper of this series since this requires a more careful analysis of the infrared spectrum. In this way, we will be able to draw a scenario to explain the different observations derived for relics and normal-size galaxies. At the time of writing, the first period of observations (ESO-P104) has concluded, and 42 hours of observations have been carried out on 17 galaxies.

The first INSPIRE Data Release, and the associated dedicated paper, comprising all the spectroscopic observations carried on during the first semester (running from October 2019 until March 2020) is planned for the first half of 2021. Meanwhile, to test our observations, data reduction, and analysis strategy, we present here the initial results for a pilot study focusing on three systems, representative of the whole sample, and for which the spectroscopic observations have already been completed: KiDS J0224−3143, KiDS J0314−3215, and KiDS J0847+0112.

Another necessary step, which is the subject of ongoing work, is the correction of the telluric atmospheric absorption lines, which severely affect the red part of the VIS and the NIR arms. Therefore, for the moment, we limit the analysis to the UVB and blue part of the VIS (from ∼350 to ∼650 nm), which is sufficient to derive the line-of-sight velocity distribution (LOSVD)⁶ and to break the stellar age-metallicity degeneracy, and thus constrain the stellar age and metallicity from full spectral fitting and α abundance from line indices with uncertainties smaller than 10%.

3. The INSPIRE pilot study: KiDS J0224−3143, KiDS J0314−3215, and KiDS J0847+0112

Observations for KiDS J0224−3143, KiDS J0314−3215, and KiDS J0847+0112, to which we gave priority to test our strategy, data reduction, and analysis routines, were completed in 2019. We highlight that these three objects are representative of the whole sample as they span a relatively wide range in redshift, stellar masses, sizes and magnitudes. In particular, KiDS J0224−3143 is the most massive system in our sample (M_⋆ = 2.7 × 10¹¹ M_⊙), and among the galaxies with the highest redshift (z = 0.3841 ± 0.0001, corresponding to a Universe that is 9.5 Gyr old⁷) KiDS J0314−3215 is slightly less massive and closer (z = 0.2874 ± 0.0001, when the Universe was 10.4 Gyr old), but it is one of the smallest systems in the current sample, with an effective radius of only R_e = 0.65 kpc, corresponding roughly to 0.16″ on the sky (a single pixel in XHS). Finally, KiDS J0847+0112 is brighter than the other two (mag_r = 18.41-pagination mag), and it is one of the INSPIRE galaxies with the lowest redshifts (z = 0.1764 ± 0.0001, when the Universe was 11.5 Gyr old).

The main morpho-photometric characteristics of these three objects are listed in Table 1. The first block (separated by vertical lines) gives the KiDS ID and the coordinates of the objects; the second block lists a number of properties derived from the r-band KiDS images (shown in Fig. 3), which is the band with the best resolution and seeing from KiDS. In particular, we list the total magnitude (MAG_AUTO, in the KiDS catalogue, which is obtained from SEXTRACTOR, Bertin & Arnouts 1996) and the average surface brightness within the effective radius in the r-band. We also give the effective radius, R_e, Sérsic index, n, and axis ratio, q, always in the same band and calculated fitting a PSF convolved Sérsic profile to the images with the code 2DPHOT (La Barbera et al. 2008). These latter two quantities demonstrate that the objects are bulge-dominated but elongated, thus hinting at a disk-like shape (see also T16; T18; S20). Finally, in the last block, we provide the spectroscopic redshifts, the median (of the values for g, r and i bands) effective radii converted into kpc, and the associated stellar masses as derived in S20.

Fig. 3. Cutouts of 1′  ×  1′ of the KiDS r-band images for the three relic candidates that we study in this paper. Spatial scale and image orientation are common to all the cutouts and are shown in the first and third panels, respectively.

Quantities that refer to the XSH observations are given in Table 2. These include the position angle, PA, which is always aligned with the photometric major axis, the final exposure time on target, the median seeing of the observations, and the spectroscopic redshift directly inferred from the XSH final spectra. The redshifts we derive from these new spectra are qualitatively consistent (within 5%) with those derived in T18 (and reported again in S20) for the first two galaxies, but more reliable thanks to the better quality of the new spectra. For J0847+0112, the redshift estimate we compute is consistent with the one inferred by the GAMA Survey (z = 0.17636).

3.1. Data reduction and extraction of the 1D spectra

The data has been reduced with the ESO XSH pipeline (v3.3.0) under the EsoReflex workflow (Freudling et al. 2013) until the production of the 2D flux-calibrated spectra. The extraction of the 1D ones, which is described in more detail below, was instead performed making use of IDL and Python routines specifically developed for the purpose.

As mentioned before, the spatial resolution of XSH is not good enough to resolve such compact objects. In fact, in all the cases, the effective radii in arcseconds (0.25″, 0.15″, and 0.47″ in the r-band for KIDS J0224−3143, KiDS J0314−3215, and KiDS J0847+0112, respectively) are much smaller than the seeing (which varied from ∼0.8″ to ∼1.2″ for these three systems). In the extreme case of KiDS J0314−3215, the half-light radius is comparable to a single detector pixel. This means that the spectra are completely seeing dominated. Thus, the most meaningful approach to compare the galaxies with each other is to try to extract the spectra from an aperture that contains more or less the same fraction of light for the three objects, given the seeing during observations and their surface brightness profiles. We therefore used the following procedure: (i) we extracted the surface brightness profiles of the three objects directly from the 2D spectra (collapsing the 2D frame in the dispersion direction and assuming spherical symmetry); (ii) we integrated these profiles up to the radius where the flux reaches the null value, which we consider as the radius containing the total light; (iii) we computed the same integral for different apertures (i.e. integrating only up to smaller and smaller radii); and (iv) we took the ratio between these integrals and the total one to obtain the fraction of light as a function of the aperture radius.

We find that extracting from a fixed aperture of r = 0.4″, corresponding to the five innermost central pixels, roughly corresponds to 30% of the total light for all three objects⁸. With this aperture, we maximise the S/Ns of the extracted spectra and make sure to probe a region dominated by baryons. This strategy will be revised in the forthcoming papers of this series, where we will also test an optimal extraction procedure, which should allow us to further increase the final S/N of the 1D spectrum (Naylor 1998) by ∼10%.

Before extracting the spectra, we applied a simple clipping procedure to identify and interpolate the bad pixels, cosmic ray, and sky residuals directly on the 2D frames of each single exposure for each arm. On average, we only interpolate ∼4% of the total pixels. We then performed the 1D extraction for each OB, in each arm separately, over an aperture of ±0.40″ centred on the photometric centre of the object (the pixel with the largest flux), which was chosen to maximise the S/N of the final spectra. We also ran the code MOLECFIT (Smette et al. 2015; Kausch et al. 2015) on the VIS arm⁹. Then, we obtained a final spectrum for each system in each arm, summing the spectra relative to different OBs, weighted for their corresponding errors. Finally, after extracting and summing, we joined together the UVB and VIS final spectra for each galaxy, as described in the next section.

3.2. UVB and VIS combination

To correctly infer the stellar population parameters, it is crucial to analyse a sufficiently large wavelength range. This is the only way to break the age-metallicity (Worthey 1994) and metallicity-[α/Fe] degeneracy (e.g., Spiniello et al. 2014) and study how these stellar population parameters correlate with galaxy masses (Spinrad 1962; Frogel et al. 1978; Faber & French 1980; Cenarro et al. 2003) or other morpho-photometric characteristics and environments (Trager et al. 2000, 2008; Trager 2006). For this reason, we need to combine the UVB and VIS arms, in order to simultaneously fit a larger number of stellar absorption features arising from different chemical species (e.g., CaK, CaH, G4300, Mg, Fe, Na, Balmer lines, TiO). However, since the two spectra from the two different arms of XSH have different nominal resolutions (R_VIS = 5000, R_UVB = 3200), a spectral convolution needs to be performed in order to bring the VIS and UVB spectra to the same final resolution, before combining the two sections. For this purpose, we first shifted the spectra to the rest-frame wavelength and then convolved each spectrum with a Gaussian function with a variable sigma (following the prescription of Cappellari 2017). The kernel of the Gaussian broadening function is computed as follows:

$$\begin{array}{c}\hfill \mathrm{kernel}(\lambda )=exp\left({\displaystyle \frac{-{\lambda }^2}{2\eta {(\lambda )}^2}}\right),\end{array}$$(1)

where λ is the wavelength measured in Å and η is defined as

$$\begin{array}{c}\hfill \eta (\lambda )={\displaystyle \frac{\sqrt{{\mathit{FHWM}}_{\mathrm{fin}}^2-{\mathit{FHWM}}_{\mathrm{XSH}}{(\lambda )}^2}}{2.3548}},\end{array}$$(2)

with a full width half maximum FWHM_XSH = λ/R_arm for each arm, and FWHM_fin = 2.51 Å at the restframe wavelength. This latter final value was chosen to be equal to that of the stellar templates and the single stellar population models used in the next section. The final combined 1D spectra are plotted in Fig. 4 in units of 10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹ and at rest-frame wavelength. In the figure, red dashed lines highlight the stellar absorption lines, typical of ETGs, which are used to derive stellar kinematics and population parameters in the next sections. Grey vertical shaded boxes show the regions where the UVB and VIS spectra were joined, while purple ones show the region contaminated by telluric lines that we corrected with MOLECFIT. To avoid possible problems of flux calibration at the boundaries of the two arms, we manually assigned a very large error to these pixels so that they would not be considered in the stellar template fitting procedure described in the next sections.

Fig. 4. Final 1D spectra extracted from an aperture of five pixels (±0.4″ from the photometric centre) for the three systems presented in this pilot paper, from the one at highest redshift (KiDS J0224−3143, top panel) to the one at the lowest redshift (KiDS J0847+0112, bottom panel). The spectra are plotted in rest-framed wavelength and fluxes (in units of 10−17 erg s−1 cm−2 Å−1). Red dashed lines highlight the stellar absorption lines, typical of ETGs, that are used to derive stellar kinematics and population parameters. Grey vertical shaded boxes show the regions where the UVB and VIS spectra were joined for each system. The purple boxes show regions contaminated by telluric lines that we corrected with MOLECFIT.

We infer the final S/N in the central part of the final spectra, from a region of 500 Å, centred around the Mg_b and Fe5270, Fe5335 absorption lines (5000–5500 Å rest frame). The resulting values are reported in Table 3. We note that, for two of the objects, the resulting S/Ns are slightly lower than expected. For the next semester of XSH observations (P105), we revised the integration times on the basis of the results obtained from these three spectra.

3.3. Stellar kinematics

The stellar velocity and velocity dispersion were derived using the Penalised Pixel-fitting software (PPXF; Cappellari & Emsellem 2004; Cappellari 2017). In particular, we used the latest python version of PPXF, v7.2.0¹⁰.

As suggested by the author of the code, to fit the kinematics we used an additive Legendre polynomial but did not use any multiplicative ones. This approach helps to correct the continuum shape during the fit, which is only performed to obtain the two first moments of the Gauss-Hermite series (velocity and velocity dispersion). The final chosen grade for the polynomial additive degree is 20, but we investigate the impact of a different choice in the Appendix A.

The fit was performed between 3500–6500 Å, where strong stellar absorption lines are present (CaK, CaH, Hγ, G-band, Hβ, Fe, Ca, Mg, NaD, see Fig. 4). During the fit, we carefully set the good-pixel region by hand, manually masking out the region where the UVB and VIS spectra were joined (i.e. the grey regions in Fig. 4) and assigning larger errors to the regions where the deepest telluric line in the optical range was removed by MOLECFIT (i.e. purple regions in Fig. 4). In addition, we also used a sigma clipping method similar to that described in Cappellari et al. (2002). In practice, we ran the fit for the first time and selected all the galaxy pixels that deviated more than 3 sigma, in terms of the RMS¹¹, from the best fitting template spectrum. We then repeated the fit a second time with the same parameters as before but masking out all the pixels selected during the previous step.

To fit the LOSVD, we did not apply any linear regularisation to the weights of the stellar templates (i.e. we set the corresponding keyword in PPXF REGUL = 0), which are instead tested to recover the stellar populations, as described in Sect. 3.4, and, in more detail, in Appendix C.

As stellar templates for fitting the velocity and velocity dispersion, we tested both single empirical stars and single stellar population (SSP) models. We thus verified that the kinematical results do not depend on the template sets used for the fit. In particular, as empirical stars, we employed the full MILES stellar library, which consists of 985 stars spanning a large range in atmospheric parameters (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011). The spectra cover the 3525–7500 Å range and have a spectral resolution of FWHM = 2.51 Ź². As SSPs, we instead used the MILES models by Vazdekis et al. (2015) with a bimodal stellar IMF with a fixed slope of Γ = 1.3 and BaSTI theoretical isochrones¹³. The stellar population parameters of the models, which are also used to study the stellar populations of the three galaxies, are described in more detail in the next section.

The result of the kinematics fit obtained with the SSP models is listed, for each galaxy, in Table 3. We note that, since the spectra are in the rest-frame wavelength, the velocity that we quote, which is the systemic velocity relative to the nominal redshift of the system, should be consistent with zero. In the same table, we also list the S/N calculated around the Mg_b, summing up the fluxes between 5000 and 5500 Å and dividing them by the sum of the errors over the same range and the reduced χ² of the fit. We do not list the results obtained with stars, but these are always consistent within a 1-sigma error with those obtained using SSP models.

Uncertainties on the estimated values for the LOS kinematics due to: i) a variation of the additive polynomial; ii) the use of different templates; and iii) re-scaling of the noise are taken into account, on top of the random errors given in output by a single fit. For this purpose, we developed a semi-automatic routine that repeats the fit multiple times, following a classical Monte Carlo perturbation. In particular, the routine simulates the observed spectrum 256 times, pixel by pixel, each time adding a random noise generated using a Gaussian distribution with the observed flux at that pixel as central value and the RMS noise as standard deviation. Then, we ran the fit on these 256 new spectra, varied the grade of the additive polynomial between 16 and 24 (±4 from the optimal value), and rescaled the noise by a constant factor in order to maintain the reduced χ′² between 0.9 and 1.1¹⁴. Hence, in this way, we obtained a distribution of velocity and velocity dispersion measurements with associated uncertainties, which are shown in Fig. 5. The final quoted values of the uncertainties listed in Table 3 were then computed as standard deviations. The central values that we report are instead the values obtained from the original observed spectrum using the best-fit configuration (minimum χ², ADEGREE = 20).

Fig. 5. Distribution of velocities and velocity dispersions computed repeating the kinematical fit 256 times, changing the parameters and setup, as described in the text, to estimate the uncertainties associated with velocity and velocity dispersion. The values reported in Table 3 are instead those obtained with the best-fit configuration on the original observed spectrum.

An estimate of the stellar velocity dispersion at the effective radius for two of these three objects was already presented in S20. Given the much lower S/N of the spectra used in that paper, the resulting uncertainties on the values were, of course, much larger than the ones inferred here. Nevertheless, the values are very close to each other for KiDS J0314−3215, and within 1.5σ error for KiDS J0224−3143. KiDS J0847+0112 has been targeted by the GAMA Survey. We downloaded the GAMA spectrum and calculated the velocity dispersion using PPXF with the same configuration. Also in this case, the result is in fair agreement and consistent within 2σ errors. The values of the velocity dispersion presented in S20 (and GAMA) from the low-quality spectra are listed in Table 3, along with the S/N and the ID of the objects.

3.4. Stellar populations

As already stated, our working hypothesis is that relics are galaxies that have been formed via a short and intense star formation ‘burst’ during the first phase of the formation scenario proposed for massive ETGs, and then have evolved passively and undisturbed ever since. As such, they should ideally not have any signs of recent star formation episodes and be well represented by a single, very old stellar population with super-solar α abundance and metallicity. The large [α/Fe] originates because the quenching of the star formation occurs before Type Ia supernova explosions can pollute the interstellar medium with iron (e.g., Gallazzi et al. 2006, 2014; Thomas et al. 2005).

The super-solar metallicity is expected, because these objects are relatively massive, and thus, following the mass-metallicity relation (e.g., Gallazzi et al. 2005), one would expect that they are populated mainly by metal-rich stellar populations.

Furthermore, PPXF allows one to extract distributions only in the 2D parameter space defined by age and metallicity, thus, for an estimate of [α/Fe], we took a different approach. We first inferred the [α/Fe] via line-index measurements and then ran the full-spectral fitting code using only models with the ‘right’ value of [α/Fe] for each system. In detail, we used the SPINDEX code from Trager et al. (2008) to calculate the line-index strengths of the Mg_b and the many iron (Fe) lines present in the wavelength range at our disposal. We then produced the classical Mg_b-Fe index-index plot, which allowed us to obtain an estimate of the [α/Fe] and a first guess of the total metallicity for the three objects. To minimize the error bars and obtain a constraint that was as secure as possible, and as independent as possible on other elemental abundances (e.g., [Ti/Fe], as for Fe5015), we took the mean of 24 different iron lines,¹⁵ creating the Fe₂₄ index that is plotted on the x-axis in Fig. 6. We overplotted the line strengths obtained from MILES SSP models with varying total metallicity and [α/Fe] and two different fixed values of the age (10 and 5 Gyr). Both galaxies and SSP spectra were smoothed to a common resolution of σ = 270 km s⁻¹ before the indices were computed.

Fig. 6. Line-index inference on the [α/Fe] abundance of the relic candidates. In this index-index plot, the classical Mgb index (Trager et al. 2008) is plotted against the Fe24 index, which is the mean of the 24 iron lines present in the spectral range we considered here (3500−6500). The grid represents the strengths of MILES SSP models with varying metallicity (dotted) and [α/Fe] (solid; we colour-coded – in black – the original models and – in green – the interpolated ones), as indicated in the plot but with a fixed age of 10 Gyr. To show the effect of a varying age, which does not affect the inference on [α/Fe], we also plot younger models with an age of 5 Gyr with less thickness. Black arrows in the top left of each panel indicate the direction in which the model grid will move if the corresponding parameters increase. Both galaxies and SSP line strengths are calculated after convolving the spectra to a final, common resolution of σ = 270 km s−1.

Given that the publicly available version of the MILES models only allows for two different [α/Fe] abundances (0.0 and 0.4), we linearly interpolated the models between these values (for each different age and metallicity values) to build a finer grid (Δ = 0.1). The interpolated models are plotted in green, while the original ones are plotted in black. We simply assumed the one of the closest models as the [α/Fe] value and assigned an error bar equal to the step between the models (0.1).

By making use of PPXF, we constrained, at this point, the mass-weighted stellar age and metallicity for the three relic candidates, after fixing the [α/Fe] to the value read from the index-index plot. We ran the fit using only SSP models with [α/Fe] = 0.2 for KiDS J0847+0112 and of 0.3 for the other two systems as templates, but, given the error bars from Fig. 6, we tested the impact that changing the [α/Fe] has on the kinematics and on the stellar populations in Appendix B.

This time, again following the prescriptions given within the code documentation, we remove the additive polynomial (ADEGREE = −1), which could affect the line strength of the spectral features, and use a multiplicative one instead (MDEGREE = 10, but we tested the impact of this assumption on the results in Appendix A).

As starting guess for velocity and velocity dispersion, we use the values obtained in Sect. 3.3, which are also listed in Table 3. The fitting results for KiDS J0224−3143, KiDS J0314−3215, and KiDS J0847+0112 are plotted in Fig. 7. The galaxy 1D spectrum is plotted (in rest-frame wavelength) in black, and the best-fit template is overplotted in red. The residuals of the fit are shown in green; blue lines are the pixels that have been masked out via our sigma-clipping routine, while grey shaded regions have been excluded from the fit by hand because they correspond to tabulated emission or sky lines (as highlighted on the top x-axis of the panel).

Fig. 7. Final results of stellar population PPXF fits on KiDS J0224−3143 (top row), KiDS J0314−3215 (middle) and KiDS J0847+0112 (bottom). Each panel shows the galaxy spectrum (black) with the best-fit stellar template (red) overplotted on it. The residuals for the good pixels are shown in green; blue lines are the pixels that have been masked out via our sigma-clipping routine, while grey shaded regions have been excluded from the fit by hand.

We note that the velocity dispersion values obtained with this new PPXF setting are slightly lower, but overall consistent (within 2σ uncertainties), with those given in the previous section (given in Table 3). The impact of performing the fit with α-enhanced SSP models is described in the Appendix B.

Furthermore, using the same semi-automatic routine already employed to estimate the errors on the kinematic parameters, we also obtain a realistic estimate of the standard deviation on the stellar population parameters. We repeated the fit 256 additional times on each spectrum (i.e. for each system) at every run, randomly simulating a new spectrum according to the 1σ RMS noise level obtained during the first fit, changing the MDEGREE from 6 to 14 (±4 from the chosen value) and the regularization from 0 to MAX_REGUL. The distribution of the different meaningful values of the age and [M/H] for each galaxy are plotted in Fig. 8. The standard deviation of each distribution is assumed as uncertainties on the values listed in Table 4. We note that we only used this procedure to compute the standard deviation on the central values of the mass-weighted stellar population parameters, which were instead calculated in the best possible configuration (producing the best χ′² and with MDEGREE = 10).

Fig. 8. Histograms of the distribution of mass-weighted mean ages and metallicities inferred when repeating the fit 256 times; changing the parameters, the set-up, and the regularization, as described in the text. These distributions were only used to estimate the uncertainties on the values reported in Table 4. The mean ages ([M/H]) of the distributions shown here are not identical to the values found with the best possible setup (listed in the table) but do not differ by more than 1 Gyr (0.04 dex).

We performed the stellar population fit in two different ways: once fixing the regularisation parameter to null (REGUL = 0), and a second time considering instead the maximum regularisation still consistent with the galaxy spectrum and its S/N (MAX_REGUL). Hence, with MAX_REGUL, the derived star formation history (SFH) is the smoothest (minimum curvature or minimum variation) one still consistent with the spectral observations. On the other hand, with REGUL = 0, the derived SFH gives a much larger weight to a much smaller number of templates, thus possibly discarding secondary star formation episodes. Considering these two extreme cases, is therefore important to interpret the inferred SFH and draw conclusions on the ‘relicness’ of the three galaxies, which is discussed in Sect. 4. Finally, in Appendix C, we describe tests on more values for the REGUL parameter and study the impact that this has on the derived stellar population parameters and on the SFH. To set the MAX_REGUL, we used a total number of SSP models ranging from 90 to 110, depending on the galaxy, with ages ranging from 1 to the age of the Universe at the redshift of each object (9.5, 10.5, and 11.5 for KiDS J0224−3143, J0314−3215, and J0847+0112, respectively) with steps of Δt = 0.5 Gyr and five different total metallicity bins (in terms of the star’s total metal abundance, [M/H] = { − 0.66, −0.25, +0.06, +0.26, +0.04}. For each galaxy, we only used models with the [α/Fe] value we estimated from Fig. 6. The MAX_REGUL found for each galaxy is reported in Table 4. We note that the spectrum of KiDS J0314−3215 allows for a much smoother solution (MAX_REGUL = 600) and thus requires a more extended SFH, with respect to the other two objects. We comment further on this in the next section.

4. Results

The final results of the PPXF stellar population fitting are reported in Table 4, both for the unregularised (REGUL = 0) case and that with the maximum regularisation (REGUL = MAX_REGUL). In the table, for each system we list the mass-weighted stellar population parameters with associated uncertainties and also the formation redshifts (z_form) for both cases. Formation redshifts were obtained with ΛCDM general cosmology in the following way. First, we calculated the age of the Universe at the redshift of each galaxy, t_uni(z_gal), then we subtracted the age of the Universe at the formation redshift, t_uni(z_form), and made it equal to the mass-weighted age estimated from the SSP fitting (also listed in the table):

$$\begin{array}{c}\hfill t_{\mathrm{uni}}(z_{\mathrm{gal}})-t_{\mathrm{uni}}(z_{\mathrm{form}})=t_{\mathrm{SSP}}.\end{array}$$(3)

Finally, we solve for z_form.

The stellar population parameters derived for the three relic candidates analysed here confirm that these objects are indeed mainly populated by old stars. This is also fully consistent with our previous findings, based on photometric colours (see Fig. 1) and low-resolution spectroscopy (T18; S20)¹⁶.

4.1. Relic confirmation

In order to assess whether there systems are relics or not, in Fig. 9, we plot the cumulative stellar mass formed at each time, starting from the Big Bang (BB). We do so for the main setups and assumptions we made in the paper (i.e. for the two extreme regularisation values, assuming the [α/Fe] as computed via line-index analysis or taking the two extreme [α/Fe] values permitted by the SSP models.

Fig. 9. Star formation histories obtained with PPXF from the spectra of the three objects extracted within an aperture of radius r = 0.4″, corresponding to five XSH pixels. The time x-axis corresponds to the time since the BB, and the vertical solid line in each panel shows the age of the Universe at the redshift of the corresponding object. Black lines show the cumulative stellar mass fraction inferred from the unregularised mass-weighted PPXF fit with the best value of [α/Fe] estimated from line indices (Fig. 6). Green lines instead show the results of the fit performed with the regularization set to the MAX_REGUL. Finally, coloured dashed lines show the impact of performing the fit (with REGUL = 0) using models with different [α/Fe].

From the figure, it is clear that KiDS J0847+0112 (bottom panel) has the most ‘extreme’ star formation history, when considering all the possible sources of uncertainties due to the different configurations and setups. It formed the totality of its stellar mass at most 2 Gyr after the BB, most likely even before and over very short timescales (∼0.5–1.0 Gyr), and then it had no further star formation episodes.

The star formation history of KiDS J0224−3143 (top panel) is not that different from that of KiDS J0847+0112 for the best configuration. However, it is fair to highlight that considering [α/Fe] = 0, which is however very unlikely given the estimate obtained from the index-index plot, only ∼70% of the stars were assembled during the first formation episode.

Finally, for KiDS J0314−3215 (middle panel), in the case of an unregularised fit, one main star formation episode formed 80% of the stars ≤0.5 Gyr after the BB, and then the assembly was completed with another fast (1 Gyr long) episode that started 2 Gyr after the BB. Considering the unlikely configuration with [α/Fe] = 0, the second star formation burst occurred instead only 4 Gyr after the BB. However, this galaxy has the highest MAX_REGUL (600), and thus a much smoother SFH is compatible with its spectrum too. In the MAX_REGUL case, the galaxy completed its assembly (i.e. ≥95% of the stellar mass in place) only ∼5 Gyr after the BB, and therefore we cannot classify it as a relic. The formation redshift of this object is also lower than that of the other two, in both the unregularised (z_{form, J0314} = 4.65) and the regularised (z_{form, J0314} = 2.64) configurations. However, we note that the majority of the stellar mass of this galaxy (80%) is characterised by an old and metal rich population. Thus, we speculate that although this object cannot be classified as a relic (i.e. having only one short star formation burst early on in cosmic time), it did not go through a numerous series of interactions during the second phase of its formation scenario. This also explains why it is so compact. Finally, for KiDS J0314−3215, we derived a velocity dispersion value lower than that inferred for the other two objects. We conclude that this system might therefore be less massive than we thought, and thus it might even not classify as a UCMG any longer (not passing the mass threshold).

Considering the most plausible setups and assumptions, KiDS J0224−3143 and KiDS J0847+0112 are characterised, for at least 95% in mass, by stellar populations that are only 1–2 Gyr younger than the age of the Universe at their redshifts (at z = 0.3841, t_Uni ∼ 9.5 Gyr and at z = 0.1764, t_Uni ∼ 11.5 Gyr). They formed through only one star formation burst at very high redshifts (z_{form, J0224} = 10.65; z_{form, J0847} = 9.20, inferred from the mass-weighted median ages calculated by PPXF for the unregularised case and z_{form, J0224} = 6.02, z_{form, J0847} = 6.04 for the regularised one) and have a large stellar velocity dispersion, a super-solar metallicity and super-solar α abundances. KiDS J0314−3215, although also mainly characterised by an old stellar population, has experienced at least two star formation episodes and reached 95% of the stellar mass ∼4.5 Gyr after the BB.

We will provide a more rigorous definition of ‘relicness’ when we analyse the full INSPIRE dataset, in future publications. In summary, we can conclude that two out of three galaxies are classified as relics in the sense that they only had one SF episode, very early on in cosmic time, and they had already assembled 95% of their total stellar mass 2 Gyr after the BB.

4.2. Comparison with the three local relics of F17

Here we briefly qualitatively compare our findings on the first two confirmed relics in the nearby Universe (0.15 <  z <  0.4) with the properties of the three local systems confirmed and studied by F17. As one can see from Fig. 2, in a stellar mass-size plane, the INSPIRE UCMGs (green points) and the three local relics of F17 (yellow triangles) lie on the same region, although our objects are on average smaller but slightly less massive.

In terms of morphology, the two new confirmed relics show a compact disk-like shape (also in line with the axis ratio inferred from Roy et al. 2018, q = 0.39, 0.38, 0.25 respectively for KiDS J0224−3143, J0314−3212, and J0847+0112), very similar to that of massive red nuggets at high-z (z ∼ 2, Buitrago et al. 2008; van der Wel et al. 2011) and of the local relics. KiDS J0224−3143, however, looks slightly more roundish from the KiDS high spatial resolution r-band image shown in Fig. 3. This is also the galaxy with the larger Sérsic index (n = 6.5), although all three objects have very peaked profiles with n >  4. This galaxy, according to the stellar mass of M_⋆ = 2.7 × 10¹¹M_⊙, estimated in T18, is the most massive relic ever found.

The velocity dispersion values we inferred for the new nearby relics are lower than the ones reported in F17 for the three local objects. We speculate that this might be due to the fact that, given the smaller apparent sizes on sky of these objects, we are integrating over larger physical radii. If the velocity dispersion profiles are centrally peaked and then decrease with radius, as, for example, in the case of NGC1277 (Trujillo et al. 2014), which is the most extreme relic in F17, we expect to measure a smaller dispersion since we cover a larger spatial aperture for the galaxies. In fact, looking at Fig. 2 in F17, one can see that for r ∼ 0.5−1 R_e and beyond, the stellar velocity dispersion values are measured to be around σ_⋆ ∼ 250 km s⁻¹, thus perfectly consistent with the values we find for KiDS J0224−3143 and KiDS J0847+0112. KiDS J0847+0112, for which we infer a velocity dispersion ∼50 km s⁻¹ lower, might still be consistent if it has a very peaked sigma profile similar to that of NGC1277, but it could also simply be a lower mass galaxy.

Finally, looking at the derived SFH, plotted in Fig. 9, we seem to observe a ‘degree of relicness’, which is in line with what was already reported by F17. The SFH of KiDS J0847+0112 is, in fact, more ‘extreme’, especially when considering all the possible setups and uncertainties in the [α/Fe] estimate.

At the same time, in this sense, even KiDS J0314−3215 might be considered as a less extreme relic, since it is nevertheless dominated by very old stars (80% of the stellar mass rather than 95%). This means that this galaxy might have experienced little interaction after the first mass assembly episode, but this ‘accreted’ component comprises at most 20% of its stellar mass¹⁷. Further investigation is required to understand the implication of this on the size growth and time evolution of this object.

In conclusion, as in F17, we cannot make any secure statement on the basis of only a few objects (five in total). We therefore stress the need for a larger galaxy sample, which is one of the main goals of the INSPIRE survey. Hence, we will further investigate the existence of the ‘degree of relicness’, especially in connection with the environment and other structure parameters (e.g., size or colours) and assess, in this way, the path that red nuggets undergo to become normal ETGs.

5. Conclusions and future perspectives

Relic galaxies, the untouched survivors of the massive and compact high-z red nuggets, provide a unique opportunity to track the formation of the in situ component formed at high-z during the first phase of the formation scenario, which is mixed with the accreted one in normal massive ETGs. Since studying in detail the stellar population of a statistically large sample of high-z red nuggets would require very long integration times, relics are ideal laboratories that allow us to investigate the physical processes that shaped the mass assembly of massive galaxies in the high-z Universe. Hence, the conditio-sine-qua-non to build a complete and comprehensive model for the star formation and cosmic evolution of massive galaxies in the Universe is enlarging the sample of relic galaxies for which stellar population properties, including the IMF, can be studied in great detail. This is one of the main goals of the INSPIRE project, which we introduced in this paper.

Further to introducing the new survey, as well as its current status and the observing strategy, we presented here results based on three systems, representative of the whole sample. KiDS J0224−3143 is the most massive UCMG in our dataset, KiDS J0314−3215 is one of the most compact objects, and KiDS J0847+0112 is one of the objects at the lowest redshift.

In particular, in this paper we have:

• (i)

  Presented the combined final XSH UVB+VIS spectrum for each object, also correcting for the different resolutions of the two arms¹⁸;

• (ii)

  Calculated the values of the integrated (from a five-pixel aperture) stellar velocity and velocity dispersion with the PPXF code, also testing the impact of the fitting setup on the final results;

• (iii)

  Inferred the [α/Fe] abundance by comparing line-index strengths of Mg and Fe lines in the galaxies with those of Vazdekis et al. (2015) SSP models with a Kroupa IMF, a range of ages, metallicities, and [α/Fe] abundances (interpolated to cover a finer grid);

• (iv)

  Inferred mass-weighted mean age and [M/H], performing full-spectral fitting of the spectra in the 3500−6500 Å wavelength range and presented results for the two possible extreme cases in terms of smoothness of the inferred SFH;

• (v)

  Estimated uncertainties on the kinematics and stellar population parameters repeating the fit multiple times, and changing at each run the noise level, the degree of the polynomial used to adjust the continuum, and the good-pixel region;

• (vi)

  Confirmed two of the three galaxies as relics, proving that they formed the majority (> 95%) of their mass at very high redshifts (z ≥ 6), and, in addition, they are also characterised by a super-solar metallicity and α abundance, similarly to the three confirmed local relics of Ferré-Mateu et al. (2017).

In conclusion, we confirmed two of the three galaxies analysed here as relics. If we extrapolate this result to the whole INSPIRE sample, this means that we might confirm ∼27 UCMGs as relics and thus increase by a factor of ∼8 the number of known relic galaxies, also extending the sample to relatively higher redshifts.

Once the full INSPIRE sample is fully reduced and analysed, we will be able to compute precise number density of relics in the redshift window 0.1 <  z <  0.5 and compare our findings with the predictions from state-of-the-art hydrodynamical cosmological simulations. This will allow us to study, in detail, the mechanisms responsible for the dramatic size growth measured from z ∼ 1 to the present-day Universe and thus disentangle the possible merger scenarios. We will also investigate the formation scenario and mass assembly of these extremely compact and massive passive objects, for both relics and younger UCMGs. In particular, the presence of a conspicuous number of young very compact and very massive galaxies will seriously challenge current hydrodynamical simulations that predict the formation of such objects only in gas-rich environments of the young, high-redshift Universe (e.g., Zolotov et al. 2015). Thus, a measurement of the relative abundance of post-star-forming and old compact galaxies and its evolution with redshift will provide a key test for any models of massive galaxy formation.

In addition, by inferring the stellar population parameters, including the IMF slope at the low-mass end (at least for the objects with highest S/N spectra) of relics, we will investigate the mechanisms that regulated the star formation at high-z, when the density and temperature of the Universe were much higher. According to theoretical work (e.g., Chabrier et al. 2014), high density, temperature, and turbulence of the gas are key parameters that drive the fragmentation of molecular clouds. Higher density and a higher temperature make fragmentation easier, forming more dwarf stars, hence a bottom-heavier IMF.

Finally, measuring the stellar velocity dispersion from spectroscopic data and stellar masses via stellar population analysis will allow us to compute the ratio between stellar and dynamical mass (M_⋆/M_dyn) for individual objects, a potentially valuable tracer of the likely mechanism by which galaxies grow (e.g., Hopkins et al. 2009; Hilz et al. 2013). Specifically, if galaxies experience merger-driven growth, the ratio measured within the effective radius should decrease with time (Tortora et al. 2014, 2018a); this decrease would be stronger in the case of minor mergers (van de Sande et al. 2013). Hence, under the hypothesis that relics experienced very little or no mergers at all, they should all show a very high M_⋆/M_dyn, and thus a small amount of dark matter within the effective radius. However, for one for the confirmed relics in the local Universe, Mrk 1216, Buote & Barth (2019) found the opposite from Chandra X-ray observations: the internal (within R_e) DM fraction of this object is two times higher than the typical fractions found in nearby massive normal-sized ETGs.

The only way to investigate further is to enlarge the sample of confirmed relics for which these kinds of studies can be performed, also pushing the redshift boundaries and finding systems living in different environments. This is indeed what we will achieve thanks to INSPIRE.

Acknowledgments

The authors thank the referee for a very constructive report which helped improving the final quality of the manuscript. CS is supported by a Hintze Fellowship at the Oxford Centre for Astrophysical Surveys, which is funded through generous support from the Hintze Family Charitable Foundation. C.S. is also very grateful to Ortwin Gerhard and his ‘Dynamics Group’ at the Max-Planck-Institut für Extraterrestrische Physik (MPE, Garching by Munich) for interesting and constructive discussions. C.T., A.G., L.H. and S.Z. acknowledge funding from the INAF PRIN-SKA 2017 programme 1.05.01.88.04. G.D. acknowledges support from CONICYT project Basal AFB-170002. AFM has received financial support through the Postdoctoral Junior Leader Fellowship Programme from La Caixa Banking Foundation (LCF/BQ/LI18/11630007). N.R.N. acknowledges financial support from the One hundred top talent programme of Sun Yat-sen University, Grant N. 71000-18841229. D.S. is a member of the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne.

References

Appendix A: Testing the effect of changing the grade of the polynomials in the fits

As highlighted in the main body of the paper, we performed two different PPXF fits, one for the kinematics and one for the stellar populations, with different configurations of the program. In particular, we used an additive Legendre polynomial to correct for the continuum shape when inferring the velocity and velocity dispersion, while we used a multiplicative polynomial when calculating age and metallicity. In this Appendix, we motivate the choice of the polynomial degree that we use in the main text in both cases, and we also demonstrate the robustness of the presented results against the use of different values, within a wide range of possibilities.

For the grade of the additive polynomial (ADEGREE), we test values from 4 to 29, with steps of 1 and finally chose ADEGREE = 20, which stabilises the results against changes of other parameters (e.g., stellar templates, good-pixel region, see the text for more details) while contemporarily minimising the reduced χ′². A demonstrative plot of the effect of the additive polynomial degree on the velocity dispersion is provided in Fig. A.1.

Fig. A.1. Velocity dispersion obtained as function of the additive polynomial degree in PPXF. We show here the results obtained using α-enhanced SSP, and we refer the reader to Appendix B to see the effect that changing [α/Fe] has on the resulting velocity dispersion values. The results stabilise, for all the galaxies around ADEGREE ∼ 20, which is the final choice we adopted for the results given in Table 3.

For the multiplicative polynomial (MDEGREE) used to infer the stellar population parameters in Sect. 3.4, we tested degrees from 4 to 25 and finally chose MDEGREE = 10. The mass-weighted mean ages and metallicities derived by varying the grade of the multiplicative polynomial are shown in Fig. A.2. For all the objects, the results are on average stable for 8 < MDEGREE < 20. Above this degree, the mean age of KiDS J0224−3143 becomes ∼0.5 Gyr younger, which is within the errors. We note that we report the results of this test for the REGUL = 0 case, but we conducted the same investigation for REGUL = MAX_REGUL, finding fully consistent results.

Fig. A.2. Mass-weighted mean metallicity (top) and age (bottom) as function of the multiplicative polynomial degree in PPXF for each galaxy. The results stabilise for MDEGREE ∼ 10 (until ∼20), which is the final choice we adopted for the results given in Table 4.

Appendix B: Testing the effect of changing the [α/Fe] values in the full-spectral fit

In the main body of this paper, we used two different approaches to derive the stellar population parameters. First, we inferred the [α/Fe] from line-index measurements (see Fig. 6), and then we obtained age and metallicity by performing a full spectral fitting with SSP models with the previously inferred α abundances, as templates. The line strength inference, although solid, might suffer from degeneracy between variation of [α/Fe] and variation of single elemental abundances (e.g., Ti/Fe, which contaminates many optical iron lines). Hence, here in this Appendix, we test the impact of using models with different [α/Fe] on the results presented in the paper. For this purpose, for each galaxy we repeated the PPXF fit in 2D (varying age and metallicity) five times, each time using models with different α abundances (from 0.0 to 0.4 with a step of 0.1, see main body). We did that for both the REGUL = 0 and REGUL = MAX_REGUL cases. The impact of a varying [α/Fe] on the stellar velocity dispersion is shown in Fig. B.1 (the regularisation does not influence the σ_⋆ values), while the impact on the stellar age and metallicity can be visualised in Fig. B.2. Increasing [α/Fe], the velocity dispersion and the metallicity decrease, while the inferred population gets older. However, this effect is much less evident for KiDS J0847+0112, the most extreme relic, and it is much larger for KiDS J0314−3215, which is the object with the most extended SFH. For KiDS J0224−3143, if one assumes a [α/Fe] = 0, then the inferred mass-weighted age becomes younger (∼6 Gyr), confirming what was already found in Fig. 9. However, as already stated, this result would contradict the results inferred from the index-index plot (Fig. 6). This test confirms our result that KiDS J0847+0112 and KiDS J0224−3143 are indeed relic galaxies.

Fig. B.1. Impact on stellar velocity dispersion of varying the [α/Fe] abundances of the SSP models used for the PPXF fit. Different colours indicate the results for the three different objects, as highlighted by the legend.

Fig. B.2. Impact on stellar population parameters of varying the [α/Fe] abundances of the SSP models used for the PPXF fit. The top panel shows the variation on the inferred (mass-weighted) metallicity, whereas the bottom panel shows the variation in the (mass-weighted) stellar age. Filled symbols refer to the fits performed with REGUL = 0, while empty symbols refer to the regularised fits. Different colours indicate the results for the three different objects, as highlighted by the legend.

Appendix C: Testing the effect of different regularisation values on the stellar populations

In order to give a meaningful interpretation of the template weights, in terms of the SFH and metallicity distributions of individual galaxies, a regularised fit should be tested. However, it is not trivial to properly set the REGUL parameter (see e.g., Sect. 4.2 in Shetty & Cappellari 2015). For this reason, in this Appendix we undertake an investigation into the impact of different choices for the regularisation parameter on the results. In particular, following the prescriptions described in the PPXF documentation, we:

– Performed a number of unregularised fits (REGUL = 0), each time rescaling the input noise spectrum until obtaining a reduced χ′² (χ² divided by the number of fitted pixels, N_{good pixels}) as close as possible to 1. All the other parameters for the fit were kept fixed (fitted region: 3500−6500, ADGREE = −1);

– Performed a number of regularised fits, with the noise level set before, each time increasing the REGUL parameter until the new χ² increases by $\mathrm{\Delta }{\chi }^2=\sqrt{(2\times N_{\mathrm{good}\mathrm{pixels}}}$). In the remainder of this Appendix, we refer to the value of the REGUL keyword that optimises the χ′² as MAX_REGUL since, as already mentioned in Sect. 3.4, it produces the smoothest possible star formation history.

For each galaxy, we tested five different values of REGUL, between 0 and MAX_REGUL (specifically i × MAX_REGUL, with i = [0.00, 0.25, 0.50, 0.75, 1.00]). As expected, when increasing the regularisation, the fit uses an increasing number of SSP templates, and the resulting solution, in terms of age and metallicity, is smoothed. However, this only causes very minor changes for the two newly confirmed relics, whereas it affects the results much more in terms of inferred mean ages and metallicities for KiDS J0314−3215. This can be visualised from Fig. C.1, where we plot, for each of the galaxies (each column) in the form of a 2D density panel, the distribution of the weights (dark blue = 0, light yellow = maximum weight) that the PPXF attributes to the SSP templates for the five cases with increasing regularisation (top row REGUL = 0, bottom row REGUL = MAX_REGUL). No matter which regularisation value is used for KiDS J0847+0112, the fit uses only very old and metal-rich models. The mean values are almost identical for all cases. For KiDS J0224−3143, a young and metal-poor population (age ∼ 1.5 Gyr) pops up when increasing the REGUL parameter. However, this contributes very little to the mass budget of the galaxy, as the inferred mean age is always consistently old. For KiDS J0314−3215, instead, increasing the smoothness of the fit, the distribution of the weights gets much more broad, and younger models with lower metallicities are used in the fit too. The inferred age is, in this case, roughly 1.5 Gyr younger.

Fig. C.1. 2D density plot showing the weights attributed to the SSP models by PPXF during the fit, in a log (age), [M/H] space (for simplicity, we also plot ages in linear scale on the top x-axis). Each column shows results for a different galaxy, while each row corresponds to the different results obtaining increasing the REGUL from 0 (top) to MAX_REGUL (bottom). The weight fractions are reported on the vertical colour-bars, where yellow indicates the maximum weight and dark blue indicates zero weight. The upper limit of the age axis is slightly higher (∼10%) than the age of the Universe at the redshift of the objects, in order to take into account possible uncertainties in the value computed from the stellar population analysis.

All Tables

All Figures

Fig. 1. g − i colour plotted versus redshift for the spectroscopically confirmed UCMGs in KiDS (∼100, T18; S20). The relic candidates, targets of the XSH INSPIRE LP, are plotted as green points. The three objects, which are the focus of this study, are plotted as magenta diamonds. Blue points are UCMGs from T18 and S20 consistent with younger stellar ages. Two out of the three confirmed relics of F17 are plotted as yellow triangles (colours are taken from the SDSS DR16). Lines show single-burst old stellar population models (Bruzual & Charlot 2003) with solar (Z = 0.02, solid lines), sub-solar (Z = 0.008, dotted), and super-solar (Z = 0.05, dashed) metallicity and ages of 12 Gyr (black) and 8 Gyr (red).

In the text

Fig. 2. Stellar mass-size relation for high-z nuggets (black points, from Szomoru et al. 2012), the F17 local relics (yellow triangles), and the relic candidates targeted by INSPIRE (green points and magenta diamonds), compared to the relation for normal SDSS ETGs (grey line, from Shen et al. 2003). By selection, all UCMGs are clear outliers of the local stellar mass-size relation, and high- and low-z systems lie on the same region of the diagram.

In the text

	Fig. 3. Cutouts of 1′  ×  1′ of the KiDS r-band images for the three relic candidates that we study in this paper. Spatial scale and image orientation are common to all the cutouts and are shown in the first and third panels, respectively.
In the text

Fig. 4. Final 1D spectra extracted from an aperture of five pixels (±0.4″ from the photometric centre) for the three systems presented in this pilot paper, from the one at highest redshift (KiDS J0224−3143, top panel) to the one at the lowest redshift (KiDS J0847+0112, bottom panel). The spectra are plotted in rest-framed wavelength and fluxes (in units of 10−17 erg s−1 cm−2 Å−1). Red dashed lines highlight the stellar absorption lines, typical of ETGs, that are used to derive stellar kinematics and population parameters. Grey vertical shaded boxes show the regions where the UVB and VIS spectra were joined for each system. The purple boxes show regions contaminated by telluric lines that we corrected with MOLECFIT.

In the text

	Fig. 5. Distribution of velocities and velocity dispersions computed repeating the kinematical fit 256 times, changing the parameters and setup, as described in the text, to estimate the uncertainties associated with velocity and velocity dispersion. The values reported in Table 3 are instead those obtained with the best-fit configuration on the original observed spectrum.
In the text

Fig. 6. Line-index inference on the [α/Fe] abundance of the relic candidates. In this index-index plot, the classical Mgb index (Trager et al. 2008) is plotted against the Fe24 index, which is the mean of the 24 iron lines present in the spectral range we considered here (3500−6500). The grid represents the strengths of MILES SSP models with varying metallicity (dotted) and [α/Fe] (solid; we colour-coded – in black – the original models and – in green – the interpolated ones), as indicated in the plot but with a fixed age of 10 Gyr. To show the effect of a varying age, which does not affect the inference on [α/Fe], we also plot younger models with an age of 5 Gyr with less thickness. Black arrows in the top left of each panel indicate the direction in which the model grid will move if the corresponding parameters increase. Both galaxies and SSP line strengths are calculated after convolving the spectra to a final, common resolution of σ = 270 km s−1.

In the text

Fig. 7. Final results of stellar population PPXF fits on KiDS J0224−3143 (top row), KiDS J0314−3215 (middle) and KiDS J0847+0112 (bottom). Each panel shows the galaxy spectrum (black) with the best-fit stellar template (red) overplotted on it. The residuals for the good pixels are shown in green; blue lines are the pixels that have been masked out via our sigma-clipping routine, while grey shaded regions have been excluded from the fit by hand.

In the text

Fig. 8. Histograms of the distribution of mass-weighted mean ages and metallicities inferred when repeating the fit 256 times; changing the parameters, the set-up, and the regularization, as described in the text. These distributions were only used to estimate the uncertainties on the values reported in Table 4. The mean ages ([M/H]) of the distributions shown here are not identical to the values found with the best possible setup (listed in the table) but do not differ by more than 1 Gyr (0.04 dex).

In the text

Fig. 9. Star formation histories obtained with PPXF from the spectra of the three objects extracted within an aperture of radius r = 0.4″, corresponding to five XSH pixels. The time x-axis corresponds to the time since the BB, and the vertical solid line in each panel shows the age of the Universe at the redshift of the corresponding object. Black lines show the cumulative stellar mass fraction inferred from the unregularised mass-weighted PPXF fit with the best value of [α/Fe] estimated from line indices (Fig. 6). Green lines instead show the results of the fit performed with the regularization set to the MAX_REGUL. Finally, coloured dashed lines show the impact of performing the fit (with REGUL = 0) using models with different [α/Fe].

In the text

Fig. A.1. Velocity dispersion obtained as function of the additive polynomial degree in PPXF. We show here the results obtained using α-enhanced SSP, and we refer the reader to Appendix B to see the effect that changing [α/Fe] has on the resulting velocity dispersion values. The results stabilise, for all the galaxies around ADEGREE ∼ 20, which is the final choice we adopted for the results given in Table 3.

In the text

	Fig. A.2. Mass-weighted mean metallicity (top) and age (bottom) as function of the multiplicative polynomial degree in PPXF for each galaxy. The results stabilise for MDEGREE ∼ 10 (until ∼20), which is the final choice we adopted for the results given in Table 4.
In the text

	Fig. B.1. Impact on stellar velocity dispersion of varying the [α/Fe] abundances of the SSP models used for the PPXF fit. Different colours indicate the results for the three different objects, as highlighted by the legend.
In the text

Fig. B.2. Impact on stellar population parameters of varying the [α/Fe] abundances of the SSP models used for the PPXF fit. The top panel shows the variation on the inferred (mass-weighted) metallicity, whereas the bottom panel shows the variation in the (mass-weighted) stellar age. Filled symbols refer to the fits performed with REGUL = 0, while empty symbols refer to the regularised fits. Different colours indicate the results for the three different objects, as highlighted by the legend.

In the text

Fig. C.1. 2D density plot showing the weights attributed to the SSP models by PPXF during the fit, in a log (age), [M/H] space (for simplicity, we also plot ages in linear scale on the top x-axis). Each column shows results for a different galaxy, while each row corresponds to the different results obtaining increasing the REGUL from 0 (top) to MAX_REGUL (bottom). The weight fractions are reported on the vertical colour-bars, where yellow indicates the maximum weight and dark blue indicates zero weight. The upper limit of the age axis is slightly higher (∼10%) than the age of the Universe at the redshift of the objects, in order to take into account possible uncertainties in the value computed from the stellar population analysis.

In the text