© The Authors 2026

1. Introduction

The James Webb Space Telescope (JWST) has enabled the discovery of a particularly enigmatic class of galaxies: a large population of red and compact objects referred to as little red dots (LRDs, Matthee et al. 2024; Kocevski et al. 2023, 2025; Labbe et al. 2025). Initially presumed to be purely active galactic nuclei (AGNs) due to their broad Balmer lines, with FWHM ≈ 2000 − 4000 km s⁻¹, associated with the broad line region (BLR) (Kocevski et al. 2023; Greene et al. 2024; Kocevski et al. 2025; Taylor et al. 2025a; Wang et al. 2024), hints of a non-negligible stellar component in the observed flux began to challenge our understanding of these sources. Indeed, the spectral energy distribution (SED) flattens in the rest-frame near infrared, and a stellar component has been suggested to better explain the rest-frame-observed SED instead of a pure AGNs torus (Williams et al. 2024; Pérez-González et al. 2024). Several definitions of LRDs have been employed since then, each focusing on red and compact objects. In the literature, they can be selected based on their broad Balmer lines (Kokorev et al. 2023; Kocevski et al. 2023; Matthee et al. 2024; Furtak et al. 2024) or through color-color selections. Among these, the criterion that selects a rest-frame UV excess and a red optical slope, forming the V shape is used (Kocevski et al. 2025; Leung et al. 2025; Euclid Collaboration: Bisigello et al. 2026), as well as more direct color-color diagrams that identify objects with red and blue colors (Labbé et al. 2023; Kokorev et al. 2024; Akins et al. 2025b; Barro et al. 2024a; Xiao et al. 2025; Ma et al. 2025). All selection criteria share a common requirement that the compactness of the galaxies, typically defined in the observed 4.4 μm band as F_F444W(0.4″)/F_F444W(0.2″) < 1.7 (Labbe et al. 2025; Akins et al. 2025b; Kokorev et al. 2024; Greene et al. 2024; Barro et al. 2024a; Xiao et al. 2025). The half-light radius can also serve as a compactness criterion: R_e < 1.5 r_h,star, where r_h,star is the effective half-light radius of stars (Kocevski et al. 2025).

When the V shape defined on the UV and optical slope by Kocevski et al. (2025) is observed, the inflection point appears to occur at the same rest-frame wavelength as the Balmer break, which is known to be associated with old stellar populations (Setton et al. 2024). Moreover, if LRDs are assumed to be pure AGN, the accretion rate is expected to vary over time, leading to variable flux (Ulrich et al. 1997; MacLeod et al. 2012), which has not yet been observed in the UV or X-ray bands (Kokubo & Harikane 2025; Zhang et al. 2025; Tee et al. 2025). Additionally, the broadening of the lines can also be explained by the compact kinematics of a star-forming galaxy (Pérez-González et al. 2024; Akins et al. 2025b; Baggen et al. 2024). However, the first instance of time variability in an LRD was detected in A2744-QSO1 (Furtak et al. 2025) by measuring the equivalent width (EW) of H_α and H_β. Approximately 20% variability was observed over a 2.4-year period, confirming that the flux is dominated by an AGN, at least for this particular LRD. On the one hand, Xiao et al. (2025) attempted to detect [CII] lines to differentiate between an AGNs and dusty star formation for two V shaped LRDs at z_spec > 7. The non-detection of [CII] disfavors the only-stellar model. This is further supported by the lack of detection in sub millimeter observations with ALMA at 1.2 mm (Akins et al. 2025b; Casey et al. 2024; Labbe et al. 2025), as well as by the upper limit on the dust mass of 10⁶ M_⊙ estimated by Casey et al. (2025) in a sample of 60 LRDs.

Due to the variety of definitions, not all LRDs exhibit a clear V shape, broad lines, or a Balmer break. Approximately 60% of LRDs display broad lines (Kocevski et al. 2025; Greene et al. 2024; Loeb 2024), and > 70% of V shaped LRDs show this feature, suggesting a heterogeneous population governed by diverse physical processes. The degeneracy between the two interpretations – pure AGNs and pure stellar model – is well described in Akins et al. (2025b); they studied 434 LRDs as part of the COSMOS-Web survey (Casey et al. 2023) between z ≈ 5 − 9. If they are interpreted as purely AGN, the inferred accretion rate from the LRDs alone competes with that of all AGNs at cosmic noon, i.e., z ≈ 2, approximately 10⁻⁴ M_⊙ yr⁻¹ Mpc⁻³ (Yang et al. 2023a). The black hole mass inferred would also be overmassive with respect to the stellar host component (Pacucci et al. 2023; Maiolino et al. 2024; Durodola et al. 2025). If they are interpreted as purely stellar, inferred stellar masses are very important, above 10¹⁰ M_⊙, and JWST/MIRI is key to inferred robust masses consistent with Λ cold dark matter (CDM). The inferred density profile is broadly in agreement with the densest core observed in the universe up to z ≈ 3 by Hopkins et al. (2010), though slightly higher. The stellar density can even reach the threshold required for runaway stellar collisions to occur Guia et al. (2024).

In both cases, the physical properties of LRDs reach extreme values, suggesting that a combination of both models might be necessary. Significant Balmer break strengths have been observed in several LRDs (Naidu et al. 2025; de Graaff et al. 2025; Taylor et al. 2025b), approximately four times higher than those in massive quiescent galaxies. This strength challenges the stellar origin of the Balmer break in LRDs. As the usual point of view of AGNs is also challenged by previous studies, a new explanation of LRDs has been suggested. In this model, an LRD is interpreted as an AGNs surrounded by extremely dense gas (n_H ≈ 10¹¹ cm⁻³, N_H ≈ 10²⁵ cm⁻²), a scenario first proposed by Ji et al. (2025), Inayoshi & Maiolino (2025) to account for the Hβ absorption lines observed in several LRDs. This configuration also enables the presence of an extremely strong Balmer break, as developed in Ji et al. (2025), Naidu et al. (2025), Rusakov et al. (2025), de Graaff et al. (2025), Taylor et al. (2025b). In the BH^* model (so named by Naidu et al. (2025)), the UV excess is attributed to the host galaxy, which consists of a young stellar population with M_* ≲ 10⁹ M_⊙. This could explain why variations in the UV have not been detected, as they would not be related to AGNs activity. Additionally, dust is no longer needed to explain the red part of the SED, which is consistent with the non-detection in ALMA observations.

Another question that remains open to this date is the evolution of LRDs. Selection criteria used in the literature predominantly identify objects in the range 4 < z < 8, regardless of whether the selection method is based on rest-frame properties or not (Kocevski et al. 2025; Matthee et al. 2024; Barro et al. 2024a; Akins et al. 2025b; Kokorev et al. 2024; Labbe et al. 2025). Because the redshift distribution exhibits a drop at z ≈ 4 (Kocevski et al. 2025; Ma et al. 2025; Zhuang et al. 2025), even in rest-frame selections such as the V shape criterion used in Kocevski et al. (2025) and Euclid Collaboration: Bisigello et al. (2026), at least one of the selection criteria must vary with redshift. Consequently, known LRDs must evolve as a function of redshift. This implies that LRDs evolve into a different state that no longer meets the selection criteria. If this assumption holds, there must be a transition phase in which these galaxies still exhibit LRD features while simultaneously developing new characteristics.

Several scenarios can explain the evolutionary path of LRDs (see Section 6). However, recent observations suggest that some LRDs exhibit a resolved young stellar component, forming an extended structure beyond their initial compact core (Rinaldi et al. 2025; Chen et al. 2025; Torralba et al. 2026; Zhuang et al. 2025). Since this stellar component may evolve over cosmic time, we investigate in this paper a specific evolutionary pathway in which LRDs acquire a young stellar component–observed as a blue star-forming periphery–while retaining their compact red inner core. To explore this scenario, we select a sample of 55 galaxies at a median redshift of z_med = 3.5, analyze their properties, and compare them to the known LRD population. We examine their potential evolution across redshift, aiming to shed light on the observed decline in the typical LRD selection around z ≈ 4.

This paper is organized as follows: In Section 2, we present the JWST/NIRCam and HST/ACS data and the CEERS field, along with our sample selection from the catalog of Merlin et al. (2024), intended to be related to LRDs at lower redshifts. In Section 3, we present the methodology used to infer physical properties. The characteristics of the selected galaxies, including common features and differences with LRDs are presented in Section 4 and Section 5. We discuss different evolutionary scenarios for LRDs in Section 6 and conclude our study in Section 7. Throughout this work, we adopt an initial mass function (IMF) by Chabrier (2003) and the cosmological parameters from Planck Collaboration VI (2020). Magnitude used are in the AB system (Oke 1974).

2. Data and sample selection

2.1. CEERS data

In this work, we use JWST data from the Cosmic Evolution Early Release Science (CEERS) survey (Finkelstein et al. 2017, 2025), in the Extended Groth Strip (EGS) field. We specifically used the 10 JWST/NIRCam (Rieke et al. 2003, 2005; Beichman et al. 2012; Rieke et al. 2023) pointings that covered approximately 97 arcmin² of the sky. Seven filters are available: F115W, F150W, F200W, F277W, F356W, F410M, and F444W with a depth at 5σ for point sources of 29.15, 29.00, 29.17, 29.19, 29.17, 28.38, and 28.58 mag, respectively. We complete our short-wavelength coverage with the ACS/WFC instrument of the Hubble Space Telescope (HST) using the F606W and F814W filters (Davis et al. 2007; Ford et al. 1998), (Grogin et al. 2011; Koekemoer et al. 2011), with depths of 28.73 and 28.50 mag, respectively.

Concerning the observations and data reduction process, we refer the reader to Bagley et al. (2023). Briefly, the CEERS1, CEERS2, CEERS3, and CEERS6 pointings have been calibrated using the JWST Calibration Pipeline version 1.7.2 and CRDS pmap 0989. The remaining six pointings (CEERS4, CEERS5, CEERS7, CEERS8, CEERS9, and CEERS10) were calibrated using version 1.8.5 of the pipeline and CRDS pmap 1023. Then, all the pointings were reduced using stages 1 and 2 of the pipeline with additional steps such as removal of 1/f noise, wisps, and distortion correction. The final mosaics were aligned through stage 3 using astrometry data from Gaia-EDR3 (Gaia Collaboration 2021). Each pointing is background-subtracted and has a pixel scale of 0.03 arcsec pix⁻¹. All 10 pointings were then PSF-matched to the F444W-band image (0.16″ FWHM), see Gómez-Guijarro et al. (2023) for details. The total area of JWST/MIRI in the CEERS survey is ≈9 arcmin² (Yang et al. 2023b; Kirkpatrick et al. 2023). Since we target a population with a surface density of ≈0.6 − 1 arcmin² (Kocevski et al. 2025; Pérez-González et al. 2024), we choose to focus only on the JWST/NIRCam coverage to increase our statistics and have an homogeneous sample.

2.2. CEERS catalog and preselection

In this work, we use the catalog constructed by Merlin et al. (2024), which includes data from multiple deep JWST extragalactic fields, including the CEERS field. Briefly summarizing the detection and flux measurement process, this photometric catalog was designed to maximize object detection, particularly for very faint galaxies across all JWST fields, using SExtractor v2.8.6 (Bertin & Arnouts 1996). The detection image combines F356W and F444W bands, with an initial signal-to-noise (S/N) threshold of ≈2. Photometry was extracted using Kron elliptical apertures (Kron 1980) defined in the F356W+F444W detection image. The measured fluxes in all bands are then corrected using the aperture factor, determined from the detection image. A total of 82 547 objects were detected in the CEERS field, with available fluxes spanning the F606W to F444W filters and basic morphological parameters, including half-light radii (R_e), ellipticity, and position angle. Photometric redshifts were derived through SED fitting using zphot for one estimate (Fontana et al. 2000), and EAzY (Brammer et al. 2008) across three independent runs. For further details, see Merlin et al. (2024). Each source thus has four independent photometric redshift estimates.

From this catalog, we first applied some prefilters. First, all bands must be available: F606W, F814W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W (i.e., the error in each filter must be different from −99), which excludes sources at the very edges of the maps that are not measured in each filter. This removes 16 070 galaxies out of 82 547. A second threshold was applied to the signal-to-noise ratio, S/N_F444W > 10 and S/N_F356W > 10, to select well-detected objects (which retains 22 376 sources). A final criterion is applied to select bright objects in the reddest band (F444W < 27 mag), targeting massive objects, similarly to LRDs (Barro et al. 2024a). At this stage, our sample contains 18 218 objects. Throughout this work, the redshift used for each object corresponds to the median of all photometric redshifts reported in the catalog of Merlin et al. (2024) when the spectroscopic redshift is not available.

2.3. Sample selection

The number density of LRDs drops at z < 4 (Kocevski et al. 2025; Inayoshi 2025; Ma et al. 2025; Zhuang et al. 2025; Pacucci & Loeb 2025), indicating that detections of these objects become rare at lower redshift. This implies that the selection criteria used to identify them may no longer apply at lower redshifts, suggesting that some of the known properties of these galaxies are likely to evolve over time. As these properties change, new characteristics may emerge that could hinder their detection using current methods. Here, we aim to identify galaxies that exhibit a compact red inner core alongside a blue, star-forming periphery. If such an evolutionary pathway is plausible for the LRD population–where LRDs acquire a blue, star-forming outskirts–there should exist a transitional phase in which the overall flux of the galaxy remains dominated by the characteristic red color of LRDs. There are several color selections for the red part of LRDs in the literature: F277W-F444W > 1.5 (Akins et al. 2025b; Pérez-González et al. 2024) and F200W – F356W > 1.0 (Kokorev et al. 2024; Labbe et al. 2025) to select LRDs with z < 6. F200W-F444W is also used to probe the optical slope from z = 4.5 to z = 7 in Barro et al. (2024a) and Barro et al. (2024b).

Therefore, the bands typically used to identify the red color in LRDs span from F200W to F444W. To select lower-redshift objects exhibiting similar rest-frame red colors, we adopt the analogous color criterion F150W – F356W > 1.5. Indeed, at 4 < z < 5, where the number density of LRDs has already started to decrease in most samples (Kocevski et al. 2025; Barro et al. 2024a; Akins et al. 2025b; Pérez-González et al. 2024; Ma et al. 2025), F444W probes rest-frame 0.74 − 0.88 μm. To trace similar wavelengths at 3 < z < 4, we use F356W that covers 0.71 − 0.89 μm versus F444W’s 0.89 − 1.1 μm range. The same logic applies when substituting F200W with F150W: F200W probes 0.3 − 0.4 μm rest-frame at 4 < z < 5, and is substituted by F150W that covers 0.3–0.375 μm at 3 < z < 4. Applying this first color selection selects 491 galaxies with z_med = 2.7.

The redshift distribution of LRDs shows a rapid decline below z ≈ 4, indicative of significant evolutionary changes. To properly study this transition phase, we must focus on systems sufficiently close to z ≈ 4, as more distant (too low redshift) selections may yield galaxies that no longer share key characteristics with LRDs. A redshift cut is then applied, selecting 263 objects with z > 2.5 as shown in Figure 1, forming the parent sample. This criterion also serves to constrain the redshift range, preventing the filters from probing significantly different rest-frame wavelengths. Following visual inspection to remove artifacts and overly blended sources, 226 objects then remain with z_med = 3.6. By removing blended objects, merger phases may be excluded from this study, as such systems often exhibit complex morphologies, as discussed in Section 6.

Fig. 1. Color-magnitude diagram illustrating the parent sample: in purple we show the sources that respect the color-magnitude criterion, and black edges denote the parent sample, respecting z > 2.5. In red, we show the color-magnitude selection for the red color of LRDs F277W – F444W > 1.5 used in Akins et al. (2025b), having a zmedian ≈ 6.5. Red and purple contours represent all CEERS sources at 25%, 50% and 80% levels shown for F150W-F356W.

To probe the LRD population, a commonly adopted compactness criterion requires that the 4.4 μm flux be predominantly concentrated toward the center. This criterion is slightly different among the literature and can be: F_F444W(0.5″)/F_F444W(0.2″) < 1.5 (Barro et al. 2024a), F_F444W(0.4″)/F_F444W(0.2″) < 1.7 (Greene et al. 2024; Kokorev et al. 2024; Xiao et al. 2025) or F_F444W(0.2″)/F_F444W(0.5″) > 0.5 (Akins et al. 2025b), all leading to mostly unresolved objects. A population directly descendant from LRDs should similarly exhibit a compact core. To this end, we apply the condition F_{F444W,  inner}/F_{F444W,  total} > 0.4, where F_{X,  inner} denotes the flux measured within an elliptical aperture with a semi-major axis of 2 × HWHM = 0.16″ (half-width at half-maximum) centered on the object, and F_X, total is the total flux measured within the segmentation map, in the filter X. To account for galaxy morphology, we compute the ellipticity as a function of radius using Autoprof (Stone et al. 2021), adopting the ellipticity measured at 0.16″. We adopt a major axis of 0.32″, which is larger than the commonly used compactness threshold of 0.2″, in order to apply a less stringent criterion that also allows partially resolved objects to be included. The segmentation map is generated from a detection image built by stacking all available JWST bands from F115W to F444W. This approach improves the signal-to-noise ratio and incorporates flux from extended, faint outer regions, enabling robust detection of low-surface-brightness structures (see Figure 2). This selection yields a sample of 110 objects.

Fig. 2. Example illustrating the definition of the outskirts and inner regions. Left panel: stacked image combining all JWST bands, defined as STACK = F115W + F150W + F200W + F277W + F356W + F410M + F444W, using PSF-matched data. The dotted line shows the segmentation map, while the solid line traces the ellipse with a semi-major axis of 0.16″ = 2HWHM. Right panel: color composite image using F150W (blue), F277W (green), and F444W (red). The white circle indicates the PSF FWHM beam at F444W. As expected, the galaxy exhibits a centrally concentrated red emission and a bluer, more extended periphery.

Recent observations suggest that some LRDs exhibit a resolved young stellar component, forming an extended structure beyond their initial compact core (Chen et al. 2025; Torralba et al. 2026; Zhuang et al. 2025). This young component may grow over cosmic time, potentially explaining why LRDs are no longer observed as a distinct population at lower redshifts. In this study, we explore an evolutionary pathway in which LRDs may develop a blue, star-forming periphery while retaining their red inner core. To probe such systems, two additional criteria are required to identify this star-forming outskirt.

First, we impose $\frac{F_{\mathrm{F}115\mathrm{W},\mathrm{outskirt}}}{F_{\mathrm{F}115\mathrm{W},\mathrm{total}}}>0.4$, ensuring that the rest-frame UV emission originates predominantly from the outer region, consistent with a younger stellar population relative to the core. Here, F_{X, outskirt} is defined as the total flux minus the inner flux. A larger threshold (above 0.4) would fail to trace the phase in which LRDs begin to develop an outskirt, as the outskirts would already be prominent. Moreover, a too low threshold would probe objects without any outskirts and therefore pure LRDs that are not in this potential transition phase: ≈48% of the objects at a threshold lower than 0.4 are LRDs present in the Kocevski et al. (2025) sample. This requirement reduces the sample to 72 objects.

Second, we apply a color cut of F115W_outskirt − F200W_outskirt < 2 to ensure the presence of young stars in the periphery. If such a criterion is not applied, we observe the presence of a strong Balmer break in the outer regions of some galaxies at z ≈ 3.2. Such outskirts are unlikely to be associated with LRDs, and the Balmer break is therefore more likely to originate from an old stellar population, indicating a fully quiescent periphery. This contradicts the scenario of a recently acquired, star-forming envelope that we aim to probe in this work. Thus, applying this criterion limits the strength of the Balmer break in the redshift range z = 2 − 4.5, therefore constraining the presence of evolved stellar populations in the outskirts. Indeed, the Balmer break starts to affect the F200W filter at z ≈ 4.5, meaning that a strong break can significantly boost the flux observed in F200W for galaxies at z ≲ 4.5, making this band sensitive to older stellar components. In contrast, the F115W filter consistently probes rest-frame wavelengths shortward of the Balmer break down to z = 2, ensuring a reliable tracer of young stellar populations.

Applying only the first of the two criteria $(\frac{F_{\mathrm{F}115\mathrm{W},\mathrm{outskirt}}}{F_{\mathrm{F}115\mathrm{W},\mathrm{total}}}>0.4)$ would allow fully quiescent systems with no UV emission in either region, while using only the second (F115W_outskirt − F200W_outskirt < 2) would include galaxies whose cores are as blue as their outskirts, contradicting the scenario of a recently acquired star-forming envelope.

After applying this selection, our final sample consists of 55 galaxies. A summary of the full selection process is provided below:

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \text{F150W}-\text{F356W}>1.5\hfill \\ \hfill & z>2.5\hfill \\ \hfill & \frac{F_{\text{F444W, inner}}}{F_{\text{F444W, total}}}>0.4\&\frac{F_{\text{F115W, outskirt}}}{F_{\text{F115W, total}}}>0.4\hfill \\ \hfill & {\text{F115W}}_{\text{outskirt}}-{\text{F200W}}_{\text{outskirt}}<2.\hfill \\ \hfill \end{array}\end{array}$$

This selection relies on a single evolutionary scenario in which LRDs develop a blue star-forming outskirt while retaining their red inner core. As a result, the population probed in this work may not encompass all potential LRD descendants, and alternative evolutionary pathways may exist, as discussed later in Section 6.

2.4. Final sample

Following the selection process, we retain a final sample of 55 objects, with z_med = 3.6 ± 1.1, each with a defined outskirt and inner part. Among these, 15 objects have spectroscopic redshifts. Six of them are already included in the LRDs sample studied in Kocevski et al. (2025), forming the high-redshift portion of our sample at z_med = 6 as visible in Figure 3. All these LRDs present the same characteristic: an extended blue outskirt emitting in the UV, potentially indicating this transitional phase. Using SED fitting and morphological analysis, we aim to study the properties of the inner parts and the outskirts in the following sections. The goal is to identify features or characteristics within our sample that may be associated with LRDs at redshifts z < 4, in order to confirm or rule out a potential evolutionary link between the two populations.

Fig. 3. Redshift distribution of the final sample in purple. Six LRDs from our sample were also in the Kocevski et al. (2025) LRDs sample, shown in red.

Since our selection targets red galaxies, the sample may contain dusty star-forming galaxies (DSFGs), such as the optically-faint galaxies (OFGs), H-dropout or red star-forming galaxies (RedSFGs, Wang et al. 2019; Gómez-Guijarro et al. 2023; Pérez-González et al. 2023; Tarrasse et al. 2025) which typically exhibit a dust attenuation of A_v > 2. Our study focuses on objects with inner cores resembling those of LRDs, located around the transitional epoch at z = 4, and thus expected to be dust-free, consistent with the properties of LRDs (Labbe et al. 2025; Akins et al. 2025b; Casey et al. 2024). SED fitting (see Section 3) yields a median attenuation in the inner regions of A_v, med = 0.84 ± 0.80. Using the FIR catalog of Henry et al. (in prep.), only three galaxies in our sample are detected in the Spitzer/MIPS 24 μm band with a S/N greater than 5, and none are detected in the Herschel/PACS bands from 100 μm to 500 μm with an S/N above 3, indicating that our sources are not significantly dust-obscured. We also examined the H-dropout selection defined in Wang et al. (2019), and found that only one source in our sample satisfies their criteria. This indicates that our sample has minimal overlap with this population.

3. Inferring physical properties

3.1. Spectral energy distribution fitting

Similar to the LRD population, the inferred physical properties in our sample can highly depend on the model we assume. Indeed, the derived properties vary significantly depending on whether the modeling assumes a purely stellar component, a stellar component with an AGN, an AGN-only or the recently proposed BH^* scenario (Pérez-González et al. 2024; Akins et al. 2025b; Leung et al. 2025; Rinaldi et al. 2025; de Graaff et al. 2025). The nature of the AGNs in LRDs remains an open question, here we adopt a stellar-only modeling approach. Consequently, all comparisons between the properties of our inner components and the characteristics of LRDs are based on this stellar-only assumption. We note that, the inferred mass of the inner region is potentially overestimated, as any AGNs contribution is not accounted for.

To measure the flux of the inner part across all bands for constructing SED_inner, we measure the flux within the elliptical aperture defined by Autoprof in all filters (see Figure 2). For SED_outskirt, we measure the remaining flux in the segmentation map mask after subtracting the inner elliptical aperture in all bands. Uncertainties were estimated by measuring the flux in 30 dark sky regions surrounding each galaxy, using the same aperture as that employed for the galaxy flux measurement. For this we used 4.2″-wide cutouts centered on the target galaxy. The final uncertainty was taken as the root mean square (RMS) of these off-source measurements. CIGALE (Boquien et al. 2019) is then used to derive the physical properties of each region, employing identical models and parameters for both components: We assume a delayed star formation history (sfhdelayed) with SFR ∝ t exp(−t/τ), a Chabrier IMF (Chabrier 2003), and bc03 stellar populations (Bruzual & Charlot 2003). We include both the nebular module and the dustatt_modified_starburst attenuation module, which combines the Calzetti et al. (2000) law with extensions from Leitherer et al. (2002) (see Table 1).

3.2. Morphological analysis

To infer morphological properties, we use Astrophot (Stone et al. 2023) in single-band mode with the F444W filter of NIRCam. This software performs 2D profile fitting and takes into account the PSF, offering several fitting methods and models. We employ the Levenberg-Marquardt method with a single Sérsic profile for each galaxy. For this, we use small cutouts of 6 arcseconds, centered on each galaxy. All galaxies in the field of view were fitted to avoid contamination from neighbors. The Sérsic index n was restricted to 0.2 < n < 10, and the half-light radius to 0.01″ < R_e < 1″. Other parameters, such as the position (x, y), the position angle θ, the ellipticity q, and the central brightness density I_e, were left free. Among the sample, 30% have R_e < HWHM, potentially leading to nonphysical characteristics.

To determine how far we can go below the PSF HWHM to study morphological features, we simulate 600 galaxies using single Sérsic profiles with Astrophot, convolving them with the PSF at F444W and injecting them into dark areas of the CEERS field of view. We fix the ellipticity at q = 0.7, corresponding to the median values of our sample with R_e > HWHM. The injected half-light radius increases from R_e = 0.01″ to R_e = 0.12″ and the Sérsic index from n = 0.5 to n = 6. Dark regions are randomly selected in the field of view, and must have S/N < 2. We also set the total magnitude of the injected sources to 26 mag, which corresponds to the faintest source among those with R_e < HWHM. Knowing the true values of R_e and n for the injected sources, we compute the total fitting error after retrieving their structural parameters using Astrophot:

$$\begin{array}{c}\hfill {\text{Err}}_{\text{tot}}=\sqrt{{\text{Err}}_{R_e}^2+{\text{Err}}_n^2}=\sqrt{{\left(\frac{R_{e,\text{obs}}-R_{e,\text{r}}}{R_{e,\text{obs}}}\right)}^2+{\left(\frac{n_{\text{obs}}-n_{\text{r}}}{n_{\text{obs}}}\right)}^2.}\end{array}$$(1)

Here, R_e, obs and n_obs denote the inferred (observed) parameters obtained from fitting the mock sources, while R_e, r and n_r correspond to the true (real) input parameters used to generate the mock data. The results are shown in Figure 4. We find that below the HWHM, the total error remains low and nearly constant until a threshold of R_e ≈ 0.04″. Below this value, more than 40% of the sources exhibit total errors exceeding 50%, with a median error approximately four times higher than that of sources with R_e > 0.04″. Adopting a lower resolution threshold would therefore risk mischaracterizing the morphology of the galaxies. However, this threshold can be considered a conservative resolution limit. Indeed, the most compact sources in our sample (R_e < HWHM) all have magnitudes brighter than 26 AB, and the effective resolution limit is expected to improve with decreasing magnitude (Akins et al. 2025b). Applying this resolution limit R_limit = 0.04″, we consider 77% of our sample to be resolved in F444W.

Fig. 4. Error in percent of the half-light radii Re and Sérsic index n as a function of the inferred half-light radii for mock sources injected into dark areas of the CEERS field of view with mag = 26. The error is defined as ${\text{Err}}_{\text{tot}}=\sqrt{{\text{Err}}_{R_e}^2+{\text{Err}}_n^2}$. The dotted line denotes the HWHM in F444W, whereas the dash-dotted line refers to the resolution limit adopted in this study.

4. Relation with LRDs

This sample has been constructed to identify potential analogs of LRDs at lower redshifts, based on their color profiles and internal light distributions. To evaluate whether this sample is truly comparable to or related to LRDs, we first compare their stellar masses, although this comparison is subject to significant uncertainties, primarily due to the strong dependence of mass estimates on the physical model adopted in previous studies. Nevertheless, meaningful comparisons remain possible, provided that the same model is consistently applied. The masses of LRDs assuming only a stellar model are notably high, approximately 10¹⁰ M_⊙ as illustrated in Figure 5. This figure compiles several studies of LRDs that consider only the stellar model. The dark purple distribution represents the mass of the inner part and the total mass of the sample. The stellar masses of LRDs inferred without MIRI data by Akins et al. (2025b) are likely overestimated at z > 5, as shown by Wang et al. (2025) and discussed in Akins et al. (2025b). The inner stellar mass of our sample is consistent with the bulk of the distributions, with a median value of 10^{9.9 ± 0.4} M_⊙, and a slightly higher total stellar mass of 10^{10.1 ± 0.4} M_⊙. This total includes the outskirt component of 10^{9.6 ± 0.6} M_⊙. The median stellar age of the outskirts is 300 ± 200 Myr, indicating that they are effectively composed of young stars formed relatively recently as the selection has been constructed (see Section 2.3). Given the median redshift of z_median = 3.6, this implies that the bulk of the outskirts began forming around z ≈ 4.2.

Fig. 5. Mass distribution for several samples of LRDs from previous works, assuming only a stellar model, from Rinaldi et al. (2025), Pérez-González et al. (2024), Akins et al. (2025b), Leung et al. (2025), Labbé et al. (2023). This sample is visible in purple for both the inner part and the total mass.

Because galaxies are known to grow in size with redshift, following the mass-size relation (van der Wel et al. 2014; Allen et al. 2025), the compactness of this sample can be evaluated with reference to this relation in addition to our selection criterion ($\frac{F_{\text{F115W, outskirt}}}{F_{\text{F115W, total}}}>0.4$). As shown in Figure 6, the mass-size relation is visible in each redshift bin described in Allen et al. (2025) in the top row. The majority of the sample falls below the relation and can thus be considered compact. We then compute the distance to the mass-size relation in the bottom of Figure 6. The median, shown in cyan, is 990 ± 330 pc below the mass-size relation. Due to the unresolved sources in F444W (for which we consider R_e = 0.04″), this can be regarded as a lower limit. In Leung et al. (2025), the mass-size relation of LRDs was studied under both stellar-only and stellar+AGNs model assumptions. Since we adopt a stellar-only model in our analysis, we show only the stellar results from their work in Figure 6, represented as red triangles. Since this sample was not resolved, the HWHM was used as R_e. The sample we constructed here exhibits a similar distance to the mass-size relation compared to LRDs, which is ≈700 pc. This indicates that this population of galaxies are as compact as LRDs relative to the mass-size relation. The median size of our overall sample is R_e,median = 520 ± 300 pc in F444W. If only the resolved sources are considered, it reaches R_{e,resolved,median} = 670 ± 230 pc. For the unresolved sources, it is R_{eunresolved,median} = 255 ± 30 pc, which is an upper limit.

Fig. 6. Top row: The stellar mass vs. Re plane in each redshift bin of the sample, with the mass-size relation from Allen et al. (2025) overlaid. The relation for quiescent galaxies (QG) from Yang et al. (2025) is shown in orange. Triangular markers represent unresolved sources (for which we consider Re = HWHM), while circular markers indicate resolved sources. Bottom: Distance to the mass-size relation for the entire sample. The cyan square represents the median of the sample, with 15% and 85% percentiles. The red triangles denote LRDs from Leung et al. (2025) which are also unresolved.

Another feature of LRDs that has been studied is their stellar surface density (Akins et al. 2025b; Leung et al. 2025; Guia et al. 2024; Khan et al. 2025a). The central density inferred for LRDs is notably high, with values reaching Σ_* ≈ 10^11.6 M_⊙ kpc⁻² and Σ_* ≈ 10¹¹ M_⊙ kpc⁻² in Akins et al. (2025b) and Leung et al. (2025), respectively. These values are broadly consistent with the maximum stellar surface density observed across a wide range of redshifts at z = 0 − 3, as reported by Hopkins et al. (2010) and Grudić et al. (2019), with Σ_max = 3 × 10¹¹ M_⊙ kpc⁻². The sample studied here has a median redshift of z_med = 3.6. Thus, with the F444W band, we are directly tracing the 1 μm rest-frame, effectively mapping the mass. Assuming that the surface brightness profile observed in the F444W band probes the surface mass profile, we use the previously inferred R_e and n from the single Sersic fit to derive the surface mass density profile. By knowing the total stellar mass M_tot = M_inner + M_outskirt, we can infer Σ_e, defined as the stellar mass surface density at R_e:

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_e& =\frac{M_{\text{tot}}}{2\pi q{R}_e^2{J}_n},\text{where}\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill J_n& ={\displaystyle {\int }_0^{\infty }}x(exp[-b_n{x}^{1/n}-1])dx,\hfill \end{array}$$(3)

where q is the ellipticity, J_n is then numerically computed. Because we are using morphological features and a non-negligible portion of our sample is unresolved in the F444W band, we compute this profile only for galaxies that are considered resolved, those with R_e > 0.04″ (see Section 3.2). Among the six LRDs in our sample, only one is resolved and therefore used to derive the surface density profile. To verify the robustness of the inferred Σ_e, we also calculated it directly using the mass-to-light ratio measured in the F444W band: ${\mathrm{\Sigma }}_e=I_e\frac{M_{\ast }}{L_{F444W}}$, I_e being the surface brightness at R_e, M_* the stellar mass inferred with CIGALE and L_F444W the total luminosity of the object in F444W.

This yields consistent results: Σ_*, median = 10^{11.6 ± 1.9} M_⊙ kpc⁻² in both cases. This method assumes a constant mass-to-light ratio across the galaxy, which is not the case in our sample. To quantify the bias of this assumption, we recompute the central mass profile using the inner part, which better describes the central region: ${\mathrm{\Sigma }}_e=I_e\frac{M_{\ast ,\text{inner}}}{L_{F444W,\text{inner}}}$. This yielded very similar results: Σ_*, median = 10^{11.5 ± 1.9} M_⊙ kpc⁻².

In Figure 7, we present the median profile of the stellar mass density in purple. The profile for LRDs, as reported by Akins et al. (2025b) and assuming only a stellar model, is shown in red. Additionally, the densest objects observed across a wide range of redshifts, as studied in Hopkins et al. (2010) and Grudić et al. (2019), are also visible. The inferred stellar density profile aligns with the LRDs profile toward the center, suggesting that our sample has very similar physical properties in the inner region. More over, both profiles match the densest objects observed in the universe. However, the profile inferred for LRDs is slightly higher than Σ_max, which may be an overestimation in the case of LRDs due to their nearly unresolved size and the presence of an AGN, as discussed in Akins et al. (2025b), which can also be a possibility in our sample. Furthermore, the physical size of our sample is larger than that of LRDs (≈520 ± 300 pc vs 100–300 pc for LRDs), indicating a more extended structure which is also visible in Figure 7.

Fig. 7. Stellar mass density profiles for our sample using only resolved galaxies in purple. In red the profile inferred by Akins et al. (2025b) assuming only a stellar model. We also add dense cores from Hopkins et al. (2010) and Grudić et al. (2019).

A selection based on a color-color diagram was developed by Barro et al. (2024a) to identify both types of LRDs: those with and without the V shape. However, at lower redshifts (z < 4), galaxies escape the selection due to the fixed bands used. To assess whether the inner cores of our objects are as red as LRDs while still satisfying the blue color criterion across a broad redshift range, we recompute the color–color diagram by shifting the filter bands for galaxies with z < 4. Figure 8 illustrates this selection for the inner regions only, by shifting the color criteria as a function of redshift: in red for z > 4, using the same bands as Barro et al. (2024a): F115W-F200W and F200W-F444W. In purple, we show our sample with z < 4, with the bands shifted accordingly to trace the same rest-frame flux at z < 4, i.e., using F814W-F150W and F150W-F356W. All the galaxies selected effectively present a red inner region, as our selection was designed to enforce. However, the blue color appears to evolve as redshift decreases. Specifically, the blue color at z < 4 (F814W-F150W) appears redder than its counterpart at z > 4 (F115W-F200W), causing approximately 50% of the galaxies at z < 4 to fall outside the LRD color selection defined by Barro et al. (2024a). This shift suggests a possible cosmic evolution of the population over time.

Fig. 8. Color–color selection of LRDs described in Barro et al. (2024a), originally defined for galaxies at z > 4, applied to the inner regions of our sample at z > 4 (shown in red). For galaxies at z < 4, the same rest-frame colors are probed by shifting the bands accordingly (shown in blue). Contours represent the 25%, 50%, and 75% density levels for each distribution. The applied color correction is indicated by a thick black line, with an average shift of ≈0.2 mag.

5. The building of an outskirt

Our sample appears to be more extended than typical LRDs (see Figure 7), yet it retains a similarly dense and red core, and becomes slightly more massive when the outskirts are taken into account, suggesting a potential growth scenario. If the outskirts formed around LRDs in a later stage, we should observe a time correlation with the outskirt mass fraction. To investigate this, we compute the outskirt mass fraction, f_outskirts = m_outskirts/(m_outskirt + m_inner), as a function of redshift, as shown in Figure 9. There is a strong correlation with redshift, with the outskirt mass fraction increasing at lower redshifts. The LRDs from Kocevski et al. (2025) present in our sample (visible as red stars) lie at the very low outskirt mass and high redshift end of the figure. This is in line with the idea of a growth scenario. The observed correlation could be affected by cosmological surface brightness dimming. Since dimming causes the surface brightness to appear fainter at higher redshift (I_obs ∝ I₀(1 + z)⁻⁴), we can only detect the brightest outskirts as redshift increases (see Appendix A.1). In other words, at higher redshift we preferentially detect the most massive outskirts, potentially missing the low-mass ones, if any. Consequently, the outskirt mass fraction at high redshift shown in Figure 9 should be considered as upper limits. Once corrected for this effect, the trend would likely appear even clearer.

Fig. 9. Outskirts’ mass fraction as a function of redshift. Gray dots represent the data, while gray contours illustrate their distribution at 0.25, 0.5 and 0.75 levels. The median is shown as a purple circle, with each point computed in bins of Δz = 1. Vertical error bars indicate the standard deviation, while horizontal bars represent the width within each bin. Red stars refer to the LRDs also in Kocevski et al. (2025).

Because of this growth, the compactness criterion may no longer select these galaxies, also correlated with an evolution of colors at lower redshift (see Figure 8).

To further investigate the spatial growth of the outskirts, we compute the median stacked image of our sample within bins of f_outskirt = [0 − 0.2, 0.2 − 0.4, 0.4 − 0.6] using the F444W filter. This stacking approach increases the S/N by at least a factor of three in each f_outskirt bin, enabling the detection of faint substructures that would otherwise remain undetectable in individual images. In each bin, we fit a single Sérsic profile using Astrophot to derive R_e as a function of f_outskirt, shown as squares in Figure 10. Since R_e is not significantly affected by cosmological dimming (Yu et al. 2023; Ormerod et al. 2024), and given that we increased our S/N to minimize its potential impact, we consider our derived R_e to be robust and physically meaningful.

Fig. 10. Half light radii Re as a function of foutskirt. The stars correspond to the median Re from stacking all the filters, and the squares correspond to only F444W. The color bar shows the mean redshift in each foutskirt bin.

As expected, as the mass fraction in the outskirts increases, the physical size of the galaxies also grows with decreasing redshift. The most compact systems, with the lowest f_outskirt, are unresolved in F444W and are located at high redshift, with an upper size limit of 260 ± 55 pc. This increases to around 285 ± 80 pc for galaxies with moderate f_outskirt, and to about 590 ± 120 pc for the most extended and massive systems. Stars in Figure 10 represent the median sizes derived across all JWST PSF-matched bands, probing not only rest-frame NIR emission but also the rest-frame UV, allowing to increase the S/N by a factor 9, allowing faint detection. Galaxies without any outskirt are still not resolved, having an upper limit of 260 pc, which is expected, and those with a moderate outskirt reached 510 ± 140 pc while the most extended systems with the largest outskirt reached 828 ± 160 pc. This suggests that the UV emission starts to be more prominent and detected when the fraction of outskirt is above 0.2, at redshift ≈4.

A common way to select LRDs is by their V shaped SED, based on continuum slopes, first used by Kocevski et al. (2025). To determine if the formation of these outskirts impacts the V shape, we reconstruct the β_UV vs. β_optical plane, a rest-frame criterion, using only the inner SED which is expected to be mostly related to LRDs. We compute β_UV,inner and β_{optical,inner} assuming f_λ ∝ λ^β in the rest-frame UV and rest-frame optical region. Each slope β_inner is compute using 3 bands blueward or redward of the Balmer break at 3645Å as in Kocevski et al. (2025), by adjusting a linear fit in log(f_λ) space. Due to the shallower depth of the HST filters (28.73 and 28.50) compared to F115W (29.15), some of our sources are undetected in either F606W or F814W, occasionally resulting in negative flux values within the central aperture due to background subtraction. In such cases (affecting 5 objects), we compute β_UV,inner using only the two available bands and conservatively set the associated uncertainty to 1.

In Figure 11, we display this plane and overlay the V shape selection from Kocevski et al. (2025) with a dashed region. Contours are also visible for different bins of f_outskirts. The red stars indicate the 6 LRDs also present in Kocevski et al. (2025). It appears that a low f_outskirts is necessary to maintain the V shape. As f_outskirts increases with decreasing redshift, galaxies depart from this V shape selection. In the top right-corner we display the evolution of β_UV,inner as a function of the redshift: as the redshift decreases, β_UV,inner increases and escapes the V shape selection (dashed region), at z ≈ 4. This departure could explain why even a rest-frame selection suffers from a drop around z = 4, due to physical evolution: as the outskirt is growing, a stellar component becomes non-negligible around this redshift for the majority of LRDs. If this stellar component erases the V shape, it could also explain why some LRDs exhibit the V shape while others do not. Broad color-color selection traces all LRDs, with or without a stellar component, and the V shape may select only those dominated by AGN, consistent with the finding that 70% of LRDs selected using the V shape present broad lines, indicative of AGNs (Kocevski et al. 2025).

Fig. 11. βUV/βopt plane of the inner regions of the sample. The fraction of the mass outskirt is represented by contours with 0.4 and 0.8 levels, increasing from red to purple and blue. Redshift is indicated by a gray color bar. The LRD selection from Kocevski et al. (2025) is shown in the light gray dashed region and LRDs in common are shown in red stars. In the top right-hand corner, we show βUV,inner as a function of redshift, the colorbar denotes the mass fraction of the outskirt foutskirt, the shaded region refers to the V shape selection.

If our evolution and growing scenario of LRDs is corect, the number density of our sample should be compatible with the one of LRDs and at lower redshifts. We compute the number density as a function of redshift in Figure 12, by dividing our sample into redshift bins and calculating the corresponding comoving volume. As we have 6 LRDs in our sample, we exclude them to compare our distribution with the one of LRDs. Our sample is visible in red and the common sources with Kocevski et al. (2025) in purple. We also compute the distribution of LRDs based on the sample of Kocevski et al. (2025) in CEERS only (visible in orange), and the ones from Akins et al. (2025b), Barro et al. (2024a), and Euclid Collaboration: Bisigello et al. (2026). Additionally, the two theoretical models of the number density developed by Pacucci & Loeb (2025) and Inayoshi (2025) are shown for comparison.

Fig. 12. Number density as a function of redshift from our sample visible in red, excluding the six common sources with the LRDs in CEERS from Kocevski et al. (2025), here shown in orange. In purple, the number density of the sources in both samples is shown. The dimming-corrected number density of our sample is shown with dark gray square symbols, interpreted as an upper limit. LRDs from Akins et al. (2025b), Barro et al. (2024a), and Euclid Collaboration: Bisigello et al. (2026) are represented by light gray squares, circles, and triangles, respectively. Horizontal error bars indicate the redshift bins used to compute the number density, and vertical error bars represent Poissonian errors. The theoretical evolution from Pacucci & Loeb (2025) is shown as the dashed black line, while the solid line refers to the model of Inayoshi (2025).

Firstly, we find that the final number density in the lowest redshift bin, z = 3 ± 0.5, is 10^{−4.15 ± 0.09} Mpc⁻³, which is consistent with the number density of LRDs at z = 6 ± 0.5 in CEERS (based on the sample of Kocevski et al. (2025)), 10^{−4.17 ± 0.11} Mpc⁻³. The fact that we recover a number density consistent with that of LRDs suggests that this evolutionary scenario may be viable and potentially the most likely pathway for the evolution of LRDs. Nevertheless, this does not exclude other possibilities, as some LRDs in a merger state have already been observed (Mérida et al. 2025; Tanaka et al. 2024), as discussed in Section 6.

Secondly, the number density of LRDs (in orange) declines at the same redshift where our distribution (in red) starts to increase, around z ≈ 5. However, since we select objects with a recently formed blue periphery, this feature may either be undetectable at higher redshifts–due to cosmological surface brightness dimming–or may not have formed yet. To account for this effect, we applied a correction to the observed number density in each redshift bin, corresponding to the fraction of galaxies likely missed due to cosmological dimming (see Appendix A). The corrected number density is shown as dark gray square symbols in Figure 12. We emphasize that this correction assumes all our sources would already possess a developed outskirts component at higher redshift. In this case, the apparent emergence of the outskirts at lower redshift would not be a physical effect, but rather an observational bias caused by cosmological dimming–implying that LRDs could host undetected peripheries at high redshift. However, stacked analyses by Akins et al. (2025b) and Delvecchio et al. (2025), who respectively combined 429 and 302 individual LRDs from different samples across filters ranging from F150W to F444W, did not detect any extended emission in their stacked images. This suggests that the majority of LRDs are intrinsically compact and do not host faint, dimmed outskirts that would remain undetectable in individual observations. Therefore, this favors the physical formation of a stellar periphery over cosmic time rather than an observational bias caused by cosmological dimming. Based on this, we conclude that the dimming-corrected number density should be considered as a conservative upper limit, while the true number density of such objects is likely closer to the observed values. It is also important to mention that it does not exclude other evolutionary scenario, as some LRDs in a merger state have already been observed (Mérida et al. 2025; Tanaka et al. 2024), as discussed in Section 6, it rather suggest that the other evolutionary scenario may not be the main channel of evolution for LRDs.

V shaped LRDs might be a pure black hole embedded in dense gas as recently suggested by Ji et al. (2025), Inayoshi (2025), Naidu et al. (2025), Rusakov et al. (2025). As the V shaped SED diminishes with the growth of the stellar component and time (see Figure 9), the influence of this component likely increases, and LRDs begin to lose their known characteristics. In this scenario, V shaped LRDs represent an earlier phase, while LRDs without the V shape represent a later one.

6. Discussion

We now propose to discuss the origin of the growing stellar outskirt. Assuming that these galaxies are initially LRDs, several scenarios could take place, allowing stars to populate the outer regions, as illustrated in Figure 13:

• a)

  Gas can be expelled from AGNs and supernova feedback from the inner region and then form stars.

• b)

  One or a succession of mergers can favor the formation of the outer region:

  • (b).1

    Wet mergers can bring gas that fuels the LRDs through cosmic time.

  • (b).2

    Dry mergers, however, are more likely to form an outer region solely through kinematic processes and tidal interactions, due to the absence of gas.

• c)

  Direct accretion from the environment by either cold streams or diffuse hot accretion may help to form the stellar outskirt.

Fig. 13. Drawing illustrating the different possible scenarios of evolution of LRDs, as discussed in Section 6: through (a) feedback, (b) mergers, or (c) accretion.

Throughout scenario (a), the gas initially present in the inner regions could be ejected by feedback from the central AGNs or supernovae, potentially forming a gaseous outskirt that may later give rise to new stars. This would result in younger metal-rich stars with higher sSFR in the outskirt and, consequently, significant opacity. In general, as discussed in Baggen et al. (2024), an expansion of the stellar distribution can occur with mass loss due to AGNs feedback or stellar winds, increasing the size of the galaxy (Fan et al. 2008) by the redistribution of the gas in the galaxy (Khan et al. 2025b) and therefore could be applied to LRDs. The molecular gas reservoir of an LRD was measured by Akins et al. (2025a) using CO(7-6) and [C_I](2 − 1), where the amount of neutral molecular gas appeared to be limited M_{neutral gas} < 10¹⁰ M_⊙. However, the presence of an extreme gas reservoir was suggested by Ji et al. (2025), de Graaff et al. (2025), Rusakov et al. (2025) and as in the Black Hole Star model (BH^*) proposed in Naidu et al. (2025), in which the gas is extremely dense (n_H ≈ 10¹¹ cm⁻³) and turbulent (v_turb = 500 km s⁻¹). Due to these extreme conditions, the gas could be spatially redistributed because of instabilities favored by either AGNs or supernovae feedback from the host galaxy, a scenario also discussed in Taylor et al. (2025b). However, the fraction of ionized gas $\chi =\frac{n_{{\text{H}}^{+}}}{n_{\text{H}}}$ remains to be studied in such systems. A large fraction of ionized gas would also require an explanation of how this energetic component becomes neutral before forming new stars. Moreover, if the velocity of the gas expelled is too important (v_gas > v_escape), it could be expelled beyond the gravitational potential of the system, thus no longer remaining bound to it and preventing star formation in the outskirts. Therefore, although plausible, this scenario remains highly speculative.

Scenario (b) could also occur, with mergers between LRDs or between a LRD and a galaxy. If the merger involves a gas-rich component (scenario (b).1), it may bring gas in the outskirt through tidal events and kinematic interactions, leading to the formation of new stars outside the core. However, mergers generally tend to concentrate the gas in the central region (Hopkins et al. 2013) and therefore also trigger star formation in the inner part, so that a higher sSFR and younger stars in the center is more likely to be observed. The stars already in place in the LRD may also be ejected into the outskirt through dry merger events (scenario (b).2), as discussed in Khan et al. (2025a), who modeled the evolution of LRDs over hundreds of millions of years, accounting for the impact of supermassive black hole (SMBH) mergers. Using initial parameters from LRDs first observed by Baggen et al. (2024) at z_spec = 6.68 and z_spec = 6.98, they found that the galaxy undergoes mergers that flatten the initial stellar mass density profiles by one order of magnitude, decreasing from 10¹² M_⊙ kpc⁻² to 10¹¹ M_⊙ kpc⁻² over a timescale of 100 − 800 Myr, which approximately corresponds to the period between z = 4.5 and z = 3. In the BH^* model of Naidu et al. (2025), the blue component of LRDs is explained by a very young stellar population that would be approximately 1 Gyr old at z ≈ 3 if formed at z ≈ 5. If these stars are then redistributed through dry merger events, the result would be an outskirt as old as the inner region with similar colors. LRDs in a merger state have already been observed by Rinaldi et al. (2025) and Chen et al. (2025), using high-resolution UV data, and a more complex system involving two LRDs merging in Mérida et al. (2025) and Tanaka et al. (2024), supporting this scenario (b) for at least some descendants of the LRD population.

In scenario (c), the gas is directly accreted by cold gas filaments or diffuse accretion of warm gas (the so called hot accretion mode). Such accretion is expected to feed a disk together with its angular momentum (Kereš et al. 2005; Dekel & Birnboim 2006; Dekel et al. 2009). Because LRDs are high redshift sources, the cold accretion mode is more likely to dominate those systems. The infalling gas would be metal-poor, leading to young stars with low metallicity in the outer region (Sánchez Almeida et al. 2014; Pichon et al. 2011). This has also been discussed in Kocevski et al. (2025), where infalling gas is suggested to bring angular momentum, promoting star formation on larger scales at later stages.

Because we are studying LRDs successor candidates over a long timescale of approximately 2 Gyr (2.5  <  z < 7.5), all those scenarios could occur. To distinguish which one is more likely, we analyze in Figure 14 ages and specific star formation rates (sSFR) of the outskirt relative to the inner part as a function of the outskirt mass fraction. Two regimes appear to emerge, depending on whether f_outskirt is below or above 0.4.

Fig. 14. Differences of the ages and sSFR between the outskirt and inner part as a function of the outskirt mass fraction foutskirt. The data are shown in gray with the median profile in pink. For each distribution, the black solid line refers to no radial changes between the inner/outksirt part.

The first regime, where f_outskirt < 0.4, 90% of galaxies exhibit a younger outskirt, favors scenarios (a) and (c). Additionally, also 90% of these galaxies present a higher sSFR in the outer region, with a significant fraction (20%) having ${\text{sSFR}}_{\text{outskirt}}-{\text{sSFR}}_{\text{inner}}>5{\text{Gyr}}^{-1}$. This excludes the merger scenario (b) when f_outskirt < 0.4, in which a flat specific star formation rate (sSFR) profile is more likely to be observed with dry mergers ((b).2), whereas a bluer inner part is more expected with wet mergers ((b).1), with a higher sSFR in the center, which is not observed in our sample. This might be a possible evolutionary path for LRDs, but our sample was not specifically constructed to trace this population. Moreover, Carranza-Escudero et al. (2025) found that the environment of LRDs tends to be less dense compared to other galaxies at similar redshifts, making the merger-driven evolutionary pathway less plausible. Finally, galaxies in this regime are very close to an LRD state, having a young, low-mass outskirt with the V shape. Therefore, all the characteristics we would expect to observe in connection with feedback from AGNs or supernovae in LRDs are most likely to be observed in this phase.

In the second regime, when f_outskirt > 0.4, the proportion of galaxies with a younger outskirt drops to 50%, the ages are more homogeneous and the V shape is no longer observed, which disfavors scenario (a) and (b).1. The sSFR profile shows the same characteristics, suggesting a homogeneous spatial repartition of both those properties, which favors the dry merger scenario of growth (b).2. Nevertheless, if the accretion rate declines with time (as it scales with (1 + z)^5/2; Dekel et al. 2013), the sSFR in the periphery would also decrease, potentially reaching values comparable to those of the inner region, thus making scenario (c) a plausible evolutionary pathway in this regime too.

If we now consider the outskirt mass fraction as a proxy for time, as suggested by Figure 9, direct accretion from filaments (scenario (c)) or AGN/supernova feedback (scenario (a)) could dominate initially even though the latter mechanism remains more uncertain (no strong evidence of gas reservoir except highly ionized gas), with the V shape in place up to f_outskirt = 0.3, from z = 6.5 to z = 4.5, corresponding to approximately 1 Gyr. In this regime, galaxies have a young, low-mass outskirt with a higher sSFR. Over time, the accretion rate declines ($\propto {(1+z)}^{\frac{5}{2}}$; Dekel et al. 2013), the sSFR and ages of the periphery respect to the inner region tend therefore to homogenies, as seen in our sample when f_outskirt > 0.4 in Figure 14. Merger events could also occur, with a rate of ≈2 Gyr⁻¹ at z = 4 (Mundy et al. 2017; Puskás et al. 2025). Through kinematics, a redistribution of stars can take place, where age and sSFR differences should no longer be observable. Because the merger mechanism alone seems insufficient (Inayoshi 2025) and even unlikely (Carranza-Escudero et al. 2025) to explain the decline in the number density of LRDs, another mechanism is needed that could be either scenario (a) or (c). This mechanism could first be a phase of expansion due to accretion (scenario (c)) or the central AGNs activity (scenario (a)), followed by a gradual decline in accretion (scenario (c)) or dry mergers (scenario (b)). At redshift 3.5 < z < 4.5, the stellar mass within the outskirts is ≈10⁹ M_⊙, consistent with the estimated accreted mass around LRDs in Chen et al. (2025). They investigated the origin of the extended emission observed in three different LRDs, attributing it to gas accretion. Then, they estimated accretion rates and the resulting stellar mass, assuming a star formation efficiency of 0.1. From redshift z = 6 to z = 4, this would yield a stellar mass of about 10⁹ M_⊙, strongly supporting the evolutionary scenario via accretion (scenario (c)).

7. Conclusion

To date, the nature, origin, and evolutionary path of LRDs remain poorly understood. Their steep decline in number density at z < 4 suggests the action of physical processes that alter their observable properties, making them increasingly difficult to identify as they evolve. In this study, we investigated the possible fate of this population by identifying a set of candidate galaxies that may represent their descendants: galaxies as red as LRDs with a blue star-forming periphery at z_med = 3.6.

The population we identified appears to meet the key conditions expected of LRD descendants. They are approximately 1 kpc below the mass-size relation, exhibit stellar masses (M_* = 10^{10.1 ± 0.4} M_⊙) and central stellar densities (Σ_* ≈ 10^{11.6 ± 1.9} M_⊙ kpc⁻²) similar to those of LRDs. Their number density at z = 3 ± 0.5 is about the same as that of the LRDs at 5 < z < 7, 10^−4.15 Mpc⁻³ and they present the same rest-frame red colors. The blue rest-frame color is also consistent for sources at redshifts z > 4. However, approximately 50% of the sources at redshifts z < 4 fall outside the blue color selection, indicating a redshift-dependent evolution.

Our interpretation is also strengthened by the finding that the fraction of mass contained in the outskirt increases with cosmic time, suggesting that we see the growth of outskirts of the LRDs taking place through cosmic time. The disappearance of the LRDs at lower z might therefore be accompanied by the acquisition of a stellar extension, which is compatible with a recently built extension. The growth scenario is strongly supported by the fact that the high-redshift sources in our sample are unresolved, with an upper limit on the effective radius of R_e ≈ 260 pc and an outskirt mass fraction of f_outskirt ≈ 0.1 at z ≈ 4.5. As redshift decreases, both f_outskirt and the effective radius R_e increase, with R_e ≈ 285 pc and f_outskirt = 0.3 at z ≈ 3.5, and R_e ≈ 590 pc and f_outskirt = 0.5 at z ≈ 3. Additionally, six LRDs previously analyzed in Kocevski et al. (2025) are found at z > 5 within our sample. All of them exhibit rest-frame UV emission in their outer regions.

Our results suggest that the decline in the number density of LRDs is driven by the acquisition of a stellar component. As this component grows, the characteristic V shaped SED fades, and the physical size of the galaxies increases. Consequently, both the V shaped SED and compactness criteria used to select LRDs may no longer effectively identify this population.

The growth of the stellar component can be explained by several mechanisms. However, accretion through cold streams or diffuse inflows is the more favored scenario, as indicated by the presence of younger outskirts when the galaxies are near the LRD phase, i.e., at f_outskirt < 0.4. When f_outskirt > 0.4, radial properties become more homogeneous, suggesting either a transition to a less significant accretion-driven growth phase or a growth phase dominated by dry merger events.

The galaxies analyzed in this study exhibit extremely dense cores that emit predominantly in the NIR, along with non-negligible rest-frame UV emission tracing their extended structures. An intriguing avenue for further investigation is the identification of potential analogs at lower redshifts, around z ≈ 2. Benton et al. (2024) report the existence of star-forming galaxies (SFGs) at z ≈ 2 with high central stellar surface densities, suggesting that dense cores can already be in place in SFGs, similar to those found in quiescent galaxies (QGs). The most similar sources in our sample are those with high outskirt mass fractions (f_outskirt > 0.4) and lower redshifts (z < 3). These galaxies may represent systems evolving into bulge-dominated SFGs, given their already established compact, dense cores and star-forming outskirts.

Further investigation is needed to constraint any AGNs counterpart through spectroscopy in our sample. An Integral Field Unit (IFU) would also provide resolved metallically measurements, offering accurate proxies for stellar ages between the inner part and the outskirt. The study of dust with ALMA as a function of the outskirt mass fraction should also bring some insight into the buildup of dust in this population.

Acknowledgments

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope (JWST). The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1345. The CEERS NIRCam Images used in this article are available on the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute via doi: 10.17909/z7p0-8481. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101148925. We would acknowledge the following open source software used in the analysis: Astropy (Astropy Collaboration 2022), photutils (Bradley et al. 2024), Numpy (Harris et al. 2020).

References

Appendix A: Outskirt detectability and surface brightness dimming

In in this work, we select objects with a young, star forming periphery. Therefore, it would be interesting to investigate wether this outskirts component remains detectable or not due to the surface brightness dimming in the universe. To this end, we first compute the outskirts surface brightness density by measuring the surface of each outskirt (number of pixels) and using the observed magnitude in the F356W filter, which is both one of the reddest (hence better tracing stellar mass) and deepest bands available (with a 5σ point source detection limit of 29.17 mag). The magnitude limit is shown as a dotted line on Figure A.1 and estimated using μ_limit = 29.17 + 2.5log₁₀(2πHWHM²), with HWHM = 0.08″. The effective magnitude limit is also visible as a dashes-dotted line, corresponding to the faintest outskirt in our sample. As the observed flux decreases with redshift following I_obs ∝ I₀(1 + z)⁻⁴, I₀ being the emitted flux, the apparent magnitude increases. Denoting by m₁ the magnitude of a source at redshift z₁, its magnitude at a different redshift z₂ can be expressed as:

$$\begin{array}{c}\hfill m_2=m_1+10{log}_{10}\left(\frac{1+z_1}{1+z_2}\right).\end{array}$$

We then compute the surface brightness μ₂ at redshift z₂ as:

$$\begin{array}{c}\hfill {\mu }_2=m_2+2.5{log}_{10}(S_2)\end{array}$$

where S₂ is the projected area of the outskirt in arcseconds squared at redshift z₂. Using this, we are able to shift our sample through redshift as visible on Figure A.1 from z = 3 to z = 7 and evaluate its detectability. As shown, outskirts remain above the detection limit (i.e., are detectable) for redshifts z < 4, indicating that the entire sample is observable in this range. However, at z = 5, approximately 50% of the outskirts fall below the detection threshold, and this fraction increases to 80% at z = 6. This implies that at z > 5, we are primarily sensitive to the most massive and extended outskirts. Consequently, a significant fraction of the population is likely missed, leading to an underestimation of the true number of such objects at high redshift.

Fig. A.1. Surface brightness density as a function of redshift for the outskirts in our sample (red dots), along with the observed surface brightness shifted to different redshifts (squares), including the dimming effect of the universe. The magnitude detection limit is shown as a dashed line, while the effective magnitude limit is indicated by a dash-dotted line.

All Tables

All Figures

Fig. 1. Color-magnitude diagram illustrating the parent sample: in purple we show the sources that respect the color-magnitude criterion, and black edges denote the parent sample, respecting z > 2.5. In red, we show the color-magnitude selection for the red color of LRDs F277W – F444W > 1.5 used in Akins et al. (2025b), having a zmedian ≈ 6.5. Red and purple contours represent all CEERS sources at 25%, 50% and 80% levels shown for F150W-F356W.

In the text

Fig. 2. Example illustrating the definition of the outskirts and inner regions. Left panel: stacked image combining all JWST bands, defined as STACK = F115W + F150W + F200W + F277W + F356W + F410M + F444W, using PSF-matched data. The dotted line shows the segmentation map, while the solid line traces the ellipse with a semi-major axis of 0.16″ = 2HWHM. Right panel: color composite image using F150W (blue), F277W (green), and F444W (red). The white circle indicates the PSF FWHM beam at F444W. As expected, the galaxy exhibits a centrally concentrated red emission and a bluer, more extended periphery.

In the text

	Fig. 3. Redshift distribution of the final sample in purple. Six LRDs from our sample were also in the Kocevski et al. (2025) LRDs sample, shown in red.
In the text

Fig. 4. Error in percent of the half-light radii Re and Sérsic index n as a function of the inferred half-light radii for mock sources injected into dark areas of the CEERS field of view with mag = 26. The error is defined as ${\text{Err}}_{\text{tot}}=\sqrt{{\text{Err}}_{R_e}^2+{\text{Err}}_n^2}$. The dotted line denotes the HWHM in F444W, whereas the dash-dotted line refers to the resolution limit adopted in this study.

In the text

	Fig. 5. Mass distribution for several samples of LRDs from previous works, assuming only a stellar model, from Rinaldi et al. (2025), Pérez-González et al. (2024), Akins et al. (2025b), Leung et al. (2025), Labbé et al. (2023). This sample is visible in purple for both the inner part and the total mass.
In the text

Fig. 6. Top row: The stellar mass vs. Re plane in each redshift bin of the sample, with the mass-size relation from Allen et al. (2025) overlaid. The relation for quiescent galaxies (QG) from Yang et al. (2025) is shown in orange. Triangular markers represent unresolved sources (for which we consider Re = HWHM), while circular markers indicate resolved sources. Bottom: Distance to the mass-size relation for the entire sample. The cyan square represents the median of the sample, with 15% and 85% percentiles. The red triangles denote LRDs from Leung et al. (2025) which are also unresolved.

In the text

	Fig. 7. Stellar mass density profiles for our sample using only resolved galaxies in purple. In red the profile inferred by Akins et al. (2025b) assuming only a stellar model. We also add dense cores from Hopkins et al. (2010) and Grudić et al. (2019).
In the text

Fig. 8. Color–color selection of LRDs described in Barro et al. (2024a), originally defined for galaxies at z > 4, applied to the inner regions of our sample at z > 4 (shown in red). For galaxies at z < 4, the same rest-frame colors are probed by shifting the bands accordingly (shown in blue). Contours represent the 25%, 50%, and 75% density levels for each distribution. The applied color correction is indicated by a thick black line, with an average shift of ≈0.2 mag.

In the text

Fig. 9. Outskirts’ mass fraction as a function of redshift. Gray dots represent the data, while gray contours illustrate their distribution at 0.25, 0.5 and 0.75 levels. The median is shown as a purple circle, with each point computed in bins of Δz = 1. Vertical error bars indicate the standard deviation, while horizontal bars represent the width within each bin. Red stars refer to the LRDs also in Kocevski et al. (2025).

In the text

	Fig. 10. Half light radii Re as a function of foutskirt. The stars correspond to the median Re from stacking all the filters, and the squares correspond to only F444W. The color bar shows the mean redshift in each foutskirt bin.
In the text

Fig. 11. βUV/βopt plane of the inner regions of the sample. The fraction of the mass outskirt is represented by contours with 0.4 and 0.8 levels, increasing from red to purple and blue. Redshift is indicated by a gray color bar. The LRD selection from Kocevski et al. (2025) is shown in the light gray dashed region and LRDs in common are shown in red stars. In the top right-hand corner, we show βUV,inner as a function of redshift, the colorbar denotes the mass fraction of the outskirt foutskirt, the shaded region refers to the V shape selection.

In the text

Fig. 12. Number density as a function of redshift from our sample visible in red, excluding the six common sources with the LRDs in CEERS from Kocevski et al. (2025), here shown in orange. In purple, the number density of the sources in both samples is shown. The dimming-corrected number density of our sample is shown with dark gray square symbols, interpreted as an upper limit. LRDs from Akins et al. (2025b), Barro et al. (2024a), and Euclid Collaboration: Bisigello et al. (2026) are represented by light gray squares, circles, and triangles, respectively. Horizontal error bars indicate the redshift bins used to compute the number density, and vertical error bars represent Poissonian errors. The theoretical evolution from Pacucci & Loeb (2025) is shown as the dashed black line, while the solid line refers to the model of Inayoshi (2025).

In the text

	Fig. 13. Drawing illustrating the different possible scenarios of evolution of LRDs, as discussed in Section 6: through (a) feedback, (b) mergers, or (c) accretion.
In the text

	Fig. 14. Differences of the ages and sSFR between the outskirt and inner part as a function of the outskirt mass fraction foutskirt. The data are shown in gray with the median profile in pink. For each distribution, the black solid line refers to no radial changes between the inner/outksirt part.
In the text

	Fig. A.1. Surface brightness density as a function of redshift for the outskirts in our sample (red dots), along with the observed surface brightness shifted to different redshifts (squares), including the dimming effect of the universe. The magnitude detection limit is shown as a dashed line, while the effective magnitude limit is indicated by a dash-dotted line.
In the text