# Evolution of the star formation rate and ΣSFR of galaxies at cosmic morning (4 < z < 10)

A. Calabrò, L. Pentericci, P. Santini, A. Ferrara, M. Llerena, S. Mascia, L.

Napolitano, L. Y. A. Yung, L. Bisigello, M. Castellano, et al.

► To cite this version:

A. Calabrò, L. Pentericci, P. Santini, A. Ferrara, M. Llerena, et al., Evolution of the star formation rate and ΣSFR of galaxies at cosmic morning (4 < z < 10). Astronomy & Astrophysics - A&A, 2024, 690, 10.1051/0004-6361/20249768. insu-04846495

HAL Id: insu-04846495

https://insu.hal.science/insu-04846495v1

Submitted on 18 Dec 2024

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
Evolution of the star formation rate and Σ<sub>SFR</sub> of galaxies at cosmic
morning (4 < z < 10)

A. Calabrò¹,*, L. Pentericci¹, P. Santini¹, A. Ferrara², M. Llerena¹, S. Mascia¹, L. Napolitano¹,
L. Y. A. Yung³, L. Bisigello⁴,⁵, M. Castellano¹, N. J. Cleri⁶,⁷, A. Dekel⁸, M. Dickinson⁹, M. Franco¹⁰,
M. Giavalisco¹¹, M. Hirschmann¹², B. W. Holwerda¹³, A. M. Koekemoer³, R. A. Lucas³, F. Pacucci¹⁴,¹⁵,
N. Pirzkal¹⁶, G. Roberts-Borsani¹⁷, L. M. Seillé¹⁸, S. Tacchella¹⁹,²⁰, S. Wilkins²¹,²², R. Amorín²³,²⁴,
P. Arrabal Haro⁹, M. B. Bagley¹⁰, S. L. Finkelstein¹⁰, J. S. Kartaltepe²⁵, and C. Papovich⁶,⁷

(Affiliations can be found after the references)

Received 27 February 2024 / Accepted 19 June 2024

ABSTRACT

The galaxy-integrated star formation rate (SFR) surface density measurement (Σ<sub>SFR</sub>) has been proposed as a valuable diagnostic of the mass accumulation in galaxies given it is more tightly related to the physics of star formation and stellar feedback than other indicators. In this work, we assembled a statistical sample of 230 galaxies observed with JWST in the GLASS and CEERS spectroscopic surveys to estimate Balmer line-based dust attenuations and SFRs (i.e., from Hα, Hβ, and Hy), along with the UV rest-frame effective radii. We studied the evolution of galaxy SFR and Σ<sub>SFR</sub> in the first 1.5 billion years of our Universe, from a redshift of z ~ 4 to z ~ 10. We found that Σ<sub>SFR</sub> is mildly increasing with redshift with a linear slope of 0.16 ± 0.06. We explored the dependence of SFR and Σ<sub>SFR</sub> on stellar mass, showing that a star-forming main sequence and a Σ<sub>SFR</sub> main sequence are present out to z = 10. This dependence exhibits a similar slope compared to the same relations at lower redshifts, but with a higher normalization. We find that the specific SFR (sSFR) and Σ<sub>SFR</sub> are correlated with the [O III] λ5007 Å/[O II] λ3727 Å ratio and with indirect estimates of the escape fraction of Lyman continuum photons; hence, they are likely to play an important role in the evolution of ionization conditions at higher redshifts and in the escape of ionizing radiation. We also searched for spectral outflow signatures in the Hα and [O III] emission lines in a subset of galaxies observed at high resolution (R= 2700) by the GLASS survey, finding an outflow incidence of 2/11 (= 20%<sup>32%</sup>/<sub>9%</sub>) at z < 6, but no evidence at z > 6 (0/6, <26%). Finally, we find a positive correlation between A<sub>V</sub> and Σ<sub>SFR</sub>, and a flat trend as a function of sSFR, indicating that there is no evidence of a drop in A<sub>V</sub> in extremely star-forming galaxies between z ~ 4 and ~10. This result might be at odds with a dust-clearing outflow scenario, which may instead take place at redshifts of z ≥ 10, as suggested by some theoretical models.

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

**1. Introduction**

The cosmic evolution of galaxy star formation rates (SFRs)
is one of the fundamental predictions of astrophysical mod-
els and cosmological simulations. It also stands as one of the
most widely studied processes observationally. Indeed, it pro-
vides essential insights into cosmic structure formation across
all scales, the accretion of gas into these structures, the effi-
ciency of conversion into stars, and, ultimately, the diffusion of
baryonic material in the intergalactic medium (IGM) through
stellar feedback (White & Rees 1978; White & Frenk 1991;
Springel & Hernquist 2003; Shapley 2011; Hopkins et al. 2012;
Behroozi et al. 2013; Madau & Dickinson 2014).

The SFR is intimately linked to other galaxy properties,
most importantly the stellar mass (M<sub>⋆</sub>). A correlation between
SFR and M<sub>⋆</sub>, known as the main sequence of star formation
(Noeske et al. 2007), has been determined across over five orders
of magnitude in M<sub>⋆</sub> at all redshifts. While the slope, m, of this
relation remains rather constant with redshift, with m ranging
between 0.5 and 1.2, depending on both selection effects and
the method to measure the SFRs (see Speagle et al. 2014 for a
review), its normalization increases at earlier cosmic epochs.
Indeed, focusing on the specific SFR (sSFR=SFR/M<sub>⋆</sub>), that is,
the SFR normalized by the total stellar mass content, it follows a
rather smooth, monotonic increase by at least one order of mag-

nitude from a redshift of 0 up to the Epoch of Reionization (EoR, Speagle et al. 2014; Pannella et al. 2015; Schreiber et al. 2015; Santini et al. 2017; Davidzon et al. 2018; Leslie et al. 2020; Popesso et al. 2023), which is supported by theoretical models and simulations (Furlong et al. 2015; Henriques et al. 2015; Behroozi et al. 2019; Donnari et al. 2019; Di Cesare et al. 2023).

In addition to the SFR and the sSFR, a quantity that is gain-
ing increasing attention (now that JWST has been equipped to
resolve galaxy sizes out to the EoR) is the surface density of star
formation (Σ<sub>SFR</sub>). This quantity represents the SFR normalized
by the surface area where it occurs and encapsulates informa-
tion about the spatial distribution of star formation. It is usually
defined by the expression Σ<sub>SFR</sub> = SFR/(2π × r<sub>e</sub><sup>2</sup>), where r<sub>e</sub> is the
half-light radius. As galaxies become increasingly more compact
at higher redshifts, since the dependence on the size is quadratic,
Σ<sub>SFR</sub> increases more rapidly with redshift compared to the sSFR.
In particular, it can rise by more than three orders of magnitudes
in typical star-forming galaxies from z ~ 0 to z ~ 6 (Wuys et al.
2011; Holwerda et al. 2015), reaching extreme conditions at the
epoch of reionization (Σ<sub>SFR</sub> ≥ 10 M<sub>☉</sub>/yr/kpc<sup>2</sup>). Therefore, it is
even more essential to characterize the peculiar conditions of SF
in the first phases of galaxy assembly.

Being more intimately related to the physics of star forma-
tion, such as the gas mass surface density through the Kennicutt-
Schmidt relation (Kennicutt 1989; Kennicutt et al. 2007) and the
effectiveness of stellar feedback, Σ<sub>SFR</sub> is thought to regulate the

<footnote>* Corresponding author; antonello.calabro@inaf.it</footnote>
redshift evolution of the sSFR and $M_\star$ (Lehnert et al. 2015). Moreover, for the same reason, Salim et al. (2023) introduced the term “$\Sigma_{SFR}$ main sequence” to indicate the $\Sigma_{SFR}$ versus $M_\star$ relation, and they claim that this is even more fundamental than the sSFR – $M_\star$ relation. Along this new “main sequence”, star-forming galaxies have remarkably similar values of $\Sigma_{SFR}$ at $z \sim 0$ (as also found independently by Schiminovich et al. 2007 and Förster Schreiber et al. 2019), with only a weak dependence (i.e., a positive correlation) on stellar mass across over three orders of magnitudes from $10^8$ to $10^{11} M_\odot$. This is indicative of a more universal and self-regulated nature of $\Sigma_{SFR}$. This relation exists up to at least $z \sim 2$ (Salim et al. 2023), even though it becomes slightly steeper compared to $z \sim 0$, which they interpret as evidence of bulge formation in the central regions of more massive systems.

The level of $\Sigma_{SFR}$ influences that of other galaxy properties. First, the increase in $\Sigma_{SFR}$ is tightly coupled to the increase of the interstellar medium (ISM) gas pressure and of the electron density, $n_e$, with redshift (Jiang et al. 2019; Reddy et al. 2023). The increase of $n_e$ (and consequently $\Sigma_{SFR}$) at fixed stellar mass may be entirely responsible for the higher ionization parameter of high redshift galaxies (Reddy et al. 2023). Thus, $\Sigma_{SFR}$ assumes the primary role that was previously attributed to metallicity. Secondly, $\Sigma_{SFR}$ is tightly related to the feedback of star formation, regulating the power of the outflows and how efficiently the gas is removed from the galaxies, thus affecting the future evolution of the entire system. For instance, previous works have shown that $\Sigma_{SFR}$ correlates with the mass loading factor of the outflow (Llerena et al. 2023) and with the maximum outflow velocity (Kornei et al. 2012; Bordoloi et al. 2014; Heckman & Borthakur 2016; Sugahara et al. 2019; Prusinski et al. 2021; Calabrò et al. 2022a). Similarly, an enhanced $\Sigma_{SFR}$ and stellar feedback reduce the gas and dust covering fraction in galaxies by creating channels that allow Ly$\alpha$ and Lyman-continuum (LyC) photons to escape into the IGM (Reddy et al. 2022). This picture is supported by observations suggesting a relation between $\Sigma_{SFR}$ and the escape fraction $f_{esc}$ of ionizing photons (Heckman 2001; Naidu et al. 2020; Flury et al. 2022), according to which galaxies with higher $\Sigma_{SFR}$ at both low and intermediate redshift have larger fractions of LyC leakers; also, LyC leakers have higher $\Sigma_{SFR}$ than the average galaxy population.

The key role of $\Sigma_{SFR}$ in driving outflows and easing the escape of ionizing radiation is also supported by models and simulations. For example, Sharma et al. (2016) and Sharma et al. (2017) implement galactic winds in the EAGLE cosmological simulations in such a way that they are capable of increasing $f_{esc}$ to more than 20% when $\Sigma_{SFR}$ reaches a value of $0.1 M_\odot/\text{yr}/\text{kpc}^2$. With this recipe, they are able to explain the redshift evolution of the average $f_{esc}$ of galaxies and the volume filling factor of ionized gas up to $z \sim 8$. Naidu et al. (2020) presented an empirical model of the evolution of $f_{esc}$ in galaxies is only dependent on $\Sigma_{SFR}$, which thus assumes a key role in reionization. With this tight connection between $f_{esc}$ and $\Sigma_{SFR}$, they claim that the bulk of reionization ($\sim 80\%$) must be driven by a small number of galaxies ($<5\%$), the so-called “oligarchs”, with extremely high $\Sigma_{SFR}$ ($\sim 10-20 M_\odot/\text{yr}/\text{kpc}^2$), compact size, and relatively high UV luminosity ($M_{\text{UV}} < -18$) and stellar mass ($M_\star > 10^8 M_\odot$).

The SFR in the first phases of galaxy assembly can also reach the super-Eddington regime, a condition in which the radiation pressure from young stars overcomes the gravitational force, leading to a galaxy-wide outflow that can clear the galaxy of its dust and gas content, pushing them away to large radii in a short timescale of the order of a few Myr (Ferrara et al. 2023;

Ferrara 2024). These extreme conditions of star formation are invoked to explain the stunning abundance of UV luminous, very blue, and massive galaxies at $z > 10$ found since the first extragalactic observations with JWST (e.g., Naidu et al. 2022; Castellano et al. 2022).

The goal of this paper is to investigate how the SFR and the SFR surface density evolve in the first 1.5 billion years of our Universe. The investigation of galaxy assembly during the reionization epoch is perfectly-timed. On the one hand, the advent of JWST has allowed us to measure more compact galaxy sizes compared to previous telescopes and to resolve internal structures at much finer levels. On the other hand, the unique sensitivity of JWST in the rest-frame optical and UV of galaxies has allowed us to perform large photometric and spectroscopic surveys, such as GLASS (Treu et al. 2022) and CEERS (Finkelstein et al., in prep.), targeting statistically representative samples of galaxies during reionization down to $M_\star$ of $10^7 M_\odot$, and constraining the low SFR levels expected for these low-mass systems. This is essential to trace the scaling relations followed by the global galaxy population since the earliest epochs of their formation. In addition to the evolution of star formation, we are interested in exploring how $\Sigma_{SFR}$ is related to other galaxy properties and compare the observed values to those required by models and simulations to trigger galaxy wide outflows and to carve channels for the escape of Lyman continuum radiation.

The paper is organized as follows. In Section 2, we describe the spectroscopic observations, sample selection, and the derivation of all the physical properties investigated in this paper. In Section 3, we show our results, including the redshift evolution of the $\Sigma_{SFR}$ from $z = 4$ to $z = 10$, the SFR and the “$\Sigma_{SFR}$ main sequence”. Finally, in Section 4, we compare our findings with theoretical model predictions, and explore the relation between $\Sigma_{SFR}$ and other galaxy properties to understand its role in the reionization process. Then we analyze the presence of outflows in NIRSpec galaxies and connect them to the star formation rate and dust attenuation. For our analysis, we adopt a Chabrier 2003 initial mass function (IMF), and we assume a standard cosmology with $H_0 = 70 \text{ km s}^{-1}\text{Mpc}^{-1}$, $\Omega_m = 0.3$, and $\Omega_\Lambda = 0.7$.

## 2. Methodology

The aim of this work is to analyze fundamental scaling relations describing the evolution of star formation in galaxies; hence, it is important to assemble a statistical sample as large as possible. Here, we use the data obtained by two JWST-based Early Release Science (ERS) programs, namely GLASS and CEERS, which represent the largest publicly available JWST survey thus far, with hundreds of high-redshift galaxies observed in the optical and near-infrared rest-frame. In the following part we analyze each survey in more detail, and explain how we derive the physical properties of galaxies from the available spectroscopic, imaging, and photometric datasets.

### 2.1. Spectroscopic observations and sample selection

We first considered spectra acquired through NIRSpec MSA observations by the GLASS-JWST ERS Program (PID 1324, PI: Treu; Treu et al. 2022) and by a JWST DDT program (PID 2756, PI: W. Chen; Roberts-Borsani et al. 2022). Both of them targeted galaxies residing over the central regions of the Frontier Field galaxy cluster Abell 2744. However, they were obtained using two different setups and different exposure times. The former were taken in high-resolution grating mode with three spectral configurations (G140H/F100LP, G235H/F170LP, and
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Redshift distribution of stellar masses (top panel) and SFRs (bottom panel) for the galaxies selected in this paper. We highlight the redshift histogram of the entire sample on top of the first diagram, while the $M_*$ and SFRs histograms are added on the right part of each panel. Galaxies coming from the CEERS and GLASS surveys are differentiated with a yellow and blue color, respectively. We note that GLASS probes galaxies with slightly lower masses and SFRs on average, thanks to the gravitational lensing magnification effect. The gray vertical bars in the top right corners represent the typical uncertainties (~0.2 dex) for the stellar masses and SFRs. The 90% mass completeness limits at $z \sim 5$ and $z \sim 9$ are highlighted, respectively, with black and lime horizontal segments.</figcaption></figure>

G395H/F290LP), covering a wavelength range from 1.5 to 5.14 µm at R ~ 2000–3000, and with a total integration time (per grating) of 4.9 hours. The latter instead were taken in low-resolution mode with the CLEAR filter + prism configuration, which provides a continuous wavelength coverage from 0.6 to 5.3 µm at R ~ 30–300. The prism resolution in the red part of the spectrum is enough to resolve and detect redshifted rest-frame optical emission lines, such as $H\beta$ and [O III] $\lambda5007$ Å. The exposure time in this case is of 1.23 hours. The details of the spectral reduction, wavelength and flux calibration are described in detail in Mascia et al. (in prep.). Since GLASS-JWST observations target sources behind a lensing cluster, we adopt the magnification factors ($\mu$) listed by Mascia et al. (2023), based on the latest lensing model of Bergamini et al. (2023), to correct line fluxes, SFRs, and galaxy sizes.

We also consider the spectra obtained by the CEERS survey through the NIRSpec MSA mode. The CEERS spectroscopic data were acquired over 3 different observing epochs and in 8 different pointings. The NIRSpec pointings named p4, p5, p7, p8, p9, and p10 were observed with the three NIRSpec medium resolution (G140M, G235M, and G395M) gratings, covering the observed-frame wavelength range 0.97–5.10 µm with a resolving power R ~ 1000. In addition, p4, p5, p8, p9, p11, and p12 were also observed with the low resolution prism. In all cases, the exposure time is of 3107s per spectrum. Three-shutter slitlets with a three-point nodding pattern were also used to enable accurate background subtraction. After the observing run, some galaxies have spectra obtained with both the prism and M-grating configuration, in which case we prioritize the latter, as the higher resolution of the M-grating improves the detection of faint emission lines with low signal-to-noise ratio (S/N). The main steps of the reduction, the optimal extraction of the final 1D spectra, including also the wavelength and flux calibration, follow those described in other works of our collaboration, namely Fujimoto et al. (2023), Kocevski et al. (2023), and Arrabal Haro et al. (2023). A more detailed description of the NIRSpec data processing and a final spectroscopic catalog of the entire survey will also be presented in Arrabal Haro et al. (in prep.).

To build a sample of high redshift galaxies for our analysis, we visually checked all the 1D spectra provided by the above surveys in order to assign a first value to the spectroscopic redshift $z_{\text{spec}}$ based on the visual identification of relatively bright, optical emission lines, including [O II] $\lambda3727$ Å, $H\beta$, [O III] $\lambda5007$ Å, and $H\alpha$. This process was guided by the previous knowledge of the photometric redshifts of all the sources, obtained with the EAZY code (Brammer et al. 2008). In the next section, we describe a more precise measurement of $z_{\text{spec}}$. We also remark that the spectroscopic redshifts were checked and confirmed independently by other team members in both the GLASS and CEERS collaborations. Part of our spectroscopic sample is already published by other works (e.g., Arrabal Haro et al. 2023; Mascia et al. 2023, 2024), in which case we just checked the consistency of our $z_{\text{spec}}$ measurements.

Among the sample observed with NIRSpec by either the GLASS or CEERS survey, we selected galaxies with a spectroscopic redshift $z_{\text{spec}} \ge 4$. The reason for this lower limit is twofold. First, we aim to study the properties of the galaxy population as we approach the epoch of reionization, where most of the scaling relations addressed in this work are still largely unexplored compared to cosmic noon ($1 < z < 3$). Secondly, targets with $z \ge 4$ were the main priority of the GLASS and CEERS spectroscopic surveys, thus, we have better statistics at our disposal in this range. This condition yields a final sample of 230 galaxies that we consider in our analysis. Among this sample, 19 have H-grating spectroscopic data (R~ 2700) in GLASS, 93 have M-grating spectroscopy from CEERS (R~ 1000), while 15 and 114 are observed in prism configuration (R~ 100) in the Abell 2744 and EGS fields, respectively. The global population spans a redshift range between $z = 4$ and 10.4, with a median value $z_{\text{med}} = 5.4$ (see Fig. 1-top). The number of sources decreases toward higher redshifts, with 11 selected galaxies residing at redshift $\ge 8$.

We can notice in Fig. 1 that the GLASS sample includes 5 galaxy members of a spectroscopically confirmed protocluster at redshift ~8 (Morishita et al. 2023). Furthermore, four galaxies in the same field have spectroscopic redshifts ranging between 4.01–4.05. All of them are more than 1' away from each other, but the whole group lies within a projected radius of 80'' (i.e.,
~0.6 pMpc at redshift 4). This is much larger than the protoclus-
ter at z ~ 8, but more similar to the size of some galaxy overden-
sities identified at redshifts ~3 (e.g., Calabrò et al. 2022b).

Given the variety of observational programs considered, it is
useful to better investigate the properties of our selected sample
to assess potential biases on our analysis. The GLASS obser-
vations with NIRSpec cover approximately 50 arcmin² in the
sky. Given the relatively long exposure times, they are among
the deepest available so far, as we can detect emission lines with
fluxes down to a few 10⁻¹⁸ erg/s/cm². In addition, the lensing
magnification is ideal for improving the detection of the intrisici-
cally faint, low-mass galaxy population at the epoch of reioniza-
tion. The primary science targets of the NIRSpec observations
were placed into open Microshutters (MSA), and were chosen
based on previous spectrophotometric catalogs available from
the GLASS team and from the literature. Many of them already
had a confirmed spectroscopic redshift, while the remaining sys-
tems were photometrically selected at z > 4 using previous HST
observations in the region of Abell 2744 (see Treu et al. 2022
for more details). The CEERS spectroscopic observations cover
instead a larger area in the sky of about 90 arcmin², but they
are shallower compared to GLASS. Some of the spectroscopic
targets are HST-selected in the UV rest-frame (among which
we find those galaxies falling outside of the NIRCam cover-
age), while additional targets benefitted from NIRCam obser-
vations undertaken in the first observing run (prior to NIRSpec
observations), and therefore are optically selected, even though
they might still be biased to UV-bright LBG-type objects as they
are chosen based on their photometric redshift. While a more
detailed description will be included in Finkelstein et al. (in
prep.), we note that it is difficult to precisely assess the selection
function even within the same survey. It is thus a better approach
to carry out an a posteriori analyis of the galaxy properties as a
whole.

Even though GLASS probes slightly lower masses and SFRs
on average compared to CEERS, due to the gravitational lensing
effect, the two subsets are similarly distributed in redshift, as
shown in the top histogram of Fig. 1 (top). Also, they are located
along the same M\* versus SFR relation. Moreover, in the red-
shift evolution of Σ\*SFR, and in all the other diagrams analyzed
in this paper, removing the small GLASS subset does not sig-
nificantly change the results and the best-fit relations, meaning
that they satisfy the same physical trends. We also note from
Fig. 1 that both subsets may suffer from slight incompleteness
at redshift ≥8, where we may lose a fraction of galaxies at the
lowest masses and SFRs (respectively log M\*/M☉ < 8 and log
SFR/(M☉/yr) < 0. However, as explained in Section 3.3, our
galaxies can be fairly considered as representative of the star-
forming main sequence for most part of our redshift range.

Finally, even though the standard STScI reduction pipeline
adopted for our spectra already includes a slit loss correction,
we further checked the absolute flux calibration of the output 1D
spectra by comparing them to the available broadband photome-
try (see Section 2.3 for the specific filters used). For all galaxies,
we considered total fluxes estimated in the CEERS and GLASS
collaborations. In the first survey, these fluxes are obtained by
first performing the photometry in elliptical Kron apertures and
then applying aperture corrections estimated from simulations
by matching the image PSFs between different filters, as detailed
in Finkelstein et al. (2024). A similar procedure for obtaining
total fluxes was followed for the GLASS survey by Paris et al.
(2023), using the code A-PHOT (Merlin et al. 2019). Despite the
fact that the calibration corrections for individual photometric
bands do not show a significant wavelength dependence, as also

found by other works in the collaboration (e.g., Napolitano et al.
2024), we limited the comparison to the redder spectral range
at 2 µm < λ < 5 µm; this is because the S/N of the continuum
rapidly drops at lower wavelengths, making this check and com-
parison more difficult. A second reason for this choice is that we
are interested only in optical rest-frame lines, which mostly lie
in the red part of the spectrum. For the correction, we applied
a constant multiplicative factor to the spectrum; this factor was
obtained as the median value of the corrections derived indepen-
dently for all the photometric bands whose central wavelengths
reside in the above spectral range.

## 2.2. Line measurements and derivation of star formation rates

We measured the emission line properties with a similar method-
ology to that adopted in Calabrò et al. (2019) and Calabrò et al.
(2023). In brief, we used a python version of the χ² minimization
routine MPFIT (Markwardt 2009) to fit the following emission
lines: [O II] λ3727 Å, Hγ, [O III] λ5007 Å, Hβ, and Hα, together
with their underlying continuum in a spectral window of ±5000
km/s around the line. We assumed a Gaussian function for the
emission lines, while the continuum was modeled with a first-
order polynomial. We assume for the lines a common redshift
and line velocity width σ, with a tolerance of 100 (1000) km/s on
the line centroid and σ for the H- or M-grating and prism spec-
tra, respectively. The fit yields the exact redshift of the sources,
the line widths, the line fluxes, and their corresponding uncer-
tainties. The quality of the fit is controlled through the reduced
χ² and visual inspection for each galaxy. We consider the lines as
detected when they are above a S/N of 3, while for non-detected
lines, we adopted 3σ upper limits, which were then propa-
gated across all the physical quantities derived. As discussed
in Calabrò et al. (2023), the underlying stellar absorption for
Balmer lines can be neglected as it does not produce significant
variations of the emission line fluxes for galaxies in our stellar
mass range. Moreover, our galaxies typically have younger ages,
lower metallicities, and lower dust content compared to the lower
redshift CEERS sample analyzed in Calabrò et al. (2023); there-
fore, the underlying absorption is expected to be even smaller
(Groves et al. 2012). For this reason, we did not correct the emis-
sion line fluxes for stellar absorption. We also noticed that in low
resolution spectra acquired with the prism configuration, Hα is
blended with the [NII] doublet. However, our galaxies are all
very low-mass and with sub-solar metallicities, for which the
expected contamination is less than 10%, according to lower
redshift studies (Faisst et al. 2018). A negligible contribution of
[N II] to Hα for galaxies at z > 4 was also found by the pho-
toionization models described in Simmonds et al. (2023) and
supported by the results from our high resolution dataset. As a
consequence, we do not attempt to make any correction on the
Hα flux in prism spectra. Finally, two sources at redshifts 5 <
z < 6 with clear broad Hα components and narrow [O III] λ5007
Å, indicative of broad-line AGN emission, are removed from the
sample. We refer to Calabrò et al. (2023) for a more detailed dis-
cussion of the emission line measurement.

As the first step of our analysis, we computed the gas-phase
attenuation, A<sub>V</sub>, from the available Balmer lines (among Hγ,
Hβ, and Hα) assuming a dust-screen, Milky Way like extinction
law (Cardelli et al. 1989). If only one Balmer line is detected,
we adopted an A<sub>V</sub> as inferred from the SED fitting proce-
dure described later in Section 2.3. As shown in the Appendix
of Calabrò et al. (2023), hydrogen recombination lines across
a large range of wavelengths (~0.4 to ~2 µm) are overall in
agreement with the predictions of typical dust attenuation models, even when the lines reside in different gratings, indicating that their relative calibration and the inferred values of $A_V$ are generally reliable. From the fluxes of the [O II] $\lambda3727\AA$ and [O III] $\lambda5007\AA$ lines, we also derived the O32 index, defined as the logarithm of the reddening corrected [O III] $\lambda5007\AA$/$[\text{O II}] \lambda3727\AA$ ratio.

We then derived the star formation rates as follows. We first computed the observed, dust-corrected $H\alpha$ luminosity $L_{H\alpha}$. If $H\alpha$ was not detected, we considered the fluxes of $H\beta$ (or $H\gamma$ as last possibility), corrected for dust attenuation and rescaled to $H\alpha$ assuming the intrinsic ratios of 2.86 and 6.10 (respectively, for $H\beta$ and $H\gamma$), valid for a case B recombination, a temperature $T = 10\,000\,\text{K}$, and an electron density $n_e = 10^2 \, \text{cm}^{-3}$ (Osterbrock 1989). We then converted $L_{H\alpha}$ to SFR assuming the calibration of Reddy et al. (2022) as $\text{SFR} = L_{H\alpha} \times 10^{-41.67}$. This is the most suitable calibration according to the typical subsolar metallicities expected for our galaxies at $z > 4$; it also reflects the greater efficiency of ionizing photon production in metal poor stellar populations. In our sample, we obtained the SFR from $H\alpha$ in 153 galaxies, from $H\beta$ in 67 systems, while for ten galaxies, it is inferred from $H\gamma$. The resulting SFR distribution is shown in Fig. 1 (bottom), where our galaxies span values between ~0.1 and 300 $M_{\odot}/\text{yr}$, with a median SFR of ~5 $M_{\odot}/\text{yr}$.

We note some possible caveats of our analysis, affecting the relation between the SFR and $H\alpha$ luminosity. For example, the diffuse warm ionized gas (DIG) or collisional excitation in the ISM might contribute to the observed $H\alpha$ emission, leading us to slightly overestimate the SFRs. On the other hand, part of the ionizing radiation may be absorbed by dust and helium, or might escape from the galaxy, reducing the amount of Lyman continuum (LyC) photons available for hydrogen line emission. This might lead us to slightly underestimate the total SFR of the galaxy. Modeling and testing these effects is difficult, especially at high redshifts. However, despite these uncertainties, the simple analysis performed in this paper and the dust correction through the Balmer decrement should work well as a first approximation, as shown by Tacchella et al. (2022).

We also note that the dust correction is another potential source of uncertainty. This is due both to a poor theoretical understanding of dust production mechanisms and dust growth in the ISM at high redshift (see discussion in Markov et al. 2023) and to the lack of observational constraints. For this reason, different dust attenuation and extinction laws measured in the local Universe are usually adopted in the literature also at early epochs. We remark that our choice for dust correction is the same adopted by similar works on NIRSpec galaxies at the epoch of reionization (e.g., Llerena et al. 2024). Adopting a Calzetti et al. (2000) attenuation law, as done, for example, in Mascia et al. (2024), does not significantly alter the results.

## 2.3. Photometric data and estimation of stellar masses through SED fitting

We derived stellar masses for our galaxies through an SED fitting procedure using the code CIGALE (Boquien et al. 2019), version 2022.1. We fit Bruzual & Charlot (2003) stellar population models with Chabrier IMF to the available observed photometry ranging from 0.5 to 8 µm. There is a homogeneous photometric coverage for all the galaxies in our sample. In particular, all the galaxies observed by the GLASS survey are covered by JWST/NIRCam imaging from the UNCOVER program (JWST-GO-2561; Bezanson et al. 2022). For this subset, fluxes

and uncertainties are obtained in the F115W, F150W, F200W, F277W, F356W, F410M, and F444W filters, using the catalog of Paris et al. (2023). The majority of our NIRSpec sources in the EGS field are also covered by JWST/NIRCam observations from the CEERS survey in the same filters listed above (Bagley et al. 2023). For this subset, we took the fluxes and uncertainties from Finkelstein et al. (2024). For all the NIRCam sources, we used the total flux estimates mentioned in Section 2.1. However, a subset of CEERS sources (79) do not have JWST imaging. To ensure that in this case, we also have access to complete coverage of the UV and optical rest-frame emission for the SED fitting, we considered the following broadband filters: HST/ACS F606W and F814W, HST/WFC3 F125W, F140W, and F160W, CFHT/WIRCAM J, H, and K_s, Spitzer IRAC 3.6 and 4.5 µm. In particular, we used the total photometric fluxes and associated uncertainties from the multi-wavelength catalog assembled by Stefanon et al. (2017).

The SED fitting with CIGALE was performed as follows. The redshift of the galaxy was fixed to the spectroscopic value estimated in Section 2.2. We considered a main stellar population with metallicities $Z_* = 0.0004, 0.004, 0.008, \text{ or } 0.02$, and a delayed SFH parameterized by an e-folding time that can assume the following values: 10, 100, 500, and 1000 Myr, and with a grid of 18 possible stellar ages ranging from 50 Myr to 1.3 Gyr. We also included the possibility of a recent exponential burst with SFH $\propto \exp(-t/\tau_{\text{burst}})$, where $\tau$ can assume the values 1, 20, and 50 Myr, while the ages were allowed to range from 5 to 50 Myr. We also allowed the mass fraction of the late burst population to vary from 0 to 80% of the total mass formed in the entire galaxy history. We included contribution from nebular emission, where the nebular component is parametrized with an electron density $n_e$ of 100 $\text{cm}^{-3}$, subsolar gas-phase metallicity equal to the stellar metallicity, and ionization parameter between -3 and -2. Finally, the dust attenuation of the stellar and nebular components were modeled with a Milky Way extinction curve (Cardelli et al. 1989) and parametrized with a color excess $E(B-V)_{\text{stellar}}$ ranging from 0 to 3 in steps of 0.05, a total to selective extinction $R_V = 3.1$, and a stellar continuum to nebular attenuation ratio of 0.44 (Calzetti 2001). Therefore, for each galaxy, $A_V$ was computed as $R_V \times E(B-V)$. The intergalactic medium (IGM) transmission is also taken into account following the model of Meiksin (2006). After the SED-fitting, we correct the stellar masses of the GLASS sources for the effects of lensing, using the magnification factors mentioned above. As demonstrated by Furtak et al. (2021), this approach yields consistent results compared to correcting the photometry before the SED-fitting.

We find stellar masses $M_*$ ranging from $10^7$ to $10^{10.5} M_{\odot}$ (see Fig. 1-top), with a nearly symmetric distribution around a median $\log_{10} M_*/M_{\odot} = 8.7$ ($1\sigma = 0.7$). Combining these results with the SFR measurements described in the previous section, we find that the bulk of our galaxy population has sSFRs between 1 and 100 Gyr$^{-1}$, with median value of 10 Gyr$^{-1}$. For galaxies fitted with $E(B-V)= 0$, we also set an upper limit on their dust attenuation as $E(B-V)= 0.05$, corresponding to our grid step and to the average uncertainties obtained with this methodology. Comparing the values of $A_{V,\text{nebular}}$ obtained from the SED fitting to those estimated from Balmer lines, we find that they are overall consistent with the 1 : 1 relation, even though with a large scatter.

For galaxies with multiple coverage, we also checked that removing NIRCam photometry and using only the remaining bands (HST + Spitzer + ground) in the SED fitting yields $M_*$ values that are in agreement with estimates based on JWST photometry in the entire range spanned by our sample. Therefore,
we did not introduce systematic biases (related to the measurement method) on the derived parameters for the subset that is not covered by JWST imaging. Our stellar masses are consistent with those obtained using the code ZPHOT (Fontana et al. 2000) with a similar setup and SFH to our CIGALE run, indicating that this quantity is rather robust and not significantly affected by the exact fitting procedure. Finally, we also tested a nebular attenuation equal to the stellar attenuation, but this choice did not produce significant variations in the final results of this paper.

## 2.4. Galaxy sizes and SFR surface density estimation

We measured the major-axis half-light radius $r_e$ of our galaxies by fitting their rest-frame UV images with the python software GALIGHT (Ding et al. 2021). A subset of 151 galaxies is covered by JWST imaging, in which case we adopted the background-subtracted mosaics described in Bagley et al. (2023) and Paris et al. (2023), for the CEERS and GLASS subset, respectively. For the remaining subset of 79 galaxies (all in the CEERS field), we use instead HST imaging (Koekemoer et al. 2011). To homogeneously probe the same rest-frame UV window, we performed the measurements in the F115W band for galaxies at $z < 5.5$, in the F150W filter for those at slightly higher redshifts ($5.5 < z < 7.5$), and in F200W for systems at $z > 7.5$. In cases where only HST imaging is available, we adopted the F125W and F160W band for galaxies at redshifts lower and higher than 5.5, respectively.

To measure the galaxy sizes, we first produced $3'' \times 3''$ cutouts, which are given as input to the GALIGHT code. All sources in the cutout, including the central galaxy of interest, were fitted simultaneously assuming that they are well represented by a Sersic profile. We constrained the following parameters to keep the fit within physically meaningful values: the Sersic index $n$ is free to vary from 0.2 to 8.0, and the axis ratio $q$ is set in the range $0.1 < q < 1$. After running GALIGHT, sources are flagged according to the quality of the fit, with a flag = 0 assigned to reliable fits, and flag = 1 for sources that are not fitted well by the above procedure (indicative of a more complex shape), as checked visually in the observed versus best-fit luminosity profile that is output by the code. As a result, 213 galaxies with a good size quality flag were considered when analyzing the $\Sigma_{SFR}$, as described in the following sections.

The size uncertainties are assessed in accordance with Yang et al. (2022) and rescaled based on the S/N derived from the photometry. Previous studies (e.g., Kawinwanichakij et al. 2021) have shown that the outcomes from GALIGHT are robust, and consistent with estimations made through conventional softwares like Galfit (Peng et al. 2002). A subset of our galaxies are indistinguishable from point sources. In particular, we performed a subset of simulations, finding that 0.025" and 0.1" are the smallest measurable sizes in JWST and HST images, respectively. Therefore, we set these values as upper limits. In our sample, 64 galaxies of 213 (~30%) are unresolved and have upper limits on $r_e$. We also notice that for the galaxies lying in the GLASS field (which are all covered by JWST imaging), the sizes are corrected for the effect of lensing. In particular, we divided the original sizes by $\sqrt{\mu}$ assuming an isotropic magnification, which represents a reasonable approximation considering the low magnification factors (median $\mu \approx 2$) found in our sample.

With the estimated sizes, we compute the galaxy-integrated SFR surface densities for all the objects. A subset of 25 extremely compact galaxies have an upper limit on the size but also an upper limit on the SFR. We exclude these galaxies from the analysis as their $\Sigma_{SFR}$ would be totally unconstrained. As

a result, the investigation of the redshift evolution of $\Sigma_{SFR}$ and its dependence on other galaxy properties ought to be made on a sample of 188 galaxies. The distribution of $\Sigma_{SFR}$ for this galaxy subset has a median value of $\log_{10} \Sigma_{SFR} = 0.6$ (in units of $M_{\star}/\text{yr/kpc}^2$, $1\sigma = 0.8$), with individual systems spanning the logarithmic range between -1.5 and 2.5.

The UV rest-frame based $r_e$ values adopted in this paper probe the size of star formation in the last ~100 Myr (i.e., the continuum emission from young and massive stars), a slightly longer timescale than the $H\alpha$ (or $H\beta$) based SFRs, which are sensitive to the last 10–30 Myr of star formation. Ideally, more suitable emission line-based sizes could be derived through IFU spectral observations or from imaging observations in photometric bands dominated by emission lines. To check whether the latter case occurs within our sample, we have estimated for all the galaxies the emission line contribution to the observed photometry. In particular, we compared for all the available photometric bands the in-band fluxes obtained by integrating the observed spectrum across each filter transmission curve and by summing the emission lines only. With this procedure, we find that the maximum emission line contribution to the photometry is of ~20% on average, exceeding 60% in 9 galaxies (~3.5% of the sample). However, for this small subset, the optical rest-frame sizes (in bands where the emission lines dominate) are in agreement with those estimated in the rest-frame UV, suggesting that using the latter does not significantly alter the final $\Sigma_{SFR}$ on average. Moreover, recent independent studies at similar redshifts based on JWST-NIRCam observations (e.g., Ono et al. 2024; Morishita et al. 2024) have shown that rest-frame UV sizes are in agreement with rest-frame optical sizes (where we might have contribution from strong optical emission lines).

## 2.5. Estimating sample completeness and potential selection effects

In this section, we explain how we derived an approximate estimation of the mass completeness limits for our sample. The majority of our NIRSpec selected galaxies at redshifts $4 \le z \le 10$ is made of JWST-NIRCam selected sources within a magnitude limit in F277W of ~28.5 mag. Using the recent results of Cole et al. (2023) and rescaling them to account for the different limiting magnitudes in F277W, we estimated a 90% mass completeness limit of $10^{7.8} M_{\odot}$ at $z \sim 5$ and of $10^{8.4} M_{\odot}$ at $z \sim 9$ (i.e., ~0.2 dex lower than the limits of Cole et al. 2023). This redshift-dependent completeness has effects on the stellar mass and SFR distributions shown in Fig. 1 and is likely responsible for the strong rise observed at $z > 7.5$ of the minimum stellar mass and SFR probed within our sample.

We also note that no additional selection effects have been added to the spectroscopic sample by our analysis, as we have a SFR measurement (or at least an upper limit) and a size measurement from UV rest-frame imaging for all the galaxies selected in Section 2.1. The galaxies without a measurement of $\Sigma_{SFR}$ (because of an upper limit on both the SFR and the size) or without a size (because of a complex shape) are randomly distributed in redshift and in the $M_{\star}$–SFR plane, so they do not produce systematic biases in the results.

We also investigate the effect of surface brightness dimming on source identification and sample selection. To this aim, we performed a set of Monte Carlo simulations, creating mock observations of galaxies in the F277W band following a similar approach to that adopted Treu et al. (2023) (see their Section 3.2). We placed the sources at different redshifts from $z = 4$ to $z = 10$ (in steps of 0.5), assuming for simplicity a single Sersic light pro-
file for the galaxies, with $n = 1$ (consistent with observations of EoR galaxies, Yang et al. 2022), $q = 1$, and $r_e$ ranging from 0.1 to 1.5 kpc (in steps of 0.05 kpc). The background level was set to match the exposure time (5207 s) in CEERS and then 50 simulations were run for each configuration, applying each time a detection algorithm (using the photutils package in python) with an S/N detection threshold of 3. From these simulations, we study how the detection rate varies as a function of redshift and $r_e$ for different total magnitudes of the source, from 28.5 (the magnitude limit of our sample) to 27, in steps of 0.3 mag.

We find that at $z \simeq 9$, the maximum size for which a galaxy can be detected in F277W (at fixed total observed magnitude) is smaller than at $z \simeq 4$ by ~0.2–0.25 dex. However, galaxies become intrinsically smaller at higher redshifts. For example, Ward et al. (2024) have shown that the typical $r_e$ of galaxies decreases by a factor of ~2 (0.3 dex) from $z = 4$ to 9 (at fixed stellar mass); also, the mass-size relation leads to the result of smaller sizes at high-$z$. This suggests that, within our “magnitude-limited” sample, size incompleteness does not play a significant role in the observed evolution of $\Sigma_{SFR}$ at higher redshifts. For example, assuming for galaxies at redshifts $8 < z < 10$ an average F277W magnitude of ~28, and a minimum detected SFR of $5 M_{\odot}/yr$ (Fig. 1-bottom), we are able to probe $\Sigma_{SFR}$ values down to $\log_{10} \Sigma_{SFR}/(M_{\odot}/yr/kpc^2) \sim 0.4$, which is more than $2\sigma$ below the $z - \Sigma_{SFR}$ relation extrapolated without the highest redshift bin.

# 3. Results

Here, we explore the redshift evolution of $\Sigma_{SFR}$ for our selected sample and then focus on fundamental scaling relations involving $M_*$, SFR, and $\Sigma_{SFR}$ in the entire redshift range between ~4 and ~10. To characterize all the 2D distributions presented in this paper and quantify the level of correlation between two quantities, we adopted two different approaches throughout this work. Both take into account (in a different way) the presence of upper and lower limits in the data. We explain them in detail in the next subsection for the redshift versus $\Sigma_{SFR}$ diagram, but the same steps of the analysis are followed for the other diagrams.

## 3.1. Fitting methods

As a first approach to characterize the redshift evolution of $\Sigma_{SFR}$, we computed the median redshift and the median $\Sigma_{SFR}$ of galaxies falling in different, well defined bins of increasing spectroscopic redshift (4–5, 5–6, 6–7, 7–8, > 8), which are shown as black empty diamonds in Fig. 2. In each bin, we also calculated the median absolute deviation (MAD) of $\Sigma_{SFR}$, indicated with black error bars. In this calculation, we assign half of the $\Sigma_{SFR}$ value to galaxies with $3\sigma$ upper limits, or twice the value in case of $3\sigma$ lower limit on $\Sigma_{SFR}$. In data science, imputation of upper (lower) limits with half (double) the detection limit is a common approach, and simulation studies have found that this is a better choice compared to assuming the detection limit itself or a zero value (for upper limits), as it introduces the least bias into the estimates (Beal 2001). Then we are able to derive a global, best-fit linear relation and $1\sigma$ limits through Monte Carlo simulations. In particular, we performed 1000 realizations of our data varying the data points according to their x and y axis $1\sigma$ uncertainties. For each realization, we repeated the same process of median $\Sigma_{SFR}$ estimation in each of the previously defined redshift bins and fit these median values with a first order polynomial. Finally, we derived the best-fit relation as the median slope and normalization of all the different realizations, including the uncertainties as their standard deviation.

In the second analysis method, we used the Python package PyMC<fnref n="1" />. This is an advanced tool for statistical modeling featuring next-generation Markov chain Monte Carlo (MCMC) sampling algorithms such as the No-U-Turn Sampler (NUTS; Hoffman & Gelman 2024), a self-tuning variant of Hamiltonian Monte Carlo (HMC; Duane et al. 1987). HMC and NUTS can manage complex posterior distributions and fitting models, taking advantage of gradient information from the likelihood to achieve much faster convergence than traditional sampling methods. A first advantage of this approach is that it can be applied to the entire dataset without dividing the sample in multiple bins as in the previous case. A second advantage is that it implements the methods of the Bayesian survival analysis in the linear regression to deal with upper and lower limits on the y axis. Most importantly, this method does not impute the limits (i.e., does not strictly substitute them with a different value), but instead integrates them out through the likelihood. In practice, the entire dataset is fitted with a linear relation, obtaining a first guess of the line parameters. Then, pyMC is run by assuming Gaussian priors for the slope and intercept, with $\sigma = 1$ and the first guesses as centers. At this step, normal measurements, upper and lower limits are considered simultaneously in the fit. Finally, the full posterior distributions of the slope and intercept are used to calculate the mean linear relation, the $1\sigma$ uncertainty, and for showing the ‘posterior plot, that is, a stack of random draws from the posterior distribution of the slope and intercept. We report in Table 1 the best-fit parameters obtained for the main relations studied in this paper and using both the linear regression methods described above.

## 3.2. The redshift evolution of the SFR surface density

Our results for the evolution of $\Sigma_{SFR}$ are shown globally in Fig. 2, where the median $\Sigma_{SFR}$ and MAD in five increasing bins of redshift (defined as 4–5, 5–6, 6–7, 7–8, and > 8) are represented with empty diamonds and vertical error bars, while the best-fit to the bins is highlighted with a black dashed line and a gray shaded area indicating its $1\sigma$ uncertainty. In addition, we also consider the Bayesian fitting method, plotting the mean linear relation as a blue continuous line, the $1\sigma$ uncertainty as blue dotted lines, and the ‘posterior plot’ as a blue shaded region around the mean relation.

Regardless of the analysis approach, we find a mild evolution of $\Sigma_{SFR}$ from $z \simeq 4$ to $z \sim 8$, with a slight but statistically significant increase by a factor of two (0.3 dex), going from $2.5 M_{\odot}/yr/kpc^2$ to ~5 $M_{\odot}/yr/kpc^2$. We note that in the short cosmic epoch that we are exploring (around one billion years of cosmic history), the linear relation yields a good description of the increase in $\Sigma_{SFR}$. In all cases we obtain a positive correlation with best-fit slopes of $0.16 \pm 0.06$ and $0.21 \pm 0.03$, respectively from the Bayesian regression and the fit to the bins. Therefore, in all cases we can exclude a flat trend at least at $2\sigma$ level.

Fitting the median $\Sigma_{SFR}$ in the five redshift bins, the slope is slightly higher compared to the other method, as due to a sudden rise of $\Sigma_{SFR}$ in the last bin at $8 < z < 10$, where we reach a median $\Sigma_{SFR}$ of ~20 $M_{\odot}/yr/kpc^2$. However, we notice that the median $\Sigma_{SFR}$ derived in the last bin is consistent with the Bayesian fit within $1\sigma$. We have also removed the galaxies in the last redshift bin and fitted again the sample with the two methods, but we do not obtain significantly different results. Therefore the upturn of $\Sigma_{SFR}$ at $z > 8$ is not significant from our data. This trend might be due to the different sample complete-

<footnote n="1">¹ https://peerj.com/articles/cs-1516/; https://doi.org/10.5281/zenodo.10371446</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Redshift evolution of $\Sigma_{SFR}$ from $z = 4$ to 10. The median values and MAD of $\Sigma_{SFR}$ in bins of increasing redshift are shown with black empty diamonds and black vertical error bars, while the best-fit to the bins is represented with a black dashed line with a gray region indicating its $1\sigma$ uncertainty. The best-fit obtained with the Bayesian linear regression is shown with a blue continuous line, with the blue dotted lines and the blue shaded regions representing the $1\sigma$ uncertainty and the ‘posterior plot’, respectively. Theoretical model predictions and comparison observations are included as specified in the legend.</figcaption></figure>

<table><caption>**Table 1.** Best-fit linear relations of the diagrams analyzed in this work.</caption><thead><tr><th>x</th><th>y</th><th>Redshift range</th><th>Best-fit bins</th><th>Best-fit pymc</th></tr></thead><tbody><tr><td>z</td><td>log Σ<sub>SFR</sub>/(M<sub>☉</sub>/yr/kpc<sup>2</sup>)</td><td>All</td><td>m = 0.21 ± 0.03<br>q = -0.7 ± 0.2</td><td>m = 0.16 ± 0.06<br>q = -0.3 ± 0.4</td></tr><tr><td>log M<sub>⋆</sub>/M<sub>☉</sub></td><td>log SFR/(M<sub>☉</sub>/yr)</td><td>Low-z</td><td>m = 0.65 ± 0.09<br>q = -4.9 ± 0.8</td><td>m = 0.63 ± 0.07<br>q = -4.8 ± 0.7</td></tr><tr><td>log M<sub>⋆</sub>/M<sub>☉</sub></td><td>log SFR/(M<sub>☉</sub>/yr)</td><td>High-z</td><td>m = 0.4 ± 0.3<br>q = -3 ± 2</td><td>m = 0.63 ± 0.10<br>q = -4.8 ± 1.0</td></tr><tr><td>log M<sub>⋆</sub>/M<sub>☉</sub></td><td>log Σ<sub>SFR</sub>/(M<sub>☉</sub>/yr/kpc<sup>2</sup>)</td><td>Low-z</td><td>m = 0.17 ± 0.13<br>q = -1.0 ± 0.9</td><td>m = 0.2 ± 0.1<br>q = -1.3 ± 0.8</td></tr><tr><td>log M<sub>⋆</sub>/M<sub>☉</sub></td><td>log Σ<sub>SFR</sub>/(M<sub>☉</sub>/yr/kpc<sup>2</sup>)</td><td>High-z</td><td>m = 0.2 ± 0.4<br>q = -1 ± 3</td><td>m = 0.16 ± 0.11<br>q = -0.3 ± 1.1</td></tr></tbody></table>

**Notes.** The first two columns indicate the x and y axis quantities, while the last two columns highlight the best fit coefficients of the linear relations ($y = m \times x + q$) obtained with two different approaches, as described in the text. The third column indicates whether the relation is derived with the subset at lower redshift ($4 < z < 6$), at higher redshift ($6 < z < 10$), or in the entire z range. Only galaxies above the 90% mass completeness limits are used for the fit.

ness limits between lower and higher redshifts, as discussed in Section 2.5, and might reflect the rise observed at $z \ge 7.5$ of the minimum stellar mass and SFR probed by our NIRSpec sample (Fig. 1). We also remark that in that redshift regime we are limited by the poor statistics, and a larger sample is required to better constrain the $\Sigma_{SFR}$ evolution at $z > 8$.

Both the normalization and the slope of our $\Sigma_{SFR}$ – redshift relation are consistent with most of the previous observations. In particular, our Bayesian best-fit line falls exactly between the observed median relations estimated by Shibuya et al. (2015) up to redshift 8 in two different bins of SFR ($1 < \text{SFR}/M_{\odot}/\text{yr} < 10$ and $10 < \text{SFR}/M_{\odot}/\text{yr} < 100$). Compared to these results, our findings indicate that the increase of $\Sigma_{SFR}$ can be extended up to at least redshift 10. The median $\Sigma_{SFR}$ between redshift 4 and 9 are overall consistent with those calculated by Morishita et al. (2024) with a sample similar to our study in terms of stellar masses and SFRs probed (fuchsia empty stars in Fig. 2). Moreover, NIRSpec galaxies ranging $8 < z < 10$ have $\Sigma_{SFR}$ comparable to those recently derived in photometrically selected systems observed by JADES at much higher redshifts (10–13) by Robertson et al. (2023). Finally, our results are slightly higher, but still marginally consistent, than the range of $\Sigma_{SFR}$ (1–20

$M_{\odot}/\text{yr/kpc}^2$) measured by Holwerda et al. (2015) for 6 galaxies at redshifts $9 < z < 10$ identified with HST in the CANDELS survey, even though these are intrinsically brighter systems compared to those analyzed in this work.

### 3.3. Main sequence of star formation at $z > 4$

We go on to focus on the fundamental scaling relations describing the evolution of the galaxy population. The most important is that involving the star formation rate and the stellar mass, which show a tight correlation at all redshifts, known as the star-forming main sequence (SFMS, Noeske et al. 2007, Daddi et al. 2007, Santini et al. 2017). This relation is fundamental to understand the process of gas conversion into stars and subsequent build-up of stellar mass across cosmic time. While the slope is relatively constant with redshift, ranging from ~0.6 to ~1, depending on the sample selection and specific tracer used for the SFR, all studies agree that the normalization increases monotonically from lower to higher redshift up to at least $z \sim 7$.

In Fig. 3, we show the $M_{\star}$ versus SFR diagram for our selected spectroscopic sample. In order to check and better appreciate variations with redshift, we divide the sample in two
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Star-forming main sequence of galaxies selected in this work in the redshift range $4 < z < 6$. The best-fit relations are derived in the mass range where we have 90% mass completeness. For comparison, the high redshift and low-mass extrapolation of the MS relations by Schreiber et al. (2015) and Speagle et al. (2014) (obtained at lower redshifts and higher $M_*$) are shown with dashed orange and dotted fuchsia lines, respectively. Bottom panel shows the same as the top, but for the redshift range $6 < z < 10$. The Bayesian best-fit relation obtained in the lower redshift bin is included as a pink continuous line for comparison.</figcaption></figure>

redshift bins, that is, $4 < z < 6$ and $6 < z < 10$ (top and bottom panels, respectively). In both cases, we further divide the sample in different bins of stellar masses from $10^7.5$ to $10^{10.5} M_\odot$ in steps of 0.5 dex (slightly larger in the first bin) when performing the first approach of the analysis explained in Section 3.1. We perform the fit considering only the galaxies above the mass completeness limits estimated in Section 2.5, that is, $\log_{10}(M_*/M_\odot)$ above 7.8 in the low redshift bin, and $> 8.4$ in the high redshift bin. With all the fitting methods, we find in both redshift ranges a tight and significant correlation between $M_*$ and the SFR, indicating that the SFMS is in place up to the highest redshifts ($z \sim 10$) explored by this work. Adopting the Bayesian regression as the reference method, the SFMS in the two redshift bins have similar slopes ($0.63 \pm 0.07$ at $z < 6$ and $0.63 \pm 0.10$ at $z > 6$), and the normalization of the relation is not significantly different between the two redshift bins in the $M_*$ range from $\sim 10^8$ to $\sim 10^{10.5} M_\odot$.

Recently, Cole et al. (2023) studied the SFMS using photometrically selected galaxies from the CEERS survey at redshifts 4.5–12, with stellar masses and SFRs estimated through SED fitting. If we focus on the stellar mass ranges where we are both mass complete at 90% (i.e., $\log_{10}(M_*/M_\odot) > 8$ at lower redshifts, and $> 8.6$ at higher redshifts), we find consistent relations with similar slopes and normalizations. This also confirms that the spectroscopic selection is not significantly biased compared to a pure photometrical selection. We have also performed a fit by including all the galaxies down to the lowest measured stellar masses, but we obtain best-fit relations that are not significantly different from those estimated above. This suggests that, even though we are incomplete at lower stellar masses, this effect might not be large.

Cole et al. (2023) also found that 10Myr-averaged SFRs are significantly more scattered around the median relation compared to 100 Myr-averaged values by ~0.3 dex (see their Fig. 7), which these authors interpret as evidence of the stochasticity of star formation. We have here the possibility to analyze the scatter of the main sequence relation using the Balmer lines as tracers of very recent star formation. To this aim, we computed the intrinsic scatter of the $M_*$ – SFR relations by subtracting in quadrature, for all the redshift and $M_*$ bins, the scatter due to the measurement uncertainty from the total observed scatter, as done in Huang et al. (2023). We find that the intrinsic scatter $\log(\sigma_{\text{intrinsic,Ha}}/(M_\odot/\text{yr}))$ ranges between 0.3 and 0.6 in both the low and high redshift bins, with a median of 0.5 across the entire stellar mass range. This result is consistent with the range of $\sigma_{10\text{Myr}}$ reported by Cole et al. (2023), and ~0.4 dex higher than $\sigma_{100\text{Myr}}$ found by the same study, thus providing a spectroscopic confirmation of the relatively larger scatter of the main sequence when SFRs are averaged over short timescales of $\ll 100$ Myr.

Both the slope and normalization of our relations are generally in agreement with the extrapolations at redshift 5 and 7.5 (and down to a stellar mass of $10^7 M_\odot$) of the SFMS relations derived for more massive and lower redshift star-forming galaxies by Speagle et al. (2014) and Schreiber et al. (2015), with a slightly better consistency with the shallower slope found by Speagle et al. (2014). We note that the evolution of SFR at fixed stellar mass expected in this redshift range (~0.15 dex) is smaller than the intrinsic scatter of the relation and too small to be appreciated with our dataset, given the uncertainties obtained in the fit.

### 3.4. The main sequence of SFR surface density

The galaxy integrated $\Sigma_{\text{SFR}}$ has been proposed as an alternative and more physical measurement of the current SF activity of galaxies, the reason being a closer connection to the gas density and the efficiency of stellar feedback compared to the specific SFR. For example, the incidence of outflows has been shown to correlate better with the $\Sigma_{\text{SFR}}$ rather than with the offset from the standard SFMS (Förster Schreiber et al. 2019). Moreover, $\Sigma_{\text{SFR}}$ also incorporates, by construction, information about the galaxy structure and morphology. For this reason, an alternative main sequence relation has been defined as the main locus occupied by star-forming galaxies in the $\Sigma_{\text{SFR}} - M_*$ plane. Compared to the sSFR, the $\Sigma_{\text{SFR}}$ properties of galaxies are more uniform at varying stellar masses, such that $\Sigma_{\text{SFR}}$ is approximately constant across the main sequence, as found by Schiminovich et al. (2007) and Salim et al. (2023) for relatively low redshift galaxies ($z \le 1$). However, there could also be deviations from this flat relation, indicating additional (or lower) star formation activity associated with different galaxy structures. For example, moderately positive correlations between $\Sigma_{\text{SFR}}$ and $M_*$ are observed at
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Σ<sub>SFR</sub> main sequence for our selected galaxies in the lower redshift bin (4 < z < 6, top) and higher redshift bin (6 < z < 10, bottom). The lines, shaded areas, and markers are the same as in Fig. 2. The average location of massive (M<sub>⋆</sub> > 10<sup>9</sup> M<sub>☉</sub>) star-forming galaxies at redshift ~0.5 presented by Salim et al. (2023) is highlighted with a cyan dashed line.</figcaption></figure>

M_★ > 10^9 M_☉ by Salim et al. (2023) for star-forming galaxies from z = 0 to z = 2, which are interpreted as the build-up of a central bulge component in those galaxies.

We now explore the Σ<sub>SFR</sub> versus M<sub>★</sub> diagram for the NIR-Spec galaxies, divided in two redshift bins (4 < z < 6 and 6 < z < 10), as done for the classical SFMS. In the fit we consider only the galaxies above the mass completeness limits defined before for each redshift. Our results are shown in Fig. 4. We find that in both redshift bins there is a very mild correlation, with Σ<sub>SFR</sub> slightly increasing in more massive galaxies: the best-fit slopes obtained with the Bayesian regression method are of 0.19 ± 0.09 and 0.30 ± 0.15, respectively, in the low and high redshift bins. This tells us that the data points are still consistent (within 3σ) with a constant Σ<sub>SFR</sub> as a function of M<sub>★</sub> in the entire redshift of our work. Considering the large range spanned by Σ<sub>SFR</sub> (by more than three orders of magnitude), it is indeed remarkable to observe a median variation of only less than 0.5 dex across more than ~3 orders of magnitude of stellar mass, from log M<sub>★</sub> = 7.5 to ~10.5. For the higher redshift sample we also find a slightly higher normalization of the Σ<sub>SFR</sub>-main sequence by ~0.5 dex (at M<sub>★</sub> = 10<sup>9</sup> M<sub>☉</sub>) compared to the subset at z < 6.

The very mild increase of Σ<sub>SFR</sub> with M<sub>★</sub> is in agreement with previous observations of star-forming galaxies at comparable masses but at lower redshifts (Förster Schreiber et al. 2019; Schiminovich et al. 2007; Salim et al. 2023). In such comparisons, we do not find clear evidence of a sudden or significant upturn of the relation in the entire mass range that we explore. However, the normalization of the Σ<sub>SFR</sub>-main sequence is significantly higher compared to lower redshift studies. The median Σ<sub>SFR</sub> of our galaxies at z ~ 5 is of ~2.5 M<sub>☉</sub>/yr/kpc<sup>2</sup> (log Σ<sub>SFR</sub> = 0.4). At fixed stellar mass, and focusing on the overlapping stellar mass region (8 < log M<sub>★</sub>/M<sub>☉</sub> < 10), our Σ<sub>SFR</sub> values are higher compared to the local Universe by approximately 2 dex at z ~ 5 and 2.5 dex at z ~ 10, indicating a strong evolution of Σ<sub>SFR</sub> with cosmic time. This is consistent with the average Σ<sub>SFR</sub> evolution found over this cosmic epoch by other observations (e.g., Shibuya et al. 2015) and simulations (e.g., Sharma et al. 2017). It can be explained as the combined effect of the higher average SFRs (by ~1 dex, Speagle et al. 2014) and more compact sizes (r<sub>e</sub> ~ 0.5 dex lower, Ward et al. 2024) of galaxies at z ~ 5 compared to the local Universe, at a fixed stellar mass of ~10<sup>9.5</sup> M<sub>☉</sub>.

**4. Discussion**

In this section, we further analyze our results to better understand how the fundamental quantity Σ<sub>SFR</sub> is related to other physical properties of the galaxies during and immediately after the epoch of reionization.

**4.1. Comparison to predictions of theoretical models**

We compare our observations to some theoretical models that have been introduced to explain galaxy evolution across cosmic epochs. Focusing on the redshift evolution of Σ<sub>SFR</sub>, we first consider the empirical model of Naidu et al. (2020). As mentioned in Section 1, the main assumption of this model is that f<sub>esc</sub> is solely dependent on Σ<sub>SFR</sub> as f<sub>esc</sub> = a × Σ<sup>b</sup><sub>SFR</sub>, with a = 1.6 ± 0.3 and b = 0.4 ± 0.1. As shown in Fig. 2, this model yields an evolution of Σ<sub>SFR</sub> up to redshift ~8 (lime curve) that is fully in agreement with our best-fit relation.

As an additional test, we compare our observations to the physical model presented by Ferrara et al. (2023) and Ferrara (2024). In detail, we take the redshift evolution of the SFR predicted by their model, and derive the Σ<sub>SFR</sub> by combining it with the empirical size evolution found by Morishita et al. (2024). At the median stellar mass of the sample plotted in Fig. 2 (log M<sub>★</sub>/M<sub>☉</sub> = 8.6), we obtain a slightly increasing trend of Σ<sub>SFR</sub>, highlighted with a dash-dotted orange line, whose slope is remarkably consistent with our data and with the model of Naidu et al. (2020). Even though it has a slightly higher normalization of ~+0.2 dex, this is well below the uncertainties of our best-fit relation. We note that the prediction is very sensitive not only to the stellar mass considered (lower M<sub>★</sub> have systematically lower Σ<sub>SFR</sub>), but especially to the effective radius. Indeed, considering the uncertainties on the best-fit redshift-size relation of Morishita et al. (2024), we would obtain an uncertainty on the normalization of the predicted Σ<sub>SFR</sub> evolution of ~0.8 dex, higher than the median absolute deviation (MAD) of our sample.

Finally, we also took the CEERS mock galaxy catalog (Yung et al. 2022) for comparison. This catalog covers an area overlapping with the observed EGS field, and contains galaxies over 0 < z < 10. The galaxies in the mock lightcone were simulated with the Santa Cruz SAM for galaxy formation (Somerville & Davé 2015; Somerville et al. 2021). The free parameters in the model were calibrated to reproduce a set of
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** $\Sigma_{SFR}$ as a function of the O32 index for our selected sample in the entire range $4 < z < 10$. The lines, shaded areas, and markers are the same as in Fig. 2. The diamonds colored in orange highlight our strong LyC leaking candidates, according to the criteria defined by Flury et al. (2022).</figcaption></figure>

galaxy properties observed at $z \sim 0$ and have been shown to well-reproduce the observed evolution in high-redshift ($z \ge 4$) one-point distribution functions of $M_{UV}$, $M_*$, and SFR (Yung et al. 2019a,b). The effective radii of these simulated galaxies are determined as described in Brennan et al. (2015) and are shown to be in good agreement with CEERS observations at $3 < z < 6$ (Kartaltepe et al. 2023).

From the mock lightcone, we extracted a random sample of galaxies with a redshift and a stellar mass distribution similar to our observed sample, and then apply a magnitude cut as $F277W < 29.15$ to match the depth of the CEERS survey (Finkelstein et al. 2024), where most of our galaxies come from. With this selection method, fitting the median $\Sigma_{SFR}$ in five bins of redshift as done for the real data, we obtained an increasing trend of $\Sigma_{SFR}$ with redshift that is very similar to the observational result. In addition, the SAM best-fit relation has a slightly lower normalization than the observed one by $\sim 0.2$ dex. As in the previous case, this is due to predictions on galaxy sizes, which are slightly larger on average compared to our values and to the size evolution obtained by Morishita et al. (2024). Despite this, the difference in $\Sigma_{SFR}$ normalization is rather small and below the $1\sigma$ uncertainties. Therefore, the SCSAM yields a good representation of the observed properties of galaxies in the entire redshift range from 4 to $\sim 10$. A similar good consistency between SCSAM model predictions and observations is obtained also for all the other diagrams studied in Section 3 and in Section 4. The models analyzed in this section also agree on the strong evolution of $\Sigma_{SFR}$ from $z=0$ to the reionization epoch mentioned in Section 3.4. To conclude, our results are consistent not only with previous observations, but also with a variety of theoretical models that are commonly adopted in the literature.

## 4.2. SFR surface density and ionization properties

The SFR surface densities of galaxies have been shown recently to be tightly related and to have a key influence on their ionizing properties (Reddy et al. 2022). We know from simulations, galaxy models, and observations, that high redshift galaxies have higher ionization parameters (at fixed stellar mass) compared to local analogs (e.g., Brinchmann et al. 2008; Nakajima et al. 2013; Steidel et al. 2014; Kewley et al. 2015; Shapley et al. 2015; Hirschmann et al. 2017, 2023; Euclid Collaboration 2024). This redshift evolution of ionization conditions is usually explained as a combined effect of lower metallicity, harder stellar ionizing spectra, and higher gas densities at $z \gg 0$. The latter increasing trend, also found with the electron density $n_e$ (Isobe et al. 2023), which is more easily measurable from optical emission lines, reflects the denser, more compact molecular clouds, and ultimately also higher values of $\Sigma_{SFR}$. Significant correlations were indeed found observationally between $n_e$ and $\Sigma_{SFR}$ (Shimakawa et al. 2015; Jiang et al. 2019; Reddy et al. 2023), confirming the close relation between the two quantities, and the important role of $\Sigma_{SFR}$ in the evolution of the ISM properties across cosmic time (see also Papovich et al. 2022).

We test here the connection between $\Sigma_{SFR}$ and the ionization conditions, comparing the measured values of $\Sigma_{SFR}$ to the O32 index defined in Section 2.2, which is usually adopted as an observational proxy for the ionization parameter. The result is shown in Fig. 5. The O32 values, represented on the y-axis, range between -0.3 and 1.3, with a median of $\sim 0.6$. We observe an increase of O32 with $\Sigma_{SFR}$, with a slope of $1.5 \pm 0.2$, indicating that the correlation is statistically significant at more than $3\sigma$ level. We notice also that the galaxies in the top right part of the diagram, with $\Sigma_{SFR} \ge 1$ and $O32 \ge 0.7$, have similarly high $\Sigma_{SFR}$ and O32 to photometrically selected extreme emission line galaxies (EELGs) identified in the CEERS field by Llerena et al. (2024). We refer to that paper for a more detailed discussion of those galaxies. Overall, these findings corroborate previous results, supporting the close relation between $\Sigma_{SFR}$ and the ionization properties of galaxies. We also note that excluding the small subset of galaxies at $z > 8$ with very high $\Sigma_{SFR}$ does not significantly change this global picture.

Finally, seven galaxies observed at high resolution in GLASS have a S/N of $[O\text{II}]\lambda3727\AA > 3$ and a physically meaningful flux ratio $[O\text{II}]_{3729}/[O\text{II}]_{3726}$ in the range 0.3–1.5 (Osterbrock 1989). From this flux ratio, we derived for those galaxies the electron density $n_e$, using the python package PYNEB (Luridiana et al. 2015) and assuming an electron temperature $T = 10000$ K, as above. We find that most of the galaxies cluster within a range of densities going from 200 to $400 \text{ cm}^{-2}$ (median $n_e = 280 \text{ cm}^{-2}$). The poor statistics does not allow us to test possible correlations between $n_e$ and $\Sigma_{SFR}$, even though we notice that the galaxy ID 100003 with the highest density in the sample ($> 10^4 \text{ cm}^{-2}$, the only one outside the range 200–400 $\text{cm}^{-2}$), also has the highest $\Sigma_{SFR}$ ($\sim 10^2 M_{\odot}/\text{yr}/\text{kpc}^2$). We defer a more detailed analysis to a future work with a larger sample.

## 4.3. SFR surface density and escape of Lyman continuum photons

Enhanced values of the O32 index and $\Sigma_{SFR}$ are usually associated with increased fractions of ionizing photons that escape from the galaxies. The connection between $\Sigma_{SFR}$ (or O32) and $f_{esc}$ is found from multiple observations and with different approaches. Lyman continuum leaking sources (LyC), identified through direct or indirect methods, are found to have $\Sigma_{SFR}$ significantly higher than the average value of the star-forming population at a given redshift (e.g., Cen 2020). In the other direction, galaxies with enhanced $\Sigma_{SFR}$ show larger fractions of LyC leakers. For example, Heckman (2001) and Heckman (2001).
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Diagram comparing $\Sigma_{SFR}$ and $f_{esc}$ for galaxies at $4 < z < 10$ with an estimate of $f_{esc}$(top). The lines, shaded areas, and markers are the same as in Fig. 2. sSFR versus $f_{esc}$ for this NIRSpec sample is shown in the bottom panel. The orange diamonds represent strong LyC leaking candidates according to Flury et al. (2022), as highlighted in Fig. 5. The two galaxies for which we detect outflow signatures in their H-grating spectra from GLASS are shown as big darkred stars. The sample representative error of 0.3 dex on $f_{esc}$ is added below the legend. We include with blue empty squares the LzLCS sample of 66 galaxies at $0.2 < z < 0.4$, with measurements (including direct estimates of $f_{esc}$) performed by Flury et al. (2022).</figcaption></figure>

(2002) claimed that galaxies with $\Sigma_{SFR}$ above $0.1 M_{\odot}/yr/kpc^2$ are capable to launch strong star-forming driven winds, which might be responsible for creating channels though which LyC photons can leak. Steidel et al. (2018), based on a sample of LBGs at $z \sim 3$, suggested that $\Sigma_{SFR}$ and $f_{esc}$ are highly correlated quantities. Similarly, the O32 index is considered as a valid indicator of LyC leaking, according to Izotov et al. (2020). Also, Flury et al. (2022) found that 40% of local galaxies with $\Sigma_{SFR} > 10 M_{\odot}/yr/kpc^2$ and O32 higher than 0.7 are LyC leakers; hence, they adopted these two conditions at higher redshifts to define a source as a strong LyC leaker candidate.

Following these observational findings, hydrodynamical simulations and models have started to implement the effects of spatially concentrated star formation in the form of turbulence and stellar-driven winds (Ma et al. 2016; Sharma et al. 2016), showing that an increased $\Sigma_{SFR}$ is responsible for launching

galaxy-wide outflows and for carving ionizing channels in the ISM where the Lyman continuum radiation of young massive stars and supernovae is free to propagate outside. According to the model of Naidu et al. (2020), $f_{esc}$ reaches unity when $\Sigma_{SFR}$ is one-third of the maximum $\Sigma_{SFR}$ (= $1000 M_{\odot}/yr/kpc^2$) that can be sustained without radiation pressure instabilities (Hopkins et al. 2010).

Triggered by this extensive observational and theoretical background, we used our dataset to explore the relation between $\Sigma_{SFR}$ and $f_{esc}$ at redshifts $4 < z < 10$. We find that at $z \ge 4$, it is almost impossible to directly detect the Lyman continuum emission due to the very low IGM transmission (Vanzella et al. 2015), and at $z > 5$ it becomes virtually impossible. For this reason, we determined $f_{esc}$ in an indirect way using the empirical relation presented in Mascia et al. (2023), which provides an estimation of the escape fraction using three galaxy parameters: the UV beta slope, the physical size $r_e$, and the O32 index; we note that these are among the properties that are better correlated with $f_{esc}$. This relation is calibrated on the low-redshift ($0.2 < z < 0.4$) Lyman Continuum survey (LzLCS), for which Flury et al. (2022) performed direct measurements of the Lyman continuum flux and $f_{esc}$ from imaging observations.

In Fig. 6 (top), we show $\Sigma_{SFR}$ as a function of $f_{esc}$. Using our two analysis approaches, we find in all cases a significant correlation between these two parameters, according to which galaxies with higher $\Sigma_{SFR}$ have systematically higher $f_{esc}$ than galaxies with lower $\Sigma_{SFR}$. In detail, $f_{esc}$ increases from log $\Sigma_{SFR} = -1$ to 3, with a best-fit slope of $0.37 \pm 0.05$ by fitting the median values in seven bins of $\Sigma_{SFR}$ defined with the following grid [-1, -0.5, 0, 0.5, 1, 1.5, 2, 3]. These results do not significantly change if we exclude from the analysis the galaxies at $z > 8$ that are responsible for the $\Sigma_{SFR}$ enhancement in the last redshift bin observed in Fig. 2. This indicates that galaxies with more concentrated SFR have the conditions that should facilitate the escape of LyC photons. A subset of ~20% of our galaxy sample have $\Sigma_{SFR}$ and O32 satisfying the conditions established by Flury et al. (2022) for being strong LyC leaker candidates, showing similar properties to the LyC leakers studied at redshift ~0.3. These systems, indicated with orange diamonds in Fig. 6 and in Fig. 5, tend to have higher $f_{esc}$ estimates compared to the average galaxy population at these redshifts.

We compare our observations to the direct measurements of $f_{esc}$ and physical properties estimated by Flury et al. (2022) for the LzLCS galaxy sample (shown as blue empty squares) in Fig. 6. We find a similar slope of the best-fit relation, even though we derive a higher normalization by ~0.5 dex. The large difference in normalization of the $\Sigma_{SFR} - f_{esc}$ relation mostly depends on the fact that, at fixed $\Sigma_{SFR}$, high-redshift galaxies have a lower UV slope on average compared to the low redshift systems analyzed in Flury et al. (2022); hence, they have a lower dust attenuation and a lower metallicity (see Calabrò et al. 2021 for the relation between attenuation and metallicity). This translates into an increased capability to form channels through the ISM for the escape of ionizing radiation. Indeed, the normalization offset almost disappears when considering a subset from the LzLCS with a median UV slope similar to the NIRSpec galaxies ($\beta_{median} \approx -2.1$). A minor part of the offset is also due to our sample probing slightly lower stellar masses (by ~0.3 dex) than the sample of Flury et al. (2022).

We also find that the slope of our correlation is consistent with the $f_{esc} \propto \Sigma_{SFR}^{0.4}$ relation predicted by the Naidu et al. (2020) model. We find a lower normalization of our relation by ~0.3 dex compared to their fiducial model with $a = 1.6$, possibly
suggesting a lower value of this parameter. For completeness, we also report in Fig. 6 (top) the model with $a = 0.8$, which is not ruled out in their work (as shown in their Fig. 6) and whose prediction is much closer to our best-fit relation. However, their $f_{esc}$ versus $\Sigma_{SFR}$ relation relies on measurements by Steidel et al. (2018), which are based on a sample of star-forming galaxies at redshift $z \sim 3$. The difference may thus originate from the different redshifts probed by our works, hence from the effects discussed above.

A big caveat to the above analysis is that the two quantities on the x and y axis are not completely independent: $\Sigma_{SFR}$ correlates with the size, while $f_{esc}$ is strongly dependent on the size by construction. While this is unavoidable, we also remark that, the $f_{esc}$ estimation is based also on other galaxy properties and not just on the $r_e$. However, in order to provide a more independent measurement, we also compare the estimated values of $f_{esc}$ to the sSFRs, which do not directly depend on the physical extension of the galaxies (Fig. 6, bottom panel). We find that the sSFR is still correlated to $f_{esc}$, although with a shallower slope and lower significance compared to the same relation with $\Sigma_{SFR}$. Also, in this case the slope of our best-fit relation is similar to what is found for the LzLCS galaxy sample at $z \sim 0.3$, suggesting that the level of star formation activity in a galaxy (regardless of the normalization method) has an impact on the escape of ionizing radiation. Overall, we are aware that all these quantities (i.e., $f_{esc}$, $\Sigma_{SFR}$, sSFR, size, UV slope, and O32) are all interrelated somehow and not fully independent, as many of them depend on the same physical processes. This tells us that the best candidates as LyC leaking galaxies have rather peculiar and similar properties, being characterized by blue UV slopes, compact sizes, high O32, and also enhanced star formation activity in the form of $\Sigma_{SFR}$ and sSFR.

### 4.4. Connecting star formation and outflows

We have seen in the previous section that $\Sigma_{SFR}$ is tightly related to the ionizing properties of galaxies and their ability to spread ionizing photons in their surrounding medium. The underlying physical mechanism through which $\Sigma_{SFR}$ impacts on $f_{esc}$ is that galaxies with more compact and enhanced star formation activity facilitate the launch of galaxy scale outflows through stellar feedback, which includes radiation from young stars, stellar winds, and supernova explosions. Indeed, several studies have found higher gas outflow velocities and mass loading factors in galaxies with higher $\Sigma_{SFR}$ (Calabrò et al. 2022a; Llerena et al. 2023), corroborating the close relation between $\Sigma_{SFR}$ and outflows. From the theoretical point of view, Sharma et al. (2016) and Sharma et al. (2017) predict that outflows are ubiquitous in galaxies with $\Sigma_{SFR} \ge 10 M_\odot/\text{yr/kpc}^2$. The final connection between outflows and $f_{esc}$ is then supported by a variety of works from the literature, both observational and theoretical (Heckman et al. 2011; Amorín et al. 2012; Borthakur et al. 2014; Sharma et al. 2017; Hogarth et al. 2020), according to which outflows can clear pathways in the surrounding neutral gas, favoring the escape of ionizing photons. Recently, Amorín et al. (2024) found a clear correlation between the velocity width of the [O III] $\lambda5007$ Å and $H\alpha$ broad-line wings and $f_{esc}$ for a sample of 20 Lyman continuum emitters at $z \sim 0.3$, with the stronger emitters also having the highest values of $\Sigma_{SFR}$. In this subsection, we investigate spectral outflow signatures in our sample and study how they relate to the host galaxy properties. We also test the impact of outflows in an indirect way through the effect that they may have on the dust attenuation.

#### 4.4.1. Exploring spectral outflow signatures in our sample

To study the presence of outflows of ionized gas, we use the GLASS high resolution grating spectra and perform a single Gaussian and a double-component (narrow+broad) Gaussian fit to the $H\alpha$ and [O III] $\lambda5007$ Å lines (when available), which are the lines with the highest S/N in this redshift range, thus the most suited for exploring low S/N outflow features. In detail, $H\alpha$ was fitted alone, as the contribution from [N II] $\lambda6583$ Å can be neglected at first order, while [O III] $\lambda5007$ Å was fitted together with $H\beta$ and [O II] $\lambda4959$ Å as explained in Section 2.2. In the double Gaussian fit, we add a second broad component, allowing the intrinsic velocity width $\sigma$ of the narrow component to vary between 10 and 250 km/s, and $\sigma$ of the broad component from 100 to 800 km/s, which is the limit expected for star-forming driven outflows (Calabrò et al. 2022a). We constrain the velocity shift between the two Gaussian centroids to <300 km/s, of the same order of the spectral resolution. We do not set any constraints on the flux ratio of the two components. In case of the $H\beta$ + [O III] triplet, a broad component is added to all the three lines. Finally, we use the Akaike Information Criterion (AIC) criterion (Akaike 1974) to decide if the double Gaussian provides a significantly better fit to the line profile compared to a single Gaussian fit. In particular, we consider the double component as the final solution if the difference ($AIC_{\text{single}} - AIC_{\text{double}}$) = $\Delta_{\text{AIC}}$ is > 10 (Burnham & Anderson 2004; Bosch et al. 2019), and if both the narrow and broad best-fit Gaussians are detected with a flux S/N > 3. The above threshold ensures that the double component model is 100 times more probable than the single component model to minimize the information loss.

With the above criteria, we find that a broad component is required to fit [O III] and $H\alpha$ of two galaxies: ID 160133 and ID 110000 (Fig. 7). Both galaxies have $\Delta_{\text{AIC}} > 15$ in [O III] (the highest S/N line), indicating a very strong case for the double component model. For the galaxy ID 110000, $H\alpha$ has a $\Delta_{\text{AIC}}$ slightly lower than 10 (~8), but the double Gaussian is still highly favored (50 times more probable than a single component). We primarily interpret these components as an outflow, with $\sigma_{\text{broad}} = 250 \pm 4$ km/s and 200 ± 10 km/s, respectively. We measure a broad to narrow line flux ratio of 0.16 (0.3) and a broad to narrow $\sigma$ ratio of 3.4 (2.8). These galaxies lie at redshifts 4.02 and 5.76, respectively, and both have a modest log $\Sigma_{SFR}/(M_\odot \text{yr}^{-1} \text{kpc}^{-2})$ of 0.6 ± 0.2 (0.4 ± 0.2), comparable to the average galaxy population at their redshifts, and log SFR/(yr$^{-1}$) of -7.6 ± 0.2 (-8.5 ± 0.3). These systems are highlighted as big darkred stars in Fig. 6. The indirect estimation of $f_{esc}$ for these two galaxies yields a value of ≈0.11 (≈0.04), which is slightly above (slightly below) the median $f_{esc}$ of galaxies at similar $\Sigma_{SFR}$. We also note that these outflow velocities would correspond to strong leakers ($f_{esc} \ge 4-10\%$) if we assume the $\sigma_{\text{broad}}-f_{esc}$ correlation of Amorín et al. (2024). Intriguingly, both galaxies with outflow signatures show low surface brightness clumpy structures around the main cores from high-resolution JWST images. This makes the scenario of gas accretion an interesting, alternative possibility that cannot be excluded with our spectroscopic data.

To better understand these results, we compare them to those obtained with a comparable dataset by Carniani et al. (2024), which includes 52 galaxies with $R= 2700$ spectra from the JWST-JADES program at redshifts, stellar masses, and SFRs similar to our work. They detect broad additional components in $H\alpha$ or [O III] in a fraction of galaxies going from 15% to 40%, which they interpret as evidence of outflows. They find that the outflow incidence increases with SFR but is rather stable
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Single Gaussian (1G) and double Gaussian (2G, narrow + broad) best fits to the [O III] λ5007 Å+Hβ and Hα lines, as indicated in each panel label, for the two galaxies (ID GLASS 160133 and 110000) in which we detect a significant outflow signature in both emission lines (upper and lower four panels, respectively). In each panel, the upper part shows the observed spectrum with the best fit line in red (green) line for the narrow (broad) component in the double Gaussian fit, the middle part represents the error spectrum, while the bottom part is the residual after subtracting the best fit result. The spectra are in the observed frame wavelength.</figcaption></figure>

(at the level of 20–30%) as a function of M* and sSFR. Limiting the comparison to the redshift range between 4 and ~6 (where most of their galaxies lie), and assuming the simple interpretation (also favored by their work) of an outflow scenario,

we derive in GLASS an outflow incidence of ~20%^{32%}_{9%} (2 out of 11 galaxies with detected [O III]), with uncertainties derived following Gehrels (1986). This result is comparable to Carniani et al. (2024). We also notice that our flux and σ ratios between the
narrow and broad Gaussian components fall within their observed range. All this suggests that we are probably probing the same physical phenomenon. Interestingly, we do not observe instead any outflow signatures in the 6 GLASS sources at $z > 6$ (upper limit of 26%). In this case, we cannot make comparisons as this redshift range is poorly sample by the other work. Within our sample, the number of high redshift sources with high resolution spectra is so low that the outflow incidence at $z > 6$ is still compatible with that derived at lower redshift. However, we can exclude from this analysis that outflows become more important at higher redshifts.

Furthermore, we do not find significant evidence of outflows in individual galaxies among the CEERS sample observed at medium resolution ($R \approx 1000$). To understand whether this is expected, and to test the effect of resolution on the outflow detectability, we perform a set of Monte Carlo simulations. In detail, we simulate spectra at different redshifts injecting an emission line at a specific wavelength (0.53 µm rest frame) assuming different line fluxes to produce a variety of S/N of the emission feature (from 3 to ~60). We also set a broad to narrow intrinsic flux ratio in the range between 0.1 and 0.7, and different $\sigma$ ratios ranging 2 to 12, with $\sigma_{narrow,intrinsic}$ fixed to 50 km/s. We vary the intrinsic broad line centroid within a resolution element. We then perform a single Gaussian and double Gaussian (narrow + broad) fit as done for the observed data, imposing the same requirements for taking the double component as the final best fit. Running these simulations (100 for each configuration) we find that the recovery rate is of the order of 2% or less for the typical range of outflow parameters expected from our works and from previous findings. This rate is not strongly dependent on the specific assumptions, except for the most extreme cases such as width and flux ratios $\ge$ 5 and $\ge$ 0.5, respectively. Therefore, the lack of significant outflow signatures in the CEERS sample is likely due to their lower spectral resolution.

Increasing the S/N of the emission lines, for example through spectral stacking, can slightly increase the outflow recovery rate at median resolution and the possibility of detecting faint broad components in case of high outflow velocities. We thus stack the entire dataset selected in this work with H- or M-grating observations. We first downgrade by convolution the GLASS high resolution spectra to $R = 1000$, convert the spectra to rest-frame using the spectroscopic redshifts determined in Section 2.2, resample them to a common wavelength grid, and normalize them to the median flux estimated around the specific lines analyzed. Then, we apply for each pixel a median stacking with $5\sigma$ clipping, and finally perform the line measurements as done for individual galaxies. As shown by our simulations, this stacking approach tends to broaden a line systematically (regardless of the S/N or intrinsic line width) by ~50 km/s, as due to the impact of the spectral resolution on the redshift estimation uncertainty when converting to rest-frame.

Stacking the spectra around the [O III] $\lambda5007$ Å line, which is the highest S/N line available for the largest number of galaxies in our sample, and fitting the $H\beta$+[O III] triplet with a single and double Gaussian as described before, we find that the double component model is not strongly favored according to the previous condition ($\Delta_{AIC} = 3$), thus there is no significant evidence of outflow. Given the lower S/N, we are not able to appreciate significant differences in the fit or to detect any outflow signatures if we divide the sample in two or more bins of redshift, sSFR, or $\Sigma_{SFR}$.

To conclude, we find outflow incidence rates in individual galaxies observed by GLASS that are consistent to those obtained by similar works in the same redshift range. As highlighted by

<figure><img src="image_9.png" /><figcaption>Fig. 8. Redshift evolution of $\Sigma_{SFR}$ (analogous to Fig. 2), color-coded by the attenuation, $A_V$, estimated as described in the text. Galaxies with $A_V$ values derived from the emission lines are encapsulated within bigger blue squares.</figcaption></figure>

Carniani et al. (2024), the relatively low incidence rate is not necessarily evidence of a minor role of outflows, as this could be related to the outflow geometry. For example, an incidence rate of ~20% can be obtained if outflows have a biconical morphology with an opening angle of ~37 deg. This also suggests that to robustly constrain the role of outflows in the $\Sigma_{SFR} - f_{esc}$ correlation (hence, to confirm in a direct way the outflow-driven LyC photon escape scenario at the EoR), we need to study how the outflow fraction varies as a function of $\Sigma_{SFR}$ and $f_{esc}$, which requires better statistics than what is available in the present analysis. This scientific question will be addressed in the future with a larger sample of galaxies with deep, high-resolution NIRSpec data.

### 4.4.2. Investigating the attenuation and the outflow incidence of extremely star-forming galaxies

Approaching the epoch of cosmic dawn, galaxies tend to have denser ISM and more concentrated star formation activity, as shown in the previous sections. For electron densities exceeding ~200 cm$^{-3}$, the radiation pressure from young stars starts to dominate over the thermal pressure sustained by supernovae explosions (Ferrara 2024). Assuming the $(1+z)^p$ redshift evolution of $n_e$, we expect this transition to occur for the average galaxy population in a redshift range between 3 and 7, depending on the value of the exponent $p$ between 1 and 2 (see Isobe et al. 2023). In these conditions of gas density and pressure, supernova explosions rapidly lose energy due to radiative losses, and the energy is converted into sound waves through shocks on a short timescale (see, e.g., Martizzi et al. 2015, Gatto et al. 2015, Kim & Ostriker 2015, Walch & Naab 2015, Fielding et al. 2017), making supernova-driven outflows highly inefficient in early galaxies. Some theoretical models predict that in compact galaxies with extreme star formation activity, a super-Eddington regime could be established, according to which radiation driven outflows can dramatically decrease the dust optical depth by pushing it to larger radii in a few Myr (Ferrara 2024). This drastic event is expected when the sSFR exceeds ~$10^{-7.6}$ yr$^{-1}$. The consequence of this process is the formation of a system that has negligible dust attenuation, hence a
very blue SED. This model was introduced to explain the excess of UV-bright and blue galaxies (compared to theoretical expectations) discovered during the first JWST observing campaigns (e.g., Naidu et al. 2022; Castellano et al. 2022). However, it is still unclear whether it is the right physical explanation to the observational findings and when exactly this proposed mechanism might start to dominate.

To test this scenario, we analyzed the outflow incidence rate in galaxies with extremely high sSFR. First, we found that the fraction of galaxies in super Eddington regime (i.e., sSFR > $10^{-7.6}$ yr$^{-1}$) in all the redshift range explored in this work ($4 < z < 10$) is consistent with the model predictions by Ferrara (2024), increasing from ~15% at redshift 4 to ~40% at redshift 10. Among these galaxies, if we consider only those five observed by GLASS at high resolution, we find an outflow incidence rate of 1/5 (~20%); otherwise, we have 1/3 if we limit the analysis to systems with a compact morphology, that is, with an effective radius of 150 pc or less. In both cases, these fractions are compatible with the incidence rate obtained from the full sample.

An alternative method to testing the effect of star formation on the dust content of galaxies, including the dust clearing scenario, is to analyze how the galaxy attenuation, $A_V$, is related to the star-forming properties. To this aim, in Fig. 8, we color code the points in the $z - \Sigma_{SFR}$ diagram as a function of $A_V$. We can see that galaxies at $z > 7$ have a large variety of $A_V$: while some galaxies have upper limits on $A_V$ consistent with 0, some others have higher values that can reach $A_V = 2$ in magnitude. Therefore, we do not see evidence for a sudden drop in $A_V$ with redshift down to the very low values ($A_V \approx 10^{-2}-10^{-3}$) expected by the dust-clearing scenario, nor that the global galaxy population at $z > 7$ is dominated by systems with $A_V \approx 0$, even if the majority of this high-$z$ subset also has sSFRs largely exceeding the super-Eddington conditions defined above.

We can better understand this result by plotting $A_V$ as a function of sSFR for the entire redshift range (top panel of Fig. 9). We can see that galaxies with higher sSFR show a large variety of $A_V$. Globally, the median $A_V$ slightly increases with sSFR but it is consistent also with a flat relation ($m = 0.07 \pm 0.09$). In the bottom panel of Fig. 9 we analyze instead the relation between $A_V$ and $\Sigma_{SFR}$. Here we find a positive correlation between the two quantities, with a best-fit slope of $0.28 \pm 0.07$ with the Bayesian method, indicating that galaxies with more concentrated star formation (i.e., higher $\Sigma_{SFR}$) have higher dust attenuations on average. This positive correlation is similar to what is found in lower redshift studies (e.g., Reddy et al. 2015; Tomičić et al. 2017), and it reflects the fact that galaxies accumulate more dust in their ISM as a consequence of their increased (and more compact) star formation activity. Dividing our sample in two redshift bins (lower and higher than 6), we find no significant differences in the global trends seen in both panels of Fig. 9.

In conclusion, we can hypothesize several different explanations to our results. First, dust clearing radiatively driven outflows (in super Eddington conditions) have been proposed to explain the excess of bright and very blue galaxies at $z \ge 10$. We should note that the two NIRSpec galaxies at $z \approx 10$ in our sample are both dust poor, consistent with this scenario. We are focusing instead on a cosmic epoch where the feedback from supernovae is likely contributing to the acceleration of winds, as incorporated in several models and simulations (e.g., Sharma et al. 2017). Therefore, supernovae feedback remains the most plausible explanation for the outflows observed in NIRSpec galaxies at $z < 6$. However, we should also be aware of the possible lesser role played by super-

<figure><img src="image_10.png" /><figcaption>Fig. 9. Diagrams comparing $A_V$ to sSFR and $\Sigma_{SFR}$ (top and bottom panels, respectively) for our selected galaxies in the redshift range of $4 < z < 10$. The lines, shaded areas, and markers are the same as in Fig. 2.</figcaption></figure>

novae as we approach the reionization epoch, following the argument explained above. Galaxies at cosmic dawn are also expected to be younger than analogs at lower redshifts, with consequently shorter times available to supernovae to impact the surrounding ISM in a stable way.

Secondly, the fraction of halo gas entrained by the outflowing shell might be conspicuous at the level of 10% or more. Consequently, the gas outflow velocities can be relatively low (see equation 15 in Ferrara 2024), to the level of 50–100 km/s that are not easily detectable even in high resolution spectral mode. Another alternative explanation is that we are observing the galaxies in an intermediate phase of dust accumulation, possibly coincident with a period of gas accretion and dust creation by the first star formation episodes, which has not yet reached the conditions of dust opacity for developing a sudden blow-out event. In this particular phase, whose duration is not well known, the galaxies might still retain most of their gas and dust content, showing attenuations significantly higher than 0, while maintaining a low-level outflowing activity.
**5. Summary**

We summarize the main findings of this paper as follows. We analyze a spectroscopic sample of galaxies in the redshift range from 4 to 10 observed with JWST-NIRSpec by the CEERS and GLASS surveys. Using the Balmer lines and UV-based effective radii, we study the evolution of the SFR and Σ<sub>SFR</sub> over this cosmic epoch.

1. We find that Σ<sub>SFR</sub> increases mildly in this redshift range (best-fit slope = 0.16 ± 0.06), rising by a factor of 2 from z = 4 to 10. This trend is consistent with the predictions from several semi-analytic models of galaxy evolution, such as the Santa Cruz models (Somerville & Davé 2015) and those from Ferrara (2024) and Naidu et al. (2020).

2. A star-forming main sequence relation holds out to a redshift of ~10, with a best-fit slope of ~0.7, consistent with the relations obtained for star-forming galaxies at lower redshifts.

3. We derive, for the first time at z > 4, the Σ<sub>SFR</sub> main sequence, finding a very mild increase (consistent with a flat relation) of Σ<sub>SFR</sub> with stellar mass. The slope of the best-fit relation is in agreement with observations at 0 < z < 2, while the normalization is about ~2.5 dex higher (at z ~ 10) compared to the Local Universe, resulting from the strong evolution of the average Σ<sub>SFR</sub> in star-forming galaxies during the last 13 billion years.

4. We show that Σ<sub>SFR</sub> is tightly related to the ionizing source properties, as probed by the O32 index. This corroborates recent studies proposing Σ<sub>SFR</sub> as a key quantity responsible for the higher ionization parameter observed in high redshift galaxies (Reddy et al. 2023).

5. We find correlations among Σ<sub>SFR</sub> and sSFR with the escape fraction, estimated in an indirect way using the prescription of Mascia et al. (2023), calibrated on direct measurements of f<sub>esc</sub> from the low-redshift Lyman continuum survey (Flury et al. 2022). The slopes of both correlations are in agreement with observations of low-redshift LyC leakers and with those predicted by the models of Naidu et al. (2020). This suggests that Σ<sub>SFR</sub> and sSFR play an important role in the escape of ionizing radiation.

6. We detect outflow signatures in Hα and [O III] λ5007 Å lines for ~20% of the galaxies at z < 6 observed in H-grating NIRSpec mode, in agreement with the outflow incidence found in galaxies at similar redshifts by Carniani et al. (2024). This relatively low incidence might be related to the outflow geometry. In particular, a biconical outflow morphology with an opening angle of ~37 deg is able to explain the observed fraction. A larger sample is instead necessary to study the outflow incidence at varying Σ<sub>SFR</sub> and f<sub>esc</sub> so that we may have a more direct confirmation of the SFR-outflow connection. In addition, this would give us a better picture of the outflow-driven LyC photon escape scenario at the EoR.

7. We find a positive correlation between A<sub>V</sub> and Σ<sub>SFR</sub>, along with a flat trend as a function of sSFR, indicating that there is no evidence of a drop in A<sub>V</sub> in extremely star-forming galaxies between z = 4 and 10. This might be at odds with a dust-clearing outflow scenario being in place at z ≥ 10.

Observing a larger number of galaxies at high resolution and testing the correlation between Σ<sub>SFR</sub> and outflow velocity predicted by supernovae feedback models can help us further constrain the properties and physical origin of the outflows detected both in our and in similar works at z ~ 4 to 10. It will also enable us to further clarify the relation between star formation and outflows in the epoch of reionization. Enlarging the sample

of spectroscopically confirmed galaxies at z > 10 is essential for testing possible transformations of the ISM properties occurring at this cosmic epoch in the future.

**Acknowledgements.** We thank the referee for the detailed and constructive comments. AC acknowledges support from the INAF Large Grant for Extra-galactic Surveys with JWST. We acknowledge support from the PRIN 2022 MUR project 2022CB3PJ3 – First Light And Galaxy Assembly (FLAGS) funded by the European Union – Next Generation EU. We thank Stefano Carniani for useful discussion on the main topic of the paper.

**References**

Akaike, H. 1974, *IEEE Trans. Autom. Contr.*, **19**, 716

Amorín, R., Vílchez, J. M., Hågele, G. F., et al. 2012, ApJ, **754**, L22

Amorín, R. O., Rodríguez-Henríquez, M., Fernández, V., et al. 2024, A&A, **682**, L25

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, ApJ, **951**, L22

Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJ, **946**, L12

Beal, S. L. 2001, J Pharmacokinet Pharmacodyn, **28**, 481

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, **770**, 57

Behroozi, P., Wechsler, R. H., Hearin, A. P., et al. 2019, MNRAS, **488**, 3143

Bergamini, P., Acebron, A., Grillo, C., et al. 2023, ApJ, **952**, 84

Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, ApJ, submitted [arXiv:2212.04026]

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, **622**, A103

Bordoloi, R., Lilly, S. J., Hardmeier, E., et al. 2014, ApJ, **794**, 130

Borthakur, S., Heckman, T. M., Leitherer, C., et al. 2014, Science, **346**, 216

Bosch, G., Hågele, G. F., Amorín, R., et al. 2019, MNRAS, **489**, 1787

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, **686**, 1503

Brennan, R., Pandya, V., Somerville, R. S., et al. 2015, MNRAS, **451**, 2933

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

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

Burnham, K. P., & Anderson, D. R. 2004, *Sociol. Methods Res.*, **33**, 261

Calabrò, A., Daddi, E., Puglisi, A., et al. 2019, A&A, **623**, A64

Calabrò, A., Castellano, M., Pentericci, L., et al. 2021, A&A, **646**, A39

Calabrò, A., Pentericci, L., Talia, M., et al. 2022a, A&A, **667**, A117

Calabrò, A., Guaita, L., Pentericci, L., et al. 2022b, A&A, **664**, A75

Calabrò, A., Pentericci, L., Feltre, A., et al. 2023, A&A, **679**, A80

Calzetti, D. 2001, PASP, **113**, 1449

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, **533**, 682

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

Carniani, S., Venturi, G., Parlanti, E., et al. 2024, A&A, **685**, A99

Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJ, **938**, L15

Cen, R. 2020, ApJ, **889**, L22

Cole, J. W., Papovich, C., Finkelstein, S. L., et al. 2023, ArXiv e-prints [arXiv:2312.10152]

Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, **670**, 156

Davé, R., Oppenheimer, B. D., & Finlator, K. 2011, MNRAS, **415**, 11

Davidzon, I., Ilbert, O., Faisst, A. L., et al. 2018, ApJ, **852**, 107

Di Cesare, C., Graziani, L., Schneider, R., et al. 2023, MNRAS, **519**, 4632

Ding, X., Birrer, S., Treu, T., et al. 2021, ArXiv e-prints [arXiv:2111.08721]

Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, **485**, 4817

Duane, S., Kennedy, A. D., Pendleton, B. J., et al. 1987, Phys. Lett. B, **195**, 216

Euclid Collaboration (Scharré, L., et al.) 2024, A&A, **689**, A276

Faisst, A. L., Masters, D., Wang, Y., et al. 2018, ApJ, **855**, 132

Ferrara, A. 2024, A&A, **684**, A207

Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, **522**, 3986

Fielding, D., Quataert, E., Martizzi, D., et al. 2017, MNRAS, **470**, L39

Finkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2024, ApJ, **969**, L25

Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022, ApJS, **260**, 1

Fontana, A., D'Odorico, S., Poli, F., et al. 2000, AJ, **120**, 2206

Förster Schreiber, N. M., Übler, H., Davies, R. L., et al. 2019, ApJ, **875**, 21

Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023, ApJ, **949**, L25

Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, **450**, 4486

Furtak, L. J., Atek, H., Lehnert, M. D., et al. 2021, MNRAS, **501**, 1568

Gatto, A., Walch, S., Low, M.-M. M., et al. 2015, MNRAS, **449**, 1057

Gehrels, N. 1986, ApJ, **303**, 336

Groves, B., Brinchmann, J., & Walcher, C. J. 2012, MNRAS, **419**, 1402

Heckman, T. M. 2001, *Gas and Galaxy Evolution*, **240**, 345

Heckman, T. M. 2002, *Extragalactic Gas at Low Redshift*, **254**, 292

Heckman, T. M., & Borthakur, S. 2016, ApJ, **822**, 9

Heckman, T. M., Borthakur, S., Overzier, R., et al. 2011, ApJ, **730**, 5

Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, **451**, 2663

Hirschmann, M., Charlot, S., Feltre, A., et al. 2017, MNRAS, **472**, 2468

Hirschmann, M., Charlot, S., Feltre, A., et al. 2023, MNRAS, **526**, 3610
Hoffman, M.D., & Gelman, A. 2024, J. Mach. Learn. Res., 15, 1593

Hogarth, L., Amorín, R., Vilchez, J. M., et al. 2020, MNRAS, 494, 3541

Holwerda, B. W., Bouwens, R., Oesch, P., et al. 2015, ApJ, 808, 6

Hopkins, P. F., Murray, N., Quataert, E., et al. 2010, MNRAS, 401, L19

Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522

Huang, R., Battisti, A. J., Grasha, K., et al. 2023, MNRAS, 520, 446

Isobe, Y., Ouchi, M., Nakajima, K., et al. 2023, ApJ, 956, 139

Izotov, Y. I., Schäerer, D., Worseck, G., et al. 2020, MNRAS, 491, 468

Jiang, T., Malhotra, S., Yang, H., et al. 2019, ApJ, 872, 146

Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJ, 946, L15

Kawinwanichakij, L., Silverman, J. D., Ding, X., et al. 2021, ApJ, 921, 38

Kennicutt, R. C. 1989, ApJ, 344, 685

Kennicutt, R. C., Calzetti, D., Walter, F., et al. 2007, ApJ, 671, 333

Kewley, L. J., Zahid, H. J., Geller, M. J., et al. 2015, ApJ, 812, L20

Kim, C.-G., & Ostriker, E. C. 2015, ApJ, 802, 99

Kocevski, D. D., Barro, G., McGrath, E. J., et al. 2023, ApJ, 946, L14

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kornei, K. A., Shapley, A. E., Martin, C. L., et al. 2012, ApJ, 758, 135

Lehnert, M. D., van Driel, W., Le Tiran, L., et al. 2015, A&A, 577, A112

Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58

Llerena, M., Amorín, R., Pentericci, L., et al. 2023, A&A, 676, A53

Llerena, M., Amorín, R., Pentericci, L., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202449904

Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42

Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Markov, V., Gallerani, S., Pallottini, A., et al. 2023, A&A, 679, A12

Markwardt, C. B. 2009, ASP Conf. Ser., 411, 251

Martizzi, D., Faucher-Giguère, C.-A., & Quataert, E. 2015, MNRAS, 450, 504

Mascia, S., Pentericci, L., Calabrò, A., et al. 2023, A&A, 672, A155

Mascia, S., Pentericci, L., Calabrò, A., et al. 2024, A&A, 685, A3

Meiksin, A. 2006, MNRAS, 365, 807

Menci, N., Gatti, M., Fiore, F., et al. 2014, A&A, 569, A37

Merlin, E., Pilo, S., Fontana, A., et al. 2019, A&A, 622, A169

Morishita, T., Roberts-Borsani, G., Treu, T., et al. 2023, ApJ, 947, L24

Morishita, T., Stiavelli, M., Chary, R.-R., et al. 2024, ApJ, 963, 9

Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109

Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJ, 940, L14

Nakajima, K., Ouchi, M., Shimasaku, K., et al. 2013, ApJ, 769, 3

Napolitano, L., Pentericci, L., Santini, P., et al. 2024, A&A, 688, A106

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43

Ono, Y., Harikane, Y., Ouchi, M., et al. 2024, PASJ, 76, 219

Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, 1989, 408

Pannella, M., Elbaz, D., Daddi, E., et al. 2015, ApJ, 807, 141

Papovich, C., Simons, R. C., Estrada-Carpenter, V., et al. 2022, ApJ, 937, 22

Paris, D., Merlin, E., Fontana, A., et al. 2023, ApJ, 952, 20

Peng, C. Y., Ho, L. C., Impey, C. D., et al. 2002, AJ, 124, 266

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Prusinski, N. Z., Erb, D. K., & Martin, C. L. 2021, AJ, 161, 212

Reddy, N. A., Kriek, M., Shapley, A. E., et al. 2015, ApJ, 806, 259

Reddy, N. A., Topping, M. W., Shapley, A. E., et al. 2022, ApJ, 926, 31

Reddy, N. A., Topping, M. W., Sanders, R. L., et al. 2023, ApJ, 952, 167

Roberts-Borsani, G., Morishita, T., Treu, T., et al. 2022, ApJ, 938, L13

Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, Nat. Astron., 7, 611

Salim, S., Tacchella, S., Osborne, C., et al. 2023, ApJ, 958, 183

Santini, P., Fontana, A., Castellano, M., et al. 2017, ApJ, 847, 76

Schiminovich, D., Wyder, T. K., Martin, D. C., et al. 2007, ApJS, 173, 315

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Shapley, A. E. 2011, ARA&A, 49, 525

Shapley, A. E., Reddy, N. A., Kriek, M., et al. 2015, ApJ, 801, 88

Sharma, M., Theuns, T., Frenk, C., et al. 2016, MNRAS, 458, L94

Sharma, M., Theuns, T., Frenk, C., et al. 2017, MNRAS, 468, 2176

Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15

Shimakawa, R., Kodama, T., Steidel, C. C., et al. 2015, MNRAS, 451, 1284

Simmonds, C., Tacchella, S., Maseda, M., et al. 2023, MNRAS, 523, 5468

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Somerville, R. S., Olsen, C., Yung, L. Y. A., et al. 2021, MNRAS, 502, 4858

Speagle, J. S., Steinhardt, C. L., Capak, P. L., et al. 2014, ApJS, 214, 15

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289

Stefanon, M., Yan, H., Mobasher, B., et al. 2017, ApJS, 229, 32

Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJS, 795, 165

Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123

Sugahara, Y., Ouchi, M., Harikane, Y., et al. 2019, ApJ, 886, 29

Tacchella, S., Smith, A., Kannan, R., et al. 2022, MNRAS, 513, 2904

Tomićić, N., Kreckel, K., Groves, B., et al. 2017, ApJ, 844, 155

Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110

Treu, T., Calabrò, A., Castellano, M., et al. 2023, ApJ, 942, L28

Vanzella, E., de Barros, S., Castellano, M., et al. 2015, A&A, 576, A116

Walch, S., & Naab, T. 2015, MNRAS, 451, 2757

Ward, E., de la Vega, A., Mobasher, B., et al. 2024, ApJ, 962, 176

White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52

White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341

Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, ApJ, 742, 96

Yang, L., Morishita, T., Leethochawalit, N., et al. 2022, ApJ, 938, L17

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., et al. 2019a, MNRAS, 483, 2983

Yung, L. Y. A., Somerville, R. S., Popping, G., et al. 2019b, MNRAS, 490, 2855

Yung, L. Y. A., Somerville, R. S., Ferguson, H. C., et al. 2022, MNRAS, 515, 5416

¹ INAF – Osservatorio Astronomico di Roma, Via Frascati 33 00078 Monte Porzio Catone, Italy

² Scuola Normale Superiore, Piazza dei Cavalieri 7, 50126 Pisa, Italy

³ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

⁴ INAF – Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy

⁵ Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, Via Marzolo 8, I-35131 Padova, Italy

⁶ Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA

⁷ George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA

⁸ Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel

⁹ NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719, USA

¹⁰ The University of Texas at Austin, 2515 Speedway Blvd Stop, C1400 Austin, TX 78712, USA

¹¹ University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA

¹² Institute of Physics, Laboratory of Galaxy Evolution, Ecole Poly-technique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland

¹³ Department of Physics and Astronomy, University of Louisville, Natural Science Building 102, 40292, KY Louisville, USA

¹⁴ Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

¹⁵ Black Hole Initiative at Harvard University, 20 Garden St., Cambridge, MA 02138, USA

¹⁶ ESA/AURA, Space Telescope Science Institute, 3800 San Martin Drive, Baltimore, MD 21218, USA

¹⁷ Department of Physics and Astronomy, University of California, Los Angeles, 430 Portola Plaza, Los Angeles 90095, CA, USA

¹⁸ Aix Marseille Université, CNRS, CNES, LAM, 38, rue Frédéric Joliot-Curie, 13388 Marseille, CEDEX 13 France

¹⁹ Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

²⁰ Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HE, UK

²¹ Astronomy Centre, University of Sussex, Falmer, Brighton BN1 9QH, UK

²² Institute of Space Sciences and Astronomy, University of Malta, Msida, MSD 2080, Malta

²³ Departamento de Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena, Chile

²⁴ ARAID Foundation. Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Unidad Asociada al CSIC, Plaza San Juan 1, E-44001 Toruñ, Spain

²⁵ Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA