# The UV-brightest Lyman continuum emitting star-forming galaxy

R Marques-Chaves, D Schaerer, J Álvarez-Márquez, L Colina, M Dessauges-Zavadsky, I Pérez-Fournon, A Saldana-Lopez, A Verhamme

► To cite this version:

R Marques-Chaves, D Schaerer, J Álvarez-Márquez, L Colina, M Dessauges-Zavadsky, et al.. The UV-brightest Lyman continuum emitting star-forming galaxy. Monthly Notices of the Royal Astronomical Society, 2021, 507 (1), pp.524-538. 10.1093/mnras/stab2187. hal-04982534

HAL Id: hal-04982534

https://hal.science/hal-04982534v1

Submitted on 21 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The UV-brightest Lyman continuum emitting star-forming galaxy

R. Marques-Chaves<sup>1,1*</sup> D. Schaerer,<sup>1,2</sup> J. Álvarez-Márquez,<sup>3</sup> L. Colina,<sup>3,4</sup> M. Dessauges-Zavadsky,<sup>1</sup>
I. Pérez-Fournon,<sup>5,6</sup> A. Saldana-Lopez<sup>1</sup> and A. Verhamme<sup>1</sup>

Accepted 2021 July 26. Received 2021 July 26; in original form 2021 July 13

## ABSTRACT

We report the discovery of J0121+0025, an extremely luminous and young star-forming galaxy ($M_{UV} = -24.11$, $\log[L_{Ly\alpha}/\text{erg s}^{-1}] = 43.8)$ at $z = 3.244$ showing copious Lyman continuum (LyC) leakage ($f_{\text{esc,abs}} \approx 40$ per cent). High signal-to-noise ratio rest-frame UV spectroscopy with the *Gran Telescopio Canarias* reveals a high significance ($7.9\sigma$) emission below the Lyman limit ($<912 \AA$), with a flux density level $f_{900} = 0.78 \pm 0.10\mu\text{Jy}$, and strong P-Cygni in wind lines of O VI 1033 Å, N V 1240 Å, and C IV 1550 Å that are indicative of a young age of the starburst ($<10$ Myr). The spectrum is rich in stellar photospheric features, for which a significant contribution of an AGN at these wavelengths is ruled out. Low-ionization interstellar medium (ISM) absorption lines are also detected, but are weak ($EW_0 \simeq 1 \AA$) and show large residual intensities, suggesting a clumpy geometry of the gas with a non-unity covering fraction or a highly ionized ISM. The contribution of a foreground and AGN contamination to the LyC signal is unlikely. Deep optical to Spitzer/IRAC 4.5 $\mu$m imaging show that the spectral energy distribution of J0121+0025 is dominated by the emission of the young starburst, with $\log(M_{\star}^{\text{burst}}/M_{\odot}) = 9.9 \pm 0.1$ and SFR = $981 \pm 232 M_{\odot} \text{yr}^{-1}$. J0121+0025 is the most powerful LyC emitter known among the star-forming galaxy population. The discovery of such luminous and young starburst leaking LyC radiation suggests that a significant fraction of LyC photons can escape in sources with a wide range of UV luminosities and are not restricted to the faintest ones as previously thought. These findings might shed further light on the role of luminous starbursts to the cosmic reionization.

**Key words:** galaxies: evolution – galaxies: formation – galaxies: high-redshift.

## 1 INTRODUCTION

Lyman-$\alpha$ emitters (LAEs) and Lyman break galaxies (LBGs) are widely studied populations of star-forming galaxies that are common in the early Universe. Deep field surveys have been used to study these star-forming galaxies, which are typically faint with magnitudes $R \sim 25$ AB at $z \sim 3$, corresponding to absolute UV magnitudes $M_{\star,\text{UV}}^*$ of about $-20$ to $-21$ at those redshifts (for LAEs and LBGs, respectively; e.g. Ouchi et al. 2008; Reddy & Steidel 2009). These studies have revealed that typical ($M_{\star,\text{UV}}^*$) LAEs and LBGs present a wide range of properties, with stellar masses $\log(M_{\star}/M_{\odot}) \sim 8.0-10.0$ (e.g. Shapley et al. 2001; Gawiser et al. 2007; Ono et al. 2010; Santos et al. 2020), star-formation rates (SFR) up to a few tens $M_{\odot} \text{yr}^{-1}$ (e.g. Shapley et al. 2003; Nakajima et al. 2012; Sobral et al. 2018), and low metallicity on average (e.g. Finkelstein et al. 2011; Nakajima et al. 2013; Kojima et al. 2017).

These star-forming galaxies are likely the dominant sources responsible for ionizing the intergalactic medium (IGM) in the early Universe, during the so-called epoch of reionization (EoR at $6 < z < 15$;

e.g. Robertson et al. 2015), due to their expected large hydrogen ionizing photon (hereafter Lyman continuum; LyC with $>13.6$ eV) escape fraction, $f_{\text{esc}}(\text{LyC})$. However, it is still unclear whether the faint and numerous, or the more luminous and rare are the main contributors (cf., Finkelstein et al. 2019; Naidu et al. 2020).

Recent progress has been made in selecting LyC emitters and understanding their properties. To date $\approx 50$ LyC emitting star-forming galaxies have been discovered and studied in detailed at low-$z$ ($z \lesssim 0.4$; e.g. Borthakur et al. 2014; Izotov et al. 2016, 2018a, b; Leitherer et al. 2016) intermediate-$z$ ($z \sim 1.4$, Saha et al. 2020), and moderately high-$z$ ($z > 2$, de Barros et al. 2016; Shapley et al. 2016; Vanzella et al. 2016; Bian et al. 2017; Steidel et al. 2018; Fletcher et al. 2019; Rivera-Thorsen et al. 2019; Ji et al. 2020) up to $z \sim 4$ (Vanzella et al. 2018), where the IGM still allows the direct observation of LyC photons. However, the relation between $f_{\text{esc}}(\text{LyC})$ and $M_{\text{UV}}$ is still a matter of debate. Some statistical studies targeting star-forming galaxies with $M_{\text{UV}}$ between $\simeq -20$ and $-21.5$ find a clear trend, where faint galaxies show larger $f_{\text{esc}}(\text{LyC})$ than luminous ones (Steidel et al. 2018; Pahl et al. 2021). However, these findings are in tension with those from other works. For example, Bian & Fan (2020) find a low $f_{\text{esc}}(\text{LyC}) < 14$ per cent ($3\sigma$) in faint LAEs ($M_{\text{UV}} \simeq -18.8$), suggesting that the LyC leakage in star-forming

<footnote n="1"><sup>1</sup>Geneva Observatory, University of Geneva, Chemin Pegasi 51, CH-1290 Versoy, Switzerland</footnote>

<footnote n="2"><sup>2</sup>CNRS, IRAP, 14 Avenue E. Belin, F-31400 Toulouse, France</footnote>

<footnote n="3"><sup>3</sup>Centro de Astrobiología (CSIC-INTA), Carretera de Ajalvir, E-28850 Torrejón de Ardoz, Madrid, Spain</footnote>

<footnote n="4"><sup>4</sup>International Associate, Cosmic Dawn Center (DAWN)</footnote>

<footnote n="5"><sup>5</sup>Instituto de Astrofísica de Canarias, C/Ía Láctea, s/n, E-38205 San Cristóbal de La Laguna, Tenerife, Spain</footnote>

<footnote n="6"><sup>6</sup>Departamento de Astrofísica, Universidad de La Laguna, E-38206 San Cristóbal de La Laguna, Tenerife, Spain</footnote>

<footnote>* E-mail: ruijmcoelho@gmail.com</footnote>
galaxies do not follow such trend. In addition, one of the strongest LyC emitters known at high redshift, with $f_{esc}(\text{LyC}) \approx 60$ per cent is also a very luminous ($M_{\text{UV}} = -22.20$) star-forming galaxy (Ion3 at $z = 4.0$; Vanzella et al. 2018).

Understanding the role of UV-luminous star-forming galaxies to the cosmic reionization requires the knowledge of two key properties: the volume density of these sources and their intrinsic properties, those governing the production (e.g. stellar population, star-formation histories), and escape of LyC photons (i.e. neutral gas and dust content and the geometry of the interstellar medium; ISM). UV-luminous star-forming galaxies are rare (e.g. Sobral et al. 2018), yet how much is still unknown, in particular at $z > 6$. Recent studies have found remarkably bright galaxies well into the EoR ($z > 7$) that are in excess compared to the generally observed Schechter component of the UV luminosity function (e.g. Bowler et al. 2014, 2015; Ono et al. 2018; Stefanon et al. 2019). For example, the most distant spectroscopically confirmed galaxy known, GN-z11 at $z \simeq 11.0$ (Oesch et al. 2016; Jiang et al. 2021), is also very luminous in the UV ($M_{\text{UV}} = -22.1$), and the inferred volume density of this source is a factor of $>15$ higher than predicted by models (Oesch et al. 2016).

On the other hand, the number of UV-luminous star-forming galaxies known is scarce, in particular those as bright as $M_{\text{UV}} < -23$ (e.g. Dessauges-Zavadsky et al. 2010; Lee et al. 2013; Harikane et al. 2020; Marques-Chaves et al. 2020a), preventing us to study their physical properties in a statistical basis. While bright sources are in principle easy and optimal targets for deep follow-up observations and detailed analysis of their properties, identifying them is extremely challenging, as it requires large area surveys and, in addition, spectroscopic follow-up to distinguish them from the much more abundant active galactic nuclei (AGN), which is time consuming.

Motivated by this, we have undertaken a search for UV-luminous star-forming galaxies at $z > 2$ in of one of the widest spectroscopic surveys ever performed, the ~9300 deg²-wide-extended Baryon Oscillation Spectroscopic Survey (eBOSS: Abolfathi et al. 2018) of the Sloan Digital Sky Survey (SDSS: Eisenstein et al. 2011). While eBOSS was specifically designed to select bright quasi stellar objects (QSOs) at these redshifts, not necessarily high-z star-forming galaxies, some QSO candidates were selected and observed based on their optical colours only (for details see Ross et al. 2012), mimicking the Lyman break technique. The combination of a wide area surveyed, the limiting *r*-band magnitude of ~22 AB, which corresponds to $M_{\text{UV}} \lesssim -23$ at $z \simeq 2.5$, and the available optical spectra for every source, makes eBOSS an alternative and efficient survey to search and explore such luminous star-forming galaxies.

In Marques-Chaves et al. (2020b) we presented the first result of this project, where an extremely luminous star-forming galaxy was analysed, BOSS-EVULG1 at $z = 2.47$. First selected as a SDSS QSO, follow-up observations revealed that the large luminosities in the UV ($M_{\text{UV}} = -24.4$) and nebular emission ($\log(L[\text{Ly}\alpha, \text{erg s}^{-1}]) = 44.0$) in BOSS-EVULG1 are powered by a vigorous starburst with SFR $\simeq 1000 M_{\odot} \text{yr}^{-1}$, a young age ($\simeq 4$ Myr) and very high sSFR (~100 Gyr$^{-1}$), with no evidence for a dominant contribution of a type-I/II AGN. Interestingly, BOSS-EVULG1 shows properties that are common in LyC leakers, such as low dust content ($E(B-V) \simeq 0.07$ and $\log(L_{IR}/L_{UV}) < -1.2$) and a highly ionized ISM ([O III] 5007 Å/[O II] 3728 Å $\simeq 18$). In addition, it presents very weak low-ionization ISM absorption, suggestive of low geometric covering fraction of the neutral gas that could be related with the powerful starburst-driven ionized outflow detected in $H\alpha$ in this galaxy (Álvarez-Márquez et al. 2021). However, no direct information on the LyC leakage is known for this luminous source.

In this work, we present the discovery and detailed analysis of a new luminous source, SDSS J012156.09+002520.3 ($\alpha, \delta$ [J2000] = 20.4837°, 0.4223°), hereafter J0121+0025, an extremely UV-luminous ($M_{\text{UV}} = -24.1$) star-forming galaxy at $z = 3.244$ showing copious LyC leakage ($f_{esc,abs} \approx 40$ per cent). The paper is structured as follows. The discovery and follow-up observations are presented in Section 2. The analysis of the rest-frame UV spectroscopic observations and imaging data are presented in Section 3. In Section 4 we discuss the LyC properties of J0121+0025 and compare them with those from other LyC emitters, and, finally, in Section 5 we present the summary of our main findings. Throughout this work, we assume a Salpeter (1955) initial mass function (IMF) and a concordance cosmology with $\Omega_m = 0.274$, $\Omega_\Lambda = 0.726$, and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. All magnitudes are given in the AB system.

## 2 DISCOVERY AND FOLLOW-UP OBSERVATIONS

J0121+0025 was discovered as part of our search for Extremely UV-Luminous Galaxies (EUVLGs; $M_{\text{UV}} < -23$) within the eBOSS survey (Abolfathi et al. 2018) of the SDSS (Eisenstein et al. 2011). J0121+0025 is part of a large sample of ~70 very luminous star-forming galaxies that were previously classified as QSOs in the Data Release 14 Quasar catalogue (Pâris et al. 2018). This sample also includes BOSS-EUVLG1 at $z = 2.469$ recently analysed in Marques-Chaves et al. (2020b) and in Álvarez-Márquez et al. (2021). The sample and the selection techniques will be presented in a separated work (R. Marques-Chaves in preparation). Briefly, the selection consists on searching for narrow $Ly\alpha$ profiles in optical SDSS spectra of $z \gtrsim 2$ sources, and blue/flat optical to mid-infrared (mid-IR) colours, i.e. properties that are not expected to be present in QSOs.

The BOSS spectrum of J0121+0025 (plate-mjd-fiberid: 4228-55484-818) shows features characteristic of an unobscured luminous star-forming galaxy, rather than an AGN. In particular, the shallow BOSS spectrum shows a narrow $Ly\alpha$ profile ($\simeq 350 \text{ km s}^{-1}$ full width half maximum; FWHM) and evidence of wind profiles in $N \text{V}$ 1240 Å and C IV 1550 Å in the form of P-Cygni that could be indicative of a young starburst. In addition, the BOSS spectrum shows emission below $\lambda_{\text{obs}} < 3870$ Å, that although detected with low significance ($\simeq 2.3\sigma$), could be an indication of LyC leakage ($\lambda_0 < 912$ Å). Given this, J0121+0025 was selected as a priority target for deep follow-up spectroscopy.

### 2.1 GTC/OSIRIS rest-frame UV spectroscopy

Optical spectroscopy was obtained with the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy instrument (OSIRIS)<fnref n="1" /> on the 10.4 m Gran Telescopio Canarias (GTC). The data were obtained in service mode over two nights, 2020 August 18 and 19 in dark Moon and sub-arcsec ($\simeq 0.8-0.9$ arcsec) conditions as part of the GTC programme GTC21-20A (PI: R. Marques-Chaves). The R1000B grism was used, with dispersion of 2.12 Å, providing a full spectral coverage of 3600–7500 Å, which corresponds to 860–1800 Å in the rest-frame at $z \simeq 3.25$. The OSIRIS 1.2 arcsec-wide slit was centered on J0121+0025 and oriented with the parallactic angle (see Fig. 1). Given this configuration, the corresponding instrumental resolution is $R \sim 800$ or $\simeq 400 \text{ km s}^{-1}$. In total 12 exposures of 900 s were acquired.

<footnote n="1">¹http://www.gtc.iac.es/instruments/osiris/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Cutout of J0121+0025 from the Subaru/HSC R-band image (left). The orientation of the GTC/OSIRIS long slit is marked in blue. The slit also includes a low-z star-forming galaxy with a photometric redshift of 1.10 ± 0.27 and located ≃6.6 arcsec north from J0121+0025. The spectroscopic 1 arcsec-radius BOSS fiber is also marked in red. J0121+0025 shows a compact morphology (≈0.55 arcsec FWHM) and it is only barely resolved in this image [with a point spread function (PSF) of ≃0.50 arcsec FWHM]. Middle and right-hand panels show the residuals from the Galfit modeling of J0121+0025 in the same band using a Sersic and PSF profiles, respectively. As shown in these panels, J0121+0025 is better modeled using a Sersic profile than a PSF, being the r.m.s of the residuals in the region encompassing J0121+0025 reduced by ≈50 per cent.</figcaption></figure>

Data were processed with standard IRAF<fnref n="2" /> tasks. Each individual 2D spectrum is bias-subtracted and flat-field corrected. The wavelength calibration is done using HgAr+Ne+Xe arc lamps data obtained in both days. Individual 2D spectra were background subtracted using sky regions around J0121+0025 (≈10 arcsec on both sides). The slit includes a faint (R = 23.6, z<sub>phot</sub> = 1.10 ± 0.27) star-forming galaxy ≃6.6 arcsec north from J0121+0025 (see Fig. 1), so this region was excluded for the background statistics. The continuum of this source is detected, but not any emission line covered in the OSIRIS spectral range. Individual 1D spectra are extracted, stacked, and corrected for the instrumental response using observations of the standard star Ross 640 observed in both nights. The reddening effect for the extinction in the Galaxy was taken into account adopting the extinction curve of Cardelli, Clayton & Mathis (1989) and using the extinction map of Schlafly & Finkbeiner (2011). Finally, the flux of the spectrum is matched to that obtained from photometry in the R-band to account for slit-losses, and the final spectrum is corrected for telluric absorption using the IRAF *telluric* routine. The reduced spectrum of J0121+0025 is shown in Fig. 2 and presents very high signal-to-noise ratio in the continuum, SNR ≃ 20 − 30 per spectral bin.

## 2.2 Ancillary data

J0121+0025 falls in the SDSS equatorial Stripe 82 and the rich ancillary data available in this field are used. These consist of optical imaging from the Hyper Suprime-Cam (HSC) on Subaru in the G, R, Z, and Y bands from the second data release of the HSC Subaru Strategic Program (Aihara et al. 2019) and MEGACAM image in the I band on the Canada-France-Hawaii Telescope (CFHT), processed and stacked using the MegaPipe image staking pipeline (Gwyn 2008). The optical images are of great quality, both in terms of seeing conditions (0.50–0.85 arcsec FWHM) and depth, reaching 5σ limits similar or above ≃25.5 for all bands. J0121+0025 is detected in all images. Fig. 1 shows a cutout of J0121+0025 from the Subaru/HSC R-band. J0121+0025 shows a compact morphology, barely resolved

(≈0.55 arcsec FWHM) only in the best seeing conditions image, R-band, with a PSF of 0.50 arcsec FWHM measured using several stars in the field. The lack of lensing structures, such as multiple images or arc-like morphologies, and the compact morphology make unlikely that J0121+0025 is being magnified by gravitational lensing. Using aperture photometry with a diameter of 2.5 × FWHM, J0121+0025 shows a flat spectral energy distribution (SED) in the optical with magnitudes of ≃21.60 AB in the R- to Y-bands, consistent with those from SDSS photometry. Table 1 summarizes the photometry of J0121+0025.

Near-IR imaging is also available in this field. J0121+0025 was observed with the VISTA InfraRed CAMera (VIRCAM) as part of the VISTA-CFHT Stripe 82 (VICS82) survey (Geach et al. 2017). It is detected in the J- and Ks-bands, with magnitudes of 21.87 ± 0.23 and 21.46 ± 0.18, respectively. J0121+0025 is also included in the Spitzer/HETDEX Exploratory Large-Area (SHELA) survey (Papovich et al. 2016) and is detected in the two first IRAC channels at 3.6 µm (I1) and 4.5 µm (I2) with magnitudes of 21.61 ± 0.16 and 22.01 ± 0.20, respectively. Finally, this region has been imaged in the X-ray by XMM-Newton with a total integration time of 2.6 ks. However, J0121+0025 is not detected with an X-ray flux limit of ≃6 × 10<sup>-15</sup> erg s<sup>-1</sup> cm<sup>-2</sup> (0.2–2 keV), corresponding to a luminosity 3σ limit of 1.2 × 10<sup>45</sup> erg s<sup>-1</sup> at z = 3.25 (considering a photon index Γ = 1.7).

# 3 RESULTS

## 3.1 The nature of the ionizing source: SFG or AGN?

The brightness of J0121+0025 (R = 21.6) and its corresponding UV luminosity at z = 3.25, M<sub>UV</sub> = −24.1, rival that of bright QSOs at similar redshift (e.g. Pâris et al. 2018). Therefore, it is critical to investigate first the nature of the ionizing source of J0121+0025.

The high SNR rest-frame UV spectrum (Fig. 2) is rich in absorption/emission features that are common in young starbursts, rather than in AGNs. Features associated with different components of the galaxy, such as stars, nebular emission and ISM, are clearly identified, and marked in Fig. 2 with different colours.

In particular, stellar wind P-Cygni profiles and photospheric absorption lines are detected with high significance (green and

<footnote n="2">²http://iraf.noao.edu/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** GTC/OSIRIS rest-frame UV spectrum of J0121+0025 (black) and its corresponding $1\sigma$ uncertainty (blue). Yellow dotted lines marked below the spectrum identify several photospheric absorption lines, some of them resolved and detected with high significance for which the systemic redshift was derived, $z_{\text{sys}} = 3.244 \pm 0.001$. Low-ionization ISM absorption, wind lines in the form of P-Cygni, and nebular emission are marked in blue, green, and magenta lines, respectively. The best-fitting Starburst99 model with age of 3 Myr, $Z_{*}/Z_{\odot} = 0.4$, and $E(B-V) = 0.04 \pm 0.02$ is shown in red. The S99 spectrum blueward Ly$\alpha$ has been corrected for the Lyman forest absorption, using the mean and standard deviation IGM transmission $T(\text{IGM}) = 0.60 \pm 0.19$ (see Section 3.3). The spectral region in green corresponds to the emission below $\lambda_{0} < 912$ Å related to LyC leakage.</figcaption></figure>

<table><caption>**Table 1.** Optical to mid-IR photometry of J0121+0025.</caption><thead><tr><th>Band</th><th>λ<sub>eff</sub><br>(μm)</th><th>Magnitude<br>(AB)</th><th>Telescope/instrument</th></tr></thead><tbody><tr><td>G</td><td>0.47</td><td>21.98 ± 0.08</td><td>Subaru/HSC</td></tr><tr><td>R</td><td>0.61</td><td>21.60 ± 0.05</td><td>Subaru/HSC</td></tr><tr><td>I</td><td>0.77</td><td>21.57 ± 0.06</td><td>CFHT/MEGACAM</td></tr><tr><td>Z</td><td>0.89</td><td>21.53 ± 0.08</td><td>Subaru/HSC</td></tr><tr><td>Y</td><td>1.00</td><td>21.58 ± 0.08</td><td>Subaru/HSC</td></tr><tr><td>J</td><td>1.25</td><td>21.87 ± 0.23</td><td>VISTA/VIRCAM</td></tr><tr><td>K<sub>s</sub></td><td>2.14</td><td>21.46 ± 0.18</td><td>VISTA/VIRCAM</td></tr><tr><td>I<sub>1</sub></td><td>3.56</td><td>21.61 ± 0.16</td><td>Spitzer/IRAC</td></tr><tr><td>I<sub>2</sub></td><td>4.51</td><td>22.01 ± 0.20</td><td>Spitzer/IRAC</td></tr></tbody></table>

yellow dashed lines in Fig. 2). The detection of photospheric lines in J0121+0025 indicates unambiguously that the UV luminosity is dominated by stellar emission, rather than an AGN. We identify more than ten photospheric features. Some of them are resolved and detected with high significance (e.g. C II 1324 Å, O IV 1343 Å, and S V 1501 Å; see Fig. 3). We use these to determine the systemic redshift $z_{\text{sys}} = 3.244 \pm 0.001$ of J0121+0025. Others are seen in blends from multiple transitions (e.g. Si III 1417 Å, C III 1427 Å, and Fe V 1430 Å at $\lambda_{0} \simeq 1415-1435$ Å). These stellar absorption lines are intrinsically weak in star-forming galaxies, with $EW_{0}$ typically well below 1 Å (e.g. Shapley et al. 2003; Steidel et al. 2016; Rigby et al. 2018). As they are formed in the photospheres of hot stars and are seen in absorption, the background radiation should be dominated by the starlight, otherwise they would not be detected. Even a small contribution of an AGN to the UV continuum ($\lesssim 25$ per cent), that is featureless in these spectral regions, would make these lines disappear at the SNR of our spectrum. In addition, the observed P-Cygni profiles in N V 1240 Å and C IV 1550 Å can be also well explained/modelled by stellar models with a very young age ($\simeq 3$ Myr burst; see Fig. 3 and Section 3.2 for details), similar to those seen in other very young starbursts (e.g. Rivera-Thorsen et al. 2019; Vanzella

et al.), some of them also very/extremely luminous (Vanzella et al. 2018; Marques-Chaves et al. 2020b). While some rare AGNs, such as broad or narrow absorption line QSOs (BAL/NAL QSOs), can show N V and C IV profiles mimic those of stellar P-Cygni, from the combination of a broad emission and a redshifted absorption (see for example Bentz, Osmer & Weinberg 2004; Appenzeller et al. 2005), photospheric lines are not present in the spectra of AGNs.

The rest-frame UV morphology of J0121+0025 appears compact, but there is evidence of a resolved structure. Using the best-seeing condition image (R-band from Subaru, 0.50 arcsec FWHM), J0121+0025 appears marginally resolved with a FWHM $\simeq 0.55''$, that corresponds to $\simeq 1.5-2.0$ kpc proper. Using GALFIT (Peng et al. 2002), the light distribution of J0121+0025 is better modeled using a Sersic profile instead of using a PSF model, with residuals in the region encompassing J0121+0025 reduced by $\simeq 50$ per cent (Fig. 1). This suggests that the source is spatially resolved. The best-fitting model (Sersic profile) gives an effective radius $r_{\text{eff}} = 0.6$ pix, that corresponds to $r_{\text{eff}} = 0.1$ arcsec assuming the Subaru pixel scale $0.168$ arcsec pix$^{-1}$, and $r_{\text{eff}} \sim 0.8$ kpc at $z = 3.244$. Although it has been shown that GALFIT can recover effective radius down to $\sim 0.5$ pix if the PSF is properly known and the source is bright enough (see Vanzella et al. 2017), we assume conservatively an $r_{\text{eff}} < 1$ kpc.

From another perspective, the Ly$\alpha$ line shows a narrow profile with an intrinsic FWHM $\simeq 350$ km s$^{-1}$ compatible with star-forming galaxies (see Fig. 3 and Section 3.2.3). The AGN population typically shows broad Ly$\alpha$ profiles, up to several hundreds or thousands km s$^{-1}$, even in the case of narrow line AGNs (e.g. Hainline et al. 2011). In addition, the observed Ly$\alpha$ equivalent width, flux and corresponding luminosity can be well explained by star-formation only (discussed in Section 3.2), without the need to invoke an AGN contribution. He II 1640 Å is also detected, showing a broad profile with $\simeq 2500$ km s$^{-1}$. Because He II is a non-resonant (and recombination) line, its origin is likely stellar and not nebular (from an AGN), otherwise we would expect to detect a similar
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Spectral features detected in J0121+0025. Top left: Lyα spectral profile seen in the GTC/OSIRIS (black) and BOSS (green) spectra. It shows a narrow profile with an intrinsic FWHM = 350 ± 40 km s⁻¹ with its peak redshifted by ≃120 km s⁻¹ relative to the systemic velocity. Top middle and right-hand panels show the N V and O VI wind lines in the form of P-Cygni (black) and the best-fitting S99 model (red). Bottom left: Stellar photospheric lines used to derive the systemic redshift. Bottom middle: Profiles of the low-ionization ISM lines Si II 1260 Å and C II 1334 Å (blue and red, respectively). These lines are weak and have their centroids redshifted respect to the systemic velocity by ≃−460 and ≃−510 km s⁻¹, respectively. The spectrum also shows two other absorption lines (marked with dashed lines), whose nature is still unclear, but likely not physically associated with J0121+0025 (e.g. outflows). Bottom right: Peculiar profile of the He II line with three peaks (yellow). S99 and BPASS models are also shown in red and blue, respectively, but they fail to reproduce the observed emission.</figcaption></figure>

or even broader profile in the resonant Lyα line, which is not the case. This will be discussed in more detail in Section 3.2. Nebular emission of [O III] 1666 Å is also detected (with low significance) and shows a narrow profile, not resolved in our OSIRIS spectrum (<450 km s⁻¹; Figs 2). This line is not present in the spectra of typical AGNs.

The presence of a UV-faint or an obscured type-II AGN is more difficult to exclude. Unfortunately, we cannot use rest-frame UV and optical line diagnostics (e.g. Baldwin, Phillips & Terlevich 1981; Nakajima et al. 2018) to discriminate between star-formation and AGN, as these lines are not covered by the OSIRIS spectrum. However, the mid-IR photometry of J0121+0025 disfavours the presence of such AGN contribution (top panel of Fig. 4). J0121+0025 shows a flat/blue SED from the rest-frame UV to near-IR with R − I2 = −0.41 ± 0.21 and I1 − I2 = −0.40 ± 0.31, where R, I1, and I2 bands probe rest-frame wavelengths of ≃0.16, 0.84, and 1.06 µm, respectively. As shown in Fig. 4, the optical to mid-IR colours of J0121+0025 place it far away from the locus of AGNs at similar redshift (Pâris et al. 2018). AGNs tend to have

red optical-to-mid-IR SEDs due to the rising emission at $\lambda_0 \gtrsim 1$ µm (rest) originated by the dust torus (Assef et al. 2013). We also check for possible variability from an AGN in the optical photometry. In the bottom panel of Fig. 4 we compare the observed magnitudes in the *g*, *r*, and *i* bands of J0121+0025 from Subaru, CFHT, SDSS, DECaLS, and Pan-STARRS1 that probe different epochs, from MJD 51544 to 58362 (~18 yr). No variability is detected in J0121+0025. Lastly, J0121+0025 is not detected in X-rays in the 2.6 ks XMM-Newton data, although the corresponding $3\sigma$ limit $L_{(0.2-2)\text{keV}} = 1.2 \times 10^{45}$ erg s⁻¹ is not deep enough to further explore a possible X-rays emission of an AGN. Note that a significant X-ray emission from star-formation is still expected in J0121+0025. Assuming SFR = 981M$_\odot$ yr⁻¹ (see Section 3.4) and following Grimm, Gilfanov & Sunyaev (2003), we expect an X-ray emission from star-formation $L_x(\text{SFR}) \sim 10^{43}$ erg s⁻¹.

Overall, the shape of the rest-frame UV spectrum, the detection of stellar features (photospheric absorption and wind lines), the resolved morphology, and the multiwavelength SED highly suggest that the luminosity of J0121+0025 is being powered by a vigorous starburst
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top: Comparison between the 3.4–4.6 and 0.6–4.6 µm colours of J0121+0025 (blue) and those from SDSS/BOSS QSOs (Pâris et al. 2018) at similar redshift (red). J0121+0025 shows bluer colours than typical QSOs. For comparison, the young starburst BOSS-EUVLG1 (≈4 Myr; Marques-Chaves et al. 2020b) is also shown in yellow. Bottom: Observed magnitudes of J0121+0025 in various filter bands (g, r, and i) with different telescopes and epochs. No variability is detected in J0121+0025.</figcaption></figure>

and that a significant contribution of an AGN at these wavelengths is unlikely.

## 3.2 Rest-frame UV properties

### 3.2.1 Young stellar population: age, metallicity, and attenuation

One of the most prominent features in the spectrum of J0121+0025 is the P-Cygni associated with wind lines that are much stronger than in typical LBGs (e.g. Shapley et al. 2003). These lines are produced by strong outflows of material from the most massive stars, whose strength and spectral shapes depend strongly on the age and metallicity (and on the IMF) of the stellar population, being N v and C IV by far the most sensitive features (Chisholm et al. 2019).

To infer the properties of the young stellar population in J0121+0025, the observed line profiles of N v and C IV are compared to those obtained with the spectral synthesis code STARBURST99 (S99: Leitherer et al. 1999), following the same methodology described in Marques-Chaves et al. (2018, 2020a). We use S99 instead of BPASS models (Stanway, Eldridge & Becker 2016), because the later are less able to match the details of the P-Cygni absorption of C IV (see discussion and fig. 5 in Steidel et al. 2016). Briefly, we generate high-resolution (0.4 Å) UV spectra using standard Geneva tracks with a grid of metallicities ($Z_*/Z_⊙$, where $Z_⊙ = 0.02$) of 0.05, 0.2, 0.4, and

1, and burst ages from 1 to 30 Myr. An IMF with a power slope index $\alpha = -2.35$ over the mass range $0.5 < M_*/M_⊙ < 100$ is considered. S99 outputs are redshifted to $z = 3.244$, smoothed to the spectral resolution of the OSIRIS spectrum, and rebined to the spectral bin of 2.12 Å. Dust attenuation is also taken into account, considering values $E(B-V)_*$ ranging from 0 to 0.2 and using the Calzetti et al. (2000) extinction curve and its extension to short wavelengths ($<0.15$ µm) provided by Reddy et al. (2016). Spectral windows that are free of absorption/emission features (Rix et al. 2004) and strong sky-subtracted residuals are used to offset the flux of S99 models. We then compare the observed N v and C IV profiles with those from S99, performing a $\chi^2$ minimization over the spectral range 1225–1245 Å for N v and 1528–1538 for C IV, excluding in the fit spectral regions that could be affected by interstellar absorption or nebular emission, which is particularly relevant for C IV.

Wind lines of N v and C IV are well reproduced by a $Z_*/Z_⊙ \simeq 0.4$ and $\simeq 3$ Myr burst of star formation (red in Fig. 2). A colour excess of the stellar continuum $E(B-V)_* = 0.04 \pm 0.02$ is also inferred, which is compatible with the observed UV slope, $\beta_{UV} = -2.05 \pm 0.10$. Scaling the best-fitting S99 model to $M_{UV} = -24.1$, this leads to a burst mass log($M_*/M_⊙$) = 9.8. According to S99, the number of O-type stars is $\sim 8 \times 10^5$ yielding an intrinsic ionizing photon production rate $N_{\text{int}}(\text{LyC}) \simeq 1.4 \times 10^{55}$ s$^{-1}$ and the production efficiency, $\xi = N_{\text{int}}(\text{LyC})/L_{\text{UV, int}}$, is log($\xi$) = 25.2.

Considering the full spectral range, the overall agreement between the best-fitting model and the observed spectrum is mixed. Some stellar features are well fitted, as the N v and C IV P-Cygni and some photospheric absorption, while others show poor agreement. The profile of Si IV 1393,1402 Å shows evidence of a P-Cygni contribution (in addition to the ISM absorption), but the model underpredicts it, which could suggest an additional contribution of a slightly older stellar population with age $\simeq 5$ Myr (see fig. 6 in Leitherer et al. 2001). It is also interesting to note that the OSIRIS spectrum shows evidence of a P-Cygni around the photospheric blanketing by Fe v and O v around $\simeq 1360-1375$ Å, but it appears slightly blueshifted with respect to the predicted S99 model, by $\simeq 400-500$ km s$^{-1}$. On the other hand, other stellar features are relatively well reproduced, such as the region around $\simeq 1420-1430$ Å from the blended emission of Si III, Fe v, and C III transitions (also called the ‘1425’ index), the photospheric S v 1501 Å line, and the Fe IV complex around $\sim 1600$ Å (see Fig. 2). The model also predicts a relatively strong P-Cygni in O VI 1031,1033 Å, but is still weak compared to the observed one (Fig. 3), which could indicate an even younger stellar population ($\le 2$ Mry). Note, however, that the profile of O VI is likely affected by the contribution of the $Ly\beta$ absorption and the IGM attenuation that could impact the observed profile.

Overall, the best-fitting S99 model with $Z_*/Z_⊙ = 0.4$, age of 3 Myr and $E(B-V)_* = 0.04$ fits reasonably well the observed spectrum of J0121+0025, in particular the wind lines N v and C IV, that are the most sensitive features to these parameters (Chisholm et al. 2019). Note, however, that the inferred $Z_*$ and age are model-dependent and should be treated with caution and considered as approximated values. This is particularly relevant for the metallicity, which is less constrained, as our analysis is limited to discrete models with a set of four metallicities. Consequently we can only say that the model with $Z/Z_⊙ = 0.4$ is favoured with respect to the other models ($Z/Z_⊙ = 0.05, 0.2, 1$). In addition, we are assuming a single burst model, which might not be realistic (nor a continuum star-formation rate). Nevertheless, the age is much better constrained and should be <8 Myr, otherwise the wind line of N v would appear much weaker ($\simeq 10$ Myr) or almost non-existent (>10 Myr). The same applies if a continuous star-formation history (SFH) is assumed. In
this case, the redshifted emission of N v could be well described with a continuous SFH with an age up to ≃20 Myr, but not the corresponding blueshifted absorption, that would be underestimated for ages ≳ 10 Myr due to the increasing contribution of B-type stars to the UV continuum. In addition, the blue mid-IR colour I1 − I2 = − 0 . 4 ± 0 . 31 highly supports a very young age of the stellar population that would be difficult to explain with a sightly older stellar population ( ≳ 15 − 20 Myr) with a continuous star-formation rate.

3.2.2 Broad He II emission

The He II 1640 Å emission is shown in detail in the bottom right-hand panel of Fig. 3 and presents a complex morphology characterized by a broad profile with two absorption in the central part of the line forming a triple peaked emission. Fitting a Gaussian profile and excluding the two absorption features in the fit, we measure a FWHM ≃ 2500 km s⁻¹. Because of its broadness, this emission is likely stellar in origin, otherwise we would expect a similar or even broader profile in other nebular lines like the resonant Lyα line. However, this is not the case as Lyα shows a narrow profile (≈350 km s⁻¹ FWHM; see Section 3.2.3). The nebular [O III] 1666 Å line is also detected and show an unresolved profile in the OSIRIS spectrum (<450 km s⁻¹ FWHM; see Fig. 2).

We measure a rest-frame equivalent width EW$_0$ = 3.2 ± 0.3 Å for He II (corresponding to the yellow region of Fig. 3) that is much larger than that inferred in the z ~ 3 LBG composite spectrum of Shapley et al. (2003) (EW LBG 0 = 1.3 ± 0.3). The He II emission in J0121+0025 is also stronger than the average EW$_0$ ≃ 2.5 Å found in some extreme Wolf–Rayet (WR) star clusters in the local Universe (Chandar, Leitherer & Tremonti 2004; although see the cases of NGC 3125-1 with EW$_0$ ≃ 7.4 Å, or the dwarf galaxy II Zw 40 with EW$_0$ ≃ 7.1 Å; Leitherer et al. 2018). As shown in Fig. 3, the best-fitting S99 model clearly underpredicts the strength of the observed He II profile. In fact, Brinchmann, Pettini & Charlot (2008) predict a EW$_0$(S99) = 0.3 Å for a S99 burst model with Z$_*$/Z$_⊙$ = 0.4 and 3 Myr age. Even considering an extreme case with a continuous SFH and Z$_*$/Z$_⊙$ = S99 models predict EW$_0$(S99) ≤ 2.4 Å (see fig. 2 of Brinchmann et al. 2008). For comparison, a BPASS binary model (Stanway et al. 2016) with the same age and metallicity is also shown in Fig. 3, and although it predicts stronger stellar He II emission than S99 models (see also Steidel et al. 2016), it still under-predicts the observed emission in J0121+0025.

The strength of He II in J0121+0025 now raises the question regarding the presence of more exotic stellar populations. For example, He II appears very strong in the spectra of very massive stars (VMS; > 100 M$_{⊙}$) in the central cluster R136 of the 30 Doradus star-forming region (see Crowther et al. 2016). In fact, the contribution of such massive stars has been proposed to explain the large EW$_0$'s in He II in a few local star-forming galaxies (with EW$_0$ ≃ 2.0–4.7 Å; Senchyna et al. 2021). Interestingly, the complex triple-peaked He II profile seen in the spectrum of J0121+0025 resembles that observed in two star-forming galaxies analysed by Senchyna et al. (2021), namely SB 179 and SB 191 (their fig. 6), that according to these authors could be the product of rapid rotation of rare Onfp stars (e.g. Walborn et al. 2010). Investigating in detail the nature of the He II line and WR and VMS content (and possible nebular contribution) is out of the scope of this work as it requires follow-up observations (e.g. high-spectral resolution of the He II line and near-IR spectroscopy to put strong constraints on the metallicity), and in addition updated stellar models with the inclusion of very massive stars.

3.2.3 The Lyα line

The Lyα line in J0121+0025 shows a spectrally unresolved profile in the OSIRIS spectrum (R ≃ 800), but it is slightly resolved in the higher resolution BOSS spectrum (Fig. 3). Fitting a Gaussian profile, we measure an intrinsic FWHM = 350 ± 40 km s⁻¹, after correcting for the instrumental broadening (≈150 km s⁻¹), and a rest-frame equivalent width EW$_0$(Lyα) = 14 ± 3 Å. The Lyα line has its peak closed to the systemic redshift, redshifted by v$_{peak}$ ≃ 120 ± 50 km s⁻¹. It is still not clear why the measured EW$_0$(Lyα) appears so low, compared to the intrinsic EW int 0 (Lyα) ~ 100 Å expected for a ≃3 Myr age burst and Z$_*$/Z$_⊙$ (Schaerer 2003). Roughly, this yields a Lyα escape fraction of ~ 14 per cent, which is lower than f$_{esc,abs}$(LyC) ≈ 40 per cent (see Section 3.3). Considerable fiber/slit losses, strong IGM attenuation near λ$_0$ ≃ 1215 Å, destruction of Lyα photons by dust, and/or large f$_{esc}$ of ionizing photons could in principle explain such differences. In addition, there could be more ionizing photons available than H I to be ionized, i.e. the ISM in J0121+0025 could be mostly ionized approaching to a density-bounded geometry. Such a scenario should be further investigated using, e.g. the [O II] 3727,29 Å and [O III] 5008 Å lines ([O III]/[O II] ratio). In addition, the H β line could provide constraints on the properties of the Lyα emission, both in terms of its intensity and spectral shape. These lines ([O II], [O III], and H β) are redshifted to the near-IR H- and K-bands and are accessible from the ground.

Using the BOSS/SDSS spectrum that is less affected by flux losses, we measure a total Lyα flux of F(Lyα) = (5.72 ± 0.10) × 10⁻¹⁶ erg s⁻¹ cm⁻². This corresponds to a luminosity log(L$_{Lyα}$/erg s⁻¹) = 43.8 ± 0.1 at the redshift of J0121+0025. To test whether or not this luminosity could be explained by star-formation, we compare the SFR of J0121+0025 obtained from SED fitting (SFR(SED) = 981 ± 232M$_{⊙}$ yr⁻¹; see Section 3.4) to the Lyα star-formation rate using the Kennicutt (1998) conversion. Assuming case-B recombination and the Salpeter (1955) IMF, the Lyα luminosity corresponds to a SFR(Lyα) ≃ 80 M$_{⊙}$ yr⁻¹. Even considering f$_{esc}$(LyC) > 0 (and so SFR(Lyα) ≳ 80 M$_{⊙}$ yr⁻¹), star-formation can naturally explain the observed Lyα flux and the corresponding luminosity.

3.2.4 ISM, kinematics, and covering fraction

In addition to the strong Lyα line and other stellar features, the high S/N spectrum reveals also absorption features that are associated with the ISM gas and produced by the resonance transition of several ionic species.

Low-ionization ISM lines (LIS) of Si II 1260 Å, C II 1334 Å, and Si II 1526 Å are detected in the spectrum. Absorption in O I 1302 Å and Si II 1304 Å is also detected but it is not resolved and, in addition, its profile is also contaminated by the photospheric complex around λ₀ ≃ 1294–1298 Å. Other LIS lines that are usually present in the spectra of LBGs, such as Fe II 1608 Å or Al II 1670 Å (Shapley et al. 2003), are not detected in J0121+0025.

Since LIS lines are seen against the continuum provided by starlight of J0121+0025, they are useful to probe the kinematics of the gas along the line of sight. Despite the low spectral resolution of the OSIRIS spectrum, the centroids of the LIS lines are clearly blueshifted with respect to the systemic redshift, by v$_{peak}$(LIS) ≃ −450 km s⁻¹ (Fig. 3). This could be an indication of fast gas outflows caused by supernova explosions and stellar winds, similar to those detected in the other extremely UV luminous starburst BOSS-EUVLG1 (v$_{peak}$(LIS) ≃ −400 km s⁻¹; Marques-Chaves et al. 2020b;
<figure><img src="image_7.png" /></figure>

**Figure 5.** 2D (top) and 1D (bottom) spectra of the far-UV region of J0121+0025. The 2D spectrum has been smoothed for visualization purpose. The GTC spectrum is in black and its corresponding $1\sigma$ uncertainty is in grey. Features associated with the Lyman series (from Ly$\alpha$ to the Lyman limit), ISM, and stellar absorption (photospheric and wind lines) are marked below the spectrum in green, blue, and yellow colours, respectively. The low-resolution best-fitting S99 model (3 Myr age, $Z_{\odot}/Z_{\oplus}$ = 0.4 and $E(B-V)_* = 0.04$) is plotted in red, one corrected for the IGM transmission ($<T(IGM)> = 0.60 \pm 0.19$, solid red) and other assuming $T(IGM) = 1$ (dashed red). The high-resolution S99 model is also shown in blue. Horizontal grey lines mark the spectral windows used to infer $T(IGM)$ (see text). These exclude the regions associated with the Lyman series and ISM from J0121+0025, so that we probe H I absorbers from the Lyman forest in the line of sight only. Vertical lines above the spectrum mark the position of four strong H I absorbing systems identified at $z = 3.076, 2.898, 2.816$, and 2.733 (#1 to #4, respectively). We also note that the 2D spectrum also shows a faint continuum offset by $\sim 6$ arcsec of J0121+0025 from a low-$z$ star-forming galaxy (with a $z_{phot} = 1.10$; see Fig. 1).

Álvarez-Márquez et al. 2021) or in the Sunburst LyC emitter (Rivera-Thorsen et al. 2017; Vanzella et al. 2021).

The OSIRIS spectrum also shows two additional absorption lines at $\simeq 5320$ Å and $\simeq 5640$ Å (Fig. 3), whose origin is still not clear. These lines appear close to Si II 1260 Å and C II 1334 Å of J0121+0025, by $\simeq -1700$ and $\simeq -1200$ km s$^{-1}$, respectively (Fig. 3). These lines could be associated with outflows from J0121+0025, but such scenario is unlikely. Because these lines have broadly similar ionization potentials, the outflowing gas traced by Si II and C II should be kinematically coherent (and co-spatial), which is not the case. Moreover, the profile of Si II 1526 Å does not show this secondary absorption component. It is possible that the absorption line at $\simeq 5320$ Å is associated with Si IV 1393 Å at $z = 2.816$, for which Ly$\alpha$ is also detected in absorption around $\simeq 4640$ Å (H I absorbing system #3 in Fig. 5). In such a case, the corresponding Si IV 1402 Å absorption of this intervening system is contaminating the profile of Si II 1260 Å of J0121+0025. For the other absorption at $\simeq 5640$ Å, we are not able to find any redshift solution, so the redshift of this intervening system is still not known.

Regarding the strength ($EW_0$) of LIS lines in J0121+0025, we fit Gaussian profiles to the absorption profiles, excluding the contribution of the secondary not-related absorption component mentioned before (Fig. 3). We measure rest-frame equivalent widths of $1.05 \pm 0.15$ Å, $0.59 \pm 0.08$ Å, and $0.81 \pm 0.10$ Å, for Si II 1260 Å, C II 1334 Å, and Si II 1526 Å, respectively. Note that the measured equivalent width of Si II 1260 Å should be considered

as an upper limit, because its profile is likely contaminated by the intervening metal lines Si IV 1402 Å at $z = 2.816$. These lines appear spectrally resolved with an intrinsic FWHM $\simeq 550-650$ km s$^{-1}$. For comparison, the $z \sim 3$ LBG composite spectrum of Shapley et al. (2003) shows large $EW_0$'s for the same lines, with $EW_0 \simeq 1.7$ Å.

The weakness of LIS lines in J0121+0025 could arise either by a low geometric covering fraction of the gas, $C_f$, a low ion column density, and/or a highly ionized ISM. Considering the linear part of the curve of growth (i.e. low column density), the ratios of the $EW_0(1260)/EW_0(1526)$ can be related through their oscillator strengths, for which we would expect a $EW_0(1260)/EW_0(1526) \simeq 5$ if these lines were not saturated. However, we measure $EW_0(1260)/EW_0(1526) \lesssim 1.3$, suggesting that at least one of these lines is saturated. This is not surprising at all, as these lines appear almost always saturated in the spectra of star-forming galaxies, even in damped-Ly$\alpha$ systems with sub-solar metallicities (e.g. Dessauges-Zavadsky et al. 2006).

Taking this in consideration, it is possible that the weakness of LIS lines arises from a low $C_f$. Assuming the optically thick regime and an ionization-bounded ISM with a uniform dust-screen geometry, the $C_f$ can be inferred using the residual intensity of the absorption line, $I$, so that $C_f = 1 - II_{I0}$, where $I_{0}$ is the continuum level. We measure $II_{I0} \simeq 0.8$ for Si II 1260 Å and C II 1334 Å, yielding to $C_f(\text{SiII}) \simeq 0.2$. Note that with the low spectral resolution of our data we are likely overestimating $II_{I0}$; however, this effect should not be
dominant as the lines are spectrally resolved. Following Gazagnes et al. (2018), this yields a neutral gas covering fraction $C_f(\text{HI}) \simeq 0.55$, for that a significant fraction of LyC photons could escape. In fact, using the prescriptions of Chisholm et al. (2018; see also Saldana-Lopez, in preparation), the inferred $C_f(\text{HI}) \simeq 0.55$ leads to a predicted $f_{\text{esc,abs}}^{\text{pred}}(\text{LyC}) \approx 0.25$, which is consistent with the observed value from the spectrum (next section).

## 3.3 Lyman continuum radiation

### 3.3.1 The direct detection of LyC

One of the most remarkable features observed in the OSIRIS spectrum of J0121+0025 is the detection of emission below $\lambda_{\text{obs}} \simeq 3880$ Å, i.e. $\lambda_0 < 911.8$ Å, that can be related to LyC leakage. This emission is real and it is not related with detector artifacts (e.g. cosmic rays, flat-field correction) or bad sky subtraction. It is detected in the spectra of different nights, as well as in the much shallower BOSS spectrum although with less SNR ($\simeq 2.3\sigma$).

A zoom-in into this region is shown in Fig. 5. The observed flux emission at rest-frame 880–910 Å has a total SNR of 7.9 and an average SNR per spectral bin of 0.98. The mean flux density in this spectral region is $f_{900}(\text{obs}) = 0.781 \pm 0.099\mu\text{Jy}$, which corresponds to a magnitude of $m_{900} = 24.17$ (AB). For comparison, the mean non-ionizing UV flux density, estimated from the OSIRIS spectrum over the rest-frame range 1490–1510 Å, is $f_{1500} = 8.874 \pm 0.399\mu\text{Jy}$. Combining these measurements, we find a ratio of the ionizing to non-ionizing flux density ($f_{900}/f_{1500}$)$_\text{obs} = 0.088 \pm 0.012$, which corresponds to a $\Delta m = 2.64$ (AB). It is worth noting that the ionizing emission $\lambda_{\text{obs}} < 3880$ Å suffers for an apparent absorption at $\lambda_{\text{obs}} \simeq 3813$ Å. After careful inspection, this absorption could be related with $\text{Ly}\beta$ and $\text{Ly}\delta$ absorption associated with two H I systems at $z = 2.733$ and $z = 3.076$, respectively (#1 and #4 in Fig. 5).

To infer the relative LyC photon escape fraction, we compare the observed ionizing flux density of J0121+0025 at 900 Å, $f_{900}(\text{obs})$, to that from the S99 model, $f_{900}(\text{S99})$, that best represents the shape of the rest-frame UV spectrum (Section 3.2.1), i.e. a burst with an age of 3 Myr, $Z_{\odot}/Z_{\odot} = 0.4$ and $E(B-V) = 0.04$ (Figs 3 and 5). The S99 model has been already corrected for $E(B-V)$, thus it incorporates the relative attenuation between 900 and 1500 Å due to the variation of the attenuation coefficient with the wavelength [$k_{900} \simeq 14.3$ and $k_{1500} \simeq 10.3$, assuming the Calzetti et al. (2000) curve and its extension to the UV provided by Reddy et al. (2016)].<fnref n="3" /> To probe the region below $\lambda_0 < 912$ Å, we use the low-resolution version of S99 model because it extends to the LyC region. Finally, we also take into account the contribution of the IGM transmission, $T(\text{IGM})$. The relative LyC photon escape fraction, $f_{\text{esc,rel}}(\text{LyC})$, can thus be expressed using the following formulation:

$$f_{\text{esc,rel}}(\text{LyC}) = \frac{f_{900}(\text{obs})}{f_{900}(\text{S99})} \times \frac{1}{T(\text{IGM})} \quad (1)$$

Both $f_{900}(\text{obs})$ and $f_{900}(\text{S99})$ are measured using the same spectral window, defined at 880–910 Å in the rest-frame. We find $f_{900}(\text{obs})/f_{900}(\text{S99}) = 0.34 \pm 0.04$, where the uncertainties arise from $f_{900}(\text{obs})$. We note that equation (1) is consistent to that used in the literature (e.g. Shapley et al. 2016; Vanzella et al. 2016; Steidel et al. 2018), where $f_{900}/f_{1500}$ is used instead of $f_{900}$, because in our case $f_{1500}(\text{S99})$ has been already matched to $f_{1500}(\text{obs})$ (see Fig. 2), thus $f_{1500}(\text{obs}) \simeq f_{1500}(\text{S99})$.

A precise estimate of $f_{\text{esc,rel}}(\text{LyC})$ is not possible given the stochastic nature of $T(\text{IGM})$ and the large fluctuation of the attenuation in one single line of sight (Inoue & Iwata 2008; Inoue et al. 2014). However, it is still possible to place rough constraints. Assuming that $f_{900}(\text{obs})/f_{900}(\text{S99}) = 0.34 \pm 0.04$ is well constrained,<fnref n="4" /> this implies that $T(\text{IGM})$ should be at least larger than $>0.34$ to keep a physical $f_{\text{esc,rel}}(\text{LyC}) < 1$ (e.g. Vanzella et al. 2012). On the other hand, $f_{\text{esc,rel}}(\text{LyC})$ must be $\ge 0.34$, where the most extreme value (0.34) stands for a completely transparent IGM.

To get a more quantitative estimate of $T(\text{IGM})$, we use the non-ionizing part of the OSIRIS spectrum, from 912–1215 Å, and compare it to that of S99 best-model. We exclude spectral regions that are associated with the Lyman series and ISM from J0121+0025 (marked in Fig. 5) that are not included in the S99 models, so that we probe H I absorbers from the Lyman forest in the line of sight only, and not the the CGM and ISM associated to J0121+0025. The spectral regions used to estimate $T(\text{IGM})$ are marked as horizontal grey lines in Fig. 5. Note, however, that with this comparison, we are inferring the $T(\text{IGM})$ in this spectral range, not necessarily at $\le 912$ Å, still it can serve as a first-order approximation. It is worth noting that the S99 model predicts a relatively strong break around the Lyman limit, which is not observed in the spectrum of J0121+0025 (and in other LyC leakers; see e.g. Steidel et al. 2018). We find a mean value and standard deviation $T(\text{IGM}) = 0.60 \pm 0.19$, which is compatible with the Inoue et al. (2014) model or that obtained in other LyC leakers at similar redshifts and $f_{\text{esc,rel}}(\text{LyC})$ (e.g. Shapley et al. 2016) using Monte Carlo simulations of the IGM transmission, but it is larger than those obtained by, e.g. Steidel et al. (2018) and Fletcher et al. (2019), with mean values of $T(\text{IGM}) \simeq 0.3$.

Using equation (1) we infer $f_{\text{esc,rel}}(\text{LyC}) \sim 0.56$, with possible values ranging from 0.34 to 1. Note that we are considering only the uncertainties due to the IGM. A slightly older stellar population, for example a burst with $\simeq 6$ Myr, will produce less ionizing photons than the model we are considering (3 Myr) by a factor of $\simeq 2$. However, the 6 Myr burst model would require less extinction to explain the observed $\beta_{\text{UV}} = -2.05$, so that the $f_{900}(\text{obs})/f_{900}(\text{S99})$ would be roughly similar in both cases. Other sources of uncertainty can impact the inferred $f_{\text{esc,rel}}(\text{LyC})$, but are extremely difficult to quantify. These include the uncertainties due to the flux calibration at the blue edge of the OSIRIS spectrum, which has a total efficiency $\lesssim 5$ per cent at $\lambda < 4000$ Å only,<fnref n="5" /> or differential slit-losses between $\lambda_0 < 912$ Å and $\lambda_0 \simeq 1500$ Å, due to the effect of the atmospheric dispersion and the broadening of the point-spread function at blue wavelengths. Nevertheless, the uncertainties due to the IGM should be dominant. Given this, the absolute LyC escape fraction, $f_{\text{esc,abs}}(\text{LyC})$, defined as (e.g. Leitet et al. 2013)

$$f_{\text{esc,abs}}(\text{LyC}) = f_{\text{esc,rel}}(\text{LyC}) \times 10^{-0.4[E(B-V) \times k_{1500}]} \quad (2)$$

is $f_{\text{esc,abs}}(\text{LyC}) \approx 0.39$, with possible values between $\simeq 0.23$ and 0.69, allowed by the constraints on the IGM ($0.37 < T(\text{IGM}) < 1$) and $f_{\text{esc,rel}}(\text{LyC}) (<1)$. This a useful quantity to estimate the number of ionizing photons escaping from J0121+0025, $N_{\text{esc}}(\text{LyC})$, such as

$$N_{\text{esc}}(\text{LyC}) = f_{\text{esc,abs}}(\text{LyC}) \times N_{\text{int}}(\text{LyC}), \quad (3)$$

<footnote n="3">³The inferred $E(B-V)$ using Calzetti et al. (2000) curve is already low, so using other extinction curves (e.g. SMC) will have little impact in our results.</footnote>

<footnote n="4">⁴For an age of 3 Myr the predicted ionizing flux at 900 Å using S99 and BPASS (including binaries) are roughly the same, see e.g. Chisholm et al. (2019).</footnote>

<footnote n="5">⁵http://www.gtc.iac.es/instruments/osiris/#Spectroscopic_Photon_Detection_Efficiency</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of the spatial profiles of the emission below (solid blue) and above (dashed red) the rest-frame $\lambda = 912$ Å in the 2D spectrum of J0121+0025. These profiles have been extracted from the 2D spectrum over the rest-frame range 880–910 Å and 920–950 Å, respectively, and both have consistent spatial morphologies and are co-spatial.</figcaption></figure>

where $N_{\text{int}}(\text{LyC}) \simeq 1.4 \times 10^{55}$ s$^{-1}$ is the intrinsic ionizing photon production rate from the S99 model scaled to $M_{\text{UV}} = -24.1$ for a burst of 3 Myr and $Z_{\odot}/Z_{\odot} = 0.4$ (and a Salpeter 1955; IMF). We thus obtain an ionizing photon escape $N_{\text{esc}}(\text{LyC}) \approx 6 \times 10^{54}$ s$^{-1}$ (with possible values ranging from $[3 - 10] \times 10^{54}$ s$^{-1}$).

### 3.3.2 Possibility of foreground or AGN contamination

We now discuss the possibility that the emission detected below $\lambda_{\text{obs}} \simeq 3880$ Å is due to foreground contamination leading to a false-positive detection of LyC signal of J0121+0025 or arising from an AGN.

In Fig. 6 the spatial profiles of the emission below and above $\lambda_0 = 912$ Å are compared. These profiles have been extracted from the 2D spectrum over the rest-frame range 880–910 Å and 920–950 Å, respectively. Both profiles have consistent spatial morphologies with $\text{FWHM}_{(\lambda_0 < 912)} = 0.98'' \pm 0.11''$ and $\text{FWHM}_{(\lambda_0 > 912)} = 1.02'' \pm 0.02''$, which are compatible with being unresolved. Their centroids appear co-spatial, without any evidence for a spatial offset.

From the best seeing condition image R-band ($\text{FWHM} \simeq 0.5''$) we do not see any evidence for the presence of an additional source to J0121+0025, which is only barely resolved (Fig. 1). Given the $5\sigma$ depth of $\simeq 26.5$ of this image, a possible contaminant would be easily detected if it was spatially offset from J0121+0025 by $\gtrsim 0.3$ arcsec, given the observed magnitude $m = 24.18$ measured from the spectrum at $\lambda_0 < 912$ Å. In addition, one strong absorption line with a residual intensity compatible with zero is detected at $\lambda_{\text{obs}} \simeq 4740$ Å, that is associated with $\text{Ly}\alpha$ from a H I absorption system at $z = 2.898$ (#2 in Fig. 5). The non-detection of any flux at $\lambda_{\text{obs}} \simeq 4740$ Å below a $2\sigma$ level of $\simeq 1 \times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$, where the OSIRIS spectrum is $\simeq 4$ times more sensitive than at $\lambda_{\text{obs}} \simeq 3850$ Å, makes the presence of a contaminant very unlikely, because it requires an extremely low, and maybe unrealistic $\beta_{\text{UV}} < -3.2$ ($2\sigma$) source to explain such colour.

The presence of a relatively bright ($m = 24.18$), very compact ($r_{\text{eff}} \lesssim 0.1$ arcsec), quasi co-spatial with J0121+0025 ($\lesssim 0.2$ arcsec),

and very blue ($\beta_{\text{UV}} < -3.2$) lower-z interloper is still possible and cannot be completely ruled out, but it is highly unlikely.

We now discuss the possible contribution of an AGN to the LyC emission observed in J0121+0025. In Section 3.1, we have shown that the contribution of an AGN to the UV luminosity of J0121+0025 should be small, if present. More specifically, the AGN contribution should be at least $\lesssim 25$ per cent, otherwise photospheric absorption lines, which are intrinsically very weak, would not be detected in the OSIRIS spectrum. Considering the most extreme case, i.e. a contribution of 25 per cent to the UV flux density, the AGN would have a rest-frame 1500 Å flux density $f_{1500,\text{obs,AGN}} \simeq 2.2\mu\text{Jy}$, corresponding to $M_{\text{UV}} = -22.7$. Assuming that the observed $f_{900}(\text{obs}) = 0.781 \pm 0.099\mu\text{Jy}$ arises from such AGN, we obtain ($f_{900}/f_{1500}$)$_{\text{obs,AGN}} \simeq 0.36$, which is a factor of $\simeq 3-7$ larger than the typical values observed in type-I AGNs with LyC detection at similar redshift and luminosities ($f_{900}/f_{1500}$)(AGN) $\simeq 0.05 - 0.14$; Steidel, Pettini & Adelberger 2001; Micheva, Iwata & Inoue 2017) or in other bright QSOs (e.g. Cristiani et al. 2016), possibly corresponding to a non-physical LyC escape fraction. Therefore, it is unlikely that the LyC emission arises from a type-I AGN. The presence of an obscured type-II AGN in J0121+0025 is more difficult to constrain, but its contribution to the LyC emission can be neglected, as these sources are by definition very obscured at short wavelengths.

An additional piece of evidence that the emission detected below $\lambda_{\text{obs}} \simeq 3880$ Å is related with escape of ionizing photons comes from the intrinsic properties of J0121+0025. In fact, J0121+0025 could be identified as a strong LyC leaker candidate, even excluding the direct information about the LyC detection. From the spectrum, J0121+0025 shows very weak LIS lines, both in terms of $EW_0$ and the residual intensity. It has been shown, from observations and simulations, that the residual intensity of LIS lines correlates with $f_{\text{esc}}(\text{LyC})$ (e.g. Heckman et al. 2001; Alexandroff et al. 2015; Chisholm et al. 2018; Mauerhofer et al. 2021; and Saldana-Lopez in preparation). Using the prescription given in Chisholm et al. (2018), the predicted $f_{\text{esc,abs}}^{\text{pred}}(\text{LyC}) \approx 0.25$ from the residual intensity of the Si II line (see Section 3.2.4) is compatible, within the uncertainties, to that observed/inferred from the spectrum, $f_{\text{esc,abs}}(\text{LyC}) \approx 0.4$. In addition, the $\text{Ly}\alpha$ line shows a narrow profile with an intrinsic $\text{FWHM} = 350 \pm 40$ km s$^{-1}$ and with its peak closed to the systemic redshift, $v_{\text{peak}} \simeq 120 \pm 50$ km s$^{-1}$. The $\text{Ly}\alpha$ profiles in the confirmed LyC leakers analysed in Steidel et al. (2018) and Fletcher et al. (2019) have their peaks less redshifted than the LyC non-leakers (their fig. 26), suggesting low neutral gas column density where $\text{Ly}\alpha$ photons, and likely LyC, could escape more easily (Verhamme et al. 2015). Other observational signatures shared by LyC leakers are also present in J0121+0025, such as the strong P-Cygni profiles and broad He II emission (e.g. Vanzella et al. 2018; Rivera-Thorsen et al. 2019; Vanzella et al. 2020), low dust attenuation ($E(B-V) \simeq 0.04$), compact morphology ($r_{\text{eff}} \sim 1$ kpc) but large SFR (next section), and large SFR surface density (e.g. Izotov et al. 2016), or the evidence for strong outflows ($v_{\text{peak}}(\text{LIS}) \simeq -450$ km s$^{-1}$).

## 3.4 Multiwavelength SED fitting

Turning to the multiwavelength properties, J0121+0025 shows a flat ($F_\nu$) spectral energy distribution (Fig. 7) from optical to the mid-IR, that is broadly consistent with a young starburst.

We perform SED-fitting with CIGALE code (Burgarella, Buat & Iglesias-Páramo 2005; Boquien et al. 2019) using the photometry from G to Spitzer 4.5µm (Table 1), covering a rest-frame wavelength 0.11–1.2 µm. The SFH is modeled using two components: a young starburst with age $\le 10$ Myr allowed by the analysis on the UV
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Best-fitting model of the spectral energy distribution of J0121+0025 using CIGALE (Burgarella et al. 2005). The SED of J0121+0025 is dominated by a young and intense burst of star formation with SFR=981 ± 232 M☉ yr⁻¹, an age of 7 ± 2 Myr, and a stellar mass of log(M*,M☉) = 9.9 ± 0.1 (black line and orange circles). The fit uses photometry from G-band to IRAC 4.5 µm (blue squares; see Table 1). The red shaded region represents the upper limit SED of the old stellar component with an age of 500 Myr and log(M*,M☉)<10.4.</figcaption></figure>

spectral features (see Section 3.2.1 and Fig. 2), and an exponentially declining SFH with age of 500 Myr. Using simultaneously two SFH components allows us to probe the properties of the young stellar population (i) that likely dominates the SED, and investigate the presence of an underlying old stellar component (ii). We adopted the stellar population models from Bruzual & Charlot (2003), and assume a fixed metallicity Z/Z☉ = 0.4 based on our analysis in Section 3.2.1. The Calzetti et al. (2000) dust attenuation law is considered with E(B-V)* < 0.1, and assuming the ratio of the stellar and nebular E(B-V)/E(B-V)neb = 0.44. We also adopt $f_{esc,abs}$(LyC) ≈ 0.40.

Fig. 7 shows the best-fitting model. The emission of J0121+0025 is dominated by a young stellar burst in the whole spectral range covered by the imaging data. The burst is characterized by a 10 Myr-weighted SFR = 981 ± 232 M☉ yr⁻¹, an age of 7 ± 2 Myr, and a stellar mass of log(M*,M☉) = 9.9 ± 0.1. This yields a high specific SFR (sSFR = SFR/M*) of 98 ± 32 Gyr⁻¹, which is a factor 20 times higher than 10¹⁰ M☉ main-sequence star-forming galaxies at 3 < z < 4 (e.g. Tomczak et al. 2016). Table 2 summarizes the properties of J0121+0025. Both the age and the stellar mass of the burst are broadly consistent with the results inferred from the S99 best-fitting model using the UV wind lines (age of ~3 Myr and log(M*,M☉) = 9.8; Section 3.2.1).

On the other hand, the old stellar population is not well constrained. The best fitting gives a stellar mass for the old stellar component log(M*,M☉) = 9.8 ± 0.6. Fig. 7 shows the upper limit SED of the old stellar population (red), corresponding to log(M*,M☉)<10.4. Overall, the observed blue colour around the rest-frame ≃1.0 µm from the two IRAC/Spitzer channels I1-I2 = -0.40 ± 0.31 limits the presence of a log(M*,M☉)>10.4 old stellar component.

## 4 DISCUSSION

With $M_{UV} = -24.1 \pm 0.1$ and $f_{esc,abs}$(LyC) ≈ 0.40, J0121+0025 is not only one of the most UV-luminous star-forming galaxies ever discovered, but also the brightest LyC leaker known among star-forming galaxies. We now discuss the implications of such discovery.

<table><caption>**Table 2.** Properties of J0121+0025.</caption><thead><tr><th></th><th>Value</th><th>Uncertainty</th></tr></thead><tbody><tr><td>R.A. (J2000)</td><td>01:21:56.09</td><td>0.1 arcsec</td></tr><tr><td>Dec. (J2000)</td><td>+00:25:20.30</td><td>0.1 arcsec</td></tr><tr><td>z<sub>sys</sub></td><td>3.244</td><td>0.001</td></tr><tr><td>M<sub>UV</sub> (AB)</td><td>-24.11</td><td>0.1</td></tr><tr><td>log(L[Lyα/erg s<sup>-1</sup>])</td><td>43.8</td><td>0.1</td></tr><tr><td>r<sub>eff</sub> (kpc)</td><td>&lt;1.0</td><td>—</td></tr><tr><td>Age (Myr)</td><td>3-7<sup>a)</sup></td><td>&lt;10<sup>2</sup><sup>a)</sup></td></tr><tr><td>Z<sub>*</sub>/Z<sub>☉</sub></td><td>0.4</td><td>[0.2-1.0]</td></tr><tr><td>E(B-V)*</td><td>0.04</td><td>0.02</td></tr><tr><td>SFR (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>981<sup>b)</sup></td><td>232<sup>b)</sup></td></tr><tr><td>log(M<sup>burst</sup>/M<sub>☉</sub>)</td><td>9.9</td><td>0.1</td></tr><tr><td>sSFR (Gyr<sup>-1</sup>)</td><td>98</td><td>32</td></tr><tr><td>ΣSFR (M<sub>☉</sub> yr<sup>-1</sup> kpc<sup>-2</sup>)</td><td>&gt;157<sup>c)</sup></td><td>—</td></tr><tr><td>f<sub>v</sub>(900 Å) (μJy)</td><td>0.781</td><td>0.099</td></tr><tr><td>f<sub>v</sub>(900 Å)/f<sub>v</sub>(1500 Å)</td><td>0.088</td><td>0.012</td></tr><tr><td>T(IGM)</td><td>0.60</td><td>0.19</td></tr><tr><td>f<sub>esc,rel</sub> (LyC)</td><td>0.56</td><td>[0.34-1.0]</td></tr><tr><td>f<sub>esc,abs</sub> (LyC)</td><td>0.39</td><td>[0.23-0.69]</td></tr></tbody></table>

(a) Age of the young stellar population obtained using UV wind lines (≈3 Myr) and the best-fitting model of the SED using CIGALE (7 ± 2 Myr);

(b) 10 Myr-weighted SFR obtained from the best-fitting SED; (c) considering $r_{eff} < 1$ kpc.

## 4.1 The brightest LyC emitter known

J0121+0025 meets the two necessary conditions to be a strong LyC leaker: it is very efficient at producing hydrogen-ionizing photons and the properties of its ISM are favourable for their escape.

The remarkably strong P-Cygni profiles in the wind lines OVI, Nv, and CIV seen in the spectrum of J0121+0025 (see Figs 2 and 3) indicates unambiguously a very young age of the burst (≈3 Myr) and the presence of a large number of O-type stars, the main-sequence stars hot enough to generate a significant number of ionizing photons. Strong P-Cygni profiles in these lines are ubiquitous in the spectra of LyC leakers, both at low-z (e.g. Borthakur et al. 2014; Izotov et al. 2016, 2018a, b) and moderately high-z (e.g. Vanzella et al. 2018, 2020; Rivera-Thorsen et al. 2019), but appear weak or absent in composite spectra of more typical LBGs or LAEs (e.g. Shapley et al. 2003; Du et al. 2018; Nakajima et al. 2018; Feltre et al. 2020; Marques-Chaves et al. 2020a). In contrast to the burstiness nature of J0121+0025, and likely other LyC leakers, smooth/continuous or declining star formation histories with ages from several tens to hundreds Myr could explain the weakness of these profiles in the spectra of typical LBGs/LAEs, which is in line with the SFHs usually inferred or assumed for them (e.g. Kornei et al. 2010; de Barros, Schaerer & Stark 2014; Arrabal Haro et al. 2020). As a result, the amplitude and the age of the burst in J0121+0025 produces a large number of ionizing photons, $N_{int}$(LyC) ≃ 1.4 × 10⁵⁵ s⁻¹ (see Section 3.3), which is a factor of ~10–30 larger than that expected to be produced by a $M_{UV}^* = -20.97$ galaxy (e.g. Reddy & Steidel 2009), assuming a continuous star formation with 100 Myr age and the same metallicity.

The presence of strong P-Cygni profiles could be a potential indicator of LyC leakage, as discussed in previous works (e.g. Izotov et al. 2018b; Chisholm et al. 2019), at least for moderately metal-rich galaxies (because the strength of the P-Cygni is also metallicity dependent; see Izotov et al. 2021). Note, however, that the presence of such wind lines indicates a high production efficiency of LyC photons, not necessarily their escape. Nevertheless, feedback from the strong winds of massive stars together with SN explosions
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Relation between the observed ratio of the ionizing to non-ionizing flux density (in $F_{\nu}$ units, left) and the observed flux density at $\sim 900$ Å (in $\mu$Jy, right) versus the UV absolute magnitude. J0121+0025 is represented with a blue circle. For comparison, we also show the other $z \gtrsim 3$ star-forming galaxies known with significant detection of LyC radiation (grey squares; de Barros et al. 2016; Shapley et al. 2016; Vanzella et al. 2016, 2018; Steidel et al. 2018; Fletcher et al. 2019; Ji et al. 2020). Empty diamonds mark the position of the estimates of several composites from Pahl et al. (2021) in different bins of $M_{\text{UV}}$.</figcaption></figure>

expected in the early phase of a burst could play a major role in shaping the ISM, creating cavities of ionized gas where ionizing photons could escape more efficiently (e.g. Heckman et al. 2011; Trebitsch et al. 2017). This might be the case of J0121+0025. The particular conditions of the ISM in J0121+0025 are in fact favourable for the escape of LyC radiation. LyC photons are easily absorbed by dust and neutral gas and these sources of opacity are apparently weak in J0121+0025, at least from our line of sight. The inferred $\beta_{\text{UV}} = -2.05$ from the spectrum is compatible with low attenuation by dust, corresponding to $E(B-V) = 0.04$. In addition, ISM absorption lines are very weak, with $EW_0 \lesssim 1$ Å, and show residual intensities of $I/I_0 \simeq 0.8$ for Si II 1260 Å and C II 1334 Å, even though they are likely saturated (Section 3.2.4). These findings suggest a clumpy ISM with a non-unity covering fraction ($C_f$ (H1) $\simeq 0.55$; Gazagnes et al. 2018). The detection of blueshifted profiles in LIS lines in J0121+0025 ($\simeq -450$ km s$^{-1}$; see Section 3.2.4 and Fig. 3) supports the presence of such strong outflows.

We now compare in Fig. 8 the LyC properties of J0121+0025 with those from other confirmed LyC leakers at $z \gtrsim 3$. The comparison sample consists of $\sim 20$ sources with significant detection of LyC radiation taken from de Barros et al. (2016), Shapley et al. (2016), Vanzella et al. (2016, 2018), Steidel et al. (2018), Fletcher et al. (2019), and Ji et al. (2020). The left-hand panel compares the observed ratio of the ionizing to non-ionizing flux density, $(f_{900}/f_{1500})_{\text{obs}}$, versus the UV absolute magnitude. We use $(f_{900}/f_{1500})_{\text{obs}}$ because it is model-independent and does not rely on assumptions about the properties of the underlying stellar population, such as SFHs, age, and metallicity. Furthermore we do not correct these values for the IGM absorption, because it is highly uncertainty and it is not always available. From this figure, our results indicate that a significant fraction of LyC photons can escape in sources with a wide range of UV luminosity, from UV faint ($M_{\text{UV}} \gtrsim -18.7$, Fletcher et al. 2019) to extremely UV luminous ($M_{\text{UV}} = -24.1$, this work), and therefore are not restricted to the faintest ones as previously thought (Steidel et al. 2018; Pahl et al. 2021, see also Bian & Fan 2020).

The importance of detecting LyC emission in such a luminous source like J0121+0025 is highlighted in the right-hand panel of

Fig. 8. Here, we compare the observed flux density at $\sim 900$ Å, $f_{900}$. Again, no IGM corrections have been applied nor corrections for the luminosity distance, although all but two LyC leakers (Vanzella et al. 2018; Ji et al. 2020) are roughly at the same redshift of J0121+0025. Therefore, this comparison should be treated with care, serving for illustrative purposes only. The combination of the extreme luminosity of the starburst ($M_{\text{UV}} = -24.1$) and the large $f_{\text{esc}}$ (LyC) makes J0121+0025 the brightest and the most powerful LyC emitter known among the star-forming galaxy population. The observed LyC flux in J0121+0025 is in fact comparable to the sum of the LyC flux of all star-forming galaxies with LyC leakage known at these redshifts. Only the highly magnified *Sunburst Arc* ($\mu \sim 30-100 \times$; Pignataro et al. 2021) or bright QSOs have comparable observed LyC fluxes ($\approx 1$ $\mu$Jy; Steidel et al. 2001; Lusso et al. 2015; Rivera-Thorsen et al. 2019).

The discovery of such a powerful LyC emitter raises now the question regarding the role of UV-luminous star-forming galaxies to the cosmic reionization (e.g. Sharma et al. 2016, Naidu et al. 2020). Answering this question is, however, out of scope of this work, as it requires the knowledge of two fundamental properties that are highly uncertain: the volume density of such luminous sources at $z \gtrsim 7$ and the physical properties connecting LyC leakage. Nevertheless, we can place rough constraints on the UV ionizing background at $z \sim 3$. To do so, we assume that luminous sources ($M_{\text{UV}} < -22$) share the same properties as J0121+0025, i.e. $f_{\text{esc,abs}}$ (LyC) $\sim 0.40$ and $\log(\xi_{\text{ion}}) = 25.2$. The co-moving production rate of hydrogen-ionizing photons ($\dot{N}_{\text{ion}}$) is given by

$$ \dot{N}_{\text{ion}} = f_{\text{esc,abs}}(\text{LyC}) \xi_{\text{ion}} \rho_{\text{UV}} [\text{s}^{-1}\text{Mpc}^{-3}], \quad (4) $$

where $\rho_{\text{UV}}$ is the dust-corrected UV luminosity density. For $\rho_{\text{UV}}$, we integrate the UV luminosity function of Reddy & Steidel (2009) down to $-22$ AB (or $\gtrsim 3L_{\text{UV}}^*$) and assume $E(B-V) = 0.04$. This yields a $\dot{N}_{\text{ion}} \sim 10^{49.8} \text{ s}^{-1} \text{ Mpc}^{-3}$, which is fairly lower than that provided by QSOs at $z \sim 3$ ($\dot{N}_{\text{ion}}$ (QSO) $\sim 10^{50.5} - 10^{51.0} \text{ s}^{-1} \text{ Mpc}^{-3}$; e.g. Becker & Bolton 2013, Cristiani et al. 2016). However, the situation may differ at very high-$z$. Recent studies have found remarkably bright galaxies at $z \gtrsim 7$ that are in excess compared
to the generally observed Schechter function of luminosity functions (e.g. Bowler et al. 2014; Oesch et al. 2014; Bowler et al. 2015; Ono et al. 2018; Stefanon et al. 2019). Assuming the double power-law luminosity function at $z \sim 6$ of Bowler et al. (2015), we find $\dot{N}_{\text{ion}} \sim 10^{48.9} s^{-1} M_{\odot} s^{-3}$ for star-forming galaxies brighter than $-22$ AB that is comparable to that inferred in QSOs at these redshifts ($\dot{N}_{\text{ion}}(\text{QSO}) \sim 10^{48.8} s^{-1} M_{\odot} s^{-3}$; Matsuoka et al. 2018c). On the other hand, it is not clear if the properties connecting the LyC leakage in J0121+0025, and therefore the large $f_{\text{esc}}$ (LyC), can be applied to other bright star-forming galaxies (see e.g. Harikane et al. 2020). This will be discussed in the next section.

## 4.2 Understanding UV-luminous star-forming galaxies: diverse properties and insights for LyC leakage

J0121+0025 shows intriguing properties that differ from those expected in bright star-forming galaxies. Here we discuss some of these properties and compare them with other UV-luminous star-forming galaxies, with particular emphasis to the properties that could be related with the LyC leakage.

A well-established trend relating $M_{\text{UV}}$ and the strength of Ly$\alpha$ and ISM lines has been found in previous works (e.g. Shapley et al. 2003; Vanzella et al. 2009; Trainor et al. 2015; Du et al. 2018), for which Ly$\alpha$ is found to be weak and LIS lines strong in UV-bright sources. While this trend has been established with statistical significance for galaxies with $M_{\text{UV}}$ between -18 and -21, a few other known $M_{\text{UV}} \simeq -23$ LBGs (Dessauges-Zavadsky et al. 2010; Lee et al. 2013; Marques-Chaves et al. 2018; Harikane et al. 2020) show the same trend, presenting a completely damped Ly$\alpha$ absorption and strong LIS lines ($EW_0 \simeq 2-4$ Å) in their spectra, which is compatible with a large column density of neutral gas ($N(H\text{I}) > 10^{20} \text{ cm}^{-2}$). However, such trend is not seen in J0121+0025 as it shows weak LIS lines and a relatively strong Ly$\alpha$ line. Another interesting difference is that P-Cygni in wind lines, in particular in N v that is very sensitive to the age, are weak or absent in these $M_{\text{UV}} \simeq -23$ LBGs (e.g. Dessauges-Zavadsky et al. 2010; Marques-Chaves et al. 2018; Harikane et al. 2020), which could indicate that their UV continua is not dominated by O-type stars ($\le 10$ Myr), at least compared to J0121+0025, but by older and less luminous stars (e.g. B-type stars). In addition, vigorous star-forming galaxies are also found to be more dusty, because the production of dust is tightly linked with star formation. However, this does not apply for J0121+0025. In fact, the derived SFR = $981 \pm 232 M_{\odot} \text{yr}^{-1}$ of J0121+0025 is comparable, in absolute terms, to those inferred in very dusty far-IR bright systems, but as opposed to them, only a small fraction of the total SFR of J0121+0025 is obscured ($\simeq 30$ per cent).

On the other hand, the properties of J0121+0025 are remarkably similar to those observed in a few, very young, and also luminous starbursts, like the recently discovered extremely luminous starburst galaxy BOSS-EUVLG1 at $z = 2.469$ ($M_{\text{UV}} \simeq -24.4$; Marques-Chaves et al. 2020b), or the less, but still luminous LyC leaker Ion3 galaxy at $z \simeq 4.0$ ($M_{\text{UV}} \simeq -22.2$; Vanzella et al. 2018). The spectra of these galaxies show very strong P-Cygni profiles in wind lines (O vi, N v, and C iv), intense nebular emission in rest-frame optical lines ($EW_0(H\alpha) \gtrsim 700$ Å) and SEDs that are compatible with a very young and intense starburst of a few Myr ($\le 10$ Myr). As J0121+0025, these galaxies show very weak ISM absorption lines ($EW_0 < 1$ Å), strong Ly$\alpha$ emission $EW_0 > 20$ Å, and are almost unobscured ($\beta_{\text{UV}} < -2.2$; Vanzella et al. 2018; Marques-Chaves et al. 2020b). Other interesting properties shared in these galaxies are the high specific SFR (sSFR $\simeq 90-100$ Gyr$^{-1}$), compact morphologies ($r_{\text{eff}} \sim 1$ kpc or less), and therefore high star-formation rate surface

density ($\Sigma$SFR $\gtrsim 100 M_{\odot} \text{yr}^{-1} \text{ kpc}^{-2}$), properties that are common in other LyC leakers (e.g. Izotov et al. 2016, 2018b), and are thought to play a key role in transforming the ISM structure (feedback).

Understanding such diversity in the properties of UV-luminous galaxies is challenging, and possibly premature for now given the lack of statistics. We note that only a handful of luminous sources are known as bright as $M_{\text{UV}} \sim -23$ and only two brighter than $M_{\text{UV}} \sim -24$. J0121+0025 (this work) and BOSS-EUVLG1 (Marques-Chaves et al. 2020b). Nevertheless, the differences are already notorious, in particular those thought to be closely related to the LyC leakage.

A possible explanation to these differences may be simply related to different stages in the evolution of galaxies. J0121+0025, BOSS-EUVLG1 or Ion3 could represent a vigorous starburst seen at very initial stages ($<10$ Myr), when the starlight is dominated by young stars and before dust-attenuation become efficient, enhancing the UV luminosity. The extreme feedback expected in such early, intense, and likely short-lived phase could eject the gas and dust from star-forming regions, allowing the escape of LyC photons (e.g. Sharma et al. 2017; Trebitsch et al. 2017; Arata et al. 2019). In fact, a powerful ionized gas outflow has been detected in BOSS-EUVLG1 with a $\log(M_{\text{out}}/M_{\odot}) \simeq 8$ and an outflowing velocity $v_{\text{out}} \simeq 600 \text{ km s}^{-1}$ (Álvarez-Márquez et al. 2021). On the other hand, if the SFR drops significantly at later stages, the UV luminosity will be still high for a considerable period of time ($\sim 50$ Myr) due to the contribution of B-type stars, but in such case SN and stellar feedback could be less effective in clearing sightlines, and neutral gas and dust would cover star-forming regions absorbing the LyC radiation. In both cases, the galaxy will appear bright in the UV, but some spectral features would appear dramatically different. Alternatively, the differences seen in the spectra of UV-luminous galaxies could arise by a non-homogeneous distribution of gas and dust in these UV-luminous galaxies, for which some of them would have star-forming regions cleared by dust and gas from a favourable sightline (J0121+0025, BOSS-EUVLG1 or Ion3), while others would not (e.g. Harikane et al. 2020).

Independent of the scenario invoked to explain these differences, the discovery of J0121+0025 (and Ion3; Vanzella et al. 2018) indicates that, at least some UV-luminous star-forming galaxies can be strong LyC emitters, which contradicts recent findings from Harikane et al. (2020). We note that the sample of $M_{\text{UV}} \simeq -23$ LBGs at $z \sim 6$ analysed in Harikane et al. (2020) was previously selected to have $\log(L[\text{Ly}\alpha/\text{erg s}^{-1}]) < 43.0$ (Matsuoka et al. 2016, 2018a, b, 2019) to avoid a possible contamination of an AGN. However, such threshold is conservative (see e.g. Ouchi et al. 2009; Sobral et al. 2015; Matthee et al. 2017; for other luminous LAEs) and can yield to a selection bias towards high $N(H\text{I})$ and/or declining SFHs. In fact, the composite spectrum presented in Harikane et al. (2020) shows a completely damped Ly$\alpha$ absorption, which is compatible with a large column density of neutral gas, and naturally explain the strong LIS lines observed in the spectrum. On the other hand, other UV-bright sources identified as AGNs by Matsuoka et al. (2016, 2018a, b, 2019), based only on $\log(L[\text{Ly}\alpha/\text{erg s}^{-1}]) > 43.0$, show intense and narrow Ly$\alpha$ emission (as narrow as $<230$ km s$^{-1}$), very weak LIS lines, and, interestingly, evidence of strong P-Cygni in N v (see fig. 9 in Matsuoka et al. 2019). The authors argue that such a P-Cygni profile could arise from a weak BAL-QSO. However, as shown here in this work, similar profiles could be naturally explained by a very young and hot stellar population that could enhance the UV and Ly$\alpha$ luminosities in these sources, similar to what is happening in J0121+0025 and BOSS-EUVLG1 (Marques-Chaves et al. 2020b).
In closing, it is clear that properties of UV-luminous star-forming galaxies are still not very well understood and must be investigated in more detail with a large statistical sample. In a future work, we will present a large sample of other ~70 extremely UV-luminous star-forming galaxies discovered within the eBOSS survey (R. Marques-Chaves, in preparation), hoping to answer some of these important questions.

## 5 SUMMARY AND CONCLUSION

This work reports the discovery of J0121+0025 at $z = 3.244 \pm 0.001$, extremely luminous in the UV ($M_{UV} \simeq -24.11$, AB) and Ly$\alpha$ line ($\log[L_{Ly\alpha}/\text{erg s}^{-1}] = 43.8$) star-forming galaxy with copious emission in the LyC spectral range ($\lambda_0 < 912$ Å). J0121+0025 is a compact starburst with $r_{\text{eff}} = 1 \pm 0.5$ kpc that is only barely resolved in very good seeing conditions ground-based imaging. The optical to mid-IR photometry is dominated by the emission of a vigorous starburst, with $\log(M_{\text{burst}}^*/M_{\odot}) = 9.9 \pm 0.1$ and a 10 Myr-weighted SFR = $981 \pm 232 M_{\odot} \text{yr}^{-1}$. This yields a high specific star-formation rate $\text{sSFR} = 98 \pm 32 \text{ Gyr}^{-1}$ and a SFR density $\Sigma\text{SFR} > 157 M_{\odot} \text{yr}^{-1} \text{kpc}^{-2}$ (considering $r_{\text{eff}} < 1$ kpc). The high SNR OSIRIS/GTC spectrum of J0121+0025 reveals strong P-Cygni in wind lines of O VI, N V, and C IV, which are well reproduced by a starburst model with an extremely young age of $\simeq 3$ Myr, which is roughly consistent with the age derived from the multiwavelength SED ($7 \pm 2$ Myr). The spectrum shows a rest-frame UV slope $\beta_{\text{UV}} = -2.05 \pm 0.10$, consistent with low dust attenuation $E(B-V)_* = 0.04 \pm 0.02$. It also shows other features characteristic of star-forming galaxies, such as stellar absorption originated in the photospheres of hot stars, for which a significant contribution of an AGN to the luminosity is ruled out. The Ly$\alpha$ is moderately strong ($EW_0[\text{Ly}\alpha] = 14 \pm 3$ Å) and shows a narrow profile (FWHM $\simeq 350$ km s$^{-1}$) with its peak redshifted, but close to the systemic velocity, by $\simeq 120$ km s$^{-1}$. Low-ionization ISM lines are also detected, but appear much weaker when compared to those observed in typical LBGs. Both the weakness ($EW_0[\text{LIS}] \simeq 1$ Å) and the large residual intensity ($I/(I_0) \simeq 0.8$) suggest a clumpy geometry with a non-unity covering fraction or a highly ionized ISM, for which a significant fraction of ionizing photons could escape. LyC radiation is detected with a significance of $\simeq 7.9$ in the OSIRIS spectrum, corresponding to a flux density $f_{\text{000}} = 0.781 \pm 0.099 \text{ µJy}$ and an ionizing to non-ionizing flux density ($f_{\text{000}}/f_{\text{1500}}_{\text{obs}} = 0.09 \pm 0.01$). The contribution of a foreground or AGN contamination to the LyC signal is discussed in detail, and although it cannot be completely ruled out, it is very unlikely. This makes J0121+0025 the most powerful LyC emitter known among the star-forming galaxy population. Our results indicate that at least some UV-luminous star-forming galaxies are strong LyC leakers, bringing new insights to the discussion of the role of luminous and very young starbursts to the cosmic reionization.

## ACKNOWLEDGEMENTS

The authors thank the referee, Eros Vanzella, for useful comments that greatly improved the clarity of this work. Based on observations made with the Gran Telescopio Canarias (GTC) installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. JAM, LC, and IPF acknowledge support from the Spanish State Research Agency (AEI) under grant numbers ESP2015-65597-C4-4-R, ESP2017-86852-C4-2-R, RyC-2015-18078, PGC2018-094975-B-C22, and MDM-2017-0737 Unidad de Excelencia 'María de Maeztu'- Centro de

Astrobiología (CSIC-INTA). ASL acknowledge support from Swiss National Science Foundation.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Abolfathi B. et al., 2018, ApJS, 235, 42

* Aihara H. et al., 2019, PASJ, 71, 114

* Alexandroff R. M., Heckman T. M., Borthakur S., Overzier R., Leitherer C., 2015, ApJ, 810, 104

* Álvarez-Márquez J., Marques-Chaves R., Colina L., Pérez-Fournon I., 2021, A&A, 647, A133

* Appenzeller I., Stahl O., Tapken C., Mehlert D., Noll S., 2005, A&A, 435, 465

* Arata S., Yajima H., Nagamine K., Li Y., Khochfar S., 2019, MNRAS, 488, 2629

* Arrabal Haro P. et al., 2020, MNRAS, 495, 1807

* Assef R. J. et al., 2013, ApJ, 772, 26

* Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5

* Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023

* Bentz M. C., Osmer P. S., Weinberg D. H., 2004, ApJ, 600, L19

* Bian F., Fan X., 2020, MNRAS, 493, L65

* Bian F., Fan X., McGreer I., Cai Z., Jiang L., 2017, ApJ, 837, L12

* Boquien M., Burgarella D., Roehlly Y., Buat V., Ciesla L., Corre D., Inoue A. K., Salas H., 2019, A&A, 622, A103

* Borthakur S., Heckman T. M., Leitherer C., Overzier R. A., 2014, Science, 346, 216

* Bowler R. A. A. et al., 2014, MNRAS, 440, 2810

* Bowler R. A. A. et al., 2015, MNRAS, 452, 1817

* Brinchmann J., Pettini M., Charlot S., 2008, MNRAS, 385, 769

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

* Burgarella D., Buat V., Iglesias-Páramo J., 2005, MNRAS, 360, 1413

* Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storch-Bergmann T., 2000, ApJ, 533, 682

* Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245

* Chandar R., Leitherer C., Tremonti C. A., 2004, ApJ, 604, 153

* Chisholm J. et al., 2018, A&A, 616, A30

* Chisholm J., Rigby J. R., Bayliss M., Berg D. A., Dahle H., Gladders M., Sharon K., 2019, ApJ, 882, 182

* Cristiani S., Serrano L. M., Fontanot F., Vanzella E., Monaco P., 2016, MNRAS, 462, 2478

* Crowther P. A. et al., 2016, MNRAS, 458, 624

* de Barros S. et al., 2016, A&A, 585, A51

* de Barros S., Schaerer D., Stark D. P., 2014, A&A, 563, A81

* Dessauges-Zavadsky M., D'Odorico S., Schaerer D., Modigliani A., Tapken C., Vernet J., 2010, A&A, 510, A26

* Dessauges-Zavadsky M., Prochaska J. X., D'Odorico S., Calura F., Matteucci F., 2006, A&A, 445, 93

* Du X. et al., 2018, ApJ, 860, 75

* Eisenstein D. J. et al., 2011, AJ, 142, 72

* Feltre A. et al., 2020, A&A, 641, A118

* Finkelstein S. L. et al., 2011, ApJ, 729, 140

* Finkelstein S. L. et al., 2019, ApJ, 879, 36

* Fletcher T. J., Tang M., Robertson B. E., Nakajima K., Ellis R. S., Stark D. P., Inoue A., 2019, ApJ, 878, 87

* Gawiser E. et al., 2007, ApJ, 671, 278

* Gazagnes S., Chisholm J., Schaerer D., Verhamme A., Rigby J. R., Bayliss M., 2018, A&A, 616, A29

* Geach J. E. et al., 2017, ApJS, 231, 7

* Grimm H. J., Gilfanov M., Sunyaev R., 2003, MNRAS, 339, 793

* Gwyn S. D. J., 2008, PASP, 120, 212

* Hainline K. N., Shapley A. E., Greene J. E., Steidel C. C., 2011, ApJ, 733, 31

* Harikane Y., Laporte N., Ellis R. S., Matsuoka Y., 2020, ApJ, 902, 117
Heckman T. M. et al., 2011, ApJ, 730, 5

Heckman T. M., Sembach K. R., Meurer G. R., Leitherer C., Calzetti D.,
Martin C. L., 2001, ApJ, 558, 56

Inoue A. K., Iwata I., 2008, MNRAS, 387, 1681

Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805

Izotov Y. I., Schaerer D., Thuan T. X., Worseck G., Guseva N. G., Orlitová I., Verhamme A., 2016, MNRAS, 461, 3683

Izotov Y. I., Schaerer D., Worseck G., Guseva N. G., Thuan T. X., Verhamme A., Orlitová I., Fricke K. J., 2018a, MNRAS, 474, 4514

Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Chisholm J., Thuan T. X., Fricke K. J., Verhamme A., 2021, MNRAS, 503, 1734

Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Thuan T. X., Fricke Verhamme A., Orlitová I., 2018b, MNRAS, 478, 4851

Ji Z. et al., 2020, ApJ, 888, 109

Jiang L. et al., 2021, Nat. Astron., 5, 256

Kennicutt R. C. Jr, 1998, ARA&A, 36, 189

Kojima T., Ouchi M., Nakajima K., Shibuya T., Harikane Y., Ono Y., 2017, PASJ, 69, 44

Kornei K. A., Shapley A. E., Erb D. K., Steidel C. C., Reddy N. A., Pettini M., Bogosavljević M., 2010, ApJ, 711, 693

Lee K.-S., Dey A., Cooper M. C., Reddy N., Jannuzi B. T., 2013, ApJ, 771, 25

Leitet E., Bergvall N., Hayes M., Linné S., Zackrisson E., 2013, A&A, 553, A106

Leitherer C. et al., 1999, ApJS, 123, 3

Leitherer C., Byler N., Lee J. C., Levesque E. M., 2018, ApJ, 865, 55

Leitherer C., Hernandez S., Lee J. C., Oey M. S., 2016, ApJ, 823, 64

Leitherer C., Leão J. R. S., Heckman T. M., Lennon D. J., Pettini M., Robert C., 2001, ApJ, 550, 724

Lusso E., Worseck G., Hennawi J. F., Prochaska J. X., Vignali C., Stern J., O'Meara J. M., 2015, MNRAS, 449, 4204

Marques-Chaves R. et al., 2018, ApJ, 854, 151

Marques-Chaves R. et al., 2020a, MNRAS, 492, 1257

Marques-Chaves R. et al., 2020b, MNRAS, 499, L105

Matsuoka Y. et al., 2016, ApJ, 828, 26

Matsuoka Y. et al., 2018a, PASJ, 70, S35

Matsuoka Y. et al., 2018b, ApJS, 237, 5

Matsuoka Y. et al., 2018c, ApJ, 869, 150

Matsuoka Y. et al., 2019, ApJ, 883, 183

Matthee J., Sobral D., Darvish B., Santos S., Mobasher B., Paulino-Afonso A., Röttgering H., Alegre L., 2017, MNRAS, 472, 772

Mauerhofer V., Verhamme A., Blaizot J., Garel T., Kimm T., Michel-Dansac L., Rosdahl J., 2021, A&A, 646, A80

Micheva G., Iwata I., Inoue A. K., 2017, MNRAS, 465, 302

Naidu R. P., Tacchella S., Mason C. A., Bose S., Oesch P. A., Conroy C., 2020, ApJ, 892, 109

Nakajima K. et al., 2012, ApJ, 745, 12

Nakajima K., Fletcher T., Ellis R. S., Robertson B. E., Iwata I., 2018, MNRAS, 477, 2098

Nakajima K., Ouchi M., Shimasaku K., Hashimoto T., Ono Y., Lee J. C., 2013, ApJ, 769, 3

Oesch P. A. et al., 2014, ApJ, 786, 108

Oesch P. A. et al., 2016, ApJ, 819, 129

Ono Y. et al., 2010, MNRAS, 402, 1580

Ono Y. et al., 2018, PASJ, 70, S10

Ouchi M. et al., 2008, ApJS, 176, 301

Ouchi M. et al., 2009, ApJ, 696, 1164

Pahl A. J., Shapley A., Steidel C. C., Chen Y., Reddy N. A., 2021, MNRAS, 505, 2447

Papovich C. et al., 2016, ApJS, 224, 28

Pâris I. et al., 2018, A&A, 613, A51

Peng C. Y., Ho L. C., Impey C. D., Rix H.-W., 2002, AJ, 124, 266

Pignataro G. V. et al., 2021, preprint (arXiv:2106.10286)

Reddy N. A., Steidel C. C., 2009, ApJ, 692, 778

Reddy N. A., Steidel C. C., Pettini M., Bogosavljević M., 2016, ApJ, 828, 107

Rigby J. R. et al., 2018, ApJ, 853, 87

Rivera-Thorsen T. E. et al., 2017, A&A, 608, L4

Rivera-Thorsen T. E. et al., 2019, Science, 366, 738

Rix S. A., Pettini M., Leitherer C., Bresolin F., Kudritzki R.-P., Steidel C. C., 2004, ApJ, 615, 98

Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ, 802, L19

Ross N. P. et al., 2012, ApJS, 199, 3

Saha K. et al., 2020, Nat. Astron., 4, 1185

Salpeter E. E., 1955, ApJ, 121, 161

Santos S. et al., 2020, MNRAS, 493, 141

Schaerer D., 2003, A&A, 397, 527

Schlafly E. F., Finkbeiner D. P., 2011, ApJ, 737, 103

Senchyna P., Stark D. P., Charlot S., Chevallard J., Bruzual G., Vidal-García A., 2021, MNRAS, 503, 6112

Shapley A. E., Steidel C. C., Adelberger K. L., Dickinson M., Giavalisco M., Pettini M., 2001, ApJ, 562, 95

Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., 2003, ApJ, 588, 65

Shapley A. E., Steidel C. C., Strom A. L., Bogosavljević M., Reddy N. A., Siana B., Mostardi R. E., Rudie G. C., 2016, ApJ, 826, L24

Sharma M., Theuns T., Frenk C., Bower R., Crain R., Schaller M., Schaye J., 2016, MNRAS, 458, L94

Sharma M., Theuns T., Frenk C., Bower R., Crain R., Schaller M., Schaye J., 2017, MNRAS, 468, L2176

Sobral D. et al., 2018, MNRAS, 477, 2817

Sobral D., Matthee J., Darvish B., Schaerer D., Mobasher B., Röttgering H. J. A., Santos S., Hemmati S., 2015, ApJ, 808, 139

Stanway E. R., Eldridge J. J., Becker G. D., 2016, MNRAS, 456, 485

Stefanon M. et al., 2019, ApJ, 883, 99

Steidel C. C., Bogosavljević M., Shapley A. E., Reddy N. A., Rudie G. C., Pettini M., Trainor R. F., Strom A. L., 2018, ApJ, 869, 123

Steidel C. C., Pettini M., Adelberger K. L., 2001, ApJ, 546, 665

Steidel C. C., Strom A. L., Pettini M., Rudie G. C., Reddy N. A., Trainor R. F., 2016, ApJ, 826, 159

Tomczak A. R. et al., 2016, ApJ, 817, 118

Trainor R. F., Steidel C. C., Strom A. L., Rudie G. C., 2015, ApJ, 809, 89

Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224

Vanzella E. et al., 2009, ApJ, 695, 1163

Vanzella E. et al., 2012, ApJ, 751, 70

Vanzella E. et al., 2016, ApJ, 825, 41

Vanzella E. et al., 2017, MNRAS, 467, 4304

Vanzella E. et al., 2018, MNRAS, 476, L15

Vanzella E. et al., 2020, MNRAS, 491, 1093

Vanzella E. et al., 2021, preprint (arXiv:2106.10280)

Verhamme A., Orlitová I., Schaerer D., Hayes M., 2015, A&A, 578, A7

Walborn N. R. et al., 2010, AJ, 139, 1283

This paper has been typeset from a `TEX/LATEX` file prepared by the author.