© The Authors 2025

1. Introduction

One of the most pressing questions in extragalactic astronomy currently is explaining the surprisingly high abundances of bright and massive galaxies in the early Universe revealed by the James Webb Space Telescope (JWST). Statistical functions such as the ultraviolet luminosity function (UVLF) and the stellar mass function (SMF) are key observational measurements that inform models of galaxy formation and evolution in the early Universe. Measurements of both the UVLF and SMF from JWST data (e.g., Donnan et al. 2023; Harikane et al. 2023; Pérez-González et al. 2023; Bouwens et al. 2023; Finkelstein et al. 2024; Donnan et al. 2024; McLeod et al. 2024; Adams et al. 2024; Robertson et al. 2024 for the UVLF; Weibel et al. 2024; Shuntov et al. 2025; Harvey et al. 2025; Wang et al. 2024 for the SMF) have revealed galaxy abundances in excess of theoretical models and simulations calibrated on the “pre-JWST Universe” (e.g., Boylan-Kolchin 2023).

To reconcile the theory with observations, several physical mechanisms have been proposed. These can be broadly separated in two classes: (1) mechanisms that increase the star formation efficiency (SFE) in the early Universe, such as feedback-free starbursts (Torrey et al. 2017; Grudić et al. 2018; Dekel et al. 2023; Li et al. 2024; Renzini 2023), and (2) mechanisms that do not need increased SFE and can instead produce abundant bright galaxies due to increased stochasticity in the UV luminosities (Mason et al. 2023; Pallottini & Ferrara 2023; Shen et al. 2023; Gelli et al. 2024) via limited-to-no dust attenuation (Ferrara et al. 2023) or a top-heavy initial mass function (IMF, Cueto et al. 2024; Hutter et al. 2025), to name a few. The increased stochasticity manifests in the dispersion in the relation between galaxy UV magnitude (M_UV) and halo mass (M_h), quantified by σ_UV. Bursty star formation histories naturally produce this stochasticity and have been identified both observationally via the scatter in the star forming main sequence (e.g., Ciesla et al. 2024; Clarke et al. 2024; Cole et al. 2025) and Hα-to-UV luminosity ratios (e.g., Endsley et al. 2024) as well as in simulations (Pallottini & Ferrara 2023; Sun et al. 2023; Muratov et al. 2015; Sparre et al. 2017), although halo assembly histories and varying dust attenuation can also contribute to σ_UV (Shen et al. 2023).

Most of these models have been developed and calibrated to reproduce one-point statistical observables such as the UVLF and the UV luminosity density. However, the galaxy-halo connection is complex (Wechsler & Tinker 2018), and there are several different mechanisms that can produce the same one-point observable, therefore making it intrinsically degenerate (Muñoz et al. 2023). Crucial information lies in how galaxies are spatially distributed, not only in how many there are. The galaxy spatial distribution, or clustering, is described by two-point statistical functions, such as the two-point correlation function (2PCF). A successful galaxy formation model has to explain both the abundances and clustering of galaxies, and therefore combining such complementary observables is key in breaking degeneracies between models (Mirocha 2020; Muñoz et al. 2023; Gelli et al. 2024).

On the observational side, using clustering as a probe alongside galaxy abundances has been done out to z ∼ 7 using Lyα emitters and Lyman-break galaxies (LBG) from Subary Suprime-Cam and Hyper Suprime-Cam (HSC) observations (Ouchi et al. 2010, 2018; Harikane et al. 2018, 2022). JWST has allowed measurements of the 2PCF out to the highest redshifts, z ∼ 10, from photometrically selected samples (Dalmasso et al. 2024a; Paquereau et al. 2025). However, JWST’s NIRCam/grism capabilities have for the first time enabled complete spectroscopically selected samples out to z ∼ 9 (e.g., Oesch et al. 2023; Kashino et al. 2023) to be compiled and their clustering (as demonstrated by Eilers et al. 2024, for a limited sample of grism-selected [O III]-emitting galaxies at 5 ≲ z ≲ 7 in a quasar field) to be measured.

On the theory side, the degeneracy issue has been explored by Mirocha (2020), Muñoz et al. (2023), who have shown that the two scenarios of high SFE versus high stochasticity that result in the same UVLF have a very distinct impact on the galaxy bias. The latter is typically measured by comparing the observed 2PCF of galaxies to the theoretical one of dark matter halos, and it describes how galaxies are preferentially distributed in the highest peaks of the underlying matter distribution. However, a comprehensive framework that models both the UVLF and 2PCF and that can be used to derive empirical constraints of the theoretical models is still lacking.

In this paper, we describe the development a theoretical framework whose purpose is to provide observational constraints on the early Universe galaxy-halo connection through the combined measurements of the UVLF and galaxy clustering via the 2PCF. This framework is based on an adaptation of the halo occupation distribution (HOD) in combination with the conditional luminosity function (CLF) and a parametric form of the SFE-M_h relation. Fitting the observed UVLF and 2PCF in an independent binning scheme allows us to provide observational constraints on the M_UV–M_h and SFE-M_h relations, stochasticity, host halo masses, and galaxy bias. To implement this framework, we measured the 2PCF of spectroscopically selected samples of Hα and [O III] emitters at 3.8 < z < 9 from the JWST NIRCam/grism surveys FRESCO (Oesch et al. 2023) and CONGRESS (Egami et al. 2023).

This paper is organized as follows. Section 2 presents in detail the theoretical framework. Sect. 3 briefly presents the NIRCam/grism data and methods used to compile the line emitting samples. In Sect. 4 we describe the methods we used to measure the 2PCF and UVLF. Sect. 5 presents the results of this work. In Sect. 6 we discuss the physical implications of our results, and in Sect. 7 we summarize and conclude the paper.

We adopted a standard ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m, 0 = 0.3, where Ω_b, 0 = 0.04, Ω_Λ, 0 = 0.7, and σ₈ = 0.82. All magnitudes are expressed in the AB system (Oke & Gunn 1983). Wherever relevant, we assumed a Salpeter (1955) IMF and rescaled the literature data points to the same IMF. The dependence of various quantities on the reduced Hubble parameter h is retained implicitly. Dark matter halo masses, scale as h⁻¹. Absolute magnitude scale as −5 log(h) in addition. When comparing to the literature, we rescale all the measurements to the cosmology adopted for this paper.

2. Theoretical framework

We develop a theoretical framework that is based on physically motivated principles and the halo model and from which we can derive fundamental statistical observables of galaxies such as the UVLF and the galaxy correlation function. At the heart of this formalism is a parametric model of the instantaneous SFE that relates the halo accretion rate to star formation, coupled with an HOD framework to predict galaxy clustering and abundances as a function of M_UV. This framework is analogous to the one developed to predict clustering and the galaxy SMF by parametrizing the stellar-to-halo mass relation (e.g., Moster et al. 2010; Leauthaud et al. 2011, 2012; Shuntov et al. 2022). The novelty of this work is the extension to the UVLF, which also sets it apart from other works that use the 2PCF in combination to the total number density of M_UV-threshold selected galaxies (e.g., Harikane et al. 2018, 2022).

2.1. The ultraviolet luminosity to halo mass relationship

Galaxy star formation and mass growth is directly connected to the growth rate of dark matter halos since the gas available for star formation is well mixed with the dark matter with a ratio of f_b ≈ 0.16 (e.g., Moster et al. 2018; Tacchella et al. 2018). This relation can be written as

$$\begin{array}{c}\hfill {\dot{M}}_{\star }=ϵf_{\mathrm{b}}{\dot{M}}_{\mathrm{h}}.\end{array}$$(1)

The halo accretion rate, Ṁ_h, is well understood from ΛCDM, and we took the Dekel et al. (2013) form as a function of halo mass and redshift. The star formation rate Ṁ_⋆ is moderated by the instantaneous SFE, ϵ. This function encapsulates all baryonic physics responsible for regulating star formation in halos and is essentially a summary of our galaxy formation model. Typically, it depends on both redshift and halo mass ϵ(z, M_h), and we adopted the double power law parametrization as a function of M_h (Moster et al. 2010, 2018):

$$\begin{array}{c}\hfill ϵ(M_{\mathrm{h}})=2{ϵ}_0{[{\left(\frac{M_{\mathrm{h}}}{M_c}\right)}^{-\beta }+{\left(\frac{M_{\mathrm{h}}}{M_c}\right)}^{\gamma }]}^{-1}.\end{array}$$(2)

This relation has four free parameters: ϵ₀, the normalization or the peak SFE; M_c, a characteristic halo mass at which the efficiency peaks; and β and γ, respectively the low- and high-mass slopes. In our implementation, we model these parameters as a function of redshift and fit for their redshift dependence (following e.g., Moster et al. 2018; Muñoz et al. 2023, more details in Sect. 4.3).

We note here that Eq. (1) is distinct from the integrated SFE, ϵ_⋆, that has been used frequently in recent high-z literature (e.g., Behroozi & Silk 2018; Boylan-Kolchin 2023). The latter quantifies the ratio between the stellar and halo mass assembled throughout the history of a halo, ${ϵ}_{\star }=M_{\star }{M}_{\mathrm{h}}^{-1}{f}_{\mathrm{b}}^{-1}$, also known as the stellar-to-halo mass relationship (SHMR, Wechsler & Tinker 2018). The two are related by a time integral and do not necessarily have the same normalization and shape (Behroozi et al. 2013; Moster et al. 2018), unless ϵ is Universal and time independent. In this paper, we mainly work with the instantaneous SFE (ϵ), and refer to it whenever we discuss the SFE, unless we state otherwise (e.g., if we refer to the integrated SFE)

To arrive at the UV luminosity, we used the fact that the SFR and UV luminosity are related with the luminosity-to-SFR conversion factor κ_UV:

$$\begin{array}{c}\hfill L_{\mathrm{UV}}=\frac{{\dot{M}}_{\star }}{{\kappa }_{\mathrm{UV}}}.\end{array}$$(3)

Here, κ_UV is fully degenerate with ϵ₀ (Muñoz et al. 2023), but in a first and simplistic approach we assumed its fiducial value ${\kappa }_{\mathrm{UV}}=1.15\times {10}^{-28}(M_{\odot }{\mathrm{yr}}^{-1})/(\mathrm{erg}{\mathrm{s}}^{-1})$ for a Salpeter IMF (Madau & Dickinson 2014). From the luminosity, we then derived the absolute magnitude M_UV in the AB system in the following way:

$$\begin{array}{c}\hfill M_{\mathrm{UV}}=51.63-2.5\log \left(L_{\mathrm{UV}}/\mathrm{erg}{\mathrm{s}}^{-1}\right).\end{array}$$(4)

Finally, we added dust by using the A_UV − β relation from Meurer et al. (1999) and β − M_UV from Bouwens et al. (2014) and following the implementation in Mason et al. (2015), Sabti et al. (2022) and references therein. This assumes that β follows a Gaussian distribution at each M_UV with a dispersion of σ_β = 0.34. The average dust extinction is then ⟨A_UV⟩ = 4.43 + 0.79 ln(10) σ_β + 1.99 ⟨β⟩.

Equations (1)–(4) relate the UV luminosity, M_UV, to the halo mass, M_h, resulting in a M_UV–M_h relationship, or UVHMR, which we denote as f_UVHMR(M_h). M_UV has a monotonically increasing relation with M_h that depends on redshift (e.g., Mason et al. 2015). The exact dependence is shaped by the change of both halo growth and the instantaneous SFE with time.

2.2. The conditional ultraviolet luminosity function

The M_UV–M_h is not a one-to-one relation and there is a scatter in M_UV at a given M_h. This can be modeled by a CLF, which gives the average number of galaxies with M_UV ± dM_UV/2 in a halo of mass M_h (Yang et al. 2003, 2009). This can be modeled as a Gaussian distribution with a scatter σ_UV:

$$\begin{array}{c}\hfill P(M_{\mathrm{UV}}|M_{\mathrm{h}})=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{UV}}}\exp [-{\left(\frac{M_{\mathrm{UV}}-f_{\mathrm{UVHMR}}(M_{\mathrm{h}})}{\sqrt{2}{\sigma }_{\mathrm{UV}}}\right)}^2],\end{array}$$(5)

where f_UVHMR(M_h) gives the average M_UV as a function of halo mass M_h. The scatter σ_UV is not only a result of intrinsic astrophysical processes, σ_UV, intr but also of observational uncertainties in measuring M_UV, σ_UV, meas (Leauthaud et al. 2011), which can be written as σ_UV² = σ_UV, intr² + σ_UV, meas². Additionally, σ_UV can be a function of halo mass (e.g., Sun et al. 2023; Gelli et al. 2024), however in our implementation in this paper we simplify it and neglect this dependence on M_h.

2.3. The halo occupation distribution

The CLF, crucially, provides a link to the HOD framework, which describes the statistical occupation of galaxies in dark matter halos and provides a formalism to model galaxy clustering via the two-point correlation function (e.g., Cooray & Sheth 2002; Berlind & Weinberg 2002). The HOD is a prescription of the mean number of galaxies residing in a halo of mass M_h, which is given as

$$\begin{array}{c}\hfill \langle N_{\mathrm{tot}}\left(M_{\mathrm{h}}\right)\rangle =\langle N_{\mathrm{cent}}\left(M_{\mathrm{h}}\right)\rangle \times [1+\langle N_{\mathrm{sat}}\left(M_{\mathrm{h}}\right)\rangle ].\end{array}$$(6)

It is made of two components describing the mean occupation numbers of central and satellite galaxies. The occupation function for central galaxies brighter than some threshold luminosity, $M_{\mathrm{UV}}^{\mathrm{th}}$, can be derived by integrating the CLF

$$\begin{array}{c}\hfill {\displaystyle \langle N_{\mathrm{cent}}\left(M_{\mathrm{h}}|M_{\mathrm{UV}}^{\mathrm{th}}\right)\rangle ={\int }_{\infty }^{M_{\mathrm{UV}}^{\mathrm{th}}}\mathrm{d}{M}_{\mathrm{UV}}P(M_{\mathrm{UV}}|M_{\mathrm{h}}).}\end{array}$$(7)

Carrying out this integration and accounting for the fact that M_UV is negative, we arrived at

$$\begin{array}{c}\hfill \langle N_{\mathrm{cent}}\left(M_{\mathrm{h}}|M_{\mathrm{UV}}^{\mathrm{th}}\right)\rangle =\frac{1}{2}[1-\text{erf}\left(\frac{f_{\mathrm{UVHMR}}(M_{\mathrm{h}})-M_{\mathrm{UV}}^{\mathrm{th}}}{\sqrt{2}{\sigma }_{\mathrm{UV}}}\right)].\end{array}$$(8)

We can obtain this functional form by analytically integrating Eq. (7), thanks to the Gaussian form of P(M_UV|M_h) and because we assumed a constant σ_UV. In principle, it is possible to allow a functional form of σ_UV with M_h and to carry out the integration of Eq. (7) numerically to obtain ⟨N_cent⟩.

The occupation of halos by satellites can be modeled by a power-law for which we adopted the Zheng et al. (2005) form

$$\begin{array}{c}\hfill \langle N_{\mathrm{sat}}\left(M_h|M_{\mathrm{UV}}^{\mathrm{th}}\right)\rangle ={\left(\frac{M_{\mathrm{h}}-M_{\mathrm{cut}}}{M_{\mathrm{sat}}}\right)}^{{\alpha }_{\mathrm{sat}}},\end{array}$$(9)

where M_cut and M_sat are the minimum and characteristic halo masses necessary to host satellite galaxies, and α_sat is the power-law slope. This is different from other implementations of this formalism, where the satellite occupation is modeled by a combination of power- and exponential-law (Leauthaud et al. 2011; Coupon et al. 2015; Shuntov et al. 2022). We opted for this simpler description in order to reduce the number of free parameters.

The average number of galaxies in a UV luminosity bin of $M_{\mathrm{UV}}^{\mathrm{th}1}<M_{\mathrm{UV}}<M_{\mathrm{UV}}^{\mathrm{th}2}$ can also be computed by simply taking the difference of the occupation numbers for two M_UV thresholds

$$\begin{array}{c}\hfill \begin{array}{cc}& \langle N_{\mathrm{cent}/\mathrm{sat}}(M_h|M_{\mathrm{UV}}^{\mathrm{th}1},M_{\mathrm{UV}}^{\mathrm{th}2})\rangle =\hfill \\ \hfill & \langle N_{\mathrm{cent}/\mathrm{sat}}\left(M_h|>M_{\mathrm{UV}}^{\mathrm{th}1}\right)\rangle -\langle N_{\mathrm{cent}/\mathrm{sat}}\left(M_h|>M_{\mathrm{UV}}^{\mathrm{th}2}\right)\rangle .\hfill \end{array}\end{array}$$(10)

2.4. Model of the galaxy ultraviolet luminosity function

Having a description of the mean occupation numbers of galaxies as a function of their M_UV and of their host halo mass M_h, we can derive an analytical form for the UVLF by simply integrating over the halo mass function in the following way

$$\begin{array}{c}\hfill {\displaystyle \mathrm{\Phi }(M_{\mathrm{UV}}^{\mathrm{th}1},M_{\mathrm{UV}}^{\mathrm{th}2})={\int }_0^{\infty }\mathrm{d}{M}_{\mathrm{h}}\langle N_{\mathrm{tot}}(M_{\mathrm{h}}|M_{\mathrm{UV}}^{\mathrm{th}1},M_{\mathrm{UV}}^{\mathrm{th}2})\rangle \frac{\mathrm{d}n}{\mathrm{d}{M}_{\mathrm{h}}}.}\end{array}$$(11)

Here, $\mathrm{\Phi }(M_{\mathrm{UV}}^{\mathrm{th}1},M_{\mathrm{UV}}^{\mathrm{th}2})$ indicates that the UVLF is computed for a M_UV bin defined by the two thresholds $M_{\mathrm{UV}}^{\mathrm{th}1}$ and $M_{\mathrm{UV}}^{\mathrm{th}2}$. For the integration limits, in practice we adopt the [10⁶ − 10¹⁵] M_⊙ range. dn/dM_h is the halo mass function (HMF), for which we use the code COLOSSUS (Diemer 2018), with the Watson et al. (2013) HMF and virial overdensity halo mass definition.

2.5. Model of the two-point correlation function

The model of the two-point angular correlation function follows closely the usual prescriptions (see, e.g., Cooray & Sheth 2002). For completeness, we detail the principal equations in Appendix D. For the computation of the 2PCF, we rely on the HALOMOD¹ code (Murray et al. 2013, 2021). HALOMOD is a library of routines to calculate a range of cosmological observables. The main ingredients that enter the modeling along with the prescriptions and assumptions that we adopt are given in Table 1.

One of the ingredients in computing the 2PCF is the halo bias, which describes how halos trace the underlying dark matter distribution. Analytical forms of the halo bias are typically calibrated in the linear regime and on large scales using N-body simulations. However, not accounting for non-linear and scale-dependent effects fails to accurately describe the ‘quasi-linear’ regime of the one-halo to two-halo transition (r ∼ 10 − 500 kpc) and can underpredict galaxy clustering by up to 30% (Reed et al. 2009; Jose et al. 2013; Mead et al. 2015; Jose et al. 2016; Mead & Verde 2021). As demonstrated in more detail in Harikane et al. (2018), Paquereau et al. (2025), accounting for non-linear and scale-dependent bias becomes important in accurately modeling the observed 2PCF at z ≳ 3. To deal with this, we use the same implementation as in Paquereau et al. (2025), based on the Jose et al. (2016) non-linear scale-dependent corrections. Finally, we also apply the halo-exclusion correction that accounts for the fact that halos cannot overlap and supresses the two-halo power at scales smaller than the halo size (Zheng 2004; Tinker et al. 2005).

2.6. Derived parameters

Using the mean halo occupation numbers and the halo mass function, we can derive various parameters related to the DM halos. The number density of galaxies predicted by the HOD model is computed as

$$\begin{array}{c}\hfill {\displaystyle n_{\mathrm{gal}}(z)=\int \mathrm{d}{M}_{\mathrm{h}}\frac{\mathrm{d}n}{\mathrm{d}{M}_{\mathrm{h}}}⟨N_{\mathrm{tot}}(M_{\mathrm{h}})⟩.}\end{array}$$(12)

The average mass of the DM halo of galaxies at a given redshift is

$$\begin{array}{c}\hfill {\displaystyle ⟨M_{\mathrm{h}}⟩(z)=\frac{1}{n_{\mathrm{gal}}(z)}\int \mathrm{d}{M}_{\mathrm{h}}{M}_{\mathrm{h}}\frac{\mathrm{d}n}{\mathrm{d}{M}_{\mathrm{h}}}⟨N_{\mathrm{tot}}(M_{\mathrm{h}})⟩.}\end{array}$$(13)

The mean galaxy bias at a given redshift is

$$\begin{array}{c}\hfill {\displaystyle b_{\mathrm{g}}(z)=\frac{1}{n_{\mathrm{gal}}(z)}\int \mathrm{d}{M}_{\mathrm{h}}{b}_{\mathrm{h}}(M_{\mathrm{h}},z)\frac{\mathrm{d}n}{\mathrm{d}{M}_{\mathrm{h}}}⟨N_{\mathrm{tot}}(M_{\mathrm{h}})⟩,}\end{array}$$(14)

where b_h(M_h, z) is the Tinker et al. (2010) linear halo bias.

2.7. Ultraviolet luminosity function, two-point correlation function, and galaxy bias dependence on model parameters

In Fig. A.1 we show how the UVLF, 2PCF and galaxy bias change with different model parameter values. This demonstrates which measurements are most constraining for the model parameters. For this we adopt a set of fiducial parameter values σ_UV = 0.76,  ϵ₀ = 0.22,  log M_c = 11.56,  β = 0.69,  γ = 0.74,  log M_cut = 8.99,  log M_sat = 12.69, α_sat = 0.76, and we vary each parameter at a time, shown in the different rows in Fig. A.1.

The UVLF is most sensitive to the parameters that govern instantaneous SFE ϵ(M_h) and σ_UV. The scatter, σ_UV, mostly affects the bright end, such that a higher scatter flattens the bright-end slope, resulting in higher abundances of bright (M_UV < −21) galaxies. Correspondingly, this reduces the clustering amplitude and galaxy bias. Increasing the maximum SFE, ϵ₀ effectively shifts the UVLF toward the brighter end while decreasing the 2PCF amplitude and bias. This happens because higher efficiencies mean that lower mass halos, that are less clustered, can host brighter galaxies. The peak halo mass M_c governs the location of the ‘knee’ of the UVLF with higher values moving the knee toward the brighter end, while increasing the 2PCF amplitude. The slopes β and γ affect the faint and bright end respectively, where increasing β, increases the 2PCF, while γ has very little impact on the 2PCF. The parameters governing the satellite occupation (M_sat,  M_cut,  α_sat) have negligible impact on the UVLF and bias, but affect mostly the smaller scales of the 2PCF (one-halo term).

2.8. Advantages and caveats

The theoretical framework we developed for this work is highly flexible in modeling galaxy observables under a common parametrization of the instantaneous SFE, the M_UV − M_h scatter, and the satellite occupation numbers. These model parameters can in turn be parametrized with redshift. This is a powerful approach that allows modeling of the UVLF at any redshift bin under the same model parameters and consistently incorporates different z-bins under the same likelihood. Furthermore, this method allows for an independent binning scheme for the 2PCF and UVLF measurements. In fact, this is how we implemented this framework in this paper (see Sect. 4.2).

Furthermore, this framework is highly modular – the same parametrization can be used to model additional statistical observables such as the SMF. This observable is complementary to the UVLF and subject to different systematics and therefore can be a powerful addition to help break degeneracies and better constrain the model. We will explore this in a future work.

As detailed in Sect. 2.1, the modeling includes a number of assumptions that can introduce systematic biases. For example, assumptions of the UV luminosity-to-SFR conversion factor κ_UV and the dust attenuation relations can be degenerate with the model parameters (we discuss this further in Sect. 6). Furthermore, there is also a degeneracy with the duty cycle that gives the fraction of galaxies that are UV-bright at a given moment (e.g., Mirocha 2020; Muñoz et al. 2023). The duty cycle can be incorporated in the model by a multiplying the total occupation number N_tot f_duty (e.g., Harikane et al. 2016) and would result in scaling the 2PCF and UVLF differently (Mirocha 2020; Muñoz et al. 2023, provide details on how exactly it affects the UVLF and galaxy bias). For example, lowering f_duty would require increasing ϵ₀ to match the UVLF, but would decrease the 2PCF amplitude. We assume a duty cycle of unity, consistent with indications from previous work (Harikane et al. 2018, 2022; Muñoz et al. 2023). Future work including larger datasets would allow us to model the duty cycle, which can even be a function of halo mass (e.g., Ouchi et al. 2001; Mirocha 2020).

3. Data

This work is anchored on the most comprehensive and complete JWST NIRCam/grism spectroscopic surveys – FRESCO (Oesch et al. 2023) in the F444W filter and CONGRESS (Sun, Egami et al. in prep., Egami et al. 2023) in the F356W. These two surveys cover the 3.1 − 5.0 μm range over 62 arcmin² in GOODS-North. In addition, FRESCO covers another 62 arcmin² in GOODS-South at 3.8 − 5.0 μm. Combined, these enable the compilation of complete samples of Hα emitters at 3.7 < z < 6.7 (4.9 < z < 6.7) in GOODS-North (GOODS-South), and [O III] emitters at 6.8 < z < 9.

We used the samples compiled by the works of Covelo-Paz et al. (2025, for Hα emitters) and Meyer et al. (2024, for [O III] emitters). We describe them briefly here but refer the reader to these references for more details. Photometric catalogs in the GOODS fields are extracted using all publicly available HST (e.g., Giavalisco et al. 2004; Grogin et al. 2011; Koekemoer et al. 2011, in ten bands) and JWST (e.g., Eisenstein et al. 2023a,b; Williams et al. 2023, in 11 bands (19 for GOODS-S). Detection is carried out on an inverse-variance stack of the F210M and F444W bands from FRESCO to ensure homogenous depth. Source extraction is done with SEXTRACTOR (Bertin & Arnouts 1996) on images that were PSF-matched to the F444W band. Aperture corrections are applied to scale the fluxes to those measured in Kron apertures on a PSF-matched version of the detection image, and then to total by dividing by the encircled energy of the Kron aperture on the F444W PSF. The detection completeness is estimated using the GLAC IAR2 software (Leethochawalit et al. 2022) by measuring the fraction of injected-recovered sources as a function of F444W magnitude. The completeness is constant around 90% out to 27 AB mag and declines to zero out to 30.5 AB mag. The NIRCam WFSS data is reduced and processed using the GRIZLI software using the standard CRDS grism dispersion files from pmap 1123 (Oesch et al. 2023; Meyer et al. 2024). In the following, we briefly describe the selection of Hα and [O III] emitters.

3.1. Hα sample

The Hα sample is constructed by first selecting all sources with Hα S/N > 3 from the GRIZLI catalog. Multi-component sources are identified upon visual inspection and components within < 0.6 arcsec separation and identical grism redshifts are merged into one galaxy and assigned the same Hα line. Quality flags are assigned via visual inspection by four team members independently. We use quality flag q ≥ 1.75, which means that the source has a clear Hα line and a matching morphology between the direct image and the 2D-spectrum and the line map, or a lower S/N but a clear morphology map. Additionally, a high quality flag is assigned for FRESCO Hα emitter at z > 5.25 that show a [O III] emission from the CONGRESS data (Covelo-Paz et al. 2025).

3.2. [O III] sample

The [O III] sample is constructed by requiring S/N ≥ 3 in F444W FRESCO band and no detection in bands blueward of Lyman-α to remove unambiguous low-z contaminants. Then, for all these sources, the 1D spectra are filtered with a Gaussian with FWHM = 50, 100, 200 km s⁻¹ and candidates are retained if they satisfy the following conditions: (1) Two lines match the [O III]λλ5008,4960 Å separation at 6.8 < z < 9.0, with a tolerance of Δν < 100 km s⁻¹ on the doublet separation, (2) S/N> 4 for the strongest line of the doublet, (3) The observed doublet ratio 1 <[O III]5008/[O III]4960 < 10. Visual inspection is carried out by multiple team members independently to remove contaminants and assign quality flags (q), similarly to the Hα catalog. We use objects with q ≥ 1.5, which is slightly less conservative than for the Hα sample; we opt for this lower value to maximize the number of objects for this intrinsically smaller sample of [O III] emitters (Meyer et al. 2024).

3.3. Selection function

The selection described above yields emission line-selected samples in three redshift bins – two Hα emitter samples, the first at 3.7 < z < 5.1 from CONGRESS, the second at 4.9 < z < 6.7 from FRESCO, and a third sample of [O III] emitters at 6.8 < z < 9. We apply a cut in the UV absolute magnitude to select only sources with M_UV < −19.1. This threshold is chosen as the brightest magnitude that selects enough sources to allow measurement of the 2PCF w(θ). Practically, this requires that there is at least one pair of objects at the angular separations θ at which w(θ) is evaluated using the Landy & Szalay (1993) estimator (Sect. 4.1, Eq. (15)).

This selection function results in 398 galaxies for the 3.7 < z < 5.1 Hα CONGRESS sample, 387 galaxies for the 4.9 < z < 6.7 for the Hα FRESCO sample, and 108 galaxies for the 6.8 < z < 9 [O III] sample. Figure 1 shows the spatial and redshift distributions of the three samples in the GOODS-North and GOODS-South fields.

Fig. 1. Spatial distribution of the three line-emitting samples (Hα and [O III] from FRESCO and Hα from CONGRESS) in the GOODS-North and GOODS-South footprints (left and middle panels) that we consider in this analysis. The right panel shows their redshift distribution. The density histogram gray background in the left and middle panels shows the footprint of the random catalogs.

4. Measurements

4.1. Galaxy clustering via two-point angular correlation function

We measure galaxy clustering via the angular 2PCF, for the three emission line-selected samples, brighter than M_UV = −19.1 mag, in the three redshift bins, using the Landy & Szalay (1993) estimator². We measured the 2PCF using Eq. (15) for GOODS-North and GOODS-South separately and combined them together using a number density weighting scheme following Durkalec et al. (2015)

$$\begin{array}{c}\hfill w(\theta )=\frac{{\displaystyle \sum _{i=1}^{n_{\mathrm{field}}}{w}_i(D{D}_i-2D{R}_i+R{R}_i)}}{{\displaystyle \sum _{i=1}^{n_{\mathrm{field}}}{w}_{i}R{R}_i}},\end{array}$$(15)

where DD_i, RR_i, and DR_i are the number of data-data, random-random, and data-random pairs in a given angular separation bin [θ, θ + δθ] normalized by the total number of galaxies and random objects for each field i = 1, 2. In addition, w_i = (N_i/V_i)² is the weight defined as the total number of galaxies divided by the volume of the corresponding field.

We constructed random catalogs for the two fields and the two samples by Monte-Carlo sampling of the LF (Hα and [O III] accordingly) to assign a line flux to each random source. To account for the fact that the spatial selection function can be wavelength dependent we use the root-mean-square (RMS) cube given in Meyer et al. (2024) along with the completeness functions for the [O III] and Hα samples (Meyer et al. 2024; Covelo-Paz et al. 2025). We compute the S/N for each source as a function of position and wavelength and retain the random draw with a probability equal to the completeness as a function of S/N. This ensures that the random catalog accounts for any spatial selection effects in the data.

Given the fact that this is a sample with highly precise redshifts from spectroscopy, it is in principle possible to compute the projected correlation as a function of comoving perpendicular separation w_p(r_p). However, our samples are too small to carry out the measurement of w_p(r_p) that is not too noisy for model fitting, and therefore we use the projected angular correlation function w(θ) for our analysis. This is also consistent with the predictions by Endsley et al. (2020) using mock galaxy samples based on the UNIVERSEMACHINE model.

Uncertainties in the 2PCF measurements come from Poisson pair-counting statistics and from cosmic variance (Norberg et al. 2009). The latter can have important impact on the measurements, especially for small field-of-view surveys like FRESCO. One way to estimate the cosmic-variance uncertainties is via jackknife, bootstrapping or Monte Carlo methods (Norberg et al. 2009, 2011). However, the small number of sources and footprint render jackknife and bootstrapping methods unsuitable, while Monte Carlo methods would require cosmological simulations and synthetic galaxy catalogs, which is expensive and out of the scope of this work. Furthermore, using such simulations, Leauthaud et al. (2011) showed that cosmic variance has an impact on the covariance matrix on large scales (i.e., in the two-halo regime) – it increases the correlation in the data at large scales, which is poorly constrained with our measurements due to the relatively small field.

We adopt a simplistic approach by considering only Poisson uncertainties, boosted by 30% in quadrature. We estimated the 30% cosmic variance uncertainty using the galcv calculator (Trapp et al. 2022) for the survey volume apparent magnitude of our sample. Additionally, we tested a more conservative boost of 10% and verified that it does not impact the results significantly. Because of the small sample, it is also difficult to derive well-defined covariance matrices that are not too noisy for likelihood computation. For this reason, we use only the diagonal elements, and therefore not consider the covariance between different θ bins (commonly adopted in the literature in case of small samples, e.g., Zehavi et al. 2011).

Due to the small volume probed in the GOODS-North and GOODS-South fields, the integral constraint (IC) suppresses/cuts w(θ) at scales near and beyond the survey size. We took this into account by adjusting the model with a correction factor that can be estimated from the double integration of the true correlation function over the survey area:

$$\begin{array}{c}\hfill {\displaystyle w_{\mathrm{IC}}=\frac{1}{{\mathrm{\Omega }}^2}\int w_{\mathrm{true}}\mathrm{d}{\mathrm{\Omega }}_1\mathrm{d}{\mathrm{\Omega }}_2.}\end{array}$$(16)

This integration can be carried out using the random-random pairs from the random catalog following Roche & Eales (1999)

$$\begin{array}{c}\hfill w_{\mathrm{IC}}=\frac{\sum w_{\mathrm{true}}(\theta )RR(\theta )}{\sum RR(\theta )},\end{array}$$(17)

where w_true(θ) comes from the HOD-predicted model. Finally, the model that we fit against the data is simply w(θ) = w_true(θ)−w_IC.

4.2. Ultraviolet luminosity functions

Since our Hα and [O III] samples are identical to those used in Covelo-Paz et al. (2025) and Meyer et al. (2024), we take the UVLF measurements presented in the respective papers. The completeness-corrected UVLFs are measured for the Hα CONGRESS at z ∼ 4.3, Hα FRESCO in two bins at z ∼ 5.3 and z ∼ 5.6, and the [O III] sample in two bins at z ∼ 7.1 and z ∼ 7.2. The advantage of these UVLFs is that they are measured from secure spectroscopic samples, therefore alleviating systematic biases that can arise from typical photometrically selected samples. These are shown as yellow/red symbols in Fig. 3.

However, these samples are limited at both the bright and faint end due the survey volume and sensitivity limits. In order to gain significant constraining power from the UVLFs, we also adopt the measurements from Bouwens et al. (2021) that use photometrically selected LBG over 1136 arcmin². We use the UVLFs in five bins at z ∼ 4, 5, 6, 7, 8 shown in gray symbols in Fig. 3. As described in Sect. 2.4, our highly flexible modeling allows us to simultaneously model all of these ten UVLFs consistently and include them under the same likelihood. Finally, we boost the UVLF errorbars to account for about 25% of cosmic variance uncertainty (e.g., Trapp et al. 2022; Sabti et al. 2022).

4.3. Fitting procedure

We fit the measurements in all redshift bins simultaneously by assuming a parametric redshift dependence of the model parameters. In this way, we model a smooth redshift evolution of ϵ(M_h, z). We adopted a linear dependence with redshift for all parameters following Muñoz et al. (2023)

$$\begin{array}{c}\hfill X=\frac{\mathrm{d}X}{\mathrm{d}z}z+C(X),\end{array}$$(18)

where X = [ϵ₀, log M_c, β, γ, σ_UV, log M_cut, log M_sat,α_sat].

We fit the models of the w(θ) and the UVLFs to our measurements using a Markov chain Monte Carlo (MCMC) approach while minimising χ²:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\chi }^2=& {\displaystyle \sum _i^{N_w}{({\mathit{w}}_i-{\mathit{w}}_i)}^T{C}^{-1}({\mathit{w}}_i-{\mathit{w}}_i)+}\hfill \\ \hfill & {\displaystyle \sum _i^{N_{\mathrm{\Phi }}}{\left(\frac{\mathrm{\Phi }(M_{\mathrm{UV},i})-\mathrm{\Phi }(M_{\mathrm{UV},i})}{{\sigma }_{\mathrm{\Phi }}}\right)}^2,}\hfill \end{array}\end{array}$$(19)

where the w terms are the measurement vectors containing w at θ, and $\stackrel{\sim }{w}$ and $\stackrel{\sim }{\mathrm{\Phi }}$ are the models for a given set of parameter values. The first term of Eq. (19) corresponds to the clustering likelihood and the second term to UVLF likelihood. The sums run over the different measurements in the three redshift bins for the clustering and ten in total UVLFs at different redshifts (five for the line-emitter samples and five for the LBG).

We carried out the fitting using the EMCEE code (Foreman-Mackey et al. 2013), that implements an affine-invariant ensemble sampler. We used 50 walkers for our five parameters and relied on the auto-correlation time τ to assess the convergence of the chain. To consider the chains converged, we require that the auto-correlation time is at least 50 times the length of the chain and that the change in τ is less than 1%. We discard the first 2 × max(τ) points of the chain as the burn-in phase and thin the resulting chain by 0.5 × min(τ). We impose flat priors on all parameters; for the mass parameters, the flat priors are on the log quantities.

For best-fit parameter values, we take the medians of the resulting posterior distribution, with the 16th and 84th percentiles giving the lower and upper uncertainty estimates. We compute the confidence intervals on all HOD model-derived quantities by computing it from 500 randomly drawn samples from the posterior, and taking the 1σ percentiles. The posterior median and 1σ uncertainties of the main model parameters X are given in Table B.1, for the three redshift bins of the Hα and [O III] samples. We also present the values for the parameters describing the linear with redshift dependence (Eq. (18)).

5. Results

5.1. Clustering of Hα and [O III] emitters

We measured galaxy clustering via the angular 2PCF for spectroscopically selected galaxies out to z ∼ 4.3 and z ∼ 5.4 for Hα emitters, and z ∼ 7.3 for [O III] emitters brighter than M_UV = −19.1 mag. These are shown in Fig. 2. The high resolution in NIRCam allowed us to probe very small scales deep into the one-halo regime, down to about 10 kpc at z ∼ 6; however, the limited survey area prevents us from probing scales larger than ∼1 Mpc at z ∼ 6.

Fig. 2. Angular 2PCF measured from the Hα CONGRESS and FRESCO samples at z ∼ 4.3 and z ∼ 5.4 and [O III] sample at z ∼ 7.3 selected at MUV < −19.1 mag. Comparison with literature measurements includes Harikane et al. (2022), Dalmasso et al. (2024a), and Paquereau et al. (2025) for samples selected at similar redshifts and MUV as our work.

The amplitude of the 2PCF for the Hα emitters increases mildly from z ∼ 4.3 to z ∼ 5.4, and continues to increase for the [O III] emitters at z ∼ 7.3. The shape of the 2PCF for the three samples resembles a power law, in contrast to the characteristic two-component behavior of the one- and two-halo terms. This is not surprising, given the fact that this is the regime where the non-linear scale-dependent bias can amplify the power at intermediate ∼100 − 500 kpc scales (Jose et al. 2017). Furthermore, relatively low-number statistics prevent us from accurately revealing the shape of the 2PCF for these samples. Finally, cosmic variance systematics (e.g., overdensities in the field) can have an influence on both the shape and amplitude of the 2PCF. Indeed, recent spectroscopic redshift searches for overdensities in the GOODS fields have revealed several significant overdensities at z ∼ 6 − 7 (Helton et al. 2024; Meyer et al. 2024; Covelo-Paz et al. 2025; Herard-Demanche et al. 2025).

In Fig. 2, we also show a comparison with recent measurements from the literature. Dalmasso et al. (2024a) measure clustering of LBG in the JADES survey, and as such their sample overlaps with ours. We compare with their measurements at z ∼ 5 and z ∼ 7, but for fainter samples than ours with mean UV magnitudes of ⟨M_UV⟩∼ − 19.2, compared to ours ⟨M_UV⟩∼ − 20. In general, there is good agreement, with differences that are likely due to the redshift and M_UV selection; for example, at z ∼ 7 our 2PCF has expectedly larger amplitude because of brighter sample. Harikane et al. (2022) measure the 2PCF of photometrically selected LBG in the HSC Subaru Strategic Program (SSP) survey over ∼300 deg² at 4 < z < 7. There is relatively good agreement with slightly higher amplitude of their brighter, M_UV ≲ −20 sample, compared to our M_UV < −19.1 mag sample. Finally, we find excellent agreement with Paquereau et al. (2025) measurements in the ∼0.5 deg² COSMOS-Web survey for stellar mass selected samples of log(M_⋆/M_⊙) > 8.5. This is comparable to our M_UV threshold of −19.1 mag assuming the Song et al. (2016) mass-to-light ratio. Slight differences likely arise from the redshift selection.

This consistency of the 2PCF with the literature is indeed reassuring and showcases that photo-z selected samples can yield robust and unbiased 2PCF, with an advantage of compiling significantly larger samples than what is possible with spec-z. This means that future implementation of our modeling framework will be possible using 2PCF measured from photo-z selection, therefore reducing errorbars and extending the analysis to z > 9.

5.2. UVHMR × HOD model fitting

We implemented the theoretical framework presented in Sect. 2 to fit simultaneously the measurements of the 2PCF in the three redshift bins (Sect. 4.1, Sect. 5.1) and the ten UVLFs across the 3.8 < z < 9.0 range (Sect. 4.2). In Fig. 3 we show the resulting median (solid lines) and 1σ uncertainty (envelope) computed from sampling the posterior distribution of the model parameters. The top row shows the best-fit models of the 2PCF, while the bottom row shows those of the UVLF at the mean redshift of the corresponding bin (e.g., ⟨z⟩∼4.3,  5.4 and 7.3). We show only one UVLF model for clarity, while in Fig. C.1 we show the individual fits for each measurement.

Fig. 3. Best-fit models of the 2PCF and UVLF. Top row: Angular 2PCF measured from the Hα CONGRESS and FRESCO samples at z ∼ 4.3 and z ∼ 5.4 and [O III] sample at z ∼ 7.3 (orange points). Bottom row: UVLF for the same line emitter samples from Covelo-Paz et al. (2025), Meyer et al. (2024) along with the photometrically selected samples from Bouwens et al. (2021). The blue curves and envelopes show the median models and 1σ uncertainty computed from the posterior. In the case of the UVLF, for simplicity we show only the models in the same redshift bins as the 2PCF measurements. We show the individual fits for each measurement in Fig. C.1.

The models show good fits of the data within the uncertainties for the 2PCF and UVLFs for all redshifts. In the case of the 2PCF, at the intermediate scales around the one-halo to two-halo transition (θ ∼ 10″) the data show slightly higher amplitudes than the data, likely due to imperfect modeling of the non-linear and scale-dependent halo bias (see Sect. 2.5). The overall good fit showcases that our framework can successfully model both the UVLF and 2PCF and derive physical insights, which we present in the following.

5.3. Mean halo mass and galaxy bias

The HOD modeling of the 2PCF and UVLF allowed us to infer the mean mass of the DM halo that hosts our line emitter samples and their bias (Sect. 2.6). Figures 4 and 5 show the mean halo mass and galaxy bias for the Hα samples at z ∼ 4.3 and 5.4 and [O III] at z ∼ 7.3.

Fig. 4. Mean halo mass for the three line emitter samples: Hα samples at z ∼ 4.3 and 5.4 and [O III] at z ∼ 7.3 with MUV < −19.1 mag. We compare with measurements from Harikane et al. (2016) and Harikane et al. (2022) for photometrically selected LBG. The colored regions are derived from the EVS (Lovell et al. 2023) formalism and indicate the confidence intervals to observe the most massive halo in the GOODS volume (2 × 60 arcmin2) within the ΛCDM model. The dotted line marks the median of the EVS distribution of the maximum plausible halo mass.

This analysis shows that the line emitters samples from our work are hosted in halos of log(M_h/M${}_{\odot })=11.{46}_{-0.04}^{+0.04}$, $11.{24}_{-0.04}^{+0.03}$, $10.{99}_{-0.06}^{+0.03}$ for the Hα and [O III] samples at z ∼ 4.3, 5.4 and 7.3 respectively. These halos are less massive by ∼0.2 dex compared to those hosting similarly bright LBG from the Harikane et al. (2016) analysis.

We interpreted the derived halo masses with the extreme value statistics (EVS; Lovell et al. 2023) formalism. EVS is a probabilistic approach to estimating the PDF$(M_{\mathrm{h}}^{\max })$ of observing the most massive halo at a given redshift and given cosmological volume. In Fig. 4 we show the confidence intervals of the PDF$(M_{\mathrm{h}}^{\max })$ derived from EVS for our GOODS area (2 × 60 arcmin²). The masses of the halos hosting our line emitter samples are lower than the most massive halo expected in GOODS by ∼0.7 dex. This is unsurprising, given the fact that we measure the mean halo mass hosting a sample of normal, star-forming galaxies with M_UV < −19.1 mag. We note that the Harikane et al. measurements are made in a much larger volume of the HSC-SSP and are not comparable with these EVS intervals derived for GOODS.

In Fig. 5 we show the galaxy bias as a function of redshift and M_UV. The blue empty circles correspond to the bias for the three line emitter samples at the corresponding median redshift. We measure a galaxy bias of $b_{\mathrm{g}}=3.{98}_{-0.15}^{+0.10}$, $5.{00}_{-0.18}^{+0.14}$, $7.{58}_{-0.36}^{+0.22}$, for the Hα and [O III] samples at z ∼ 4.3, 5.4 and 7.3 respectively. Since we fit the continuous redshift dependence of our model, we can derive the redshift evolution of the galaxy bias as a function of M_UV. This is shown in the blue lines in Fig. 5; however, we note that this is an extrapolation beyond the z and M_UV regime we probed in our work. Our results show that the bias increases with both redshift and luminosity, as expected from the theory (e.g., Kaiser 1984). Our measurements are comparable those of similarly bright LBGs from Dalmasso et al. (2024a). The stellar mass selected sample of Paquereau et al. (2025) also shows comparable values for the bias, although slightly elevated (by ≲1) at 5 < z < 8. This can be due to differences in the redshift and M_UV selection as well as in the modeling.

Fig. 5. Galaxy bias as a function of redshift and MUV. The empty circles with error bars show the bias for the three line emitter samples Hα at z ∼ 4.3 and 5.4, and [O III] at z ∼ 7.3 with MUV < −19.1 mag. The solid blue lines show the bias as a function of redshift for four MUV thresholds computed from our model. We compare with measurements from Harikane et al. (2016), Dalmasso et al. (2024b) and Dalmasso et al. (2024a) for photometrically selected LBG, as well as Paquereau et al. (2025) for stellar mass selected, normal galaxies.

5.4. Ultraviolet magnitude-halo mass relation

Our framework allows us to measure the M_UV − M_h relationship (UVHMR) and its evolution with redshift. In Fig. 6 we show our results on the UVHMR at the median redshifts of our three line emitters samples. However, in principle, since we fit for the redshift evolution of our model parameters, we have information of the continuous evolution of the UVHMR over the redshift range that we probe in this work, as well as extrapolations beyond. In this work, we only probe a limited range in M_UV and correspondingly in M_h with the 2PCF, which is marked by the solid lines and bold envelopes in Fig. 6. The transparent envelopes extend to the minimum M_UV that is probed by the UVLF. We show the dust-attenuated M_UV.

Fig. 6. Relationship of MUV − Mh (or UVHMR) for the redshift bins of our three line emitter samples. The solid lines and bold envelopes mark the MUV range that we probe with the 2PCF, while the transparent envelopes mark the range probed by the UVLF. We compare with the model by Mason et al. (2015, 2023) that includes dust attenuation in dashed lines and the observational measurements from Harikane et al. (2022) in dots.

The UVHMR shows the typical monotonic increase of UV luminosity with increasing M_h. As the halo mass increases, the slope of the UVHMR decreases, with a pivot mass at M_h ∼ 3 × 10¹¹ M_⊙. This is a similar pivot mass scale as the galaxy SHMR (Wechsler & Tinker 2018). The UVHMR evolves with redshift such that halos of fixed mass host brighter galaxies at earlier epochs. The high mass and luminosity end becomes shallower at lower redshifts due to the more important effect of dust attenuation. The relatively large uncertainties prevent us from inferring any potential mass or luminosity dependence of the redshift evolution.

We compare our measurements with the theoretical model from Mason et al. (2015, 2023) and the observational measurements from Harikane et al. (2022), which also include dust-attenuation. In the M_UV and M_h regime probed by our work, our UVHMR shows shallower slopes for all redshift bins, indicating luminosities that increase faster with halo mass. At fainter luminosities, the constraints from UVLF-only are in better agreement with the Mason et al. (2015, 2023) model.

5.5. The instantaneous star formation efficiency

The central feature of our modeling is the parametric form of the instantaneous SFE, and its dependence on halo mass and redshift, ϵ(z, M_h) (Eq. (2)), that we constrain from observations of the UVLF and 2PCF. As a reminder, ϵ is defined as the ratio between the SFR and the halo accretion rate times the universal baryon fraction $ϵ={\dot{M}}_{\star }{\dot{M}}_{\mathrm{h}}^{-1}{f}_{\mathrm{b}}^{-1}$ (Eq. (1)), and describes the efficiency of converting gas to stars.

Figure 7 shows our results on the SFE as a function of halo mass for the three redshift bins of our three line emitter samples. The solid lines and bold envelopes indicate the M_UV (and correspondingly M_h) range probed with the 2PCF, while the transparent envelopes mark the range probed by the UVLF. We note that the upper M_h range is derived from the M_UV–M_h relation, and it can extend to halo masses beyond the maximum mass expected in the GOODS volume, as indicated in Fig. 4. Therefore, we caution the interpretation of these results at the highest masses that can be considered as extrapolation.

Fig. 7. Star formation efficiency as a function of halo mass for the redshift bins of our three line emitter samples. The solid lines and bold envelopes mark the MUV range that we probe with the 2PCF, while the transparent envelopes mark the range probed by the UVLF. We compare with the observational measurements using 2PCF and HOD from Harikane et al. (2022) in dotted lines, empirical model by Tacchella et al. (2018) in dashed line, and results from cosmological simulations by Ceverino et al. (2024) and Feldmann et al. (2025) in solid colored lines.

The SFE shows the characteristic dependence with halo mass – increasing up to a peak halo mass scale of M_h ∼ 3 × 10¹¹ M_⊙ and decreasing for halos more massive than this peak mass. This characteristic shape of the SFE-M_h relation is typically explained by the effect of stellar and AGN feedback suppressing star-formation at the low- and high-mass end (e.g., Silk & Mamon 2012).

As a function of redshift, our results show very little evolution in the 4 ≲ z ≲ 7 range. At the low-mass end (M_h < 2 × 10¹¹ M_⊙) the SFE shows a decrease by about 0.3 dex from z ∼ 4.3 to z ∼ 7.3. However, at the high-mass end (M_h > 10¹² M_⊙) there are tentative indications that the trend is reversing, and the SFE increases with redshift, albeit with large uncertainties. This trend becomes clearer when we look at the redshift dependence of the best-fit parameters in Fig. 8. The normalization of the SFE, ϵ₀, shows a mild increase with redshift with a positive slope $\mathrm{d}{ϵ}_0/\mathrm{d}z=0.{02}_{-0.03}^{+0.06}$. However, the lower uncertainty limit is also consistent with zero slope, therefore making it difficult to draw a robust conclusion. The peak mass scale, M_c, also increases with redshift with a slope $\mathrm{d}\log M_{\mathrm{c}}/\mathrm{d}z=0.{14}_{-0.10}^{+0.12}$, meaning that the SFE shifts toward high masses, which creates the decrease of the SFE with redshift at the low-mass end. The high-mass end slope γ decreases ($\mathrm{d}\gamma /\mathrm{d}z=-0.{06}_{-0.12}^{+0.10}$), causing the SFE to increase with redshift at a fixed halo mass larger than M_c. This suggests increasing efficiencies of more massive halos in the early Universe. However, this remains a tentative interpretation, given the fact that the 1 σ upper limit on dγ/dz is also consistent with an increasing slope. Additionally, the increase of M_c with redshift is primarily driven by the measurements in the z ∼ 7.3 bin, where the 2PCF is relatively noisy. As shown in Fig. A.1, lower clustering amplitude would lower M_c. These uncertainties highlight the need for larger galaxy samples for more precise clustering measurements in order to conclusively infer the redshift evolution (or lack thereof) of the SFE.

Fig. 8. Resulting values of the model parameters. These values were obtained by drawing 1000 samples from the posterior and evaluating the models at a given z using Eq. (18). The orange points and errorbars mark the redshift bins of the 2PCF measurements. The dashed lines mark the z-parametrized functions (Eq. (18)) for our model parameters, evaluated at the median posterior.

We compare our results with the observational measurements using the 2PCF and HOD modeling from Harikane et al. (2022), the empirical model by Tacchella et al. (2018) and the results from the cosmological simulations FIRSTLIGHT (Ceverino et al. 2024) and FIREBOX^HR (Feldmann et al. 2025). Our analysis is closest to that of Harikane et al. (2022) and the SFE is in relatively good agreement in a 1 dex range around the peak halo mass. At lower (< 10¹¹ M_⊙) and higher (> 10¹² M_⊙) halo masses, our measurements indicate higher and lower SFE, correspondingly. The Tacchella et al. (2018) empirical model assumes a universal SFE constant with redshift and has a lower peak halo mass by about 0.2 − 0.3 dex and higher efficiencies at < 3 × 10¹¹ M_⊙. Compared to the cosmological simulations, Feldmann et al. (2025) find a non-evolving SFE in FIREBOX^HR, although they can only probe the very low mass end due to the limited simulation volume. Their results are in close agreement to ours at z ∼ 4.3, with the difference that our results indicate a mild redshift evolution. Additionally, the slope of their SFE−M_h relation is considerably shallower and can have important implications for derived galaxy statistical quantities such as the UVLF and the SFR density, which we discuss in Sect. 6. On the other hand, Ceverino et al. (2024) find a SFE that evolves with redshift in FIRSTLIGHT. Their SFE in the only overlapping z = 6 − 7 bin shows elevated efficiencies compared to ours by about 0.2 − 0.4 dex, with an indication of a peak occurring at a similar halo mass scale.

6. Discussion

6.1. Stochasticity versus star formation efficiency

Crucially, our model allows us to break the degeneracy between the high SFE versus high stochasticity scenarios by quantifying both ϵ and σ_UV observationally. In the following, we discuss the implications of our measurements on these two scenarios.

Our results indicate that there is some level of stochasticity, described by modest values of the scatter σ_UV ∼ 0.65. From our measurements in 4 ≲ z ≲ 8 range, we find that σ_UV remains roughly constant with redshift with a slope $\mathrm{d}{\sigma }_{\mathrm{UV}}/\mathrm{d}z=-0.{03}_{-0.08}^{+0.07}$. In Fig. 9 we show the redshift evolution of σ_UV inferred from our work and compared with the literature. We show in blue points with errorbars our measurements at the mean redshifts of the three line emitter samples. The blue solid line and envelope show the best-fit and 1σ uncertainty as a linear function of redshift. We note that beyond z ≳ 8 these are extrapolations since we do not include any observational (UVLF) constraints.

Compared to the literature, there is a relatively good agreement out to z ∼ 9 from independent approaches. Muñoz et al. (2023) fit the UVLF with a similar modeling and find consistent values with ours. Their linear with redshift parametrization of σ_UV is also consistent with our measurement. However, when they fit for independent z-bins they find that σ_UV needs to increase sharply at z > 10 to fit the observed JWST UVLFs, a regime currently unconstrained by our framework.

Indeed, there is observational evidence that the transition from secular to stochastic star formation occurs at z ≳ 6 − 9 (Ciesla et al. 2024; Cole et al. 2025; Endsley et al. 2024; Looser et al. 2025; Dressler et al. 2023; Langeroodi & Hjorth 2024). Ciesla et al. (2024), for instance, study SF histories and show that at z ∼ 9 about 87% of massive galaxies have evidence for a stochastic star formation in the last 100 Myr, while only 15% at z < 7. However, this is not reflected in the redshift independent σ_UV ∼ 1.2 that they infer (filled squares in Fig. 9). This is likely due to the different definition of σ_UV which is simply the M_UV dispersion computed from SED fitting of their sample. Additionally, they estimate σ_UV in an alternative way by decomposing it into different components and using the SFR-M_⋆ dispersion, and find σ_UV ∼ 0.68, in agreement with our results (empty squares in Fig. 9). However, we note that σ_UV (i.e., the scatter in UV magnitude at a given halo mass) that we analyze in our work is not directly comparable with empirically measured values of SF burstiness, usually derived from ratios of SFR from emission lines and UV continuum and/or non-parametric SFH modeling.

Fig. 9. Scatter in the MUV−Mh relation, represented as σUV, as a function of redshift. The blue points with errorbars mark σUV at the median redshift of our three line emitter samples, while the blue line and envelope mark the best-fit and 1σ uncertainty as a linear function of redshift. From the literature compilation, the colored boxes at the right axis show the σUV for a range of Mh independent of redshift from the FIRE and FIREBOXHR simulations (Sun et al. 2023; Feldmann et al. 2025). For Muñoz et al. (2023), we show the 1σ contours for σUV obtained for independent z-bin fits; while the dashed purple line shows their best-fit linear function with redshift. For Ciesla et al. (2024) we show results from two different methods in the filled and empty squares.

In the cosmological hydrodynamical simulations FIRE-2 (Sun et al. 2023) and FIREBOX^HR (Feldmann et al. 2025), σ_UV is a decreasing function of M_h in the range ∼1.40 − 0.35 at 9 < log(M_h/M_⊙) < 12 for FIRE-2, and ∼1.55 − 0.70 at 9 < log(M_h/M_⊙) < 11 for FIREBOX^HR, both independent of redshift. These are indicated in the colored boxes on the right axis of Fig. 9. Our results align well with these simulations, as our sample primarily resides in halos with 11 < log(M_h/M_⊙) < 11.5 (see Sect. 5.3), a halo mass regime that would correspond to similar σ_UV values in the simulations. Finally, in the SERRA simulation Pallottini & Ferrara (2023) find a σ_UV at z ∼ 7.7 consistent with ours (diamond symbol in Fig. 9).

Gelli et al. (2024) discuss the physical motivations behind a halo mass and redshift dependent σ_UV. Under the assumption that the stochasticity is driven by the ability of halos to retain the gas despite the feedback processes, then it would scale inversely to the escape velocity, i.e., ${\sigma }_{\mathrm{UV}}\sim {\mathit{v}}_{\mathrm{esc}}^{-1}\sim M_{\mathrm{h}}^{-1/3}{(1+z)}^{-1/2}$. This would mean σ_UV decreases with redshift, which is qualitatively consistent with our measurement of $\mathrm{d}{\sigma }_{\mathrm{UV}}/\mathrm{d}z=-0.{03}_{-0.08}^{+0.07}$.

Our work adds further observational evidence that stochasticity alone cannot explain the observed UVLFs and galaxy clustering measured via the 2PCF out to z ∼ 9. However, the regime in which there is some disagreement with the literature remains at z > 9, where the transition to stochastic SF is identified. For example, at the highest redshifts (z > 10 − 16) Shen et al. (2023) and Kravtsov & Belokurov (2024) find that values as high as σ_UV ∼ 2 are required to explain the observed UVLF. This is also found by Muñoz et al. (2023) as mentioned earlier. This redshift regime is beyond the scope of this work, but future application of our framework to larger datasets extending to z > 9 promises to provide valuable observational constraints.

We carried out a simple exercise using our theoretical framework to investigate the high stochasticity versus SFE scenarios at z > 9. We use the model parameters at z = 8 which corresponds to the upper redshift limit of our data, and compute the UVLF at z > 8. Since the parameters fix the HOD at z = 8, this means that the redshift evolution is primarily driven by the HMF evolution. Figure 10 shows the UVLF predictions from our model in solid lines, and compared to measurements from the literature (symbols). Unsurprisingly, at z ≥ 9 our model consistently underpredicts the observed UVLF with the difference increasing with redshift.

Fig. 10. Predictions of UVLF out to z = 13 from our model constrained on z < 8 data and compared to measurements from the literature. The solid lines correspond to the model evaluated at different redshifts from the best-fit linear redshift dependence of the model parameters fit on 4 ≲ z ≲ 8 data. The dashed and dotted lines illustrate the two extreme scenarios of increasing SFE versus increasing stochasticity and show the model where we tune by hand ϵ0 and σUV to approximately match the observations, while keeping the other parameters nominal. The filled area shows the effect of dust attenuation on the UVLF, with the upper thin lines marking the unatenuated UVLF.

We note that in simply increasing ϵ₀, our empirical modeling remains agnostic about other physical mechanisms that can be triggered by increasing the SFE. These include stellar feedback, gas ejection, accretion modulation, increase in metallicity and dust content, which can act to suppress the SFE and impact the resulting UVLF. In this sense, the SFE that we assume by increasing ϵ₀ should be considered as the effective SFE that encapsulates all the effects triggered by the increase in ϵ₀. However, since our modeling is agnostic to these, the derived UVLF are more consistent with upper limits.

We illustrate how two extreme cases would fit the UVLF. One case involves increasing only the maximum SFE ϵ₀, and the second involves increasing only the stochasticity via σ_UV. We tuned both parameters by hand at each z to roughly match the UVLF while keeping the other parameters to their nominal values (e.g., Fig. 8). At z = 10, ϵ₀ ∼ 35% (compared to the nominal value of ∼25%) would be required to match the UVLF in the first case, while the second case requires σ_UV = 1.1, which is almost double the nominal. At the highest redshifts, z ≥ 12, both cases require extreme values, increasing ϵ₀ to ∼95% and σ_UV to ∼1.8.

This exercise shows that neither of the two scenarios alone is likely to explain the observed UVLFs, since they require extreme values for ϵ₀ and σ_UV. Such values are typically not predicted by current simulations, although some models (e.g., feedback-free starbursts; Torrey et al. 2017; Grudić et al. 2018; Dekel et al. 2023; Li et al. 2024; Renzini 2023) can produce such high efficiencies. However, taken together, the impact of different SFEs, and plausible increases in stochasticity, may well result in those high abundances of luminous galaxies. Therefore, the solution is likely a contribution from both the stochasticity and the SFE along with its dependence on halo mass (e.g., Gelli et al. 2024). A mass dependence of σ_UV that we neglect in our work can have important contributions, but since our framework allows us to implement a halo mass dependence, we will explore this in future work.

To highlight the importance of the SFE dependence on halo mass, we can compare with the results by Feldmann et al. (2025) based on the FIREBOX^HR simulation. They find that a non-evolving, weakly mass dependent SFE (shown in Fig. 7) can explain the observed UV luminosity density at z > 10. The reason for this is the shallow slope of their SFE-M_h relation at 9 < log(M_h/M_⊙) < 11 that increases the contribution from the more numerous lower mass halos at higher redshifts – this can integrate the UVLF to the observed values. This suggests that the global SFE in low-to-intermediate mass halos are higher than previously determined empirically (e.g., Tacchella et al. 2018; Harikane et al. 2022) or in simulations such as FIRSTLIGHT (Ceverino et al. 2024). Compared to these models, our empirical results also show shallower SFE-M_h slopes at the low-mass end, but slightly steeper than Feldmann et al. (2025), which given the steep increase of the HMF at low masses, could be sufficient to explain why our SFE cannot reproduce the observed UVLF, and correspondingly the UV luminosity density at z > 10.

Higher SFE due to a shallow SFE-M_h slope at low-to-intermediate masses is certainly one viable explanation for the observed abundances of bright galaxies at z > 10. However, it remains to be tested against other galaxy observables such as the 2PCF and SMF and extended at log(M_h/M_⊙) > 11. Indeed, the high-mass regime can be crucial in revealing the full picture. For example, in our case the mild decrease of the SFE with z at the low-mass end is balanced by an increase at the high-mass end and this kind of evolution explains both the 2PCF and UVLF at 4 < z < 8. The importance of the bright and massive end in constraining the models and distinguishing between physical mechanisms at z > 10 is illustrated in Fig. A.1 in the context of the SFE model parameters and σ_UV. Figure A.1 shows that at z ∼ 6, varying ϵ₀ = 0.2 ± 75% changes Φ [Mpc⁻³ mag⁻¹] by ∼2 dex at M_UV = −22 mag, but < 1 dex at M_UV = −19 mag. Similarly, variation of σ_UV = 0.6 ± 75% results in ΔΦ [Mpc⁻³ mag⁻¹]∼2.5 dex at M_UV = −22 mag, but only ∼0.2 dex at M_UV = −19 mag.

We can conclude from this that it is important to determine the SFE in a large M_h range, spanning several orders of magnitude and beyond the characteristic halo mass of ∼10¹² M_⊙ at different epochs in order to understand its (non)evolution with redshift. This can be made possible by implementing our framework to wider and deeper surveys to probe the bright and faint end accordingly, in a wedding-cake approach that can include the whole JWST extragalactic surveys archive. This is necessary in order to provide suffficiently large samples for clustering measurements as well as robust measurements of the bright end of the UVLF.

6.2. Model degeneracies and how to break them

The framework that we developed in this work is highly predictive, but relies on several physical assumptions that can be degenerate. However, our framework is flexible enough to be extended to model and constrain these against additional observables.

As mentioned in Sect. 2, the luminosity-to-SFR conversion factor κ_UV is fully degenerate with the maximum SFE, ϵ₀ (Inayoshi et al. 2022; Muñoz et al. 2023; Shen et al. 2023). This means that it is effectivelly the ϵ₀/κ_UV ratio that determines the inferred SFE, when constrained against the UVLF. However, κ_UV might not remain universal at all epochs. For example, different IMF assumption can change κ_UV from 1.15 × 10⁻²⁸ (M_⊙ yr⁻¹)/(erg s⁻¹), for a Salpeter IMF (adopted in this work) to κ_UV to 0.72 × 10⁻²⁸ (M_⊙ yr⁻¹)/(erg s⁻¹), for a Chabrier IMF (Madau & Dickinson 2014). Additionally, lower metallicities and/or top-heavy initial mass functions, as expected for Population III stars, can result in an even lower κ_UV (e.g., Bromm et al. 2002; Inayoshi et al. 2022), that would drive ϵ₀ down in order to reproduce the same UVLF. Finally, κ_UV depends also on the duration of the previous star formation (Wilkins et al. 2019) and on the age of the stellar population, which can also change at earlier times.

This has been implicitly explored in Donnan et al. (2025) who assume an integrated SFE³ constant with time, tuned to reproduce the SMF at z ∼ 7, and fit for a M_⋆ − M_UV relationship that fits the UVLF at 6 < z < 13. The changing mass-to-light ratio evokes younger and thus brighter stellar populations at earlier times which effectively tweaks κ_UV toward lower values (κ_UV ∝ 1/L_UV), while forcing ϵ₀ constant. This showcases the ϵ₀/κ_UV degeneracy.

Stellar masses and SMFs provide complementary probes that can help break this degeneracy because they do not depend on κ_UV. This is shown by the fact that the stellar mass is a result of the integrated SFR times a mass loss function⁴, $M_{\star }=\int \mathrm{d}t{\dot{M}}_{\star }(1-f_{\mathrm{r}})$. This also means that the SMF can be modeled under the same parametrization as the UVLF and the 2PCF. M_⋆ being independent of κ_UV means that the SMF can break the ϵ₀/κ_UV degeneracy. However, observational stellar mass estimates depend on the IMF assumption during SED fitting (e.g., Conroy et al. 2009) which can complicate this application.

The SMF is an important probe to consider because accurate measurements of the instantaneous SFE need to be able to reproduce the observed SMF. In the context of the SMF, it is the integrated SFE, ${ϵ}_{\star }=M_{\star }{M}_{\mathrm{h}}^{-1}{f}_{\mathrm{b}}^{-1}$, that describes the relationship between the galaxy stellar mass and host halo mass (SHMR). The SHMR has been found to show some level of evolution with redshift by numerous works in the past (e.g., Conroy & Wechsler 2009; Behroozi et al. 2010; Moster et al. 2010, 2013; Behroozi et al. 2019; Girelli et al. 2020; Shuntov et al. 2022). Additionally, recent measurements from JWST at z ≳ 7 using abundance matching and HOD analysis show an important evolution toward high efficiencies at early times (Shuntov et al. 2025; Paquereau et al. 2025). This seems in contrast with the indications of a very mild evolution of the instantaneous SFE. This suggests that additional modeling components and observational probes (such as the SMF) need to be incorporated in order to accurately inform our models and break degeneracies such as ϵ₀/κ_UV and σ_UV. Therefore, we identify that the way forward in unveiling a complex picture of intertwined processes is constraining physical models on several complementary probes including the 2PCF, UVLF, SMF and potentially others.

Finally, negligible dust attenuation in the early Universe can also be responsible for the observed abundances of bright galaxies (Ferrara et al. 2023; Mason et al. 2023). In our work we used the canonical A_UV − β (Meurer et al. 1999) and β − M_UV (Bouwens et al. 2014) relations to attenuate the M_UV when fitting the observed UVLFs, and extrapolating to predict the UVLF at z > 9. However, this is uncertain because there is evidence from JWST that the β − M_UV relationship ma changes toward bluer β slopes at z > 8 (e.g., Cullen et al. 2023; Topping et al. 2024), indicating negligible dust attenuation, as well as indications of the opposite (e.g., Saxena et al. 2024). We showcase the effect of dust attenuation on the UVLF in Fig. 10 with the filled areas, where the upper thin lines mark the unatenuated UVLF. This makes it clear that the level of dust attenuation adds additional degeneracy in the models. Interestingly, in the case of very high stochasticity, the attenuation in the UVLF appears to be negligible. This is unsurprising, because low luminosity galaxies scattered in high-mass halos will indeed have little attenuation. Incorporating the SMF in our formalism would allow for relaxing of the assumptions of dust attenuation and for them to be fit directly. However, the difficulty of measuring accurate stellar masses at high redshift is likely to limit that application.

7. Conclusions

In this paper, we have presented a theoretical framework developed to combine measurements of UVLF and galaxy clustering via the 2PCF into a probe of the galaxy-halo connection. We based this powerful framework on the CLF in order to model the UV luminosity–halo mass relationship (which we called UVHMR) in combination with the HOD. It is a highly flexible adaptation of the HOD framework that allows for the use of independent redshift binning schemes for each probe (UVLF and 2PCF) and parametrizes with redshift the model parameters. As such, when fit on the UVLF and 2PCF, it provides empirical constraints on the redshift evolution of the SFE-M_h and M_UV–M_h relationships, stochasticity, host halo masses, and galaxy bias. Importantly, this formalism allows for the incorporation of independent probes to test and break degeneracies between different physical models proposed to explain the surprising abundances of bright and early galaxies seen by JWST.

To implement this framework, we measured the 2PCF of spectroscopically selected Hα and [O III] line emitters from FRESCO and CONGRESS JWST NIRCam/grism surveys in 124 arcmin² in the GOODS-North and GOODS-South fields. Our UVHMR×HOD modeling successfully fit the 2PCF and UVLF measurements in three and ten redshift bins at 3.8 < z < 9, respectively. Our main findings are summarized as follows:

• The clustering of M_UV < −19.1 mag Hα and [O III] emitters shows 2PCF amplitudes in agreement with similarly bright LBG and normal galaxy samples from the literature. The 2PCF amplitude increases mildly from z ∼ 4.3 to z ∼ 7.3, and the shape resembles a power law, thus requiring modeling the non-linear scale-dependent halo bias in order to boost the amplitude at quasi-linear scales (∼100 − 500 kpc).

• The dark matter halos that host the Hα and [O III] samples at z ∼ 4.3, 5.4, and 7.3 have mean masses of log$(M_{\mathrm{h}}/M_{\odot })=11.{46}_{-0.04}^{+0.04}$, $11.{24}_{-0.04}^{+0.03}$, and $10.{99}_{-0.06}^{+0.03}$, respectively. The galaxy bias increases with redshift with values of $b_{\mathrm{g}}=3.{98}_{-0.15}^{+0.10}$, $5.{00}_{-0.18}^{+0.14}$, $7.{58}_{-0.36}^{+0.22}$ for the three samples, showing that these galaxies reside in high peaks of the underlying matter density field. However, EVS analysis revealed that their host halos do not correspond to extreme overdense environments for their respective epochs.

• The SFE-M_h relation rises gradually with halo mass, peaks at about 20% at M_h ∼ 3 × 10¹¹ M_⊙, and declines at higher halo masses. The SFE-M_h shows only a mild evolution with redshift and tentative indications that low-mass halos decrease but the high-mass halos increase in efficiency with increasing redshift. However, uncertainties are large and consistent with no evolution scenario.

• The scatter in the M_UV–M_h relationship, quantified by σ_UV, suggests a modest stochasticity of ∼0.7. σ_UV is constant with redshift and in agreement with the literature out to z ∼ 9.

• Extrapolating our model at z > 9 shows that a constant SFE−M_h fixed at z = 8 cannot reproduce the observed UVLF. A high maximum SFE or high stochasticity alone cannot explain the high abundances of luminous galaxies seen by JWST since they require extreme and unlikely values (e.g., ϵ₀ ∼ 0.95, σ_UV ∼ 1.8). Additionally, negligible dust attenuation can result in an up to 0.5 − 1 dex increase of the abundance of luminous galaxies. Therefore, it is not implausible that the combination of high SFE, high stochasticity, and low dust attenuation could conspire to create galaxy formation events that differ from what is seen at later times, resulting in an unexpected abundance of highly luminous galaxies at early times.

Extending this analysis to z > 9 and including additional and complementary probes will be crucial in constraining the galaxy-halo connection and revealing the dominant physical mechanisms in early star formation. Judging from the consistency of our spec-z selected 2PCF measurements with the photo-z selected measurements from the literature, it is promising that photo-z selected samples throughout all JWST surveys will provide sufficiently large samples for 2PCF and UVLF measurements out to the highest redshifts (consistent with the predictions from mocks by Endsley et al. 2020). This will be possible because photometric selection allows for the selection of fainter samples, which are more numerous. However, survey width will also be crucial in constraining the UVLF and clustering of bright galaxies for which the model parameters are more sensitive and can thus shed light on the underlying physical mechanisms.

Data availability

To make our results transparent and facilitate comparison we provide all our measurements which can reproduce the figures in this paper in tabulated form at https://github.com/mShuntov/Clustering_and_HOD_of_Ha-and-OIII_emitters. and https://zenodo.org/records/15386639

Acknowledgments

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The raw data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. Some of the data products presented herein were retrieved from the Dawn JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center, which is funded by the Danish National Research Foundation under grant No. 140 (DNRF140). This work was made possible by utilizing the CANDIDE cluster at the Institut d’Astrophysique de Paris, which was funded through grants from the PNCG, CNES, DIM-ACAV, and the Cosmic Dawn Center and maintained by S. Rouberol. We thank Giorgos Nikopoulos for kindly sharing a neat UVLF literature compilation. We thank Anne Hutter for the valuable discussions. This work has received funding from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00072, as well as from the Swiss National Science Foundation (SNSF) through project grant 200020_207349. This work was further supported by funding from JWST-GO-01895, provided through a grant from the STScI under NASA contract NAS5-03127. This work made use of the following Python packages: numpy (Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), astropy (Astropy Collaboration 2022), HALOMOD (Murray et al. 2021)

References

Appendix A: UVLF, 2PCF, and bias dependence on model parameters

In Fig. A.1 we show how the UVLF, 2PCF and galaxy bias at z = 6 change with changing values of each of the model parameters. We vary one parameter at a time while keeping the others at their fiducial values (shown in the top right panel).

Fig. A.1. Dependence of the UVLF, 2PCF and galaxy bias at z = 6 on model parameters. Different rows show the UVLF (left), 2PCF (middle) an bias (right column) for different values of the seven model parameters σUV,  ϵ0,  log Mc,  β,  γ,  log Mcut,  log Msat, αsat. One parameter is varied at a time while keeping the others at their fiducial values (shown in the top right panel).

Appendix B: Best-fit parameter values

Appendix C: Fits of the ultraviolet luminosity function

In Fig. C.1 we show the UVLF fits for each individual measurement.

Fig. C.1. Fits of UVLF for each individual measurement, color coded accordingly. We plot the models only in the MUV range of the UVLF measurement. These show a good fit between the model and measurements in all z-bins.

Appendix D: Details on the halo occupation distribution-model derivation of the correlation functions

D.1. Clustering correlation function

Described briefly, a central component of the model is the galaxy-galaxy power spectrum. It can be separated into contributions from the clustering of galaxies within the same halo (one-halo term) and from galaxies between different halos (two-halo term):

$$\begin{array}{c}\hfill P_{\mathrm{gg}}(k,z)=P_{\mathrm{gg}}^{1\mathrm{h}}(k,z)+P_{\mathrm{gg}}^{2\mathrm{h}}(k,z).\end{array}$$(D.1)

These two components of the power spectrum of galaxies can be modeled under the HOD framework (Cooray & Sheth 2002; Berlind & Weinberg 2002) where the one-halo and two-halo terms are given by

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle P_{\mathrm{gg}}^{1\mathrm{h}}(k,z)=\frac{1}{{\overline{n}}_g^2}\int \mathrm{d}{M}_h\frac{\mathrm{d}n}{\mathrm{d}{M}_h}[{⟨N_s⟩}^2{u}_s^2(k)+2⟨N_s⟩u_s(k)],}\hfill \\ \hfill & P_{\mathrm{gg}}^{2\mathrm{h}}(k,z)=\hfill \\ \hfill & \frac{1}{{\overline{n}}_g^2}{\left[{\displaystyle \int \mathrm{d}{M}_h\frac{\mathrm{d}n}{\mathrm{d}{M}_h}{b}_h(M_h,z)[⟨N_c⟩+⟨N_s⟩u_s(k)]}\right]}^2{P}_{\mathrm{lin}}(k,z).\hfill \end{array}\end{array}$$(D.2)

In these equations, ${\overline{n}}_g=\int \mathrm{d}{M}_h\frac{\mathrm{d}n}{\mathrm{d}{M}_h}(⟨N_c⟩+⟨N_s⟩)$ is the mean number density of galaxies, b_h(M_h, z) is the large-scale halo bias which we chose the one given by Tinker et al. (2010). u_s(k) is the Fourier transform of the over-density profile of satellite galaxies, for which we assume that it follows the NFW profile Navarro et al. (1997) with a mass-concentration relation as calibrated by Duffy et al. (2008). The important thing to note here is that the power spectrum is defined in terms of the occupation distributions of centrals and satellites ⟨N_c⟩ and ⟨N_s⟩ specified by Eq. 8 & 9 and is where the HOD parametrization enters the model. The linear power spectrum P_lin(k, z) enters the two-halo term and dominates large scales. Details on deriving the spatial galaxy correlation function ξ(r) and the angular correlation function w(θ) using Limber (1953) equation are given in Blake et al. (2008).

All Tables

All Figures

	Fig. 1. Spatial distribution of the three line-emitting samples (Hα and [O III] from FRESCO and Hα from CONGRESS) in the GOODS-North and GOODS-South footprints (left and middle panels) that we consider in this analysis. The right panel shows their redshift distribution. The density histogram gray background in the left and middle panels shows the footprint of the random catalogs.
In the text

	Fig. 2. Angular 2PCF measured from the Hα CONGRESS and FRESCO samples at z ∼ 4.3 and z ∼ 5.4 and [O III] sample at z ∼ 7.3 selected at MUV < −19.1 mag. Comparison with literature measurements includes Harikane et al. (2022), Dalmasso et al. (2024a), and Paquereau et al. (2025) for samples selected at similar redshifts and MUV as our work.
In the text

Fig. 3. Best-fit models of the 2PCF and UVLF. Top row: Angular 2PCF measured from the Hα CONGRESS and FRESCO samples at z ∼ 4.3 and z ∼ 5.4 and [O III] sample at z ∼ 7.3 (orange points). Bottom row: UVLF for the same line emitter samples from Covelo-Paz et al. (2025), Meyer et al. (2024) along with the photometrically selected samples from Bouwens et al. (2021). The blue curves and envelopes show the median models and 1σ uncertainty computed from the posterior. In the case of the UVLF, for simplicity we show only the models in the same redshift bins as the 2PCF measurements. We show the individual fits for each measurement in Fig. C.1.

In the text

Fig. 4. Mean halo mass for the three line emitter samples: Hα samples at z ∼ 4.3 and 5.4 and [O III] at z ∼ 7.3 with MUV < −19.1 mag. We compare with measurements from Harikane et al. (2016) and Harikane et al. (2022) for photometrically selected LBG. The colored regions are derived from the EVS (Lovell et al. 2023) formalism and indicate the confidence intervals to observe the most massive halo in the GOODS volume (2 × 60 arcmin2) within the ΛCDM model. The dotted line marks the median of the EVS distribution of the maximum plausible halo mass.

In the text

Fig. 5. Galaxy bias as a function of redshift and MUV. The empty circles with error bars show the bias for the three line emitter samples Hα at z ∼ 4.3 and 5.4, and [O III] at z ∼ 7.3 with MUV < −19.1 mag. The solid blue lines show the bias as a function of redshift for four MUV thresholds computed from our model. We compare with measurements from Harikane et al. (2016), Dalmasso et al. (2024b) and Dalmasso et al. (2024a) for photometrically selected LBG, as well as Paquereau et al. (2025) for stellar mass selected, normal galaxies.

In the text

Fig. 6. Relationship of MUV − Mh (or UVHMR) for the redshift bins of our three line emitter samples. The solid lines and bold envelopes mark the MUV range that we probe with the 2PCF, while the transparent envelopes mark the range probed by the UVLF. We compare with the model by Mason et al. (2015, 2023) that includes dust attenuation in dashed lines and the observational measurements from Harikane et al. (2022) in dots.

In the text

Fig. 7. Star formation efficiency as a function of halo mass for the redshift bins of our three line emitter samples. The solid lines and bold envelopes mark the MUV range that we probe with the 2PCF, while the transparent envelopes mark the range probed by the UVLF. We compare with the observational measurements using 2PCF and HOD from Harikane et al. (2022) in dotted lines, empirical model by Tacchella et al. (2018) in dashed line, and results from cosmological simulations by Ceverino et al. (2024) and Feldmann et al. (2025) in solid colored lines.

In the text

	Fig. 8. Resulting values of the model parameters. These values were obtained by drawing 1000 samples from the posterior and evaluating the models at a given z using Eq. (18). The orange points and errorbars mark the redshift bins of the 2PCF measurements. The dashed lines mark the z-parametrized functions (Eq. (18)) for our model parameters, evaluated at the median posterior.
In the text

Fig. 9. Scatter in the MUV−Mh relation, represented as σUV, as a function of redshift. The blue points with errorbars mark σUV at the median redshift of our three line emitter samples, while the blue line and envelope mark the best-fit and 1σ uncertainty as a linear function of redshift. From the literature compilation, the colored boxes at the right axis show the σUV for a range of Mh independent of redshift from the FIRE and FIREBOXHR simulations (Sun et al. 2023; Feldmann et al. 2025). For Muñoz et al. (2023), we show the 1σ contours for σUV obtained for independent z-bin fits; while the dashed purple line shows their best-fit linear function with redshift. For Ciesla et al. (2024) we show results from two different methods in the filled and empty squares.

In the text

Fig. 10. Predictions of UVLF out to z = 13 from our model constrained on z < 8 data and compared to measurements from the literature. The solid lines correspond to the model evaluated at different redshifts from the best-fit linear redshift dependence of the model parameters fit on 4 ≲ z ≲ 8 data. The dashed and dotted lines illustrate the two extreme scenarios of increasing SFE versus increasing stochasticity and show the model where we tune by hand ϵ0 and σUV to approximately match the observations, while keeping the other parameters nominal. The filled area shows the effect of dust attenuation on the UVLF, with the upper thin lines marking the unatenuated UVLF.

In the text

	Fig. A.1. Dependence of the UVLF, 2PCF and galaxy bias at z = 6 on model parameters. Different rows show the UVLF (left), 2PCF (middle) an bias (right column) for different values of the seven model parameters σUV,  ϵ0,  log Mc,  β,  γ,  log Mcut,  log Msat, αsat. One parameter is varied at a time while keeping the others at their fiducial values (shown in the top right panel).
In the text

	Fig. C.1. Fits of UVLF for each individual measurement, color coded accordingly. We plot the models only in the MUV range of the UVLF measurement. These show a good fit between the model and measurements in all z-bins.
In the text