A Little Red Dot emerging from a nitrogen-rich UV-luminous galaxy with Very Massive Stars

Introduction

JWST has already delivered transformative results for our understanding of the earliest galaxies and black holes (for a review, see e.g., Matthee [132]; Übler [213]; Stark et al. [199]). In particular, JWST has unveiled an abundant population of luminous galaxies ($M_{\mathrm{UV}} \lesssim-20$) up to $z \approx14$ (e.g., [26]; [151]; [62]; [49]; [225]; [23]; [150]), exceeding the expectations based on pre-JWST predictions (e.g., [14, 15]; [160]). What physical mechanisms are responsible for the prevalence of such UV-luminous galaxies at early times is a subject of debate. Some possible explanations include an increased star-formation efficiency in low-mass halos (e.g., [73]; [62]; [45]); stochastic star-formation at early times, which introduces significant scatter biasing the bright end of the UV luminosity function (e.g., [203]; [196]; [144]; Mason et al. [131]); or the presence of exotic AGN and/or stellar populations with possibly top-heavy initial mass functions (IMFs; e.g., [212]; [233]; [62]; [70]).

Deep spectroscopy of the most luminous galaxies at high redshift has revealed unexpected properties, such as strong UV nitrogen emission, suggesting that many UV luminous galaxies are highly N-enhanced (e.g., [190]; Ji et al. 2024; [27]; [150]; [31]; [153]). In standard chemical evolution models, moderate nitrogen enrichment occurs through core-collapse supernovae (CCSNe). However, significant N-enhancement can only happen in timescales of $\sim100$ Myr, from the winds of asymptotic giant branch stars (AGB; e.g., [218]; [91]). The fact that some of the most N-enhanced systems are observed at only a few hundreds of Myr after the Big Bang demands alternative, faster routes of enrichment (e.g., [76]; [194]; [221]). Some processes that could explain the observed chemical abundances include a high incidence of populations of Wolf-Rayet stars (e.g., Berg et al. [9]; Kobayashi & Ferrara [90]; Senchyna et al. [194]), tidal disruption events (TDEs) in dense star clusters (e.g., Isobe et al. [76]; Cameron et al. [20]; Watanabe et al. [221]), winds of Very Massive Stars (VMSs, $\gtrsim100\,M_\odot$; e.g., Vink [219]; Shi et al. [197]), or supermassive stars (SMS, $\gtrsim10^{3}\,M_\odot$; e.g., Denissenkov & Hartwick 2014; Charbonnel et al. [30]; Nagele & Umeda [146]). Very recent works have reported robust detection of VMSs populations in luminous high-z galaxies (Marques-Chaves et al. [121]; Yang et al. [232]; Chen et al. [31]; Nakane & Ouchi [154]), identified through prominent stellar P Cygni features and broad He ii in the UV spectra (e.g., Upadhyaya et al. [215]). Some of these galaxies are also strong N iv] emitters with high implied N/O enhancement, for example, CEERS-1019 (Larson et al. [104]) at $z = 8.7$ or GN-z11 at $z = 10.6$ (Oesch et al. [161]; Bunker et al. [18]). Therefore, these galaxies could provide an observational link between nitrogen enhancement in early galaxies and their featured stellar populations.

Another surprising and extensively discussed JWST result is the discovery of an abundant population of red compact objects, the Little Red Dots (LRDs; e.g., Labbé et al. [102]; Kocevski et al. 2023; Matthee et al. [133]; Greene et al. [57]). The origin and nature of the LRDs’ central engines is highly debated, with current hypotheses ranging from super-Eddington accreting supermassive black holes (SMBHs) obscured by nuclear dust (e.g., Madau & Maiolino [115]; Scholtz et al. [192]; Billand et al. [11]), to intermediate-mass black holes ($M_{\rm BH} \sim10^{3}$–$10^{5}\,M_\odot$; IMBHs) enshrouded in dense gas that reprocesses their radiation (e.g., Ji et al. 2025; Naidu et al. [148]; de Graaff et al. [43]; Liu et al. [108]; Naidu et al. [147]; Sneppen et al. [198]; Kido et al. [88]), to supermassive stars (e.g., Nandal & Loeb [156]; Zwick et al. 2026; Chisholm et al. [32]; Kokorev et al. [97]; Lu [113]), or quasi-stars (QS; Begelman & Dexter [6]; Martins et al. [127]; Cantiello et al. [21]; Santarelli et al. [187]). Regardless of the physical interpretation, the LRDs seem to represent a common and distinct phase of massive black hole growth at high-z, that becomes increasingly rare towards the nearby Universe (e.g., Lin et al. [106]; Park et al. [168]; Kapoor et al. [84]; Loiacono et al. [112]).

Some LRD systems seem to have very low metallicities, suggestive that they were formed in near-pristine environments (e.g., Maiolino et al. [118]; Ivey et al. [77]), although metal lines are often present in their spectra, such as Fe ii, indicating some metal enrichment (e.g., Lin et al. [106]; Torralba et al. [210]; D’Eugenio et al. [47]; Kokorev et al. [96]). Remarkably, some luminous LRDs show prominent UV lines, including strong nitrogen emission, suggesting that the UV source of some LRDs could share the enrichment patterns of some of the most luminous high-z galaxies (e.g., Labbe et al. [101]; Isobe et al. [75]; Torralba et al. [211]; Tang et al. [207]; Morel et al. [145]; Isobe et al. [75]). Recent works have reported that the chemical abundances of LRDs and high-z broad-line sources in general are compatible with globular clusters in formation (e.g., Isobe et al. [75]; Ji et al. 2026; Chisholm et al. [32]; Kokorev et al. [97]). Constraining the global metal content of LRDs (M/H) further allows to rule out formation pathways. For example, IMBHs with masses up to $\sim10^{5}\,M_\odot$ can form through direct collapse of near-pristine gas, if a nearby luminous source of UV photons is present that is capable of dissociating H$_2$, suppressing molecular cooling and preventing fragmentation (e.g., Omukai [164]; Oh & Haiman 2002). For a massive black hole to be formed through direct collapse and rapid accretion, the metallicity of the gas must not exceed $\sim10^{-3}$–$10^{-4}\,Z_\odot$ (e.g., Chon & Omukai [36]; Chon et al. [35]). This process could lead to the formation of an SMS as an intermediate step (e.g., Bromm & Loeb [17]; Agarwal et al. [1]; Latif et al. [105]); and/or a dense nuclear star cluster (e.g., Omukai et al. [165]; Devecchi & Volonteri [48]; Gieles et al. [56]) which would ultimately lead to “seeding” of a massive black hole. The direct-collapse scenario has received attention in the LRD context since many of these objects show nearby companions that could have acted as the UV sources needed to supress molecular cooling (Bagger et al. 2026; Barger & Cowie [5]; Pacucci & Urry [166]).

Observations and data reduction

COLA1: A double-peaked Ly$\alpha$ emitter at $z = 6.6$

Here we present a spectroscopic analysis of COLA1 (Hu et al. 2016; Matthee et al. [136]), a luminous galaxy at $z = 6.6$ that we found to be a strong N-emitter, presents signs of containing VMSs, and shows broad lines in the optical spectrum (Meyer et al. [142]) similar to LRDs. The new JWST/NIRSpec data of this source unifies three of the most puzzling recent findings from JWST into a single source. The structure of this paper is organized as follows. In Sect. 2 we describe the observations and data reduction. In Sect. 3 we present our empirical results, including stellar population fitting, emission line modeling and chemical abundance estimation. In Sect. 4 we discuss physical interpretations of our results. Finally, in Sect. 5 we summarize our main findings. Throughout this paper we adopt a ΛCDM cosmology as described by Planck Collaboration et al. [173]; $\Omega_\Lambda= 0.69$, $\Omega_m = 0.31$, $H_0 = 67.7\ {\rm km\,s^{-1}\,Mpc^{-1}}$. All photometric magnitudes are given in the AB system (Oke & Gunn [163]).

JWST/NIRSpec

COLA1 was observed as primary target by the JWST Cycle 4 program #7722 (PIs: Torralba & Matthee), using NIR-Spec/MOS. Spectroscopy was taken using the grating/filter combinations G140M/F070LP and a G235M/F170LP, with a total exposure time of 25 209.6 s per grating ($\sim7$ h). A 3-point nodding pattern was used. The MSA data were reduced using the latest version of the software msaexp$^{1}$ (v0.9.17; Heintz et al. [65]; de Graaff et al. [41]), which builds upon the jwst pipeline$^{2}$ (v1.16.1). We use the CRDS context file jwst_1464.pmap. We use the extended-wavelength reduction, which extracts the spectra from the full coverage beyond the nominal range of the gratings, extrapolating the flat-field and photometric reference files so that the trace can be extracted across the full detector. The flat-field was undone at the calibration stage and applied at extraction, and the empirical read-noise rescaling was turned off. The 1D spectra were extracted with optimal extraction ([68]) using a fitted profile, with the background obtained by differencing the nod positions and the trace position fixed from the MSA source position and nod offsets. The MSA mask was designed so the trace of COLA1 does not overlap with that of any other object, and so COLA1 is placed in an optimal position on the detector to retrieve up to the observed H$\alpha$ wavelength on G235M. We show the spectra in Figure 1, highlighting the most relevant features.

By combining the two gratings observed by the program GO #7722 we obtain a complete coverage of the rest-UV-to-optical spectrum of COLA1; from $\sim1000$ to $\sim7000$ Å rest wavelength. However, the detector gap falls on top of the critical H$\gamma$ and [O iii] $\lambda4363$ emission lines. We complement our analysis with the JWST NIRSpec/IFU data taken by GO #3767 (PI: Meyer); 40 min on-source exposure time with the grating G395M. Those data provide a sufficient S/N of the H$\gamma$+[O iii] $\lambda4363$ emission lines. We retrieved the IFU datacube from the DAWN JWST Archive⁴. The 1D spectrum was obtained using the optimal extraction method ([68]), with a 2D Gaussian extraction kernel, fitted to the morphology of COLA1 in the whitelight image of the IFU cube. The line fluxes of the overlapping MSA (G235M) and IFU (G395M) generally agree within uncertainties. For consistency, we rescale the flux of the IFU 1D spectrum to match the [O iii] $\lambda\lambda4960, 5008$ flux expected from slit-loss correction uncertainties.

JWST/NIRCam imaging and morphological analysis

We also use photometric images from JWST/NIRCam to measure the UV size of COLA1. In particular, we use imaging in the F115W, F150W, and F200W filters from the JWST GO program #1933 (PIs: Matthee & Naidu; presented in Torralba et al. [209]). We fit the morphology of COLA1 with a Sérsic profile, and we use a second Sérsic to model the faint foreground contaminant (see Torralba et al. [209]). We fit using the pysersic software (Pasha & Miller [169]) and the synthetic PSFs from webbpsf (Perrin et al. [171]). We find effective radii of $R_{\rm eff} = 209 \pm18$ pc and $R_{\rm eff} = 193 \pm15$ pc in the F115W and F150W bands, respectively; the size is completely unresolved in F200W and F356W. Therefore, we adopt a fiducial UV size of $\approx200$ pc. Note that this size is slightly smaller than the PSF half-radii ($\approx230$ pc and 250 pc for F115W and F150W, respectively), so COLA1 is only marginally resolved in the rest-frame UV.

ALMA far-IR upper limits

ALMA observations of COLA1 were conducted in Band 6 during Cycle 7 (project 2019.1.00741.S, PI: Matthee), using compact-to-intermediate 12m array configurations with 27–46 antennas across six execution blocks and a total on-source time of $\sim266$ min. The spectral setup consisted of two sidebands, each composed of two spectral windows with total bandwidths of 3.75 GHz and channel width resolutions of 7.8 MHz, and the main one aimed at the redshifted [C ii] $\lambda158\,\mu$m line. The observations have a synthesized beam of $0.374^{\prime\prime}$. The data were calibrated using the automatic pipeline provided by the Common Astronomy Software Applications package (CASA; McMullin et al. [140]), version 6.2.1.7.

Property	Value	Units
$M_{\mathrm{UV}}$	$-21.36 \pm0.01$	–
$\beta_{\mathrm{UV}}$	$-2.5 \pm0.1$	–
$\beta_{\mathrm{opt}}$	$-1.1 \pm0.7$	–
$r_{\mathrm{eff,UV}}$	$0.21 \pm0.02$	kpc
$L_{\mathrm{H}\alpha}$ (Total)	$9.9^{+0.8}_{-0.3} \times10^{42}$	erg s$^{-1}$
$L_{\mathrm{H}\alpha}$ (Narrow)	$2.10^{+0.28}_{-0.14} \times10^{42}$	erg s$^{-1}$
$f_{\mathrm{esc}}^{\mathrm{Ly}\alpha}$	$0.67 \pm0.3$	–
$\xi_{\mathrm{ion}}/(1-f_{\mathrm{esc}}^{\mathrm{LyC}})$	$10^{25.56\pm0.02}$	Hz erg$^{-1}$
$12+\log_{10}(\mathrm{O}/\mathrm{H})$	$7.62 \pm0.01$	–

From the non-detection in the integrated continuum band, we estimate a continuum upper $3\sigma$ limit of $<24.6$ µJy beam$^{-1}$ at observed 242 GHz. The [C ii] $\lambda158$ µm is also undetected. We estimate an upper $3\sigma$ limit of $L_{\mathrm{[C\,II]}} < 2 \times10^{7}\,L_{\odot}$, assuming a top-hat line with FWHM $=125$ km s$^{-1}$. This upper limit corresponds to an upper limit of $\log_{10}(\mathrm{SFR}/M_{\odot}\,\mathrm{yr}^{-1}) \lesssim0.2$ using extrapolations of the standard calibrations (e.g., De Looze et al. [44]; Schaerer et al. [189]), but see Harikane et al. [63]; Matthee et al. [135] for discussions whether these could be applicable to LAEs. For comparison, the narrow H$\alpha$ emission implies a somewhat higher $\log_{10}(\mathrm{SFR}/M_{\odot}\,\mathrm{yr}^{-1}) \sim0.9$ using the calibrations from Kramarenko et al. [99].

Results

Flux measurements of UV and optical lines

The data shown in Figure 1 reveals strong UV emission lines including N iv] $\lambda\lambda1486$, 1483, N iii] $\lambda1750$ complex, C iii] $\lambda\lambda1907$, 1909, broad He ii $\lambda1640$ (Figure 2), P Cygni features in N v $\lambda\lambda1239$, 1243, Si iv] $\lambda\lambda1394$, 1403, and C iv $\lambda\lambda1548$, 1551, and a very blue continuum slope ($\beta_{\mathrm{UV}}=-2.5\pm0.1$). In the rest-optical, it shows broad Balmer lines, H$\alpha$ absorption, and a redder continuum ($\beta_{\mathrm{opt}}=-1.1\pm0.7$).

We use the Python package unite (Hviding [71]) to measure the emission lines of the G140M and G235M spectra of COLA1. We first performed a preliminary fit to the most luminous narrow emission lines in the rest-optical (i.e., [O iii] $\lambda\lambda4960$, 5008, [Ne iii] $\lambda3869$, $\lambda3968$) to pinpoint COLA1’s systemic [O iii] redshift at $z=6.59054\pm0.00002$. Then, we fit every emission line as a single Gaussian, adding a second broad Gaussian component when needed. The absorptions of the lines showing P Cygni profiles are fitted as Gaussians with negative fluxes and blueshifted velocities with respect to the emission counterparts. Broad components are added after a visual inspection, to fit the lines N iv] $\lambda\lambda1483$, 1486, C iv $\lambda\lambda1548$, 1559, He ii $\lambda1640$, [Ne iii] $\lambda3869$ and [O iii] $\lambda\lambda4960$, 5008. In addition, the H$\alpha$ and H$\beta$ emission lines present broad features, but these will be treated separately in Sect. 3.4. The full list of measured emission line fluxes is listed in Tables A.1 and A.2. The fitted components for all the emission lines are shown in Figs. A.1 and A.2.

Gas-phase O/H metallicity

From the emission line measurements we estimate the single-zone gas-phase metallicity of COLA1 using the direct method (e.g., Peimbert [170]). We obtain a mean electron density of $\log_{10}(n_{\mathrm{e}}/\mathrm{cm}^{-3})=3.4^{+0.3}_{-0.4}$ from the [O ii] $\lambda\lambda3727$, 3729 doublet.

Another density tracer is the UV doublet C iii] $\lambda\lambda1906$, 1909, however the measured ratio is close to the theoretical maximum, so it lies in a regime where it is hardly constraining (see Sect. 3.5). We estimate a direct-method temperature of $(2.2\pm0.2)\times10^{4}$ K from the [O iii] $\lambda4363$ and [O iii] $\lambda\lambda5008$ fluxes using PyNeb (Luridiana et al. [114]). Finally, we estimated a gas-phase metallicity of $12+\log_{10}(\mathrm{O}/\mathrm{H})=7.62\pm0.01$. From this result, we adopt a fiducial metallicity of $Z=0.1\,Z_{\odot}$ throughout this paper. However, note that this value should be interpreted with caution, given the evidence for a multiphase ISM with heterogeneous enrichment within COLA1 (see Sect. 3.5). None of our key interpretations rely on a precise measurement of the oxygen abundance.

Stellar population synthesis fitting

We model the UV spectrum of COLA1 with a combination of single stellar population models (SSP) including VMSs of up to 300 $M_{\odot}$ (Martins et al. [125]). We use the stellar population fitting code FICUS (Chisholm et al. [34]; Saldana-Lopez et al. [185]), following the procedure described by Marques-Chaves et al. [121]. Briefly, FICUS models the UV spectrum as a linear combination of SSPs, together with the nebular continuum computed with PyNeb (Luridiana et al. [114]). We consider two sets of stellar population models: BPASS models with a maximum stellar mass of 100 $M_{\odot}$, and the models from Martins et al. [125], which include VMSs with masses up to 300 $M_{\odot}$. All models adopt an extended Salpeter IMF with mass exponent $\alpha=2.35$.

For the VMS models, we explore two prescriptions for stellar mass loss following Martins et al. [125]: a prescription (VMS_GR), corresponding to the mass-loss rates inferred for LMC stars, and a metallicity-dependent prescription (VMS_GRZ), in which the mass-loss rates scale linearly with metallicity relative to the LMC values. We consider combination of SSPs, considering 10 models with ages ranging from 0 to 40 Myr. The nebular continuum is computed with PyNeb assuming an electron density of $n_{\mathrm{e}}=2500$ cm$^{-3}$ and an electron temperature of $T=2.2\times10^{4}$ K. We adopt a Reddy et al. [178] dust attenuation law. We also use the models in Martins et al. [125] that include the late stages of evolution of the VMSs. These models explicitly include the advanced stages ($<10\%$ of the star lifetimes) where VMSs are predicted to emit increased high-ionization photons, particularly He ii ionizing continuum that is relevant to assess the nebular emission of He ii $\lambda1640$ and N v 1240, among others [124, 125]. These models are explored in more detail in Sect. 3.6.

The observed spectrum is fitted in the range 1200–2000 Å assuming a metallicity (of both stars and ISM) of $0.1Z_{\odot}$, matching the inferred gas-phase metallicity of COLA1. Since COLA1 is a known strong leaker of hydrogen ionizing photons [136, 209] we adopt an escape fraction of $f_{\mathrm{esc}}^{\mathrm{LyC}}=0.6$, which translates to downscaling the nebular continuum that FiCUS computes as 60% assuming a picket-fence escape of ionizing photons. This choice for $f_{\mathrm{esc}}^{\mathrm{LyC}}$ is motivated in Sect. 4.7. Note that a significantly lower $f_{\mathrm{esc}}^{\mathrm{LyC}}$ would produce a strong Balmer jump that is not observed. Our results are robust against small differences in the assumption of the gas metallicity or $f_{\mathrm{esc}}^{\mathrm{LyC}}$; these would only marginally affect the global $\chi^2$ of all the models due to different infilling of the P Cygni absorptions from nebular continuum.

We find a best-fit stellar model with mean age $2.43 \pm0.75$ Myr, and stellar mass $\approx10^8\,M_{\odot}$ (assuming an IMF with slope $\alpha=2.35$; Salpeter [186]); but we note that $>50\%$ of the UV light is emitted by $\sim3 \times10^4$ stars with masses $>100\,M_{\odot}$. The best-fit models for all the populations tested are shown in Figure 3. The best fit is given by the VMS model with GRZ tracks, followed by the VMS models with GR tracks ($\Delta\mathrm{BIC}=22$–26), and performing significantly better than the BPASS models with stars up to only $100\,M_{\odot}$ ($\Delta\mathrm{BIC}=53$). The models including VMS with a low metallicity, $0.01Z_{\odot}$, perform fairly worse

(ΔBIC = 119) due to the lack of absorption features in the stellar spectra. Therefore, the presence of strong P Cygni profiles suggests that the metallicity of the stars in COLA1 cannot be much lower than 10% the solar metallicity.

Balmer line modeling

The Hα and Hβ lines of COLA1 show clear broad component emission (Meyer et al. [142]). In addition, Hα shows a blueshifted Balmer absorption (see Figure 5), feature that is characteristic of a large fraction of LRDs (e.g., Matthee et al. [133, 137]; Juodžbalis et al. [82]; Yanagisawa et al. [230]; Davis et al. [40]). We fit the Hα and Hβ lines with the models described in Torralba et al. [210] and Matthee et al. [137], tailored to LRD’s Balmer lines. Briefly, our model consists of a narrow component at the line center, an intermediate-width Gaussian, and broad wings described by the intermediate component convolved with a symmetric exponential kernel simulating electron scattering. We also add a Gaussian absorber that affects the broad component and the continuum (see Matthee et al. [137], for a more detailed description). Alternatively, we also fit the same model but replacing the broad exponential wings with a broad Gaussian, and a Lorentzian (Fig. B.1). An exponential broad component is preferred over a Lorentzian (ΔBIC = 13) or

!Hα, Hβ, and [O iii] profiles of COLA1 bbox=0.06,0.07,0.49,0.43

a broad Gaussian (ΔBIC = 43), as found in many LRDs (e.g., Rusakov et al. [182]; Scholtz et al. [192]).

Independent of the model choice, we measure the broad and narrow component fluxes of the Balmer lines with the procedure described in Sect. 3.1. We fix the intrinsic width of the narrow components of the balmer series from Hα to H(10-2), and we add an arbitrary number (3) of broad Gaussian components to describe the broad Hα and Hβ emission, with the only aim of measuring robust integrated fluxes (reported in Table A.2). Note that for higher-order Balmer transitions no broad component is detected due to sensitivity limitations. In Figure 6 we show the Balmer decrements of the narrow components with respect to Hβ. The Balmer decrements of the narrow lines across the series are compatible with Case B recombination with negligible dust attenuation (see also Figure 5), compared to the expected decrements using a Cardelli et al. [22] dust attenuation law. In turn, the decrement of the broad components is Hα/Hβ(broad) ≳ 9, far from the 2.86 Case B value. The broad components of Hγ and successive lines in the Balmer series undetected in our data, possibly due to their low S/N.

The complex broad Balmer line profiles of LRDs could be arising from radiative-transfer effects in a thick wind (e.g., Rusakov et al. [182]; Naidu et al. [148, 147]; Sneppen et al. [198]; Torralba et al. [210]; Chang et al. [29], but see Madau et al. [116] for an alternative interpretation). We test this hypothesis in order to show that the high Balmer decrement of COLA1 can arise from column density effects in a thick wind rather than localized dust attenuation. For that, we run a grid of Sirocco (Matthee et al. [132]) models (see also Sneppen et al. [198]; Kaaz et al. [83]; Mascia et al. [128]). For simplicity, we simulated a homologous wind, analogous to a spherically expanding layer from a Type IIn supernova (e.g., Naidu et al. [147]), with a central blackbody source, and wind metallicity $Z = 0.1\mathrm{Z}_{\odot}$. We ran a grid with parameters spanning $\log_{10}(T_{\mathrm{BB,in}}/\mathrm{K}) \in[4.7, 4.8, 5.0]$, $\log_{10} U \in[0, 0.25, 0.5, 0.75, 1]$, base density $\log_{10}(n_{\mathrm{base}}/\mathrm{cm}^{-3}) \in[6, 7, 8]$, wind maximum velocity $v_{\max} \in[200, 250, 300, 350]~\mathrm{km\,s}^{-1}$, and column density $N_{\mathrm{H}} \in[24.5, 24.8, 25.0, 25.2, 25.5]\ \mathrm{cm}^{-2}$. The expanding wind has a density law $n \propto r^{-2}$. We fit to the broad wings and absorption of the continuum-subtracted Hα line, after masking the central narrow emission (masked $|v-v_{\mathrm{H}\alpha,0}|<400\ \mathrm{km\ s}^{-1}$). For each grid model, the only free parameter for the fit is an arbitrary flux scale $A$.

As shown in Figure 7, the best-fit model is capable of reproducing the broad Hα wings and blueshifted absorption. That same configuration predicts a Hβ profile with very weak wings as compared to the Case B recombination expectation, matching the observations. Moreover, this configuration also produces significant narrow emission from the wind with a high decrement, which we interpret as the reason why Hα/Hβ is higher than expected without the need of any dust attenuation. Indeed, integrating the residuals after subtracting the best-fit model to the observed COLA1 spectrum, the flux ratio is Hα/Hβ $\sim3$, consistent with the extra emission coming from a different region governed by Case B recombination and no dust attenuation. While the Hα absorption is clearly detected, Hβ does not show it clearly. However, this is expected due to the predicted wind component of Hβ being weaker, the signal-to-noise is lower, and the absorption could be totally infilled by the strong narrow component. After finding a best fit model, we ran a finer grid in column density and $\log_{10} U$ to obtain a best-fit model with $\log_{10}(T_{\mathrm{BB,in}}/\mathrm{K})=5$, $\log_{10} U=0.2$, $\log_{10}(n_{\mathrm{base}}/\mathrm{cm}^{-3})=7$, $v_{\max}=300\ \mathrm{km\ s}^{-1}$, and $\log_{10}(N_{\mathrm{H}}/\mathrm{cm}^{-2})=24.8$ (shown in Figure 7).

Nitrogen and Carbon abundance

In the following sub-sections we estimate the relative abundances of oxygen, nitrogen, and carbon in the ISM of COLA1. We use two complementary approaches: a single-zone method based on PyNeb line emissivities (Sect. 3.5.1), and a two-zone density photoionization model using Cloudy (Sect. 3.5.2).

Single-zone CNO abundances

The low ratio N IV] 1483/1486 $\approx0.09$ implies an electron density of $n_{\mathrm{e}}\gtrsim10^{6.5}\ \mathrm{cm}^{-3}$ (e.g., Keenan et al. [86]; Curti et al. [38]), while the high C III] 1906/1909 $\approx1.7$ needs $n_{\mathrm{e}}\lesssim10^{3.5}$ (e.g., Schmidt et al. [191]). This suggests that the ISM of COLA1 is composed of a two-zone density system, which could have different chemical abundances, as found for other young starbursts in the literature (e.g., Martinez et al. [123]; Arellano-Córdova et al. [2]; Chen et al. [31]).

We first measure the total abundances using a single-zone density model for direct comparison with other sources in literature. We estimate the nitrogen and carbon abundances from the UV emission line fluxes by comparing to the emissivities from PyNeb (Luridiana et al. [114]), adopting $T_{\mathrm{e}}=2.2\times10^{4}\ \mathrm{K}$, and electron density $n_{\mathrm{e}}=10^{4}\ \mathrm{cm}^{-3}$ as estimated from the optical oxygen lines (see Sect. 3.2); the results depend weakly on these parameters. We use the narrow N IV] 1486 flux as a tracer of N$^{+3}$, the total flux of the N III] $\lambda\lambda1750$ quintuplet for N$^{+2}$, C IV $\lambda\lambda1548,1551$ for C$^{+3}$, and C III] $\lambda\lambda1906,1909$ for C$^{+2}$. C IV has a stellar component to the flux, but subtracting the P Cygni profile of the best-fit stellar population from the observed data only affects the measurements by $\sim0.01$ dex, well within the $1\sigma$ errors. From the observed positive line fluxes we obtain

hence

If we use instead a higher density of $10^{6.5}\ \mathrm{cm}^{-3}$ as the N IV] doublet suggests, the measured abundances change by only roughly 0.01 dex. We have assumed an ICF of O$^{+2}$/O $=1$ (see e.g., Martinez et al. [123]). The nitrogen ionic population of N$^{+3}$ dominates over N$^{+2}$, although we have disregarded N$^{+4}$ and higher states, assuming that lines are powered only by collisional excitation of the hot photoionized gas. We detect some nebular N V $\lambda1240$ emission (probably photoionized by $>77\ \mathrm{eV}$ photons), therefore our nitrogen abundance might be conservative (see Sect. 3.6).

Two-zone CNO abundance photoionization model

Now turning to a two-zone approach, we fit the relative fluxes of all the C, N, and O emission lines across the UV-to-optical spectrum to emissivities using the photoionization code Cloudy (last described in Gunasekera et al. [60]). We ran a Cloudy grid of slabs with metallicity [M/H] = $-1$ and default abundances, density $\log_{10}(n_{\mathrm{H}}/\mathrm{cm}^{-3})$ from 2.5 to 8 in 0.5 dex steps, ionization parameter ranging $\log_{10} U \in[-3, -2.5, -2, -1.5, -1, -0.5, 0]$. We also vary the relative N and C abundances in a grid of multiplicative factors $\times\mathrm{N} \in[1, 2, 4, 6, 8, 12, 16, 20, 25]$ and $\times\mathrm{C} \in[0.03, 0.05, 0.1, 0.2, 0.4, 0.7, 1, 2]$. This is motivated by multiple recent works finding anomalous N/C abundance patterns regularly in high-$z$ galaxies (e.g., Rai & Roberts-Borsani [175]; Isobe et al. [75]; Shi et al. [197]; Berg et al. [9]). We run Cloudy with our best-fit stellar population SED fitted to the UV as an input (Sect. 3.3). We use the stopping condition stop efrac $-1$, to simulate a radiation-bounded cloud. We fit the relative fluxes of the main CNO lines of COLA1, measured as described in Sect. 3.1, using only the narrow components of each line. We use the total positive flux of C iv without splitting between broad and narrow components due to its uncertain radiative transfer. We fit the relative line fluxes as a linear combination of any two Cloudy models of the grid, fitting for respective weights $0 < w < 1$. We also use a Calzetti et al. [19] dust attenuation curve, fitting for $E(B-V) \in[0, 0.4]$. We find the best-fit using a differential evolution algorithm (Storn & Price [201]) implemented in the Python package lmfit (Newville et al. [158]). In order to estimate the uncertainty of the best-fit parameters, we perform a $\chi^2$ profiling by fixing each parameter in an array of values, and re-fitting for the rest of them, finding the best minimum $\chi^2$ for a fixed parameter. We show the $\chi^2$ profiles for each parameter, and the fitted emissivities per density zone for the modeled lines in Figure 8 and Figure 9, respectively.

The model again converges to two different zones with density $\log_{10}(n/\mathrm{cm}^{-3}) = 7.5^{+0.43}_{-0.5}$ and $\log_{10}(n/\mathrm{cm}^{-3}) = 3.18^{+0.32}_{-0.19}$. We find a best-fit dust attenuation parameter $E(B-V) = 0.14^{+0.11}_{-0.05}$, corresponding to a very low dust obscuration. While no dust attenuation is preferred by the narrow Balmer decrements (Figure 6), this value could be accomodated especially considering a complex dust geometry around the line-emitting regions. We estimate a nitrogen abundance of $\log_{10}(\mathrm{N/O})_{\mathrm{dense}} = 0.00^{+0.13}_{-0.10}$ for the dense zone. This zone requires a very low $\log_{10}(\mathrm{C/O})$ in order to not overproduce the C iii] $\lambda1909$ component, so we can only estimate an upper limit to this quantity $\log_{10}(\mathrm{C/O})_{\mathrm{dense}} < -1.4$ ($2\sigma$), yielding $\log_{10}(\mathrm{N/C})_{\mathrm{dense}} > 1.3$ ($2\sigma$; see Figure 8). The C iv doublet should contain non-negligible contribution from stellar origin (see Figure 3), and due to its resonant nature, the RT effects are uncertain for this line. However our key results do not strongly depend on a precise measurement of C iv, since our models predict that most of the C iv emissivity is produced in the diffuse zone (see Figure 9). In turn, for the diffuse zone the nitrogen abundance is unconstrained, given that the N iv] $\lambda1483$ component is undetected and blended with the broad component, hence we measure $\log_{10}(\mathrm{N/O})_{\mathrm{diff}} < -0.45$ ($2\sigma$), while the carbon abundance is well defined by the C iv and C iii] lines; thus $\log_{10}(\mathrm{C/O})_{\mathrm{diff}} = -0.28^{+0.21}_{-0.09}$ and $\log_{10}(\mathrm{N/C})_{\mathrm{diff}} < 0.02$ ($2\sigma$). All the measured single- and two-zone abundances are listed in Table 2.

Ratio	Dense zone	Diffuse zone	Integrated
$\log_{10}(\mathrm{N/C})$	$> +1.3$	$< +0.02$	$+0.58^{+0.12}_{-0.10}$
$\log_{10}(\mathrm{N/O})$	$0.00^{+0.13}_{-0.10}$	$< -0.45$	$-0.07^{+0.04}_{-0.05}$
$\log_{10}(\mathrm{C/O})$	$< -1.4$	$-0.28^{+0.21}_{-0.09}$	$-0.64^{+0.11}_{-0.09}$

Notes. Upper or lower limits are given for a $2\sigma$ confidence level. Dense and diffuse zone abundances are obtained through the two-zone cloudy emissivity model, and the integrated independently from PyNeb emissivities (see Sect. 3.5).

UV emission line photoionization modeling

In Sect. 3.3 we fitted the UV continuum using stellar population synthesis spectra, however we do not include nebular emission lines in our FiCUS models. We fit the UV lines of COLA1 separately. We then use these values as the starting point for an MCMC analysis, in which the model grid is interpolated to evaluate the likelihood at arbitrary parameter values. This allows us to explore the parameter space continuously and obtain a more precise estimate of the best-fitting parameters, together with robust uncertainties.

The spectral fit is performed in the range of rest-frame wavelengths 1200–1950 Å. We mask the C iv and Si iv emission lines in this process, as they present P Cygni stellar wind features that could also be impacted by interstellar absorption, the latter processes not captured by these Cloudy models. As shown in the first panel of Fig. C.1, the best-fit using our fiducial stellar population as an input for the Cloudy runs is capable of reproducing the main UV nebular lines (N iv, O iii], N iii], C iii]). The two-zone model converges to a high-density component with a very elevated $\log_{10}(\mathrm{N/C})$ and a more diffuse zone with less extreme abundances, to reproduce the observed N iv] and C iii] ratios. The broad He ii $\lambda1640$ component arises from the stellar SED, but there is a lack of He ii-ionizing photons from this stellar population model to power the relatively strong narrow He ii (IP 54.42 eV), which we interpret as nebular. Moreover, the narrow component of the N v $\lambda\lambda1240$ doublet (IP 77.47 eV) is also under-estimated, requiring an ionizing spectrum with a large amount of high-energy photons.

We also also consider the possibility that there is an additional source of ionizing photons such as an AGN accretion disk. We run an additional Cloudy grid including an AGN SED⁵ mixed with the fiducial stellar model, at fractions of [0, 12.5, 25, 50, 75]% of the UV luminosity at 1500 Å. The best-fit model is found with an AGN fraction of 12.5% (fifth panel of Fig. C.1). Including an AGN overshoots He ii without producing enough N v flux. These results therefore point towards a very hard ionizing far-UV spectrum, but with a cutoff not exceeding by much the production IP of N v ($\lesssim80$ eV). No other high ionization lines are found across COLA1’s spectrum. For example, [Ne v] $\lambda3426$ (IP 97.11 eV) is undetected with EW $< 1.2$ Å at $3\sigma$ (compare to e.g., Chisholm et al. [33]; McKaig et al. [139]; Reiss et al. [180]). The Ne iv $\lambda2424$ is also undetected with EW $< 1.4$ Å at $3\sigma$, although its IP of 63.45 eV is lower than that of N v, this line is expected to be luminous in typical AGN-dominated sources (e.g., Castellano et al. [27]; Scholtz et al. [193]).

It is well established that stars with lower metallicity produce harder ionizing spectra (e.g., Schaerer [188]). We repeat the fitting procedure with models using the 0.01$Z_{\odot}$ stellar population model (see Sect. 3.3). In this case, the nebular He II is well reproduced, but the fit to this line is noticeably worse, overshooting the broad component. Additionally, this SED is still incapable of producing significant N V emission, and the overall continuum fit is worse. The models including the advanced phases of VMS evolution of Martins et al. [125] produce enough high energy photons to power both nebular He II and N V, despite the over-

all fit being statistically disfavored ($\Delta$BIC = 122). These results highlight the uncertainties of the stellar models of very massive stars, but they show that under certain assumptions these stars could have the ionizing properties needed to power the observed high-ionization UV emission lines (Martins et al. [125]).

Discussion

A Little Red Dot outshone by a population of VMSs

The UV of COLA1 shows all the classical signatures of a young population of massive stars, while the optical continuum and the luminous broad Balmer lines are typically associated with an AGN or an LRD. The strong N iv] emission indicates the presence of highly nitrogen-enriched gas (e.g., [234, 31, 153, 150]). The absence of low-ionization absorption lines (e.g., Si ii, O i, C ii, Al ii), the strong double-peaked Ly$\alpha$ emission and weak Balmer jump are indications of a high ionizing continuum escape fraction (e.g., [78, 87, 153]).

The optical continuum of COLA1 cannot be described with the same SSP model that reproduces the UV, requiring an extra component. In Figure 10 we show the complete G140M+G235M spectra, degraded to the resolution of NIRSpec/PRISM ($R \sim100$). We compare to the NIRSpec/PRISM stack of LRDs from Matthee et al. [137], and to the prototypical “Little Blue Dot” GS_3073 ([214, 16]). COLA1 has a redder optical slope than the LRD stack, and the overall UV-to-optical continuum is more similar to that of GS_3073. However, a critical difference is that COLA1 shows Balmer absorption in the H$\alpha$ line, which LBDs lack. Other notable differences can be seen between the general UV shapes of COLA1 and GS_3073, with the latter showing relatively stronger UV emission lines (see Ji et al. (2024)). We also compare to a median stack of star forming galaxies (SFGs), selected from the JWST Dawn archive by imposing a secure visual inspection redshift $z > 5$, and blue UV slopes $\beta< -2$. As one can see in Figure 10, the optical continuum of typical blue SFGs remains blue, while COLA1 shows a slight turnover, that can be interpreted as an extra component in addition to the young stellar light.

The residual of subtracting the best-fit stellar model (fit only to the UV) from COLA1’s spectrum reveals a continuum component that dominates the rest-optical (first panel of Figure 11), and has a turnover near the Balmer limit (H$_{n=2\to\infty}$; 3645 Å). This component could in principle be identified as an older stellar population with a Balmer break, or a dust-obscured AGN. However, we find a more compelling interpretation that the residuals trace a red, compact component, such as the optical continuum of an LRD. We favor this interpretation because we also detect broad Balmer emission lines with extreme Balmer decrements, as well as blueshifted H$\alpha$ absorption—the latter being one of the defining properties of LRDs rarely seen in any population of AGN (e.g., [137, 230, 195] (Juodžbalytė et al. 2026)). Some works have shown that the characteristic UV-to-optical “V-shape” of many LRDs can be decomposed into a blue star-forming galaxy component and a red source similar to the red component we find here dominating the optical light (e.g., [148, 205, 42, 222] (Ji et al. 2025)). We propose that the spectrum of COLA1 can be well described by the same scenario, except that here we catch the stellar population at an exceptionally young age.

Evolution of the stellar population – COLA1 as a rare snapshot of an emerging LRD

Due to the steep stellar mass-luminosity relation (e.g., [183, 66, 50]), the UV spectrum of COLA1 is dominated by the light of the most massive stars, which have a short lifetime ([125]). If we assume no significant star formation occurs after the observed burst, the UV will fade as the most massive stars end their lives in a timescale of 2–5 Myr. We test this directly by evolving our best-fit stellar population. In the top right panel of Figure 11 we show the stellar model evolved by $\Delta t = 5$ Myr into the future. In order to build the evolved spectrum, we used the SSP models that correspond to the ages of the two starburst of our best fit (see Figure 1) plus 5 Myr, with the same initial mass, weighted by the observed light fractions at the present epoch, essentially preserving the star-formation history. That is, we replace the population at 3 Myr by one at 8 Myr, and since 7 Myr is not available in our library, we linearly interpolate the models with 5 Myr and 8 Myr. We also add nebular continuum and emission lines by giving the evolved SSPs as a input to Cloudy, with the same parameters as the fiducial $\Delta t = 0$ model, and re-scaling the ionization parameter by the decrease of ionizing continuum ($Q_{\rm H^+}$) in the integrated stellar population model SEDs. The bottom panel of Figure 11 shows the $M_{\rm UV}$ evolution of our best-fit SSP model up to 25 Myr.

The UV luminosity of the aged stellar population fades as the most massive stars disappear. After $\sim5$ Myr, our evolutionary toy model will leave a UV-to-optical spectrum that fully resembles an LRD with a “V-shape” ($\beta_{\rm UV,aged} = -1.7 \pm0.1$ and $\beta_{\rm opt,aged} = +0.1 \pm0.1$), and would pass standard LRD selection criteria (e.g., [42, 72] (Kocevski et al. 2026)). Here we assume that the lifetime of the red component is longer than that of the VMSs and that it does not evolve significantly during this process, although in a realistic scenario some evolution might be expected (see Sect. 4.5 for a more detailed discussion about timescales).

The appearance of dust can make the $M_{\rm UV}$ decline over time even steeper. We illustrate this with a very simple toy model of dust production by SNe. We assume that every $> 8\,M_\odot$ star undergoes a supernova at the end of its life producing 0.02 $M_\odot$ of net dust mass, after reverse-shock destruction ([10]). In each time step we add this dust mass to a spherical shell with radius 200 pc (roughly the UV size of COLA1). The UV attenuation is estimated using a mass absorption coefficient of $\kappa_{1500} = 5 \times10^4\ {\rm cm}^2\,{\rm g}^{-1}$ (for a Milky Way attenuation law, $R_V = 3.1$; [223]). We show the effect of applying this dust model evolution as a dashed line in the bottom panel of Figure 11.

The UV fading also depends heavily on the IMF. The SSP models that we employ assume a slope of $\alpha= 2.35$ ([186, 124, 125]). To illustrate the effect of a top-heavy IMF we adopt a parametrised description to modify the available SSP. In this toy model, we assume that the starburst fades as the time-dependent main-sequence turnoff mass $M_{\rm TO}(t)$ decreases. Therefore the UV luminosity can be approximated as

with $\alpha$, the IMF slope and $k$ the effective luminosity–mass exponent, $L_\star\propto M^k$. We use as main-sequence turnoff masses the H-burning end ages from Ekström et al. [51].

The exponent $k$ is calibrated requiring the expression above to reproduce the measured fading of the BPASS models for the utilized slope $\alpha= 2.35$ yielding $k = 2.5$. Then we compute the evolution of $M_{\rm UV}$ from Eq. (1), normalizing to the observed value at $t = 0$. The evolution of $M_{\rm UV}$ of this toy model is shown in the bottom panel of Figure 11 as a gold solid line, assuming a top-heavy IMF with $\alpha= 1.85$ (e.g., [216]).

For this assumption, the starburst reaches $M_{\mathrm{UV}} = -20$ in only $\sim2$ Myr, much faster than the 5 Myr of the fiducial SSP assuming a Salpeter-like IMF. Note that hydrodynamical simulations of very dense star clusters predict that, due to converging mass flows, the most massive stars can grow efficiently, and extreme IMFs can be achieved up to log-flat shapes ($\alpha= 1$; Chon et al. [35]). Generally, the fading process to a typical LRD UV magnitude of $\sim-19$ can happen faster than the duty cycle of $\sim10$–$20$ Myr, inferred from the incidence of LRDs amongst $\sim10^{8}\ M_{\odot}$ galaxies (Sun et al. [205]), suggesting that the red engine may outlive the VMS population.

The origin of the high N/C abundance

In Sect. 3.5 we have reported a very high single-zone N/O abundance, only comparable to the most enriched systems at high-z (e.g., Ji et al. 2026; Naidu et al. [150]; Morel et al. [145]). Nitrogen is a bottleneck element in the stellar CNO cycle. Convection and rotation of massive stars can efficiently cause mixing in the stellar envelope, bringing the N-rich, C-poor material of the H-burning core to the surface, allowing it to be ejected through stellar winds (e.g., Maeder & Meynet [117]). For this reason, the N/C abundance ratio has the potential to be a strong indicator of enrichment by stellar yields. Below we explore models of yields by rotating VMSs’ winds, massive Pop III stars (up to $\sim6000\ M_{\odot}$), and SMS nucleosynthesis ($\gtrsim10^{4}\ M_{\odot}$).

VMS winds

We computed the $\log_{10}(\mathrm{N/C})$ abundances of the stellar wind yields of the PARSEC VMS models, which include stellar tracks for masses up to 600 $M_{\odot}$ (Costa et al. [37]; Nguyen et al. [159]). In Figure 12 we show the abundance as a function of an increasing IMF high-mass cutoff. The total yields of the stellar population depends on age, so we plot the yields for our fiducial best-fit population of 2–3 Myr, weighted by total mass, but note that these ages are close to what would yield the maximum N/C achievable by a VMS population (see e.g., Shi et al. [197]). We also consider the ejecta of fast-rotating 0.14$Z_{\odot}$ VMSs ($v/v_{\mathrm{crit}} = 0.4$; Higgins et al. [67]).

In Figure 12 we compare the $\log_{10}(\mathrm{N/C})$ abundances of COLA1 (dense zone, diffuse zone, and integrated) with the yields of VMS populations. Achieving the measured single-zone $\log_{10}(\mathrm{N/C}) = +0.58$ with such stellar populations would require an IMF cutoff of at least $\approx300$ $M_{\odot}$ for non-rotating stars, and $\approx100$ $M_{\odot}$ for the case of $v/v_{\mathrm{crit}} = 0.4$ rotating stars. However, note that the green curves plotted in Figure 12 assume zero dilution with the ISM metals; only net stellar yields are considered. This would correspond to an extreme case where the original intracluster gas has been blown away by a previous generation of stars (see e.g., Marques-Chaves et al. [122]). Therefore, we conclude that it is unlikely that the measured single-zone N/C abundance can be produced by the winds of VMS only.

Turning to the two-zone abundance estimates, the dense-zone $\log_{10}(\mathrm{N/C}) > 1.3$ simply cannot be reached with this models with an IMF cutoff of up to $\gtrsim600$ $M_{\odot}$, even with fast-rotating stars. In order to explain the observed $\log_{10}(\mathrm{N/C})$ in our toy model we would need $\gtrsim50\%$ of the stars rotating at an extreme $v/v_{\mathrm{crit}} = 0.7$ (Meynet et al. [143]), which renders infeasible compared to the observed rotation of individual massive stars in the local Universe (e.g., Harries et al. [64]; Ramírez-Agudelo et al. [176]).

PopIII stars

Alternative mechanisms need to be invoked to explain the extreme $\log_{10}(\mathrm{N/C})$ values in COLA1. Nandal et al. [157] predicted that the ejecta of Pop III stars of up to $\sim6000$ $M_{\odot}$ would be composed mainly of H, He and N. In Figure 12 we show the $\log_{10}(\mathrm{N/C})$ prediction by Nandal et al. [157] as a result of the Pop III ejecta in the embryonic halos of the Regan et al. [179] hydrodynamical simulation. The abundances of the Pop III ejecta of this model broadly match the integrated value of COLA1. The upper and lower limits of the yellow box correspond to the assumptions that the 10% and 100% of the material above the CO core in the Pop III stars is ejected, respectively. Note that the integrated ejecta depends strongly on the assumption of how the IMF is sampled; higher mass stars produce higher N/C in general. As discussed in Sect. 3.3, the best-fit stellar population to the UV continuum needs a metallicity not significantly lower than 0.1 $Z_{\odot}$ to produce the observed metal absorptions. However, it could be possible that the observed N/C abundances were originated from Pop III nucleosynthesis, and the currently observed stellar population is a second generation of stars. In this case, the metals traced by the stellar absorption features and the nebular lines (e.g., [O iii]) must have been produced by an earlier population, most likely through core-collapse supernovae, whose C- and O-rich ejecta would dilute the N/C ratio of the Pop III material. The young age of the observed population ($\sim2.5$ Myr) implies that it has not yet had time to alter the gas through its own AGB stars ($\gtrsim30$ Myr; Karakas [85]) or core-collapse supernovae of the current population (Pipino & Matteucci [172]). In addition, the nitrogen excess appears to be concentrated in the dense zone of our two-zone model (Sect. 3.5; see also Figure 8), which favours a recent, localized enrichment over a pre-existing one that would have mixed throughout the ISM.

Supermassive stars

Supermassive stars (SMS; $\sim10^{4}$–$10^{5}$ $M_{\odot}$) have been proposed as candidates responsible for the elevated N/O enrichment through hot CNO nucleosynthesis in dense stellar clusters and high-z galaxies (Nagele & Umeda [146]; Charbonnel et al. [30]; Marques-Chaves et al. [122]), or the anomalous chemical enrichment in globular clusters (Carretta et al. [25]; Denissenkov & Hartwick 2014). Charbonnel et al. [30] modeled the nucleosynthesis of a fast growing SMS reaching $10^{4}$ $M_{\odot}$ using the stellar evolution code MESA and following the prescriptions of Fragoso et al. (2023). In Figure 12 we show the N/C abundance of the [Fe/H] $= -1.14$ model of Charbonnel et al. [30] without mixing with the intracluster gas as an upper limit. The very high maximum $\log_{10}(\mathrm{N/C}) \approx1.5$ of this model is compatible with the dense-zone abundance of COLA1, albeit with the limiting case of assuming no mixing with the cluster gas or pollution from other metal-enrichment sources. Nagele & Umeda [146] presented a model using the code HOSHI and a Takahashi et al. [206] mass loss prescription. We also show in Figure 12 their $10^{5}$ $M_{\odot}$ (close to the theoretical limit due to general-relativistic instabilities; Chandrasekhar [28]) as an upper limit assuming no mixing. As shown in Figure 12, SMSs of $10^{4}$–$10^{5}$ $M_{\odot}$ could easily explain the N/C enrichment after including some amount of dilution with the original intracluster gas. Additionally, the yields of SMSs could explain the extreme N/C abundance of COLA1 in its dense zone.

Tidal disruption events

When a star passes too close to a massive black hole, its structure can be torn apart by the gravitational field causing a TDE (e.g., Lacy et al. [103]; Evans & Kochanek [52]). The TDEs yield the contents of the stellar interiors as debris, polluting the ISM with the unprocessed CNO-cycle materials, and could also explain the elevated N/C abundance measured (e.g., Kochanek [95]; Cameron et al. [20]). Extreme N/C values have been reported for TDEs in the local universe as high as $\log_{10}(\mathrm{N}/\mathrm{C}) > 0.9$ (Yang et al. [231]). In a massive, very dense star cluster containing massive stars, where a IMBH or SMBH is probably also present (see discussion in Sect. 4.5 below), an elevated number of TDEs is an expected consequence (e.g., Sakurai et al. [184]; Isobe et al. [75]). The intracluster gas could be originally polluted by the winds of the VMS detected in the UV spectrum of COLA1 in first instance, thence TDEs could possibly further boost the N/C abundance. The material disrupted by the TDEs could also be related with to the formation of a dense gaseous envelope around a massive black hole leading to the LRD phenomenon (e.g., Loeb & Ulmer [111]; Metzger [141]; Bellovary [8]), although this interpretation may be difficult to reconcile with the weak observed variability on a population level among LRDs (e.g., Kokubo & Harikane [98]; Liu et al. [110]; Stone et al. [200]).

Broad Niv] emission and dense outflows

COLA1 presents a broad component in multiple emission lines (see Appendix A). In particular, we detect broad Gaussian components in the N IV] $\lambda\lambda1483$, 1486 doublet (FWHM = 1800 $\pm$ 200 km s$^{-1}$, $\Delta v = +200 \pm200$ km s$^{-1}$) and C IV $\lambda1550$ (FWHM = 1040 $\pm$ 50 km s$^{-1}$, $\Delta v = -86 \pm200$ km s$^{-1}$)⁶ in the rest-UV. The C IV component is consistent with the systemic redshift, although very uncertain due to the degeneracy with the P Cygni absorber, and could also have a component due to resonant scattering. The N IV] doublet broad components in turn are redshifted (moving away from the line of sight), although compatible with being at the systemic velocity within 1$\sigma$ (see Fig. A.1). Similar outflows were detected in other systems with hints of VMS and similar spectra to COLA1, such as GNz11 (Chen et al. [31]; Nakane & Ouchi [154]). A key difference of our results with respect to GNz11 is that we report broad emission associated with both components of the N IV] $\lambda\lambda1483$, 1486 doublet, with the $\lambda1483$ component being stronger. This is made evident in Fig. A.1, where one can see a broad tail extending to the blue side of the doublet, not captured by a narrow component. This measurement could in fact be affected by model degeneracies, but at face value it suggests that the N IV] outflow has a density $n_{\mathrm{e}} \lesssim10^{6.5}$ cm$^{-3}$. In the rest-optical, we also detect outflow components in the [O III] $\lambda\lambda4960$, 5008 doublet (FWHM = 431 $\pm$ 7 km s$^{-1}$, $\Delta v = 53 \pm3$ km s$^{-1}$) and [Ne III] $\lambda3860$ (FWHM = 590 $\pm$ 70, $\Delta v = 10 \pm20$). The [Ne III] outflow is again consistent with systemic velocity, but the broad [O III] $\lambda5008$ emission is clearly redshifted. The optical outflows are much narrower than the strong and broad N IV], but the fact that we observe them in several lines at the same time suggests that there is outflowing gas, with density $< 10^{6}$ cm$^{-3}$. The difference of widths across the line components could reflect a structure in the ionization, and chemical enrichment pattern of the outflow (i.e., N/O and N/C; see Sect. 3.5).

Physical interpretation – Formation of a massive seed through runaway stellar mergers

The three defining properties of COLA1—a rest-UV continuum dominated by very massive stars (VMSs), an LRD-like core that dominates the rest-optical emission, and an extreme nitrogen enhancement—may have a causal relation. It has been argued that many LRDs are well described as the combination of a young starburst and a red compact engine which could be emerging from the same physical process (e.g., Inayoshi et al. [74]; Chisholm et al. [32]). As an alternative explanation, the red continuum and broad Balmer lines could be tracing a rapid (super-Eddington) accretion phase in a pre-existing massive black hole with significant dust that obscures our line of sight (e.g. Madau & Maiolino [115]), in which case our data still suggest such phenomena likely take place amidst dense star clusters.

Runaway collisions in very dense star clusters can rapidly form IMBHs, with possibly the appearance of an SMS as intermediate step (e.g., Portegies Zwart et al. [174]; Gieles et al. [56]). In this scenario, the most massive stars sink to the cluster center and merge repeatedly. Hydrodynamical simulations show that an SMS can grow to $\sim10^{4}\ M_{\odot}$ in timescales of $< 1$ Myr through runaway collisions, and can last up to $\sim5$ Myr through stellar-collision rejuvenation (e.g., Gieles et al. [56]; Vergara et al. [217]; Rantala [177]; Williams et al. [226]). In most models, the SMS lives as long as it is continuously accreting material, until it eventually collapses to an IMBH of $M_{\mathrm{BH}} \sim10^{3}$–$10^{4}\ M_{\odot}$ due to relativistic instability (e.g., Woods et al. [227]; Haemmerlé et al. [61]). In this process, a population of VMSs is expected as a natural by-product. The timescales of SMS formation are comparable to the lifetimes of the VMSs, so both could be present within similar timescales.

Our measured (barely resolved) UV effective radius of $\sim200$ pc, and compactness with a stellar surface mass density of $\Sigma_{*} \sim400\ M_{\odot}\ \mathrm{pc}^{-2}$, match individual star cluster complexes (e.g., Bouwens et al. [15]). The stellar density is also in line with hydrodynamical simulations where a handful of very dense stellar clusters can form, containing VMSs and SMSs (Rantala [177]; Chon et al. [35]; Williams et al. [226]). Our data cannot resolve individual clusters, but the compactness of COLA1 leaves room for such a configuration. Consider for example, the highly lensed system ($\mu\approx20$) reported by Yanagisawa et al. [229], where LRDs appear to form in compact off-nuclear clusters. This formation setup offers an explanation for the nitrogen enhancement through VMS stellar winds, TDEs produced by the newly formed IMBH (Sakurai et al. [184]), or the ejecta of an SMS as discussed in detail in Sect. 4.3 above.

Now we consider the possibility that the optical emission of COLA1 could be related to an SMS formed within the environment of a dense stellar cluster. SMSs are expected to have inflated, cool envelopes presenting effective temperatures as low as $\sim4000$ K, and broad permitted lines with P Cygni profiles in the optical spectrum (e.g., Gieles et al. [56]; Nandal & Loeb [156]; Zwick et al. 2026). Nevertheless, some of the current models of SMS atmospheres struggle to reproduce the red colors of LRDs’ optical spectra without invoking heavy dust attenuation (e.g., Martins et al. [126, 127]). In addition, the short lives inferred for SMSs also challenge an interpretation of the red continuum of LRDs as active SMSs directly. The red optical continuum and broad emission lines that we observe in COLA1 (and in LRDs in general) could more probably represent a stage following an SMS, where the core has collapsed due to gravitational instability as an IMBH seed (e.g., Inayoshi et al. [74]; Sun et al. [204]; Curtis et al. [39]; Liu et al. [109]; Nandal et al. [155]). Before the collapse, the SMS may rapidly lose mass through wind pulsations, providing a possible explanation for the assembly of a very dense cocoon of gas giving rise to the observed LRD features (e.g., Naidu et al. [147]). This cocoon formed from the SMS ejecta would be expected to contain N-rich material from the SMS nucleosynthesis (e.g., Nandal et al. [155]), which may explain the high N/O and N/C recently measured in the dense gas ($n_{\mathrm{e}} > 10^{6.5}\ \mathrm{cm}^{-3}$) of some LRDs (e.g., Isobe et al. [75]; Ji et al. 2026; Tang et al. [207]). In this case, accretion into the newly formed massive black hole surrounded by the dense gas is likely to be the powering source of the LRD’s bolometric luminosity. The engine may be sustaining the envelope through radiation pressure, emitting at roughly the Eddington ratio (e.g., Begelman et al. [7]; Begelman & Dexter [6]), and the lifetime of this phase could last for up to $\sim10$–$100\ \mathrm{Myr}$, mostly depending on the mass of the system and the uncertain accretion physics (e.g., Roman-Garza et al. [181]; Santarelli et al. [187]; Wang & Inayoshi [220]).

Comparison to other LRDs and N-emitters

In Figure 13 we compare the main nitrogen and carbon UV lines of COLA1 to the galaxy CEERS-1019 ($z = 8.678$; Larson et al. [104]) and the LRD SPURS-GN-2 ($z = 7.186$; Tang et al. [207]). These two sources have highly similar N iv] profiles, with a prominent $\lambda1486$ component and a broad base, indicative of outflows as discussed in Sect. 4.4. They also have high measured N/C abundances: $\log_{10}(\mathrm{N/C}) = 0.36$, and $\log_{10}(\mathrm{N/C}) \approx0.6$ (Tang et al. [207]; Marques-Chaves et al. [122]). SPURS-GN-2 presents an inverted C iii] doublet ratio with respect to COLA1, which suggest that significant carbon emission is produced in the dense gas ($\log_{10}(n_{\mathrm{e}}/\mathrm{cm}^{-3}) \gtrsim6.5$). In turn, CEERS-1019 shows C iii]$1906/1909 \approx1$; this observed ratio would be compatible with a single-zone density of $\approx10^{4}\ \mathrm{cm}^{-3}$. However in a two-zone scenario as proposed for COLA1, this would be degenerate with a low density zone producing most of the C iii] emission with some contribution from the high density zone (the same that dominates the N iv] emission). A similar case is that of GN-z11, presenting a similar spectrum as CEERS-1019, and line ratio C iii] $\sim1$ (Chen et al. [31]).

![Comparison of N iv], C iv, N iii], and C iii] emission lines](media/fig13.png) bbox=0.12,0.07,0.94,0.26

Another noteworthy difference between COLA1, CEERS-1019, and SPURS-GN-2 is the C iv profile. SPURS-GN-2 shows absorption imprinted on broad C iv, suggesting that the UV is dominated by the central engine rather than a stellar population (Figure 13). COLA1 notably shows enhanced emission near systemic velocities in C iv. This is in line with COLA1 having a very high $f_{\mathrm{esc}}^{\mathrm{LyC}}$, as clear channels in the ISM facilitate the escape of the resonant C iv photons (Naidu et al. [149]; see further discussion in Sect. 4.8).

Ionizing properties of COLA1

The new NIRSpec UV spectrum confirms the hypothesis that COLA1 is a strong producer and leaker of ionizing photons. In the first place, there is a very luminous Ly$\alpha$ line as observed by previous works (Hu et al. 2016; Matthee et al. [136]) with $\mathrm{EW}_{0} = 400^{+40}_{-30}\ \mathring{\mathrm{A}}$. In Fig. D.1 we show the Ly$\alpha$ line profile from JWST/NIRSpec data, where the increased spectral resolution of the 2nd and 3rd order spectrum reveals the double-peak shape. We estimate the Ly$\alpha$ escape fraction within the NIRSpec aperture comparing to the narrow H$\alpha$ flux, assuming Case B recombination ($\mathrm{Ly}\alpha/\mathrm{H}\alpha\approx8.6$) and negligible dust attenuation (Sect. 3.4). This calculation assumes that the Ly$\alpha$ is solely powered by the same ionizing source that powers the narrow component of H$\alpha$, which could be a combination of stars and any ionizing output from the central engine of the LRD. We obtain an escape fraction of $f_{\mathrm{esc}}^{\mathrm{Ly}\alpha} = 0.68 \pm0.03$. As discussed by Torr alba et al. [209], the high Ly$\alpha$ escape fraction, and narrow peak separation is suggestive of a high escape fraction of ionizing photons (Izotov et al. [78]; Maji et al. [119]; Kimm et al. [89]; Naidu et al. [149]). We estimate the escape fraction using the models of Mascia et al. [129], which employ Cox models to infer $f_{\mathrm{esc}}^{\mathrm{LyC}}$ using as inputs: $(E_{\mathrm{B}} - V)$, $M_{\mathrm{UV}}$, stellar mass and \log_{10}(\mathrm{H}\beta+ [\mathrm{O}\,\textsc{iii}]). With this model we obtain $f_{\mathrm{esc}}^{\mathrm{LyC}} = 0.67 \pm0.04$.

We estimate the ionizing photon production efficiency as

using a recombination coefficient $c_{\mathrm{H}\alpha} = 1.128 \times10^{-12}\ \mathrm{erg}^{-1}$ for $n_{\mathrm{H}} = 10^{3}\ \mathrm{cm}^{-3}$ and $T_{\mathrm{eff}} = 2 \times10^{4}\ \mathrm{K}$ from PyNeb (Luridiana et al. [114]). We use only the narrow H$\alpha$ flux in this calculation, implicitly assuming that the UV is only of stellar origin, and that the narrow H$\alpha$ is powered by the same stars. We obtain $\log_{10}(\xi_{\mathrm{ion},0}/\mathrm{Hz}\ \mathrm{erg}^{-1}) = 25.52 \pm0.02$, where $\xi_{\mathrm{ion},0} = \xi_{\mathrm{ion}}\left(1-f_{\mathrm{esc}}^{\mathrm{LyC}}\right)$. If we use the previously quoted number of $f_{\mathrm{esc}}^{\mathrm{LyC}} = 0.67$, the intrinsic ionizing production efficiency estimated is $\log_{10}(\xi_{\mathrm{ion}}/\mathrm{Hz}\ \mathrm{erg}^{-1}) \approx26$. Our fiducial SSP population including VMS (Sect. 3.3) yields an intrinsic $\log_{10}(\xi_{\mathrm{ion}}/\mathrm{Hz}\ \mathrm{erg}^{-1}) = 25.77$, which would imply a somewhat lower escape fraction of $f_{\mathrm{esc}}^{\mathrm{LyC}} \approx0.44$, albeit still high compared to confirmed leakers at low-$z$ (e.g., Flury et al. [54]). This discrepancy could also be explained by some amount of ionizing photons escaping from the LRD engine.

Impact on cosmic reionization

We investigate whether COLA1 is capable of carving the ionized region needed in order to explain the observed double-peaked Ly$\alpha$ profile. We explore the growth of a spherical ionized region powered by a central source as a Strömgren sphere (Strömgren [202]). We assume $\log_{10}(\dot{\xi}_{\mathrm{ion},0}/\mathrm{Hz}\ \mathrm{erg}^{-1}) = 25.52$ (Sect. 4.7) to compute the escaping ionizing photon rate as

defined as the number of photons produced by the source that are able to escape to the IGM and power an ionized bubble. We assume a mean IGM density of $\langle n \rangle= 1.88 \times10^{-7} \cdot(1+z)^3\ \mathrm{cm}^{-3}$ (Mason & Gronke [130]). In Figure 14 we show the evolution of the ionized region for different values of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ and for three fractions of the mean IGM density ($\Delta\equiv n_{\mathrm{H}}^{\mathrm{IGM}}/\langle n_{\mathrm{H}}^{\mathrm{IGM}}\rangle$), and assuming a clumping factor $C = 3$ (e.g., Finlator et al. [53]). We compare these growth curves with the expected ionized region size in the line of sight that is needed to achieve sufficient transmission of the blue Ly$\alpha$ peak. The significant Ly$\alpha$ transmission down to $-250\ \mathrm{km}\ \mathrm{s}^{-1}$ implies the existence of a proximity zone of at least $\sim0.32\ \mathrm{pMpc}$ (Matthee et al. [136]). We compute the size of the proximity zone $R_{\alpha}$ of our ionization model (Bolton & Haehnelt [12]). The size of the proximity zone depends on the IGM density, roughly as $\propto\Delta^{-1}$. As shown previously by Mason & Gronke [130], there needs to be an underdensity of $\Delta\lesssim0.5$ in order for a source like COLA1 to sustain a proximity zone of $> 0.32\ \mathrm{pMpc}$, but the total ionized region must be larger if the effects of a Ly$\alpha$ damping wing by the neutral IGM is taken into account. In particular, Mason & Gronke [130] derived a size of $R_{\alpha} \approx0.7$ Mpc, assuming that the nearest neutral patch has $\Delta= 1$. However, the $R_{\alpha}$ requirement can be relaxed towards smaller ionized regions if the density of the neutral IGM is lower. A complete IGM-ionization study is out of the scope of this work. To show consistency with our results, in Figure 14 we show that COLA1 is able to ionize past the required proximity zone of 0.32 Mpc for an underdense sightline of $\Delta\gtrsim0.5$ in a timescale much shorter than the age of the stellar population ($t \lesssim2$ Myr). However, the requirement of significant ionization

inside the bubble to sustain a Ly$\alpha$ transparent zone implies even stronger underdensities. As shown in Figure 14, an underdensity as low as $\Delta= 0.3$ is needed so COLA1 can sustain a Ly$\alpha$ transparent proximity zone to match the observed blue peak transmission of $\sim250\,\mathrm{km\,s}^{-1}$. Such strong constraints would require a very rare underdense line of sight, for example, placing COLA1 at the edge of a cosmic void as discussed by Park et al. [167].

Interestingly, the galaxies CEERS-1019 and GN-z11, which have very similar UV spectra to COLA1 (see discussion in Sect. 4.6) are also strong Ly$\alpha$ emitters at $z > 7$. Chen et al. [31] have discussed that GN-z11 must reside in an ionized bubble of $\sim0.5$ pMpc, and Whitler et al. [224] found evidence for a $\sim0.5$–1 pMpc ionized bubble at $z = 8.7$ in an overdensity containing CEERS-1019. Moreover, the Ly$\alpha$ observed in some high-$z$ LRDs suggest that some of these sources also trace large ionized regions (Tang et al. [208]). All these results suggest a connection between the LRD and N-enhanced galaxy phenomenon and the presence of large ionized IGM regions. It is possible that a previous burst of ionizing photons powered the proximity zone before turning off (Meyer et al. [142]), which should have happened in the last $\sim0.1$–$0.7$ Myr (for $\Delta= 0.5$ to 0.2), before the IGM naturally recombines (Park et al. [167]). This could suggest that there was a much more luminous episode of LRD or AGN activity in the past, or that there was recent even stronger burst of star formation. The latter case could be represented by the extremely luminous compact galaxy recently reported by Yang et al. [232], whose UV spectrum is dominated by the light of VMSs and presents very weak nebular emission lines pointing towards a very high $f_{\mathrm{esc}}^{\mathrm{LyC}}$ (see also Marques-Chaves et al. [120]; Wu et al. 2026). Given the high number densities of LRDs at $z > 6$ (e.g., Kocevski et al. 2025; Zhuang et al. [235]), if an early COLA1-like phase turns out to be common around most LRDs, these objects have the potential to have played a critical role in local reionization through short starburst episodes with high ionizing output powering Mpc-scale bubbles (e.g., Naidu et al. [152]; Matthee et al. [134]).

Summary

The spectrum of the galaxy COLA1 presents a striking combination: a UV spectrum dominated by very massive stars, a red LRD-like component in the optical, and extremely high nitrogen enhancement. Below we summarize our main findings.

1. We fit the UV spectrum of COLA1 with stellar population synthesis models, favoring those containing VMS (limited to masses up to 300 $M_\odot$). We find a best-fit stellar population with mean age $2.43 \pm0.75$ Myr, and stellar mass $\approx10^8 M_\odot$ (for an IMF with slope $\alpha= 2.35$). In particular, the VMSs are needed to reproduce the P Cygni-like profiles of N v $\lambda\lambda1239$, 1243, Si iv] $\lambda\lambda1394$, 1403, and C iv $\lambda\lambda1548$, 1551, and the broad component of He ii $\lambda1640$ (Sect. 3.3).

2. While the rest-UV spectrum has a fairly blue slope ($\beta_{\rm UV} = -2.5$), the rest-optical is significantly redder ($\beta_{\rm opt} = -1.1$), shows strong He i, and broad H$\beta$ and H$\alpha$ with a narrow Balmer absorption. The narrow components of the Balmer lines follow Case B recombination ratios with little-to-no dust attenuation, while the broad components present high decrements (H$\alpha$/H$\beta\sim9$) that could be explained by attenuation by nuclear dust, but also by RT effects in a dense wind akin an LRD envelope (Sect. 3.4).

3. We measured the nitrogen and carbon abundance of COLA1 using a single-zone density model against line emissivities from PyNeb, obtaining a supersolar $\log_{10}(\mathrm{N/O}) = -0.07 \pm0.05$, and sub-solar $\log_{10}(\mathrm{C/O}) = -0.64 \pm0.10$. The high value of $\log_{10}(\mathrm{N/C}) = +0.58 \pm0.10$ requires exotic enrichment channels, such as VMS winds, TDEs, or the yields of SMSs. Moreover, decomposing the spectrum with a two-zone model using a Cloudy grid, we find evidence for the presence of a dense zone ($n_{\mathrm{e}} \approx10^{7.5}\,\mathrm{cm}^{-3}$) with extreme $\log_{10}(\mathrm{N/C}) > 1.3$ ($2\sigma$), and a more diffuse zone with $\log_{10}(\mathrm{N/C}) < +0.02$. Such abundances can only be explained by the most extreme enrichment channels such as $>10^4 M_\odot$ SMSs (Sects. 3.6, 3.5, and 4.3).

4. We argue that the UV-to-optical spectrum of COLA1 can be well described by the combination of a dense stellar cluster of young massive stars, including VMSs ($>100 M_\odot$), dominating the UV light; and a red LRD component rising towards the rest-optical wavelength range. The lives of the VMSs that dominate the UV are predicted to be short ($\sim1$–2 Myr), therefore the UV of this system will fade by 1–2 mag in timescales of $\sim5$ Myr (depending on the IMF), then COLA1 would be observed as a typical V-shaped LRD (Sects. 4.1 and 4.2).

5. We discuss possible mechanisms leading to the single- and two-zone N/C abundances observed in COLA1. The winds of rapidly rotating VMS could account for the measured single-zone abundance only under very optimistic conditions. Producing such high $\log_{10}(\mathrm{N/C}) \approx0.6$ with stellar yields would require models of $\gtrsim1000 M_\odot$ Pop III stars, or SMSs with $10^4$–$10^5 M_\odot$. In particular, the inferred dense zone abundance of $\log_{10}(\mathrm{N/C}) > 1.3$ ($2\sigma$) requires enrichment from the yields of NSs, or TDEs of massive stars formed in a dense star cluster (Sects. 3.5).

6. We speculate that runaway collisions in a dense star cluster could have planted a massive black hole seed, that manifests in the observed LRD-like spectrum through dense-gas enshrouded accretion into a newly formed IMBH; ultimately leading to the early stages of the formation of a SMBH. The same process may have formed the VMSs observed in the UV as a by-product (Sect. 4.5).

7. COLA1 also shows extraordinary ionizing properties, a high inferred escape fraction ($f_{\mathrm{esc}}^{\mathrm{LyC}} \approx0.67$), and lives inside a $>0.3$ Mpc highly ionized bubble. We argue that if COLA1 represents a common phase in the LRD evolution, given the elevated number densities at $z > 6$, this population may have a crucial impact on cosmic reionization. The light of the observed VMSs, but possibly also past intense starbursts may have been responsible for powering COLA1’s ionized bubble (Sects. 4.6, 4.7, and 4.8).

JWST is revealing the remarkable, unprecedented properties of early stellar populations. Here, we show that these stellar populations are intimately linked to the perplexing BH populations at high-$z$ such as LRDs, and may be critical to unraveling the seeding and birth of present-day massive BHs. Recent studies are revealing that populations of VMS are common in the most luminous high-$z$ galaxies (Upadhyaya et al. [215]; Marques-Chaves et al. [121]; Nakane & Ouchi [154]; Chen et al. [31]; Yang et al. [232]). With COLA1, we may be witnessing how the same dense environment harbouring a population of VMSs is the melting pot in which a supermassive star can form and collapse, seeding the massive black holes of the nearby Universe.

Acknowledgments

AT thanks Debasish Dutta, Ruggero Valli, Sunmyon Chon, and Zoltán Haiman for useful conversations about stellar winds of massive stars, stellar evolution, supermassive star formation and massive black hole seeding. We acknowledge funding by the European Union (ERC, AGENTS, 101076224). This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. 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 programs #7722, #3767, and #1933. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.00741.S ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Some of the data products presented herein were retrieved from the Dawn JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center (DAWN), which is funded by the Danish National Research Foundation under grant DNRF140. Generative AI (Anthropic Claude Opus 5.5) assisted with coding and literature searching. The authors reviewed and approved all scientific content and take full responsibility for the manuscript. The text was entirely written, edited, and reviewed by humans.

References

Line	$\lambda_{\mathrm{rest}}$ [Å]	Flux	EW$_0$	FWHM
Ly$\alpha$	1215.67	$706^{+5}_{-6}$	$400^{+40}_{-30}$	$280^{+30}_{-30}$
N v	1238.82	$7^{+10}_{-4}$	$2^{+3}_{-1}$	$660^{+150}_{-120}$
N v (abs)	1238.82	$-11^{+3}_{-5}$	$-3.2^{+0.7}_{-1.4}$	$580^{+80}_{-160}$
N v	1242.80	$22^{+4}_{-3}$	$5.8^{+1.1}_{-1.0}$	$660^{+150}_{-120}$
N v (abs)	1242.80	$-15^{+4}_{-6}$	$-4.2^{+1.1}_{-1.6}$	$600^{+80}_{-120}$
Si iv	1393.76	$6^{+3}_{-1}$	$1.5^{+0.7}_{-0.3}$	$350^{+130}_{-130}$
Si iv (abs)	1393.76	$-7.9^{+1.7}_{-1.9}$	$-1.9^{+0.4}_{-0.4}$	$640^{+50}_{-90}$
Si iv	1402.77	$5^{+2}_{-2}$	$1.3^{+0.6}_{-0.4}$	$350^{+130}_{-130}$
Si iv (abs)	1402.77	$-3.6^{+1.3}_{-2.0}$	$-0.9^{+0.3}_{-0.4}$	$570^{+100}_{-180}$
[O iv]	1405.46	$4.2^{+1.9}_{-1.6}$	$1.0^{+0.4}_{-0.4}$	$700^{+200}_{-300}$
[N iv]	1483.32	$1.4^{+1.3}_{-1.0}$	$0.4^{+0.3}_{-0.3}$	$360^{+30}_{-30}$
[N iv] (broad)	1483.32	$24^{+16}_{-16}$	$6^{+4}_{-4}$	$1800^{+200}_{-200}$
N iv]	1486.50	$48^{+3}_{-4}$	$12.6^{+0.9}_{-0.9}$	$360^{+30}_{-30}$
N iv] (broad)	1486.50	$20^{+16}_{-16}$	$5^{+4}_{-4}$	$1800^{+200}_{-200}$
C iv	1548.19	$0.6^{+1.0}_{-0.5}$	$0.2^{+0.3}_{-0.1}$	$250^{+40}_{-90}$
C iv (abs)	1548.19	$-23^{+6}_{-9}$	$-7^{+2}_{-3}$	$680^{+10}_{-20}$
C iv (broad)	1548.19	$70^{+60}_{-70}$	$20^{+20}_{-20}$	$1040^{+50}_{-30}$
C iv	1550.77	$1.1^{+1.6}_{-0.8}$	$0.3^{+0.5}_{-0.2}$	$250^{+40}_{-90}$
C iv (abs)	1550.77	$-64^{+16}_{-20}$	$-20^{+5}_{-6}$	$680^{+10}_{-20}$
C iv (broad)	1550.77	$60^{+50}_{-50}$	$17^{+16}_{-15}$	$1040^{+50}_{-30}$
He ii	1640.42	$3.7^{+1.2}_{-1.2}$	$1.3^{+0.4}_{-0.4}$	$244^{+22}_{-24}$
He ii (broad)	1640.42	$7.4^{+1.5}_{-1.7}$	$2.5^{+0.5}_{-0.6}$	$970^{+280}_{-200}$
O iii]	1660.81	$4.8^{+0.7}_{-0.6}$	$1.7^{+0.2}_{-0.2}$	$170^{+30}_{-40}$
O iii]	1666.15	$12.7^{+0.7}_{-0.7}$	$4.5^{+0.3}_{-0.3}$	$170^{+30}_{-40}$
N iii]	1747.85	$0.8^{+1.1}_{-0.6}$	$0.3^{+0.5}_{-0.3}$	$650^{+80}_{-70}$
N iii]	1748.65	$3.3^{+1.5}_{-1.5}$	$1.4^{+0.6}_{-0.6}$	$650^{+80}_{-70}$
N iii]	1751.22	$9^{+2}_{-3}$	$3.6^{+0.8}_{-1.1}$	$650^{+80}_{-70}$
N iii]	1751.66	$2^{+2}_{-1}$	$0.7^{+1.1}_{-0.5}$	$650^{+80}_{-70}$
Si iii]	1892.03	$2.2^{+0.7}_{-0.8}$	$1.1^{+0.4}_{-0.4}$	$200^{+200}_{-100}$
[C iii]	1906.68	$16.1^{+0.9}_{-1.1}$	$7.9^{+0.5}_{-0.6}$	$310^{+40}_{-40}$
C iii]	1908.73	$9.4^{+1.0}_{-0.9}$	$4.6^{+0.5}_{-0.5}$	$310^{+40}_{-40}$
[Ne iv]	2424.26	$<1.7$	$<1.4$	$133^{+4}_{-4}$
Mg ii]	2796.35	$3.0^{+0.4}_{-0.4}$	$2.8^{+0.4}_{-0.3}$	$132^{+4}_{-4}$
Mg ii]	2803.53	$2.0^{+0.3}_{-0.3}$	$1.8^{+0.4}_{-0.2}$	$132^{+4}_{-4}$
[Fe iv]	2830.19	$1.2^{+0.3}_{-0.4}$	$1.0^{+0.4}_{-0.3}$	$132^{+4}_{-4}$
[Fe iv]	2836.57	$2.6^{+0.4}_{-0.4}$	$2.4^{+0.5}_{-0.3}$	$132^{+4}_{-4}$
[Ar iv]	2854.51	$0.6^{+0.3}_{-0.3}$	$0.4^{+0.3}_{-0.2}$	$132^{+4}_{-4}$
[Ar iv]	2869.05	$0.13^{+0.16}_{-0.10}$	$0.04^{+0.11}_{-0.03}$	$132^{+4}_{-4}$
He i	3188.67	$1.6^{+0.3}_{-0.3}$	$1.8^{+0.4}_{-0.2}$	$132^{+4}_{-4}$
[Ne v]	3426.86	$<0.8$	$<1.2$	$133^{+4}_{-4}$

Notes. Fluxes are in units of $10^{-19}$ erg s$^{-1}$ cm$^{-2}$, rest-frame EW in Å, and FWHM in km s$^{-1}$. Upper limits are given for a $3\sigma$ confidence level.

Line	$\lambda_{\mathrm{rest}}$ [Å]	Flux	EW$_0$	FWHM
[O ii]	3727.09	$3.2^{+0.4}_{-0.4}$	$5.7^{+0.8}_{-0.5}$	$132^{+4}_{-4}$
[O ii]	3729.88	$2.0^{+0.4}_{-0.4}$	$3.4^{+0.7}_{-0.5}$	$132^{+4}_{-4}$
H(10-2)	3798.98	$1.9^{+0.3}_{-0.3}$	$3.1^{+0.8}_{-0.5}$	$122^{+6}_{-6}$
H(9-2)	3836.48	$2.5^{+0.3}_{-0.4}$	$4.4^{+0.7}_{-0.7}$	$122^{+6}_{-6}$
[Ne iii]	3869.86	$19.1^{+1.1}_{-1.0}$	$38^{+2}_{-2}$	$132^{+4}_{-4}$
[Ne iii] (broad)	3869.86	$12.3^{+1.1}_{-1.1}$	$23^{+2}_{-2}$	$590^{+70}_{-60}$
He i	3889.75	$1.8^{+1.7}_{-1.3}$	$1^{+1}_{-1}$	$132^{+4}_{-4}$
H(8-2)	3890.17	$5.4^{+1.4}_{-1.7}$	$8^{+4}_{-2}$	$122^{+6}_{-6}$
[Ne iii]	3968.59	$7.5^{+0.4}_{-0.4}$	$15.1^{+0.9}_{-0.6}$	$132^{+4}_{-4}$
He $\delta$	3971.20	$7.0^{+0.4}_{-0.4}$	$14.0^{+0.9}_{-0.8}$	$122^{+6}_{-6}$
H$\delta$	4102.90	$12.0^{+0.4}_{-0.4}$	$25.8^{+1.1}_{-1.1}$	$122^{+6}_{-6}$
H$\gamma$	4341.68	$21.4^{+1.6}_{-1.6}$	$47^{+6}_{-4}$	$122^{+6}_{-6}$
[O iii]	4363.21	$13.2^{+1.6}_{-1.8}$	$28^{+5}_{-3}$	$121^{+3}_{-4}$
He ii	4687.02	$2.5^{+0.4}_{-0.3}$	$4.9^{+0.5}_{-0.5}$	$132^{+4}_{-4}$
[Ar iv]	4712.58	$1.4^{+0.3}_{-0.3}$	$2.5^{+0.8}_{-0.5}$	$132^{+4}_{-4}$
[Ar iv]	4741.45	$1.3^{+0.3}_{-0.3}$	$2.3^{+0.8}_{-0.5}$	$132^{+4}_{-4}$
H$\beta$	4862.69	$36.8^{+1.0}_{-1.0}$	$83^{+3}_{-2}$	$122^{+6}_{-6}$
H$\beta$*	4862.69	$9^{+4}_{-4}$	$13^{+9}_{-5}$	$2900^{+700}_{-110}$
H$\beta$*	4862.69	$8.1^{+1.2}_{-1.1}$	$16^{+3}_{-3}$	$330^{+20}_{-20}$
H$\beta$*	4862.69	$15^{+3}_{-3}$	$29^{+7}_{-7}$	$1460^{+100}_{-110}$
[O iii]	4960.30	$71.7^{+1.4}_{-1.1}$	$163^{+5}_{-3}$	$132^{+4}_{-4}$
[O iii] (broad)	4960.30	$40.5^{+1.1}_{-1.3}$	$91^{+4}_{-7}$	$431^{+7}_{-7}$
[O iii]	5008.24	$214^{+4}_{-4}$	$484^{+17}_{-9}$	$132^{+4}_{-4}$
[O iii] (broad)	5008.24	$120.7^{+3.3}_{-4.4}$	$270^{+12}_{-12}$	$431^{+7}_{-7}$
He i	5877.25	$8.0^{+0.6}_{-0.5}$	$24^{+2}_{-1}$	$132^{+4}_{-4}$
[N ii]	6549.86	$1.9^{+1.0}_{-1.0}$	$2^{+2}_{-1}$	$132^{+4}_{-4}$
H$\alpha$	6564.64	$121^{+5}_{-5}$	$286^{+12}_{-12}$	$122^{+6}_{-6}$
H$\alpha$*	6564.64	$64^{+6}_{-5}$	$146^{+16}_{-9}$	$330^{+20}_{-20}$
H$\alpha$*	6564.64	$92^{+16}_{-14}$	$190^{+40}_{-40}$	$1460^{+100}_{-110}$
H$\alpha$*	6564.64	$195^{+17}_{-18}$	$430^{+60}_{-40}$	$2900^{+700}_{-110}$
[N ii]	6585.28	$0.6^{+0.6}_{-0.4}$	$0.5^{+1.0}_{-0.3}$	$132^{+4}_{-4}$
He i	6680.00	$2.5^{+0.8}_{-0.7}$	$5^{+2}_{-1}$	$132^{+4}_{-4}$

Notes. Fluxes are in units of $10^{-19}$ erg s$^{-1}$ cm$^{-2}$, rest-frame EW in Å, and FWHM in km s$^{-1}$. These components are modeled as broad Gaussians to capture the general shape of the lines, a dedicated modeling of H$\alpha$ and H$\beta$ in particular is detailed in Sect. 3.4.

Emission line additional flux measurements

In Figs. A.1 and A.2 we show the UV and optical emission line fits, respectively. The models and procedure are described in Sect. 3.1. We report the flux, rest-frame EW and Gaussian FWHM of all the fitted components in Tables A.1 and A.2.

H alternative fits

In Fig. B.1 we show the fit to the H$\alpha$ emission line profile of COLA1. The model is described in Sect. 3.4, but see also Tornabene et al. (2026a) and Matthee et al. [137].

Two-zone abundance fits

Fig. C.1 shows the two-zone decomposition using Cloudy models to the UV spectrum of COLA1. This analysis is detailed in Sect. 3.6.

NIRSpec G140M Ly$\alpha$ spectrum

In Fig. D.1 we show the NIRSpec Ly$\alpha$ profile of COLA1. Due to the redshift and the choice of the grating/filter combination G140M/F070LP, we detect the Ly$\alpha$ line in the first-, second- and third-order spectral traces. The line profile is totally unresolved

in the first order, but has increased resolution in subsequent orders as shown in Fig. D.1. We compare with the VLT/X-Shooter $R \approx4000$ spectrum, that clearly shows the double-peak Ly$\alpha$ ([136, 209]).