Extreme-value modelling of the brightest galaxies at 𝑧 ≳ 9

Data from the James Webb Space Telescope have revealed an intriguing population of bright galaxies at high redshifts. In this work, we use extreme-value statistics to calculate the distribution (in UV magnitude) of the brightest galaxies in the redshift range 9 ≲ 𝑧 ≲ 16. We combine the Generalised Extreme Value (GEV) 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 parametrised by the stellar formation efficiency 𝑓 ∗ . We find that the JWST data is broadly consistent with 𝑓 ∗ of 5% − 10%, and that the brightest galaxy at 𝑧 ∼ 16 will have 𝑀 UV ≈ − 23 . 5 0 . 8 0 . 4 . If 𝑓 ∗ is dependent on halo mass, we predict 𝑀 UV ≈ − 22 . 5 0 . 5 1 . 5 for such an object. We show that extreme-value statistics not only predicts the magnitude of the brightest galaxies at high redshifts, but may also be able to distinguish between models of star formation in high-redshift galaxies.


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 et al. 2020;Robertson 2022).In pursuit of this goal, a new generation of multi-wavelength surveys has shed a new light on galaxies at unprecedentedly high redshifts (Gardner et al. 2006;Laureĳs et al. 2011;Fomalont et al. 2015;Spergel et al. 2015;DeBoer et al. 2017;Marchetti et al. 2017).
One of the essential ingredients for studying galaxy evolution is the rest-frame ultraviolet (UV) luminosity function (LF), which characterises the density of galaxies within a given volume based on their luminosity.The UV LF is well constrained within the redshift range of 2 ≲  ≲ 10 based on observations from the Hubble Space Telescope (HST) and Spitzer, predating the launch of the James Webb Space Telescope (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. 2015Bouwens et al. , 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. 2014Bowler et al. , 2020;;McLeod et al. 2024).Some authors have proposed that the bright end of the UV LF at  ∼ 4 − 7 ( UV < −24) cannot be adequately explained by the Schechter function alone, suggesting the necessity for either a ★ E-mail:cameron.heather@warwick.ac.uk double power law or a modified Schechter function to better account for observational findings (Ono et al. 2018).
The early release of JWST/NIRCam data in July 2022 has enabled us to investigate the UV LF at redshifts beyond 10, with the furthest candidate at approximately  ∼ 16 through photometry (Castellano et al. 2022(Castellano et al. , 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. 2023b,a;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 indicates a more gradual decrease, better captured by a double power law (Finkelstein et al. 2022;Naidu et al. 2022;Donnan et al. 2023Donnan et al. , 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 et al. 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 amongst 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 luminosity function from those luminous galaxies especially at the bright end of the UV LF.In other words, it is crucial to determine whether the UVbright 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  < 3 using what is now known as the Generalised Extreme-Value 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 organised as follows: in section 2, we will describe models of the galaxy luminosity function.Alongside the Schechter and the double power law functions, we will introduce the star formation efficiency model, which is parametrised by the star formation efficiency  * .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.

THE LUMINOSITY FUNCTION
The luminosity function gives the number density of galaxies per intrinsic luminosity interval [,  + dL] at a given redshift.We will use the absolute UV magnitude,  UV , as a measure of luminosity, , via the relation where  * UV is a characteristic magnitude and  * 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 luminosity function of high-redshift galaxies.The first method is to fit the luminosity function 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.
Fig. 1 shows this luminosity function (black solid lines, labelled DPL) plotted against magnitude  UV , along with the other luminosity functions that we will be exploring in this work.
Next, the Schechter model (Schechter 1976) for the luminosity function is given by where  * and α are the model parameters.Combining this with equation (1), the Schechter luminosity function can be expressed as where the parameters are given by (Harikane et al. 2023): We also plot the Schechter function against magnitude in Fig. 1 (black dashed lines).

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 stellar formation rate (SFR).A particularly simple model of the SFR (previously used in Harikane et al. (2024)) is given by: where  * is the star formation efficiency,  b ≡ Ω b /Ω m is the cosmic baryon fraction, and d h / d is the matter accretion rate for a halo of mass  h .We will use the model of the accretion rate given in Behroozi & Silk (2015).
In this section, we will take  * to be constant, taking the values 0.02, 0.05, 0.15, 0.4, 1 (section 4 discusses a more sophisticated model of  * ).We then convert the SFR into UV luminosity,  UV , via Combining equations ( 5) and ( 6), we obtain an equation for luminosity in terms of halo mass,  h .Next, luminosity is converted to UV magnitude using values taken from Oke & Gunn (1983), such that equation (1) becomes This gives us the following relation between the UV magnitude and halo mass.
Taking the derivative with respect to  h , we obtain We now introduce the halo mass function, d/ d h , defined as the number of dark matter haloes per unit mass interval Following Behroozi & Silk (2015), we will use the Tinker mass function (Tinker et al. 2008) in which where  is the variance of the smoothed linear density field, and the parameters ,  and  depend on redshift.The luminosity function can then be expressed via the chain rule as: In equation ( 13), there appears to be no direct dependence on  * .However, the  * dependence does arise through the conversion formulae in equations ( 6) and ( 8).
We plot the various luminosity functions calculated in this section at four redshifts ( = 9, 10, 12, 16) in Fig. 1, where  * is varied between 0.02 and 1 (the curves for the  * model are normalised to match those of Harikane et al. (2024)).We see that the Schechter luminosity function is comparable to the  * model with  * 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.

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 extremevalue statistics, which has previously been used in astrophysics to find the brightest galaxies in clusters by Bhavsar & Barrow (1985), quantify the abundances of the most massive Pop III stars (Chantavat et al. (2023)), predict the most massive galaxy clusters (Davis et al. 2011;Waizmann et al. 2012;Chongchitnan & Silk 2012) and estimate the abundances of extreme primordial black holes and extreme-spin primordial black holes (Chongchitnan et al. 2021;Chongchitnan & Silk 2021).
In particular, we will be working with the generalised extremevalue (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  distinct blocks, from which we identify the brightest galaxy in each block using a probability density function outlined below.Analogous to the Central Limit Theorem, these brightest galaxies will have a large- limit distribution corresponding to one of the GEV distributions (equation ( 18)).
We summarise the GEV pipeline here (more details can be found in White (1979); Davis et al. (2011)).First, we calculate the number density (<  UV ) of galaxies with UV magnitude less than  UV as Note that, by the definition of  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 [ 0 ,  1 ]: Next, to predict the brightest magnitude that galaxies could have at high redshifts, we consider the probability  0 that no galaxies in a given volume exceeds the extreme brightness  UV .In our calculation for the volume , we take  sky = 1 in order to obtain generic predictions on the extreme luminosities.The probability distribution  0 ( UV ) can be modelled as a Poisson distribution with the following cumulative distribution function (cdf): Differentiating this with respect to  UV , we obtain the probability density function (pdf): The Fisher-Tippett-Gnedenko theorem implies that in the large- limit, the cdf (equation ( 16)) approaches the GEV distribution given by: where  := ( UV − )/, with  describing the location of the peak, and  describing the scale of the pdf.The sign of the parameter  determines the GEV type, with  = 0,  > 0 and  < 0, corresponding to the Gumbel, Fréchet and Weibull distributions respectively.We can express the GEV parameter , ,  in terms of astrophysical parameters by Taylor-expanding the Poisson distribution  0 ( UV ) and the GEV distribution  ( UV ) around the peak of the pdf at  peak to cubic order.By equating the coefficients, we find While these equations provide a good approximation to the values ,  and , 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 the pdfs are well described by the Gumbel distribution, taking  ≈ 0.

Extreme-value pdf of the brightest galaxies
We plot the pdf for the distribution of the brightest galaxies in Fig. 2 for 4 redshift bins:  ∈ [9, 10], [10,11], [12, 13] and [16, 17].The solid lines are pdfs calculated using equation ( 17), whilst the dashed lines are the numerical fits assuming that the pdfs are Gumbel type ( = 0 in equation ( 18)).The parameters  and  for the Gumbel fit are given in Table A1.Note that in each bin,  is very close to the magnitude of  UV where the peak of the pdf occurs.
From Fig. 2, we see that at higher redshifts, the pdfs shift to higher values of  UV (meaning that the extreme galaxies are less bright).We also observe that increasing  * produces brighter extreme galaxies.In particular, increasing  * from 0.02 to 0.15 results in a pdf peak shift of Δ UV ≈ −2.13)).The colourful curves correspond to  * = 0.02, 0.05, 0.15, 0.4, 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).

Comparison with JWST data
We now study the redshift dependence of the extreme-value predictions more closely.In Fig. 3 and 4, we plot the peak of the Gumbel pdf for the brightest galaxies over the redshift range  = 9 − 16 (dashed lines).We also indicate the 95 th and 99 th percentiles of the pdfs (shaded dark and light colours).Fig. 3 shows the extremevalue predictions for  UV for the double power law and Schechter models.Fig. 4 shows the same predictions for the  * model with  * = 0.02, 0.05 and 0.15.
On these plots, we have included data from the JWST Early Release Science programs CEERS and GLASS.In particular, we have selected the brightest galaxies from a collection of recent studies on these galaxies (Donnan et al. (2023)

Double power law and Schechter models
In Fig. 3, we note that neither the double power law or 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  > 14.5, the model predicts extreme brightness of  UV ≲ −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  β/ in equation ( 2).We will revisit this point in the discussion section.
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 the extreme data points, and, more worryingly,  17)), assuming 3 values for the star formation efficiency:  * = 0.02 (purple), 0.05 (blue) and 0.15 (green).The dashed lines are the Gumbel distribution (equation ( 18)) with parameters ,  given in Table A1.we also see observed galaxies that are brighter than the predicted extreme brightness.

The 𝑓 * model
In Fig. 4 three of the star formation efficiency models have been plotted,  * = 0.02, 0.05, 0.15.The model with  * = 0.02 is comparable to the Schechter model discussed above.The model with  * = 0.05 (blue band) appears to be consistent with the data points, with the implication that these galaxies are expected to be the brightest we would observe.The model with  * = 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  * = 0.15, at  ∼ 16 we expect to see galaxies no brighter than  UV ∼ −23.5 0.8 0.4 .Future observations can be added to such a plot, giving us additional constraints on the star formation efficiency  * .

EXTENDING THE 𝑓 * MODEL
A natural extension to the  * 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  * has a dependence on the halo mass  h in the form of the double power law function: where  0 is a normalisation constant,  c is a characteristic mass where the efficiency is equal to  0 , and  and  are slopes that determine the behaviour at high and low masses.This equation has been calibrated at  = 4, with parameter values: Tacchella model: ( 0 ,  c , , ) = (0.22, 6.3 × 10 10 M ⊙ , 0.89, 0.4).( 21) We also see the same form in Harikane et al. (2022a) with parameter values: Harikane model: ( 0 ,  c , , ) = (3.2× 10 −2 , 10 11.5 M ⊙ , 1.2, 0.5).( 22) These two models of  * are plotted in Fig. 5.We can follow the same method as section 2.2 to find the luminosity function for this mass-dependent  * , by substituting equation ( 20) into the working.In particular, our equation for  UV ( h , ) (equation (8)) is now where  * =  h / c .We can then take the derivative with respect to  h and obtain the luminosity function: SFR/(fbMh) Comparison of two stellar efficiency functions  * as a function of halo mass  h .The double-power law form is given in equation ( 20), with parameter values in equations ( 21) and ( 22) for the Tacchella (pink) and Harikane (teal) models respectively.
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 95 th and 99 th percentile bands, varied over our redshift range, in Fig. 6.
The main difference between Fig. 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  * = 0.15 model.At  ∼ 16, the Tacchella curve predicts an extreme brightness of  UV = −22.5 +0.5 −1.5 , which is less bright than what the  * = 0.15 model predicts.The confidence interval in the Tacchella model is also wider at the high- 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 99 th percentile band.Around  ∼ 10, the extremevalue predictions in the Harikane model is comparable to that with constant  * model with  * of a few percent.However, the band curves . Extreme-value pdf profile (same as Fig. 3) plotted as a function of redshift, for the halo mass-dependent  * model (equation ( 20)).The teal (upper) band corresponds to the Harikane model (equation ( 22)) whilst the pink (lower) band is the Tacchella model (equation ( 21)).
up rapidly towards higher redshifts.Finally, like the Tacchella model, the band is also much wider than the constant  * 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  ≳ 9.

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  UV ) of the brightest galaxies at redshifts  ∼ 9 − 16.We found that such a distribution is well approximated by the Gumbel distribution.We studied a number of models of the galaxy luminosity function and derived the corresponding extreme-value distributions as a function of redshift.Our main results are shown in Fig. 4 and 6, in which our theoretical results are compared with data from the JWST.
Further summary and discussion points are given below.
• The Generalised Extreme-Value (GEV) formalism (equation ( 18)) was used to calculate pdf of the brightest galaxies.We find that the GEV parameter  (see Table A1) is so small that the extreme-value pdf for  UV is well approximated by the Gumbel distribution, as shown in Fig. 2. The pdf profile (along with the 95 th and 99 th percentiles) are shown in Fig. 3, 4 and 6, assuming different galaxy luminosity functions as shown in Fig. 1.
• The Schechter model of the luminosity function has previously been used to model the luminosity of galaxies at  ∼ 5.The parameters used in this work (equation (3)) were calibrated using observations at  ∼ 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  ≳ 10.
• 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 parametrisation 5HGVKLIW] (equation ( 2)).Indeed, we found an ad-hoc correction by switching the sign of d β/ d to negative.Donnan et al. (2024) gave new constraints on the redshift evolution of the double power law luminosity function.We use their data in place of the parametrisation (equation ( 2)), and obtained the extremevalue pdf profile shown in Fig. 7, which appears to be consistent with data.We note that in their work, d β/ d is indeed negative.Since the parameter β describes the bright-end slope of the curve, we conclude that a physical double power law parametrisation must satisfy the condition d β/ d < 0.
• The  * model of the luminosity function (as outlined in Harikane et al. (2024)) was considered as an alternative to the Schechter and DPL models.This model requires the stellar formation rate  * , 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,  * = 0.02 produces similar result to that of the Schechter model, whilst  * between 0.05 and 0.15 produces an extreme-value pdf that is consistent with JWST data.
The model predicts that at  ∼ 16, the brightest galaxies will have a UV magnitude 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 predict the halo number density accurately at low redshifts ( ∼ 1−2) and a recalibration to high- simulations is needed to confirm if the parameter values remain accurate for  ∼ 9 − 16.
• Extending the  * 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  * = 0.15 model, whilst the Harikane model is comparable to  * of a few percent, with the predictions diverging at  ∼ 16.As more high redshift galaxies are observed, the extreme-value formalism will allow us to distinguish whether  * is constant or halo massdependent.13)) with  * = 0.05 (blue) and  * = 0.15 (green), variable  * (equation ( 24)) (pink), and the doublepower law model (equation ( 2)) assuming the parameters of Donnan et al. (2024) (black).
that are consistent with JWST observation of the brightest highredshift galaxies.
• 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.

Figure 7 .
Figure7.Extreme-value pdf profile (same as Fig.3) plotted as a function of redshift, for the double power law luminosity function with parameters given inDonnan et al. (2024).

Fig. 8 Figure 8 .
Fig.8summarises 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  * and variable  * models