# Extreme-value modelling of the brightest galaxies at z9

Cameron Heather, Teeraparb Chantavat, Siri Chongchitnan, Joseph Silk

► To cite this version:

Cameron Heather, Teeraparb Chantavat, Siri Chongchitnan, Joseph Silk. Extreme-value modelling of the brightest galaxies at z9. Monthly Notices of the Royal Astronomical Society, 2024, 534 (1), pp.173-181. 10.1093/mnras/stae2051. hal-04590353

HAL Id: hal-04590353

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

Submitted on 18 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Extreme-value modelling of the brightest galaxies at $z \gtrsim 9$

Cameron Heather,<fnref n="1" />★ Teeraparb Chantavat,<fnref n="2" /> Siri Chongchitnan<fnref n="1" /> and Joseph Silk<fnref n="3" />,<fnref n="4" />,<fnref n="5" />

Accepted 2024 August 29. Received 2024 August 2; in original form 2024 April 21

## ABSTRACT

Data from the JWST have revealed an intriguing population of bright galaxies at high redshifts. In this work, we use extreme-value statistics to calculate the distribution [in ultraviolet (UV) magnitude] of the brightest galaxies in the redshift range $9 \lesssim z \lesssim 16$. We combine the generalized extreme-value approach with modelling of the galaxy luminosity function. We obtain predictions of the brightest galaxies for a suite of luminosity functions, including the Schechter and double power-law functions, as well as a model parametrized by the stellar formation efficiency $f_*$.

**Key words:** methods: statistical – galaxies: luminosity function, mass function – galaxies: star formation – galaxies: statistics.

## 1 INTRODUCTION

One of the primary objectives of contemporary astrophysics is to use the latest observational technologies to improve our understanding of the physics of galaxy formation and evolution (Bromm & Yoshida 2011; Stark 2016; Ouchi, Ono & Shibuya 2020; Robertson 2022). In pursuit of this goal, a new generation of multiwavelength surveys has shed a new light on galaxies at unprecedentedly high redshifts (Gardner et al. 2006; Laureijs et al. 2011; Fomalont et al. 2015; Spergel et al. 2015; DeBoer et al. 2017; Marchetti, Serjeant & Vaccari 2017).

One of the essential ingredients for studying galaxy evolution is the rest-frame ultraviolet (UV) luminosity function (LF), which characterizes the density of galaxies within a given volume based on their luminosity. The UV LF is well constrained within the redshift range $2 \lesssim z \lesssim 10$ based on observations from the *Hubble Space Telescope* and *Spitzer*, predating the launch of the JWST in 2021. Many studies at that time suggest a rapid decline of the UV LF at the bright end and find that it is well described by the Schechter function (Ellis et al. 2013; Madau & Dickinson 2014; Bouwens et al. 2015, 2021; Finkelstein et al. 2015; Ishigaki et al. 2018; Oesch et al. 2018). However, some evidence suggests that the decline at the bright end of the UV LF is not as steep as predicted by the Schechter function (Hathi et al. 2012; Finkelstein et al. 2013; Bowler et al. 2014, 2020; McLeod et al. 2024). Some authors have proposed that the bright end of the UV LF at $z \sim 4–7$ ($M_{UV} < -24$) cannot be adequately explained by the Schechter function alone, suggesting the necessity

for either a double power-law (DPL) or a modified Schechter function to better account for observational findings (Ono et al. 2018).

The early release of JWST/Near Infrared Camera (NIRCam) data in July 2022 has enabled us to investigate the UV LF at redshifts beyond 10, with the furthest candidate at approximately $z \sim 16$ through photometry (Castellano et al. 2022, 2023; Finkelstein et al. 2022; Naidu et al. 2022; Adams et al. 2023; Atek et al. 2023; Bradley et al. 2023; Labbé et al. 2023; Tacchella et al. 2023a, b), along with some spectroscopic redshift determinations (Schaerer et al. 2022; Arrabal Haro et al. 2023a; Bunker et al. 2023; Curti et al. 2023; Curtis-Lake et al. 2023; Heintz et al. 2023a, b; Harikane et al. 2024). Instead of the rapid decline in UV LF described by the Schechter function at the bright end, recent data from JWST indicate a more gradual decrease, better captured by a double power law (Finkelstein et al. 2022; Naidu et al. 2022; Donnan et al. 2023, 2024).

The reasons for the deviation from a Schechter function to a double power law at high redshifts are varied and speculative. Explanations include a high star formation efficiency at high redshifts (Harikane et al. 2024), the top-heavy nature of the stellar initial mass function for Population III stars (Finkelstein et al. 2023; Bovill et al. 2024; Lapi et al. 2024), the removal of interstellar dust through radiatively driven outflows during the initial stages of galaxy formation (Ferrara, Pallottini & Dayal 2023; Fiore et al. 2023; Ziparo et al. 2023), short-term variation in dust attenuation (Mirocha & Furlanetto 2023), reduced negative or even positive feedback at early times (Silk et al. 2024), or a combination of these effects.

High-redshift galaxies observed by the JWST are expected to be among the most luminous within the field of view at their respective redshifts. Consequently, it is essential to account for any potential biases that may arise when interpreting the UV LF

<footnote n="1">¹Warwick Mathematics Institute, University of Warwick, Zeeman Building, Coventry CV4 7AL, UK</footnote>

<footnote n="2">²Institute for Fundamental Study, Naresuan University, Phitsanulok 65000, Thailand</footnote>

<footnote n="3">³Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="4">⁴William H. Miller III Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="5">⁵BIPAC, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote>∗E-mail: cameron.heather@warwick.ac.uk</footnote>
from those luminous galaxies, especially at the bright end of the UV LF. In other words, it is crucial to determine whether the UV-bright galaxies observed are representative of the extreme population from the tail end of the luminosity distribution. In this study, we will address this question using extreme-value statistics and derive a semi-analytical model of the luminosity distribution for the brightest galaxies at high redshifts and compare that with JWST observation.

A pioneering work on using extreme-value statistics to study bright galaxies was Bhavsar & Barrow (1985). They modelled the distribution of the brightest galaxies in groups and clusters at a redshift range $z < 3$ using what is now known as the generalized extreme-value (GEV) distribution (which we will discuss later in Section 3). Our work naturally follows their footsteps but extends the scope to a larger population of galaxies using a more sophisticated modelling of galaxy characteristics.

This paper will be organized as follows: in Section 2, we will describe models of the galaxy LF. Alongside the Schechter and the double power-law functions, we will introduce the star formation efficiency model, which is parametrized by the star formation efficiency $f_*$ in Section 3, we will introduce a formalism for determining the distribution of the most luminous high-redshift galaxies using extreme-value statistics and compare our theoretical model with the JWST data. In Section 4, we expand our model by introducing variations in the star formation efficiency with respect to halo mass. Our conclusions and discussion will be presented in Section 5.

We will assume a flat Lambda cold dark matter cosmology with $H_0 = 68.0 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.307$, $\Omega_\Lambda = 0.693$, $\Omega_b = 0.046$, and $\sigma_8 = 0.823$.

## 2 THE LUMINOSITY FUNCTION

The LF gives the number density of galaxies per intrinsic luminosity interval $[L, L+dL]$ at a given redshift. We will use the absolute UV magnitude, $M_{UV}$, as a measure of luminosity, $L$, via the relation

$$M_{UV} - M_{UV}^{*} = -2.5 \log_{10} \left( \frac{L}{L^{*}} \right), \qquad (1)$$

where $M_{UV}^*$ is a characteristic magnitude and $L^*$ is a characteristic luminosity, which are introduced to allow the functions to be fitted to some observed data.

We now discuss two approaches to obtain the LF of high-redshift galaxies. The first method is to fit the LF to observational data using the double power-law and Schechter functions. The second method uses a halo mass function combined with a model of the stellar formation rate (SFR).

### 2.1 Double power-law and Schechter models

The double power-law function is used to model the luminosity of galaxy groups and clusters (Holmberg 1969; Tempel et al. 2009) and of radio galaxies and quasars (Dunlop & Peacock 1990). We take the form of the double power-law LF from Harikane et al. (2024):

$$\frac{dn}{dM_{UV}} = \frac{\ln 10}{2.5} \varphi^* \times \left[ 10^{0.4(\bar{\alpha}+1)(M_{UV}-M_{UV}^*)} + 10^{0.4(\bar{\beta}+1)(M_{UV}-M_{UV}^*)} \right]^{-1}, \quad (2)$$

where the parameters are given by (Harikane et al. 2023)

$$\begin{aligned} M_{UV}^{*} &= -0.09(z-9) - 19.33, \\ \log \varphi^{*} &= -0.28(z-9) - 3.50, \\ \bar{\alpha} &= -2.10, \\ \bar{\beta} &= 0.15(z-9) - 3.27. \end{aligned}$$

Fig. 1 shows this LF (black solid lines) plotted against magnitude $M_{UV}$, along with the other LFs that we will be exploring in this work.

Next, the Schechter model (Schechter 1976) for the LF is given by

$$\frac{dn}{dL} = \varphi^* \left(\frac{L}{L^*}\right)^{\bar{\alpha}} \exp\left(-L/L^*\right), \qquad (3)$$

where $\varphi^*$ and $\bar{\alpha}$ are the model parameters. Combining this with equation (1), the Schechter LF can be expressed as

$$\frac{dn}{dM_{UV}} = \frac{\ln 10}{2.5} \varphi^* \left(-10^{0.4(M_{UV}-M_{UV}^*)}(\bar{\alpha}+1)\right) \times \exp\left[-10^{0.4(M_{UV}-M_{UV}^*)}\right], \qquad (4)$$

where the parameters are given by (Harikane et al. 2023)

$$\begin{aligned} M_{UV}^{*} &= -0.32(z - 9) - 21.24, \\ \log(\varphi^{*}) &= -0.08(z - 9) - 4.83, \\ \bar{\alpha} &= -2.35. \end{aligned}$$

We also plot the Schechter function against magnitude in Fig. 1 (black dashed lines).

### 2.2 The star formation efficiency model

In this model, we start by modelling the abundance of bright galaxies using a halo mass function coupled with some model of the SFR. A particularly simple model of the SFR [previously used in Harikane et al. (2024)] is given by

$$\text{SFR} = f_* \times f_b \times \frac{dM_h}{dt}(M_h, z), \qquad (5)$$

where $f_*$ is the star formation efficiency, $f_b \equiv \Omega_b/\Omega_m$ is the cosmic baryon fraction, and $dM_h/dt$ is the matter accretion rate for a halo of mass $M_h$. We will use the model of the accretion rate given in Behroozi & Silk (2015).

In this section, we will take $f_*$ to be constant, taking the values of 0.02, 0.05, 0.15, 0.4, and 1 (Section 4 discusses a more sophisticated model of $f_*$). We then convert the SFR into UV luminosity, $L_{UV}$, via

$$L_{UV}(\text{erg s}^{-1} \text{Hz}^{-1}) = \text{SFR}(M_{\odot} \text{yr}^{-1})/(1.15 \times 10^{-28}). \qquad (6)$$

Combining equations (5) and (6), we obtain an equation for luminosity in terms of halo mass, $M_h$.

Next, luminosity is converted to UV magnitude using values taken from Oke & Gunn (1983), such that equation (1) becomes

$$M_{UV} = -2.5 \log_{10} \left( \frac{L_{UV}}{\text{erg s}^{-1} \text{Hz}^{-1}} \right) + 51.6. \qquad (7)$$

This gives us the following relation between the UV magnitude and halo mass:

$$M_{UV}(M_h, z) = -2.5 \log_{10} (\dot{M}_h) - 2.5 \log_{10} \left( \frac{f_* \times f_b}{1.15 \times 10^{-28}} \right) + 51.6. \qquad (8)$$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Comparison of the LFs $dn/dM_{UV}$ at $z = 9$ (upper left), $z = 10$ (upper right), $z = 12$ (lower left), and $z = 16$ (lower right). The black solid and dashed curves are the double power-law and Schechter LFs, respectively, calculated in Section 2.1. The colourful curves correspond to the $f_*$ model of the LF in Section 2.2 (see equation 13), taking $f_* = 0.02, 0.05, 0.15, 0.4$, and 1, with amplitudes matching those of Harikane et al. (2024). The red data points correspond to spectroscopic data from Harikane et al. (2024), and the black points correspond to photometric data from Donnan et al. (2024).</figcaption></figure>

Taking the derivative with respect to $M_h$, we obtain

$$ \frac{dM_{UV}}{dM_h} = -\frac{2.5}{\ln(10)} \frac{\dot{M}_h}{\dot{M}_h^2}. \qquad (9) $$

We now introduce the halo mass function, $dn/dM_h$, defined as the number of dark matter haloes per unit mass interval $[M_h, M_h + dM_h]$:

$$ \frac{dn}{dM_h} = f(\sigma) \frac{\bar{\rho}_m}{M_h} \frac{d \ln \sigma^{-1}}{dM_h}. \qquad (10) $$

Following Behroozi & Silk (2015), we will use the Tinker mass function (Tinker et al. 2008) in which

$$ f(\sigma) = A \left[ \left( \frac{\sigma}{b} \right)^{-a} + 1 \right] e^{-c/\sigma^2}, \qquad (11) $$

where $\sigma$ is the variance of the smoothed linear density field, and the parameters $A$, $b$, and $c$ depend on redshift. The LF can then be expressed via the chain rule as

$$ \frac{dn}{dM_{UV}} = \frac{dn}{dM_h} \left| \frac{dM_h}{dM_{UV}} \right|, \qquad (12) $$

$$ = \frac{dn}{dM_h} \times 0.4 \ln(10) \left| \frac{\dot{M}_h^2}{\dot{M}_h} \right|. \qquad (13) $$

In equation (13), there appears to be no direct dependence on $f_*$. However, the $f_*$ dependence does arise through the conversion formulae in equations (6) and (8).

We plot the various LFs calculated in this section at four redshifts ($z = 9, 10, 12$, and 16) in Fig. 1, where $f_*$ is varied between 0.02

and 1 [the curves for the $f_*$ model are normalized to match those of Harikane et al. (2024)]. We see that the Schechter LF is comparable to the $f_*$ model with $f_*$ roughly between 0.02 and 0.05. On the other hand, the double power-law function is flatter and predicts a higher number of brighter galaxies than the other models.

### 3 EXTREME-VALUE MODELLING OF THE BRIGHTEST GALAXIES

Our goal is to obtain a statistical description of the most extreme luminosities of high-redshift galaxies. We will do this using extreme-value statistics (EVS), which have previously been used in astrophysics to find the brightest galaxies in clusters by Bhavsar & Barrow (1985), quantify the abundances of the most massive Population III stars (Chantavat, Chongchitnan & Silk 2023), predict the most massive galaxy clusters (Davis et al. 2011; Chongchitnan & Silk 2012; Waizmann, Ettori & Moscardini 2012), and estimate the abundances of extreme primordial black holes and extreme-spin primordial black holes (Chongchitnan & Silk 2021; Chongchitnan, Chantavat & Zunder 2021).

In particular, we will be working with the GEV approach, also known as the block maxima method (Gumbel 1958). In this method, we divide the galaxy population in a given redshift bin into $N$ distinct blocks, from which we identify the brightest galaxy in each block using a probability density function (pdf) outlined below. Analogous to the central limit theorem, these brightest galaxies will have a large-
N limit distribution corresponding to one of the GEV distributions (equation 18).

We summarize the GEV pipeline here [more details can be found in White (1979) and Davis et al. (2011)]. First, we calculate the number density $n(<M_{UV})$ of galaxies with UV magnitude less than $M_{UV}$ as

$$n(<M_{UV}, z) = \int_{-\infty}^{M_{UV}} \frac{dn}{dM_{UV}} dM'_{UV}. \qquad (14)$$

Note that, by the definition of $M_{UV}$, equation (14) counts the brightest galaxies. The redshift dependence can then be integrated out to obtain the number count of the brightest galaxies in the redshift interval $[z_0, z_1]$:

$$n(<M_{UV}) = \int_{z_0}^{z_1} n(<M_{UV}, z) dz. \qquad (15)$$

Next, to predict the brightest magnitude that galaxies could have at high redshifts, we consider the probability $P_0$ that no galaxies in a given volume exceed the extreme brightness $M_{UV}$. In our calculation for the volume $V$, we take $f_{sky} = 1$ in order to obtain generic predictions on the extreme luminosities. The probability distribution $P_0(M_{UV})$ can be modelled as a Poisson distribution with the following cumulative distribution function (cdf):

$$P_0(M_{UV}) = \exp[-n(<M_{UV})V]. \qquad (16)$$

Differentiating this with respect to $M_{UV}$, we obtain the pdf:

$$\frac{dP_0}{dM_{UV}} = -V \frac{dn(<M_{UV})}{dM_{UV}} \exp[-n(<M_{UV})V]. \qquad (17)$$

The Fisher-Tippett-Gnedenko theorem implies that in the large-N limit, the cdf (equation 16) approaches the GEV distribution given by

$$G(M_{UV}) = \begin{cases} \exp\left[-\frac{(1 + \gamma y)^{-1/\gamma}}{\gamma} \right] & (\gamma \neq 0), \\ \exp\left[-\frac{e^{-y}}{\gamma} \right] & (\gamma = 0), \end{cases} \qquad (18)$$

where $y := (M_{UV} - \alpha)/\beta$, with $\alpha$ describing the location of the peak and $\beta$ describing the scale of the pdf. The sign of the parameter $\gamma$ determines the GEV type, with $\gamma = 0$, $\gamma > 0$, and $\gamma < 0$, corresponding to the Gumbel, Fréchet, and Weibull distributions, respectively.

We can express the GEV parameters $\alpha$, $\beta$, and $\gamma$ in terms of astrophysical parameters by Taylor-expanding the Poisson distribution $P_0(M_{UV})$ and the GEV distribution $G(M_{UV})$ around the peak of the pdf at $M_{peak}$ to cubic order. By equating the coefficients, we find

$$\gamma = n(<M_{peak})V - 1, \\ \alpha = M_{peak} - \frac{\beta}{\gamma} [(1 + \gamma)^{-\gamma} - 1], \\ \beta = \frac{(1 + \gamma)^{1+\gamma}}{\frac{dn}{dM_{UV}}|_{M_{peak}} V}. \qquad (19)$$

While these equations provide a good approximation to the values of $\gamma$, $\alpha$, and $\beta$, using a numerical approximation produces a closer fit to the Poisson pdf (equation 17), which we use to plot the GEV distribution in Fig. 2. This method assumes that the pdfs are well described by the Gumbel distribution, taking $\gamma \approx 0$.

## 3.1 Extreme-value pdf of the brightest galaxies

We plot the pdf for the distribution of the brightest galaxies in Fig. 2 for four redshift bins: $z \in [9, 10], [10, 11], [12, 13],$ and $[16, 17]$. The

solid lines are pdfs calculated using equation (17), while the dashed lines are the numerical fits assuming that the pdfs are Gumbel type ($\gamma = 0$ in equation 18). The parameters $\alpha$ and $\beta$ for the Gumbel fit are given in Table A1. Note that in each bin, $\alpha$ is very close to the magnitude of $M_{UV}$ where the peak of the pdf occurs.

From Fig. 2, we see that at higher redshifts, the pdfs shift to higher values of $M_{UV}$ (meaning that the extreme galaxies are less bright). We also observe that increasing $f_*$ produces brighter extreme galaxies. In particular, increasing $f_*$ from 0.02 to 0.15 results in a pdf peak shift of $\Delta M_{UV} \approx -2$.

## 3.2 Comparison with JWST data

We now study the redshift dependence of the extreme-value predictions more closely. In Figs 3 and 4, we plot the peak of the Gumbel pdf for the brightest galaxies over the redshift range $z = 9-16$ (dashed lines). We also indicate the 95th and 99th percentiles of the pdfs (shaded dark and light colours). Fig. 3 shows the extreme-value predictions for $M_{UV}$ for the double power-law and Schechter models. Fig. 4 shows the same predictions for the $f_*$ model with $f_* = 0.02, 0.05$, and $0.15$.

On these plots, we have included data from the JWST Early Release Science program GLASS and from Cosmic Evolution Early Release Science (CEERS) observations. In particular, we have selected the brightest galaxies from a collection of recent studies on these galaxies [Donnan et al. (2023), Bowler et al. (2020), Harikane et al. (2022b), Harikane et al. (2023), Arrabal Haro et al. (2023a), Arrabal Haro et al. (2023b), Tacchella et al. (2023b), and Harikane et al. (2024)]. Most of these studies obtained the photometric redshifts for the galaxies using the EAZY or PROSPECTOR pipeline. While the redshifts have not been spectroscopically confirmed, the papers have published redshift uncertainties, which we have displayed as error bars. We also include spectroscopically confirmed galaxies from Arrabal Haro et al. (2023a, b), Tacchella et al. (2023b), Harikane et al. (2024), Castellano et al. (2024), and Carniani et al. (2024).

### 3.2.1 Double power-law and Schechter models

In Fig. 3, we note that neither the double power-law nor Schechter model produces a satisfactory prediction for the peak magnitudes of the brightest galaxies. We discuss these separately below.

The double power-law model (left panel) overestimates the brightness of extreme galaxies. For example, for $z > 14.5$, the model predicts extreme brightness of $M_{UV} \lesssim -30$. More so, it gives an unphysical prediction that the most luminous galaxies are brighter at higher redshifts. We note that the unphysical trend can be corrected by an ad hoc change in the sign $d\tilde{\beta}/dz$ in equation (2). We will revisit this point in Section 5.

Next, using the Schechter model (right panel), we find that the width of the pdf is much smaller. As a result, the coloured band barely contains any of the extreme data points, and, more worryingly, we also see observed galaxies that are brighter than the predicted extreme brightness.

### 3.2.2 The $f_*$ model

In Fig. 4, three of the star formation efficiency models have been plotted, $f_* = 0.02, 0.05$, and $0.15$. The model with $f_* = 0.02$ is comparable to the Schechter model discussed above. The model with $f_* = 0.05$ (blue band) appears to be consistent with the data points, with the implication that these galaxies are expected to
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The pdfs for the most luminous galaxies calculated in Section 2.2 for four redshift bins: $z = [9, 10]$ (upper left), $[10, 11]$ (upper right), $[12, 13]$ (lower left), and $[16, 17]$ (lower right). The solid coloured lines are the pdfs derived from the Poisson distribution (equation 17), assuming three values for the star formation efficiency: $f_* = 0.02$ (right curve), 0.05 (middle curve), and 0.15 (left curve). The dashed lines are the Gumbel distribution (equation 18) with parameters $\alpha$ and $\beta$ given in Table A1.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The peak magnitudes $M_{\text{peak}}$ (dashed lines) of the extreme-value (Gumbel) pdf for $M_{\text{UV}}$ plotted as a function of redshift, assuming the double power-law LF from equation (2) (left) and the Schechter LF from equation (4) (right). The darker/lighter regions correspond to the 95th and 99th percentiles of the pdf. The data points are based on the most extreme samples from previous studies [Donnan et al. (2023), Bowler et al. (2020), Harikane et al. (2022b), Harikane et al. (2023), Arrabal Haro et al. (2023a), Arrabal Haro et al. (2023b), Tacchella et al. (2023b), Harikane et al. (2024), Castellano et al. (2024), and Carniani et al. (2024)].</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Extreme-value pdf profile (same as Fig. 3) plotted as a function of redshift, for the $f_*$ model of the LF in Section 2.2 (equation 13). The values for the star formation efficiency shown are $f_* = 0.02$ (upper band), 0.05 (middle band), and 0.15 (lower band).</figcaption></figure>

be the brightest we would observe. The model with $f_* = 0.15$ is also consistent with the data, in the sense that all the data points correspond to galaxies that are less bright than the extreme-value prediction.

The extreme-value modelling allows us to extrapolate to higher redshifts. For example, assuming $f_* = 0.15$, at $z \sim 16$ we expect to see galaxies no brighter than $M_{UV} \sim -23.5_{0.4}^{0.8}$. Future observations can be added to such a plot, giving us additional constraints on the star formation efficiency $f_*$.

## 4 EXTENDING THE $f_*$ MODEL

A natural extension to the $f_*$ model is to consider the case when the star formation efficiency is not a constant. Our goal is to obtain the extreme-value pdf for the brightest galaxies in this model.

One such model is that of Tacchella et al. (2018) in which $f_*$ has a dependence on the halo mass $M_h$ in the form of the double power-law function:

$$f_* = \frac{\text{SFR}}{\dot{M}_h} = 2\epsilon_0 \left[ \left( \frac{M_h}{M_c} \right)^{-\mu} + \left( \frac{M_h}{M_c} \right)^{\nu} \right]^{-1}, \quad (20)$$

where $\epsilon_0$ is a normalization constant, $M_c$ is a characteristic mass where the efficiency is equal to $\epsilon_0$, and $\mu$ and $\nu$ are slopes that determine the behaviour at high and low masses. This equation has been calibrated at $z = 4$, with parameter values

$$\text{Tacchella model: } (\epsilon_0, M_c, \mu, \nu) = (0.22, 6.3 \times 10^{10} M_{\odot}, 0.89, 0.4). \quad (21)$$

We also see the same form in Harikane et al. (2022a), with parameter values

$$\text{Harikane model: } (\epsilon_0, M_c, \mu, \nu) = (3.2 \times 10^{-2}, 10^{11.5} M_{\odot}, 1.2, 0.5). \quad (22)$$

These two models of $f_*$ are plotted in Fig. 5.

We can follow the same method as Section 2.2 to find the LF for this mass-dependent $f_*$, by substituting equation (20) into the working. In particular, our equation for $M_{UV}(M_h, z)$ (equation 8) is

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of two stellar efficiency functions $f_*$ as a function of halo mass $M_h$. The double power-law form is given in equation (20), with parameter values in equations (21) and (22) for the Tacchella (solid) and Harikane (dashed) models, respectively.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Extreme-value pdf profile (same as Fig. 3) plotted as a function of redshift, for the halo mass-dependent $f_*$ model (equation 20). The upper band corresponds to the Harikane model (equation 22), while the lower band is the Tacchella model (equation 21).</figcaption></figure>

now

$$M_{UV}(M_h, z) = -2.5 \log_{10} (\dot{M}_h) + 2.5 \log_{10} (M_*^{-\mu} + M_*^{\nu}) \\ -2.5 \log_{10} \left( \frac{2\epsilon_0}{K_{UV}} \right) + 51.6, \quad (23)$$

where $M_* = M_h/M_c$. We can then take the derivative with respect to $M_h$ and obtain the LF:

$$\frac{dn}{dM_{UV}} = \frac{dn}{dM_h} \times 0.4 \ln(10) \left| \frac{\dot{M}_h}{\dot{M}_h^2} + \frac{\mu M_*^{-\mu} - \nu M_*^{\nu}}{M_h (M_*^{-\mu} + M_*^{\nu})} \right|^{-1}. \quad (24)$$

From this, we introduce the extreme-value methodology from Section 3, where we find the peak magnitude for each redshift bin, which we then model using a GEV distribution. We plot the peak of Gumbel pdf, with the 95th and 99th percentile bands, varied over our redshift range, in Fig. 6.

The main difference between Figs 4 and 6 is the upward curvature of the bands at high redshifts. We also see that the Tacchella model gives a similar extreme-value prediction to the $f_* = 0.15$ model. At $z \sim 16$, the Tacchella curve predicts an extreme brightness of
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Extreme-value pdf profile (same as Fig. 3) plotted as a function of redshift, for the double power-law LF with parameters given in Donnan et al. (2024).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** A comparison of various EVS models in this work that are consistent with *JWST* data. Each solid line shows the peak of the Gumbel distribution (percentile bands not shown). The models for the LF $dn/dM_{UV}$ are constant $f_*$ (equation 13) with $f_* = 0.05$ (dash-dotted) and $f_* = 0.15$ (dashed), variable $f_*$ (equation 24) (solid), and the double power-law model (equation 2) assuming the parameters of Donnan et al. (2024) (dotted).</figcaption></figure>

$M_{UV} = -22.5_{-1.5}^{+0.5}$, which is less bright than what the $f_* = 0.15$ model predicts. The confidence interval in the Tacchella model is also wider at the high-$z$ end of the graph.

The Harikane model, on the other hand, predicts a more conservative set of extreme brightnesses, with some observed bright galaxies barely grazing the 99th percentile band. Around $z \sim 10$, the extreme-value prediction in the Harikane model is comparable to that with constant $f_*$ model with $f_* = 0.05$ of a few per cent. However, the band curves up rapidly towards higher redshifts. Finally, like the Tacchella model, the band is also much wider than the constant $f_*$ model.

We conclude that both the Tacchella and the Harikane models are consistent with *JWST* observations. However, considering the conservative volume coverage of *JWST*, the Harikane model will require modification to remain consistent with the observations of the brightest galaxies at $z \gtrsim 9$.

## 5 CONCLUSION AND DISCUSSION

We have applied extreme-value statistics to predict how bright the most luminous galaxies at high redshifts could be. In particular, we derived the probability distribution (in the UV magnitude, $M_{UV}$) of the brightest galaxies at redshifts $z \sim 9-16$. We found that such a distribution is well approximated by the Gumbel distribution. We studied a number of models of the galaxy LF and derived the corresponding extreme-value distributions as a function of redshift. Our main results are shown in Figs 4 and 6, in which our theoretical results are compared with data from the *JWST*.

Further summary and discussion points are given below:

(i) *The GEV formalism* (equation 18) was used to calculate pdf of the brightest galaxies. We find that the GEV parameter $\gamma$ (see Table A1) is so small that the extreme-value pdf for $M_{UV}$ is well approximated by the Gumbel distribution, as shown in Fig. 2. The pdf profiles (along with the 95th and 99th percentiles) are shown in Figs 3, 4, and 6, assuming different galaxy LFs as shown in Fig. 1.

(ii) *The Schechter model* of the LF has previously been used to model the luminosity of galaxies at $z \sim 5$. The parameters used in this work (equation 3) were calibrated using observations at $z \sim 9$. We found that the resulting extreme-value pdf (right panel in Fig. 3) appears to be in tension with the brightest *JWST* galaxies at $z \gtrsim 10$.

(iii) *The double power-law model* (equation 2), on the other hand, resulted in an unphysical extreme-value prediction (left panel in Fig. 3), where extreme galaxies appear brighter at higher redshifts. This might be due to the limited accuracy of the parametrization (equation 2). Indeed, we found an ad hoc correction by switching the sign of $d\tilde{\beta}/dz$ to negative.

Donnan et al. (2024) gave new constraints on the redshift evolution of the double power-law LF. We used their data in place of the parametrization (equation 2) and obtained the extreme-value pdf profile shown in Fig. 7, which appears to be consistent with data. We note that in their work, $d\tilde{\beta}/dz$ is indeed negative. Since the parameter $\tilde{\beta}$ describes the bright-end slope of the curve, we conclude that a physical double power-law parametrization must satisfy the condition $d\tilde{\beta}/dz < 0$.

(iv) *The $f_*$ model of the LF* [as outlined in Harikane et al. (2024)] was considered as an alternative to the Schechter and DPL models. This model requires the SFR $f_*$, which we initially assumed to be constant, and a halo mass function (taken to be the Tinker mass function). Fig. 4 shows that for constant formation rate, $f_* = 0.02$ produces a similar result to that of the Schechter model, while $f_*$ between 0.05 and 0.15 produces an extreme-value pdf that is consistent with *JWST* data.

The model predicts that at $z \sim 16$, the brightest galaxies will have a UV magnitude of around $-23$. Such a prediction depends on the choice of the halo mass function, which may be rapidly evolving at high redshifts (Donnan et al. 2024). The Tinker mass function is known to predict the halo number density accurately at low redshifts ($z \sim 1-2$) and a recalibration to high-$z$ simulations is needed to confirm whether the parameter values remain accurate for $z \sim 9-16$.

(v) *Extending the $f_*$ model* to include a dependence on the halo mass was studied in Section 4. In particular, we considered the models of Tacchella et al. (2023b) (equation 21) and Harikane et al. (2022a) (equation 22). The extreme-value predictions are shown in Fig. 6, which shows an upward curvature of the bands at high redshifts. The Tacchella model gives a similar result to the constant $f_* = 0.15$ model, while the Harikane model is comparable to $f_* = 0.05$ of a few per cent, with the predictions diverging at $z \sim 16$. As more high-redshift galaxies are observed, the extreme-value formalism
will allow us to distinguish whether $f_*$ is constant or halo mass-dependent.

Fig. 8 summarizes the main results of this work: each solid line shows the peak of the extreme-value pdf as a function of redshift, assuming the double power-law, constant $f_*$, and variable $f_*$ models that are consistent with JWST observation of the brightest high-redshift galaxies.

(vi) *Uncertainties in the data.* We have taken into account neither the effects of cosmic variance (in the sense of the statistical uncertainty in observation of rare objects at high redshifts) nor Eddington bias (Trenti & Stiavelli 2008; Teerikorpi 2015). Both of these could reduce the deduced brightness of extreme galaxies. We leave the calculation of statistical bias to a future work.

## ACKNOWLEDGEMENTS

For the purpose of open access, the authors applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission. Cameron Heather was supported by the Warwick Mathematics Institute Centre for Doctoral Training, and gratefully acknowledges funding from the University of Warwick and the UK Engineering and Physical Sciences Research Council (grant number: EP/W524645/1). The authors would like to thank Suraphong Yuma for his useful suggestions.

## DATA AVAILABILITY

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

## REFERENCES

* Adams N. J. et al., 2023, MNRAS, 518, 4755

* Arrabal Haro P. et al., 2023a, Nature, 622, 707

* Arrabal Haro P. et al., 2023b, ApJ, 951, L22

* Atek H. et al., 2023, MNRAS, 519, 1201

* Behroozi P. S., Silk J., 2015, ApJ, 799, 32

* Bhavsar S. P., Barrow J. D., 1985, MNRAS, 213, 857

* Bouwens R. J. et al., 2015, ApJ, 803, 34

* Bouwens R. J. et al., 2021, AJ, 162, 47

* Bovill M. S., Stiavelli M., Wiggins A. I., Ricotti M., Trenti M., 2024, ApJ, 962, 49

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

* Bowler R. A. A., Jarvis M. J., Dunlop J. S., McLure R. J., McLeod D. J., Adams N. J., Milvang-Jensen B., McCracken H. J., 2020, MNRAS, 493, 2059

* Bradley L. D. et al., 2023, ApJ, 955, 13

* Bromm V., Yoshida N., 2011, ARA&A, 49, 373

* Bunker A. J. et al., 2023, A&A, 677, A88

* Carniani S. et al., 2024, preprint (arXiv:2405.18485)

* Castellano M. et al., 2022, ApJ, 938, L15

* Castellano M. et al., 2023, ApJ, 948, L14

* Castellano M. et al., 2024, ApJ, 972, 143

* Chantavat T., Chongchitnan S., Silk J., 2023, MNRAS, 522, 3256

* Chongchitnan S., Silk J., 2012, Phys. Rev. D, 85, 063508

* Chongchitnan S., Silk J., 2021, Phys. Rev. D, 104, 083018

* Chongchitnan S., Chantavat T., Zunder J., 2021, Astron. Nachr., 342, 648

* Curti M. et al., 2023, MNRAS, 518, 425

* Curtis-Lake E. et al., 2023, Nat. Astron., 7, 622

* Davis O., Devriendt J., Colombi S., Silk J., Pichon C., 2011, MNRAS, 413, 2087

* DeBoer D. R. et al., 2017, PASP, 129, 045001

* Donnan C. T. et al., 2023, MNRAS, 518, 6011

* Donnan C. T. et al., 2024, MNRAS, 533, 3222

* Dunlop J. S., Peacock J. A., 1990, MNRAS, 247, 19

* Ellis R. S. et al., 2013, ApJ, 763, L7

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

* Finkelstein S. L. et al., 2013, Nature, 502, 524

* Finkelstein S. L. et al., 2015, ApJ, 810, 71

* Finkelstein S. L. et al., 2022, ApJ, 940, L55

* Finkelstein S. L. et al., 2023, ApJ, 946, L13

* Fiore F., Ferrara A., Bischetti M., Feruglio C., Travascio A., 2023, ApJ, 943, L27

* Fomalont E. B. et al., 2015, ApJ, 808, L1

* Gardner J. P. et al., 2006, Space Sci. Rev., 123, 485

* Gumbel E. J., 1958, Statistics of Extremes. Columbia University Press, New York, US

* Harikane Y. et al., 2022a, ApJS, 259, 20

* Harikane Y. et al., 2022b, ApJ, 929, 1

* Harikane Y. et al., 2023, ApJS, 265, 5

* Harikane Y., Nakajima K., Ouchi M., Umeda H., Isobe Y., Ono Y., Xu Y., Zhang Y., 2024, ApJ, 960, 56

* Hathi N. P., Mobasher B., Capak P., Wang W.-H., Ferguson H. C., 2012, ApJ, 757, 43

* Heintz K. E. et al., 2023a, Nat. Astron., 7, 1517

* Heintz K. E. et al., 2023b, ApJ, 944, L30

* Holmberg E., 1969, Ark. Astron., 5, 305

* Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, ApJ, 854, 73

* Labbé I. et al., 2023, Nature, 616, 266

* Lapi A. et al., 2024, Universe, 10, 141

* Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

* McLeod D. J. et al., 2024, MNRAS, 527, 5004

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

* Marchetti L., Serjeant S., Vaccari M., 2017, MNRAS, 470, 5007

* Mirocha J., Furlanetto S. R., 2023, MNRAS, 519, 843

* Naidu R. P. et al., 2022, ApJ, 940, L14

* Oesch P. A., Bouwens R. J., Illingworth G. D., Labbé I., Stefanon M., 2018, ApJ, 855, 105

* Oke J. B., Gunn J. E., 1983, ApJ, 266, 713

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

* Ouchi M., Ono Y., Shibuya T., 2020, ARA&A, 58, 617

* Robertson B. E., 2022, ARA&A, 60, 121

* Schaerer D., Marques-Chaves R., Barrufet L., Oesch P., Izotov Y. I., Naidu R., Guseva N. G., Brammer G., 2022, A&A, 665, L4

* Schechter P., 1976, ApJ, 203, 297

* Silk J., Begelman M. C., Norman C., Nusser A., Wyse R. F. G., 2024, ApJ, 961, L39

* Spergel D. et al., 2015, preprint (arXiv:1503.03757)

* Stark D. P., 2016, ARA&A, 54, 761

* Tacchella S., Bose S., Conroy C., Eisenstein D. J., Johnson B. D., 2018, ApJ, 868, 92

* Tacchella S. et al., 2023a, MNRAS, 522, 6236

* Tacchella S. et al., 2023b, ApJ, 952, 74

* Teerikorpi P., 2015, A&A, 576, A75

* Tempel E., Einasto J., Einasto M., Saar E., Tago E., 2009, A&A, 495, 37

* Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

* Trenti M., Stiavelli M., 2008, ApJ, 676, 767

* Waizmann J. C., Ettori S., Moscardini L., 2012, MNRAS, 420, 1754

* White S. D. M., 1979, MNRAS, 186, 145

* Ziparo F., Ferrara A., Sommovigo L., Kohandel M., 2023, MNRAS, 520, 2445

## APPENDIX A

Table A1 shows the GEV parameters ($\alpha$ and $\beta$) and peak magnitudes for each LF at various redshifts (see equation 18). The parameters have been calculated numerically using least-square fitting, assuming $\gamma = 0$ (i.e. assume that the GEV distribution is the Gumbel type). Alternatively, one can use equation (19) to approximate these values.
<table><caption>**Table A1.** The peak magnitude ($M_{UV}$) and the Gumbel parameters ($\alpha$ and $\beta$) for the LFs considered in this paper. The notation $z \sim 9, 10, 12$, and 16 refers to the redshift bins [9, 10], [10, 11], [12, 13], and [16, 17], respectively.</caption><thead><tr><th>Luminosity function</th><th>Redshift</th><th>Peak M<sub>UV</sub></th><th>α</th><th>β</th></tr></thead><tbody><tr><td rowspan="4">DPL<br>(equation 2)</td><td>z ~ 9</td><td>-27.88</td><td>27.86</td><td>0.4960</td></tr><tr><td>z ~ 10</td><td>-28.27</td><td>28.21</td><td>0.5211</td></tr><tr><td>z ~ 12</td><td>-29.06</td><td>29.09</td><td>0.6093</td></tr><tr><td>z ~ 16</td><td>-31.56</td><td>31.58</td><td>0.9298</td></tr><tr><td rowspan="4">DPL<br>(Donnan et al. 2024)</td><td>z ~ 9</td><td>-26.01</td><td>26.02</td><td>0.3679</td></tr><tr><td>z ~ 10</td><td>-25.69</td><td>25.69</td><td>0.3578</td></tr><tr><td>z ~ 12</td><td>-25.13</td><td>25.14</td><td>0.3392</td></tr><tr><td>z ~ 16</td><td>-24.30</td><td>24.31</td><td>0.3075</td></tr><tr><td rowspan="4">Schechter<br>(equation 4)</td><td>z ~ 9</td><td>-23.59</td><td>23.59</td><td>0.09729</td></tr><tr><td>z ~ 10</td><td>-23.25</td><td>23.24</td><td>0.09917</td></tr><tr><td>z ~ 12</td><td>-22.56</td><td>22.55</td><td>0.1030</td></tr><tr><td>z ~ 16</td><td>-21.17</td><td>21.12</td><td>0.1112</td></tr><tr><td rowspan="4">f<sub>*</sub> = 0.02<br>(equation 13)</td><td>z ~ 9</td><td>-23.84</td><td>23.82</td><td>0.1651</td></tr><tr><td>z ~ 10</td><td>-23.43</td><td>23.38</td><td>0.1711</td></tr><tr><td>z ~ 12</td><td>-22.53</td><td>22.52</td><td>0.1751</td></tr><tr><td>z ~ 16</td><td>-20.91</td><td>20.92</td><td>0.1907</td></tr><tr><td rowspan="4">f<sub>*</sub> = 0.05</td><td>z ~ 9</td><td>-24.85</td><td>24.81</td><td>0.1647</td></tr><tr><td>z ~ 10</td><td>-24.34</td><td>24.37</td><td>0.1706</td></tr><tr><td>z ~ 12</td><td>-23.54</td><td>23.52</td><td>0.1747</td></tr><tr><td>z ~ 16</td><td>-21.92</td><td>21.92</td><td>0.1902</td></tr><tr><td rowspan="4">f<sub>*</sub> = 0.15</td><td>z ~ 9</td><td>-25.96</td><td>26.00</td><td>0.1644</td></tr><tr><td>z ~ 10</td><td>-25.56</td><td>25.56</td><td>0.1703</td></tr><tr><td>z ~ 12</td><td>-24.75</td><td>24.71</td><td>0.1745</td></tr><tr><td>z ~ 16</td><td>-23.13</td><td>23.11</td><td>0.1895</td></tr><tr><td rowspan="4">f<sub>*</sub>(M<sub>h</sub>) Tacchella<br>(equation 21)</td><td>z ~ 9</td><td>-25.45</td><td>25.42</td><td>0.1277</td></tr><tr><td>z ~ 10</td><td>-25.15</td><td>25.18</td><td>0.1379</td></tr><tr><td>z ~ 12</td><td>-24.65</td><td>24.59</td><td>0.1835</td></tr><tr><td>z ~ 16</td><td>-22.42</td><td>22.41</td><td>0.3374</td></tr><tr><td rowspan="4">f<sub>*</sub>(M<sub>h</sub>) Harikane<br>(equation 22)</td><td>z ~ 9</td><td>-23.84</td><td>23.77</td><td>0.1451</td></tr><tr><td>z ~ 10</td><td>-23.43</td><td>23.40</td><td>0.1981</td></tr><tr><td>z ~ 12</td><td>-22.22</td><td>22.11</td><td>0.3296</td></tr><tr><td>z ~ 16</td><td>-18.26</td><td>18.20</td><td>0.4495</td></tr></tbody></table>

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