alphaFeedback

WorksPhysics and Astronomy

A combined PRIMER+MINERVA study of the evolution of \mathbf{[O\mbox{III}]+Hβ} equivalent width from z=3\mathbf{z=3} to z=8.5\mathbf{z=8.5}

AI contributions
Not recorded for this paper
View details

Abstract

We present an analysis of the evolution of [O III]+Hβ[O\,III]+\mathrm{H}β equivalent width (WλW_λ) over the redshift range z=3z=3 to z=8.5z=8.5, based on combining the JWST NIRCam PRIMER and MINERVA imaging in the UDS survey field. We find that the combined medium- and broad-band data set reduces the photometric redshift outlier rate to η=4.7%η= 4.7\% at mF444W≤28m_{\rm F444W}\leq28, a 42%42\% reduction compared with PRIMER imaging alone. We also demonstrate that the combined data set facilitates the accurate recovery of [O III]+Hβ[O\,III]+\mathrm{H}β equivalent widths in redshift bins centred on z∼3.2, 4.0, 5.0, 7.2, 8.3z\sim3.2,\, 4.0,\, 5.0,\, 7.2,\, 8.3; corresponding to the F210M, F250M, F300M, F410M and F460M filters. Based on a final sample of N=1964N=1964 galaxies with MUV≤−19.25M_{\rm UV}\leq-19.25, we first provide new constraints on the bright end of the z=7−8z=7-8 [O III]+Hβ[O\,III]+\mathrm{H}β luminosity function (log⁡10ϕ∗=−4.3−0.5+0.2\log_{10}ϕ^{*} = -4.3_{-0.5}^{+0.2}, log⁡10L∗=−43.4−0.2+0.3\log_{10}L^{*} = -43.4_{-0.2}^{+0.3}), before proceeding to investigate the redshift evolution of [O III]+Hβ[O\,III]+\mathrm{H}β equivalent width. Using SED modelling based on non-parametric continuity star-formation histories, we find that the median [O III]+Hβ[O\,III]+\mathrm{H}β equivalent width monotonically increases with redshift, rising from log10(Wλ)=2.54−0.01+0.01\mathrm{log}_{10}(W_λ)=2.54^{+0.01}_{-0.01} at z≃3z\simeq 3 to log10(Wλ)=3.05−0.03+0.03\mathrm{log}_{10}(W_λ)=3.05^{+0.03}_{-0.03} at z≃8z\simeq 8. In tandem, the distribution of log10(Wλ)\mathrm{log}_{10}(W_λ) is observed to become increasingly asymmetric with redshift, producing a fraction of extreme emission-line galaxies (EELGs, Wλ>1000 W_λ>1000\,Å) that reaches fEELG=0.50±0.10f_{\rm EELG}=0.50\pm{0.10} by z=8.3z=8.3. Our SED modelling results indicate that the physical driver of this evolution is increasingly bursty star-formation, with EELGs the result of particularly high specific star-formation rates over the 10 Myr prior to observation.

Introduction

The photometric and spectroscopic capabilities of JWST have established that galaxies at high redshift are synonymous with extremely strong nebular emission (e.g. Katz et al. 2025). Given the comparative scarcity of extreme emission-line galaxies (EELGs) in the local Universe (e.g. Amorín et al. 2015; Berg et al. 2021), understanding the evolution in nebular emission will provide direct insight into the star-formation rates, ages, metallicities and general interstellar medium (ISM) conditions of the high-redshift galaxy population (e.g. Stark 2016). Moreover, nebular emission can also provide crucial information on the production, and subsequent escape, of ionizing photons from galaxies within the epoch of reionization.

The strength of rest-frame optical emission lines such as [O iii] λλ4959,5007\lambda\lambda4959, 5007 Å reflects the production rate of massive O- and B-type stars on a timescale of ≃10\simeq10 Myr. This is in contrast to the optical continuum, which will typically reflect the contribution of the lower-mass, longer-lived stellar population. Consequently, the high equivalent widths first observed at z≳6z \gtrsim6 using a combination of Hubble Space Telescope (HST) and Spitzer photometry (e.g. Labbé et al. 2013; Smit et al. 2014; Roberts-Borsani et al. 2016), suggest that many high-redshift galaxies have recently experienced a rapid increase in their specific star-formation rate (sSFR).

This naturally leads to the conclusion that star formation in high-redshift galaxies is increasingly bursty/stochastic (Tang et al. 2019; Vanzella et al. 2020), due to a combination of low stellar masses, high gas fractions and short dynamical timescales (e.g. Faucher-Giguère 2018; McClymont et al. 2025). Indeed, bursty star formation has been widely invoked to explain the high number of UV-bright galaxies seen at z≳9z \gtrsim9 (e.g. Castellano et al. 2022; Naidu et al. 2022; Donnan et al. 2023; Shen et al. 2023; Mason et al. 2023; Gelli et al. 2024; McLeod et al. 2024) although, strictly, it appears that accounting for increasingly youthful stellar populations at earlier times is a sufficient explanation (Donnan et al. 2025).

In the bursty star-formation scenario, the nature of high-redshift galaxies that display low equivalent widths is potentially as interesting as that of extreme line emitters. Most of these galaxies are likely to be in a phase of reduced star-formation activity following a rapid burst, i.e. so-called mini-quenched galaxies (e.g. Looser et al. 2024; Trussler et al. 2025). However, the low equivalent widths displayed by at least some of these galaxies could indicate ultra-low metallicities and/or a high escape fraction of Lyman continuum photons (Endsley et al. 2025).

Prior to the advent of JWST, our ability to spectroscopically study rest-frame optical nebular emission was confined to z≤3.5z \le3.5 by the available atmospheric windows in the near-IR (e.g. Kriek et al. 2015; Holden et al. 2016; Reddy et al. 2018). However, this situation was revolutionized by the availability of NIRSpec, which quickly made the detailed study of strong line-emitting galaxies at high redshift relatively straightforward (e.g. Carnall et al. 2023; Tacchella et al. 2023; Roberts-Borsani et al. 2024).

At intermediate redshifts (z∼2z \sim2), the equivalent width (WλW_{\lambda}) of [O iii] is known to be highly correlated with ionizing photon production efficiency (ξion\xi_{\mathrm{ion}}; Tang et al. 2019). Recent NIRSpec studies have found a similar WλW_{\lambda} − ξion\xi_{\mathrm{ion}} relation at z≥6z \ge6, together with a strong correlation between [O iii] WλW_{\lambda} and β\beta, the spectral slope of the UV continuum (e.g. Begley et al. 2026; Papovich et al. 2026). Both correlations are consistent with the expectation that the nebular emission of high-redshift galaxies is dominated by a young, low-metallicity stellar population, with low dust attenuation.

Whilst the rest-frame optical nebular emission of individual high-redshift objects is best studied via spectroscopy, deep NIRCam imaging from legacy surveys such as PRIMER (Dunlop et al. 2021), CEERS (Finkelstein et al. 2025) and JADES (Eisenstein et al. 2026) has provided an opportunity to study the evolution of the equivalent width of [O iii] + Hβ\beta for large, statistically complete samples. Recent photometric studies have demonstrated that extreme emission-line galaxies (EELGs: sources with [O iii] + Hβ\beta Wλ≥1000W_{\lambda} \ge1000 Å) were much more common in the early Universe, and that the typical [O iii] + Hβ\beta WλW_{\lambda} increases from ≃100\simeq100 Å at z≃2z \simeq2 (Boyett et al. [9]) to ≃700\simeq700 Å at z≃7.5z \simeq7.5 (e.g. Begley et al. 2025; Endsley et al. 2025).

Until recently, photometric studies of [O iii] + Hβ\beta WλW_{\lambda} have been somewhat limited by the lack of suitably deep, multi-band, JWST medium-band imaging over large areas. In this context, deep medium-band imaging in multiple filters is crucial because, compared to broad-band imaging alone, it has superior emission-line sensitivity and a greater ability to constrain the underlying continuum level. The Medium band Imaging with NIRCam to Explore Revolutionary Astrophysics (MINERVA, Muzzin et al. 2025) survey aims to address this issue by observing areas of sky with existing NIRCam broad-band imaging in eight NIRCam medium-bands. The large number of medium-band filters, when combined with previous broad-band data, opens up the opportunity for accurate and consistent nebular emission-line studies over a large range in redshift. Crucially, combining the MINERVA data with existing HST and JWST imaging in the PRIMER UDS and PRIMER COSMOS fields, allows photometric catalogues to be constructed featuring 19 filters, covering the wavelength range 0.40.4–5.0 μm5.0\,\mu\mathrm{m}.

Motivated by the completion of the MINERVA medium-band imaging in the PRIMER UDS field, we have exploited the combined PRIMER+MINERVA data set to consistently study the evolution of the [O iii] + Hβ\beta equivalent width distribution and the fraction of EELGs, within a MUVM_{\mathrm{UV}}–complete galaxy sample spanning the redshift range 3<z<8.53 < z < 8.5. As part of this process, we have also investigated the improvement in photometric redshift accuracy provided by the MINERVA medium-band imaging, and provide new constraints on the bright end of the [O iii] + Hβ\beta luminosity function at z≃7z \simeq7. Based

on the results of SED fitting using non-parametric star-formation histories, we also investigate the physical properties of emission-line galaxies as a function of redshift and [O iii] + Hβ\beta equivalent width.

The paper is structured as follows. In § 2 we describe the production of the combined photometric catalogue, summarising the aperture photometry, point spread function (PSF) homogenization, dust corrections, and depth measurements. In § 3 we describe the data analysis, including details of the photometric-redshift fitting, sample selection, the extraction of accurate equivalent width information from the photometry, and an estimate of the z≃7z \simeq7 [O iii] + Hβ\beta luminosity function. We present our [O iii] + Hβ\beta emission-line results in § 4 and discuss their implications in § 5. Finally, we summarise our conclusions in § 6.

Throughout this paper we use the following cosmological parameters: H0=70 km s−1 Mpc−1H_0 = 70\,\mathrm{km\,s^{-1}\,Mpc^{-1}}, Ωm=0.3\Omega_{\mathrm{m}} = 0.3, ΩΛ=0.7\Omega_{\Lambda} = 0.7, along with a Kroupa (2001) initial mass function (IMF). All magnitudes are presented in the AB system (Oke & Gunn 1983).

Data

Survey Imaging

This work combines the JWST PRIMER and MINERVA data sets with the existing HST ACS imaging (Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey: Grogin et al. 2011; Koekemoer et al. 2011) in the UDS survey field, to produce a joint photometry catalogue featuring a total of 19 wide- and medium-band filters. The PRIMER survey (GO-1837, PI: Dunlop, Dunlop et al. 2021; Magee et al., in preparation) was a JWST Cycle-1 Treasury Program providing deep photometric imaging observations in the CANDELS UDS and COSMOS fields. The PRIMER observations consist of imaging in eight NIRCam filters (F090W, F115W, F150W, F200W, F277W, F356W, F410M, F444W), with parallel imaging in two MIRI filters (F770W and F1800W). The MINERVA survey (GO-7814, PI: Muzzin, Muzzin et al. 2025) is an on-going Cycle-4 program obtaining NIRCam medium-band imaging over an area of ≃542\simeq542 arcmin2^2 in the UDS, COSMOS, EGS and GOODS-N survey fields.

The first MINERVA field to be completed was the PRIMER UDS field, for which all observations were successfully taken by January 2026. These observations consist of NIRCam imaging in eight medium-band filters (F140M, F162M, F182M, F210M, F250M, F300M, F360M, F460M), together with MIRI parallel imaging in F1280W and F1500W. Both the PRIMER and MINERVA data sets were reduced with PENCIL: a bespoke version of the JWST pipeline v1.17.1 (Magee et al., in preparation).

The area of the UDS field covered by the combined 16 NIRCam-filter imaging of PRIMER and MINERVA is ∼234\sim234 arcmin2^2, while an area of ≃155\simeq155 arcmin2^2 is also covered by all three HST ACS filters (F435W, F606W & F814W). Below, we briefly detail the production of the combined photometry catalogue, with full details describing the production of combined PRIMER+MINERVA catalogues over both the PRIMER UDS and PRIMER COSMOS fields deferred to a future paper (Holst et al., in preparation).

PSF Homogenisation

In order to obtain accurate photometry, each image was first PSF homogenised to match the angular resolution of the F444W imaging. For each filter, PSFs were initially constructed using stacks of bright, unsaturated and isolated stars. Transfer kernels were then calculated in Fourier space to match each filter’s PSF to the F444W PSF, using varying window functions to reduce high-frequency artefacts such as fringing. Following convolution with the appropriate transfer kernel, the curve-of-growth in each filter matched that of the F444W imaging to within ∼1%\sim1\% at all radii.

Photometry

Photometry was measured using SourceExtractor ((Arnouts et al. 1999), Version 2.28.0) in dual-image mode; using the unconvolved F356W image for object detection and the PSF-homogenised images to measure the photometry. The photometry was extracted in both 0.3′′0.3^{\prime\prime} and 0.5′′0.5^{\prime\prime}-diameter circular apertures, and corrected for galactic extinction using the Gordon et al. (2023) extinction curve and AV=0.089A_V = 0.089, as appropriate for the UDS field 1. Robust, object-by-object photometry errors were calculated from the aperture-to-aperture rms of the nearest 200 blank sky apertures (using object masks based on the SourceExtractor segmentation maps). The median global depths for each of the filters in 0.3′′0.3^{\prime\prime}-diameter circular apertures (corrected to total) are presented in Table 1.

Filter5σglobal5\sigma_{\mathrm{global}}Filter5σglobal5\sigma_{\mathrm{global}}
F435W†27.39F140M27.10
F606W†28.08F162M27.16
F814W†28.09F182M27.43
F090W27.71F210M27.25
F115W27.75F250M27.33
F150W27.95F300M27.69
F200W28.12F360M27.69
F277W28.38F460M26.71
F356W28.37
F410M27.72
F444W27.83

Table 1. The 5σ5\sigma global depths in all 19 filters. The depths are calculated in a 0.3′′0.3''-diameter circular aperture and then corrected to total assuming a point-source morphology. The left-hand side lists the HST †ACS and PRIMER NIRCam filters, with the right-hand side listing the MINERVA NIRCam filters.

Analysis

Photometric Redshifts

Photometric redshifts were obtained using both eazyPy ((Brammer et al. 2008)) and LePhare ((Arnouts et al. 1999; Ilbert et al. 2006), using the python wrapper). With eazyPy we used four different template sets: two sets introduced in Larson et al. (2023) and Hainline et al. (2024) to model high-redshift galaxies, the "AGN_BLUE_SFHZ_13" 2 template set (which crucially contains the SED of a Little Red Dot from Killi et al. (2024)) and an adaptation of the default eazyPy template set with three additional high WλW_\lambda emission-line templates (Begley et al. (2025)). With LePhare, we adopted the "BC03_CHAB_MOD" Bruzual & Charlot (2003) templates, along with the pHYS emission-line prescription.

The photometric redshift estimates were performed on both the 0.3′′0.3^{\prime\prime} and 0.5′′0.5^{\prime\prime}-diameter aperture photometry. For the vast majority of sources (87% of the catalogue), the final photometric redshifts were based on the higher SNR 0.3′′0.3^{\prime\prime}-diameter photometry. However, for bright objects (mF356W<25.2m_{\mathrm{F356W}} < 25.2) we adopted the photometric redshifts derived from the 0.5′′0.5^{\prime\prime}-diameter photometry. For this subset of sources, that are typically extended, testing revealed that the larger aperture delivered better colours and superior photometric-redshift accuracy.

For each object, we adopted the median of the five different redshift runs as the final photometric redshift. To evaluate their accuracy, we compared our photometric redshifts to a compilation of N=3145N = 3145 high-quality spectroscopic redshifts within the 234 arcmin2^2 PRIMER+MINERVA footprint. The majority of these (2491) were grade=3 spectroscopic redshifts from the Dawn JWST Archive (DJA, Valentino et al. 2025; Pollock et al. 2026). The remaining 654 spectroscopic redshifts were drawn from a variety of ground-based spectroscopic surveys including: VANDELS (McLure et al. 2018; Garilli et al. 2021), UDSz (Bradshaw et al. 2013; McLure et al. 2013), MOSDEF (Kriek et al. 2015) and Galametz et al. (2018). The ground-based spectroscopic comparison sample is typical of the bright, blue samples of spectroscopic objects traditionally used to test the accuracy of photometric redshift fitting. In contrast, the spectra of DJA spectra are fainter, redder and more typical of the galaxies featuring in our PRIMER+MINERVA photometry catalogue.

To evaluate the performance of our photometric-redshift fitting, we measured: dz=(zphot−zspec)/(1+zspec)d_z = (z_{\mathrm{phot}} - z_{\mathrm{spec}})/(1 + z_{\mathrm{spec}}), the catastrophic outlier fraction, η\eta (where ∣dz∣≥0.15|d_z| \geq0.15), the standard deviation, σ=1.483×MAD(dz)\sigma= 1.483 \times\mathrm{MAD}(d_z) and the bias, b=Median(zphot−zspec)b = \mathrm{Median}(z_{\mathrm{phot}} - z_{\mathrm{spec}}). The full set of statistics are reported in Table 2 for those objects that have SNR ≥5\geq5 in the F356W filter and benefit from the full suite of relevant filters. Based on the PRIMER NIRCam plus HST ACS photometry, our photometric redshifts have an overall accuracy of σ=0.033\sigma= 0.033 and an outlier fraction of 8.1%. However, when the medium-band data from MINERVA are included, the overall accuracy improves to σ=0.015\sigma= 0.015 and the outlier fraction drops to 4.7%. The inclusion of MINERVA’s medium-bands therefore improves the photometric-redshift accuracy by more than a factor of two and reduces the catastrophic outlier rate by 42%.

(a) Total
PhotometryNη\etaσ\sigmabb
PRIMER25998.12%0.0330.008
PRIMER+MINERVA31454.71%0.0150.001
PRIMER+MINERVA σz≤0.66\sigma_z \leq0.6629722.69%0.0140.001
(b) Ground
PhotometryNη\etaσ\sigmabb
PRIMER6542.91%0.0230.004
PRIMER+MINERVA6540.92%0.012-0.002
PRIMER+MINERVA σz≤0.66\sigma_z \leq0.666500.92%0.012-0.001
(c) DJA
PhotometryNη\etaσ\sigmabb
PRIMER19459.87%0.0390.010
PRIMER+MINERVA24915.70%0.0160.002
PRIMER+MINERVA σz≤0.66\sigma_z \leq0.6623223.19%0.0150.002

Table 2. The photometric redshift statistics for: (a) the total spectroscopic sample, (b) the ground-based spectroscopic sample and (c) the DJA spectroscopic sample. These results are shown with and without MINERVA photometry and with and without the redshift run standard deviation (σz\sigma_z) cut.

One benefit of calculating our photometric redshifts from five different runs is that we are able to easily identify objects where the five best-fitting redshifts have a large standard deviation (σz\sigma_z). As might be hoped, this provides a useful method for identifying objects that are more likely to be catastrophic outliers. In the final row of Table 2, we

display the statistics for the PRIMER+MINERVA photometry with the additional criterion that σz≤0.66\sigma_z \le0.66, resulting in an outlier rate of only η=2.69%\eta= 2.69\%, down to a magnitude limit of mF444W≤28m_{\mathrm{F444W}} \le28. This criterion is therefore capable of reducing the outlier fraction by a further 43%, while excluding ≤7%\le7\% of the overall sample. Figure 1 shows the comparison between photometric and spectroscopic redshifts for the three photometric samples presented in Table 2.

Figure 1. A comparison of photometric vs spectroscopic redshift for the PRIMER photometry (left), the PRIMER+MINERVA photometry (centre), and the PRIMER+MINERVA photometry including a σz≤0.66\sigma_z \le0.66 cut (designed to remove objects with photometric redshift distributions that have a high standard deviation, see § 3.1 for discussion). Black dotted lines show the outlier limits (∣dz∣≥0.15|d_z| \ge0.15), outside of which an object is classified as a catastrophic outlier and plotted as a grey triangle. The non-outlier objects are coloured by the spatial density of points. We see a 42% reduction in the outlier fraction (η\eta) when the MINERVA photometry is added to the PRIMER photometry. We see a further 43% reduction in the outlier rate when including the σz≤0.66\sigma_z \le0.66 cut.

As expected, in all cases, the fidelity of the photometric redshifts is significantly better when tested against the bright, blue, ground-based spectroscopic sample, than the fainter, redder, DJA spectroscopic sample. That said, the improvement in photometric-redshift outlier rates for the more representative DJA spectroscopic sample is particularly significant. Using the photometric redshifts provided by the DJA (PRIMER data only) results in an outlier fraction of 12.9%, compared to our outlier fraction of 5.7% (or 3.2% including the σz\sigma_z cut). Moreover, isolating only objects from the RUBIES survey (PID-4233; PIs: de Graaff & Brammer), we measure η=3%\eta= 3\%; compared to the η=6%\eta= 6\% obtained by de Graaff et al. (2025) using only PRIMER data.

[Oiii] + Hβ\beta Equivalent Widths

In addition to improving overall photometric redshift accuracy, the inclusion of the MINERVA medium-band imaging makes it possible to study nebular line emission over a wide range in redshift. In particular, the availability of F406M imaging can constrain the [O iii] + Hβ\beta equivalent width out to 8≤z≤8.58 \le z \le8.5, while the F410M PRIMER imaging provides an opportunity to benchmark our results against previous studies at z∼7z \sim7 (e.g. Meyer et al. 2024; Begley et al. 2025).

To determine the [O iii] + Hβ\beta equivalent widths (WλW_\lambda) we modelled the 19-filter photometry data set with bagpipes (Bayesian Analysis of Galaxies for Physical Inference and Parameter Estimation, Carnall et al. 2018), closely following the methodology described in Begley et al. (2025). We tested both the Binary Population and Spectral Synthesis (BPASS) v2.2.1 (Eldridge et al. 2017; Stanway & Eldridge 2018) and BC03 (Bruzual & Charlot 2003, updated 2016 version) stellar population synthesis (SPS) models, along with a variety of different star formation histories (SFHs): delayed-τ\tau models (see e.g. Carnall et al. 2019 for parametric model analysis), constant, and an assortment of hyper-parametrisations of the non-parametric continuity SFH presented in Leja et al. (2019).

To fine tune our bagpipes setup we used spectuhs to measure the [O iii] + Hβ\beta equivalent widths for 115 galaxies drawn from our full sample that have publicly available NIRSpec PRISM spectra on the DJA. These sources have SNR ≥5\ge5 in the F356W filter and MUV≤−19.25M_{\mathrm{UV}} \le-19.25, to match the constraints we use in § 3.3 for photometric selection. We allow the redshift of the galaxies to vary with a Gaussian prior (σ=0.00333\sigma= 0.00333). To ensure accurate spectroscopic measurements of the [O iii] + Hβ\beta equivalent widths, these spectra were chosen to have SNR ≥5\ge5 over the [O iii] + Hβ\beta complex, together with significant detections of the adjacent continuum: median SNR/pix ≥1.5\ge1.5 in the continuum fitting regions.

The best recovery of the spectroscopic equivalent widths were delivered by a non-parametric continuity SFH, using the BPASS SPS models and the bagpipes default BPASS cloudy nebular emission models (Ferland et al. 2013). To make the Leja et al. (2019) SFH allow larger variations in SFR, we increased the log⁡10(SFR)\log_{10}(\mathrm{SFR}) prior scale from 0.3 to 1, and decreased the log⁡10(SFR)\log_{10}(\mathrm{SFR}) prior d.f. from 2 to 1.5. These hyper-parameters dictate σ\sigma and ν\nu, respectively, in the Student’s tt-distribution, with the chosen values similar to those used in previous studies (see e.g. Tacchella et al. 2022; Begley et al. 2025). We use age bins of 0–3, 3–10, 10–100, 100–250, and 250–500 Myr. The full set of bagpipes parameters used in our final analysis is presented in Table 3.

ParameterRangePriorHyper-parameters
Redshiftμ−3σ≤zphot≤μ+3σ\mu- 3\sigma\le z_{\mathrm{phot}} \le\mu+ 3\sigmaGaussianμ=zphot, σ=0.05†\mu= z_{\mathrm{phot}},\ \sigma= 0.05^{\dagger}
Total stellar mass formed106≤M∗/M⊙≤101210^6 \le M_\ast/M_\odot\le10^{12}log⁡10\log_{10}
Stellar Metallicity0.005≤Z∗/Z⊙≤10.005 \le Z_\ast/Z_\odot\le1Uniform
SFR Change (i→i+1i \to i+1)−10≤Δi(log⁡10(SFR))≤10-10 \le\Delta_i\left(\log_{10}(\mathrm{SFR})\right) \le10Student-tprior scale = 1, prior d.f. = 1.5
V-band attenuation0≤AV≤30 \le A_V \le3Uniform
Nebular ionisation parameter−4≤log⁡(U)≤−1-4 \le\log(U) \le-1Uniform

Table 3. The parameter ranges and priors used in the bagpipes SED fitting. † This was changed to 0.00333 for the spectroscopic sample due to the higher reliability of these redshifts.

For each object, we obtain the line and continuum fluxes directly from the bagpipes output parameter posteriors, generating a distribution of Nposterior=1000N_{\mathrm{posterior}} = 1000 equivalent widths. We take the 50th percentile of the resulting distribution to be our WλW_{\lambda} measurement, and the 16th and 84th percentiles as our 1σ1\sigma errors. Figure 2 shows the comparison of our photometric and spectroscopic equivalent widths, coloured by the filter in which the emission lines were observed. The relationship between photometric and spectroscopic equivalent widths is close to one-to-one, with a scatter of σ(Wλ)=337\sigma(W_{\lambda}) = 337 Å. To account for the minor discrepancy between our photometric measurements and their spectroscopic counterparts, we fit the data with a linear relationship, finding Wλ,spec=(1.04±0.04)Wλ,photW_{\lambda,\mathrm{spec}} = (1.04 \pm0.04)W_{\lambda,\mathrm{phot}}. All future WλW_{\lambda} measurements are corrected using this relation.

Figure 2. Photometric versus spectroscopic measurements of rest-frame [O III] + HβH\beta equivalent width for our sub-sample of N=115N = 115 galaxies with high-quality NIRSpec PRISM spectra (see Section 3.2). The dashed line shows the 1:1 relation. The data points are coloured by the filter within which the emission lines were observed (F460M in green, F410M in orange, F300M in blue, F250M in red, and F210M in pink). The upper panel shows the relative difference in WλW_{\lambda} against the spectroscopic WλW_{\lambda}. We find a relationship of Wλ,spec=(1.04±0.04)Wλ,photW_{\lambda,\mathrm{spec}} = (1.04 \pm0.04)W_{\lambda,\mathrm{phot}}.

Sample Selection

Our initial sample of galaxies comprised objects with SNR ≥5\geq5 in the F356W filter and a photometric redshift in the ranges 2.95≤z≤3.42.95 \leq z \leq3.4, 3.83≤z≤4.23.83 \leq z \leq4.2, 4.65≤z≤5.34.65 \leq z \leq5.3, 6.7≤z≤7.66.7 \leq z \leq7.6, or 8.0≤z≤8.58.0 \leq z \leq8.5, corresponding to the [O III] + HβH\beta complex being redshifted into the F210M, F250M, F300M, F410M and F460M filters. In each redshift bin, the sample was further restricted to those objects with SNR ≥3\geq3 in the corresponding medium-band filter and a careful visual inspection was performed to remove any remaining artefacts. In this study we deliberately avoid using the F360M filter to measure [O III] + HβH\beta equivalent widths in the redshift range 5.8≤z≤6.655.8 \leq z \leq6.65. This is because the F360M filter sits within the broad F356W filter, meaning that our parent photometry catalogue is biased towards high [O III] + HβH\beta equivalent width objects.

To exclude potential active galactic nuclei (AGN) we removed N=13N = 13 sources detected in X-rays within 0.5″ (using Kovcevski et al. 2018), and N=35N = 35 sources in common with the Little Red Dot (LRD) samples of Kovcevski et al. (2023) and Perger et al. (2025). At this stage, our galaxy sample consists of N=5325N = 5325 objects, over the redshift range 2.95<z<8.502.95 < z < 8.50. However, in order to mitigate the evolving MUVM_{\mathrm{UV}} limit, we restrict our final sample to those objects with MUV≤−19.25M_{\mathrm{UV}} \leq-19.25, after conservatively estimating that our sample is complete to this limit across the full redshift range. The final sample, used for all subsequent analysis, consists of N=1964N = 1964 objects, with N=633N = 633 in the 2.95≤z≤3.42.95 \leq z \leq3.4 bin, N=584N = 584 in the 3.83≤z≤4.23.83 \leq z \leq4.2 bin, N=430N = 430 in the 4.65≤z≤5.34.65 \leq z \leq5.3 bin, N=268N = 268 in the 6.7≤z≤7.66.7 \leq z \leq7.6 bin and N=49N = 49 in the 8.0≤z≤8.58.0 \leq z \leq8.5 bin.

[Oiii] + H Luminosity Function

As a final consistency check on our sample selection and SED-fitting methodology, we computed the [O III] + HβH\beta luminosity function (LF) at z≃7.3z \simeq7.3 by combining our F410M and F460M sub-samples.

LF estimate is based on the 1/Vmax⁡1/V_{\max} method (Schmidt 1968), with final errors calculated as the quadrature sum of bootstrap resampling and cosmic variance (Trenti & Stiavelli 2008) uncertainties.

In Figure 3 we show our new estimate of the bright-end of the LF at z≃7.3z \simeq7.3, together with recent literature results at z≃7.1−8.3z \simeq7.1 - 8.3 from Meyer et al. (2024), Meyer et al. (2026), Wold et al. (2025), Wold et al. (2026), and Korber et al. (2026). The Schechter (Schechter 1976) function fit to the z≃7z \simeq7 LFs by Wold et al. (2026) is also shown as a pink dotted line. Where necessary, we have corrected [O iii]λ\lambda5007 luminosities to [O iii]+Hβ\beta luminosities using a factor of 1.48 (Begley et al. 2025). Our new LF estimate is plotted at log⁡10L[O III]+Hβ≥42.6\log_{10} L_{\mathrm{[O\,III]+H\beta}} \geq42.6, which corresponds to our conservative completeness limit, based on the 5σ5\sigma limiting line flux from the bagpipes SED modelling and the turnover in the number counts at fainter luminosities.

!The [O III] + Hβ luminosity function bbox=0.53,0.07,0.90,0.34

Figure 3. A new estimate of the bright-end of the [O III] + HβH\beta luminosity function at redshift z≃7.3z \simeq7.3 (filled blue circles), based on combining our 8.0≤z<8.58.0 \leq z < 8.5 and 6.7≤z<7.66.7 \leq z < 7.6 sub-samples. We also show two LF estimates at z≃7.1z \simeq7.1 from Meyer et al. (2024) and Meyer et al. (2026), an LF estimate at z≃7z \simeq7 by Wold et al. (2025) and the two highest luminosity data points from the z≃8.3z \simeq8.3 LF estimate of Korber et al. (2026). Where necessary we have corrected from [O III]λ\lambda5007 to [O III] + HβH\beta luminosities. The dashed line shows the Schechter function fit at z≃7z \simeq7 from Wold et al. (2026). The solid black line shows a combined fit to our new bright-end data and the fainter Wold et al. (2026) data (see § 3.4 for details).

Our luminosity function estimate agrees well with the z=7z = 7 LF of Wold et al. (2026) and also follows cleanly on from the z=7z = 7 results of Wold et al. (2025). Our bright-end constraints are also consistent with, albeit marginally higher than, the data points from Meyer et al. (2024) and Meyer et al. (2026) at z≃7.1z \simeq7.1, and the Schechter function fit from Wold et al. (2026).

Due to our completeness limit, our new results provide no constraint on the faint-end slope of the LF. We therefore combine our data with that of Wold et al. (2026) (disregarding their brightest two single-object bins) and fit the combined data set at log⁡10L[O III]+Hβ≥41.4\log_{10} L_{\mathrm{[O\,III]+H\beta}} \geq41.4. Our best-fitting Schechter function is shown as the solid black line in Figure 3 and has the following parameters: log⁡10(ϕ∗/Mpc−3 dex−1)=−4.3−0.5+0.2\log_{10}(\phi^*/\mathrm{Mpc}^{-3}\,\mathrm{dex}^{-1}) = -4.3^{+0.2}_{-0.5}, log⁡10(L∗/erg s−1 dex−1)=−43.4−0.2+0.2\log_{10}(L^*/\mathrm{erg\,s}^{-1}\,\mathrm{dex}^{-1}) = -43.4^{+0.2}_{-0.2}, and α=−2.15−0.07+0.05\alpha= -2.15^{+0.05}_{-0.07}.

It can be seen that our new Schechter function is 0.5 dex higher than the Wold et al. (2026) fit at log⁡10L[O III]+Hβ≥43.75\log_{10} L_{\mathrm{[O\,III]+H\beta}} \geq43.75. This could be a result of the different observational methods, with Wold et al. (2026) using the slitless-spectroscopy of Meyer et al. (2024) to constrain the bright end of their LF. It will be possible to confirm or refute this result with the larger area that will be provided by the final complete MINERVA + PRIMER/CEERS/JADES data set.

Results

In Section 3.2 we demonstrated that detailed SED fitting of the combined PRIMER+MINERVA data set provides accurate emission-line measurements for z≃3−5z \simeq3 - 5 galaxies. In this section we investigate the evolving [O iii] + Hβ\beta equivalent width distributions and underlying physical properties of our galaxy samples at z≃3.2z \simeq3.2, z≃4.0z \simeq4.0, z≃5.0z \simeq5.0, z≃7.2z \simeq7.2, and z≃8.3z \simeq8.3.

As [O iii] + Hβ\beta equivalent width provides a direct probe of nebular emission and recent star-formation, tracing its evolution should offer insight into changes in the typical star-formation history of our sample as a function of redshift. Moreover, by focusing on the prevalence and physical properties of EELGs, it should be possible to obtain a clearer picture of the importance of bursty star formation at high redshift.

W Distributions

We present the distributions of photometric [O iii] + Hβ\beta equivalent widths for all redshift samples in Figure 4. For each object in the samples, a random WλW_{\lambda} was chosen from its posterior distribution, and a Gaussian kernel density estimate (KDE) was taken of the random draws to better visualise the underlying distribution. This method was repeated 1000 times, creating a posterior distribution of normalised KDE counts at each value of log⁡10(Wλ)\log_{10}(W_{\lambda}), from which the 16th^{\mathrm{th}}, 50th^{\mathrm{th}}, and 84th^{\mathrm{th}} percentiles were taken. The median values of log⁡10(Wλ)\log_{10}(W_{\lambda}) for the F210M, F250M, F300M, F410M, and F460M samples are listed in the Figure 4 legend. It is clear that the median value of log⁡10(Wλ)\log_{10}(W_{\lambda}) is monotonically increasing with redshift. Moreover, it can also be seen that the WλW_{\lambda} distributions become increasingly asymmetric at the highest redshifts, with longer tails towards low equivalent widths. We discuss the implications of these results, along with comparisons to previous literature studies in Section 5.1.

Figure 4. Kernel density estimates of the [O iii] + Hβ\beta equivalent width distributions for the F210M (purple, 2.95≤z≤3.42.95 \leq z \leq3.4), F250M (red, 3.83≤z≤4.23.83 \leq z \leq4.2), F300M (green, 4.65≤z≤5.34.65 \leq z \leq5.3), F410M (orange, 6.7≤z≤7.66.7 \leq z \leq7.6) and F460M (blue, 8.0≤z≤8.58.0 \leq z \leq8.5) sub-samples. The shaded regions surrounding the lines are the ±1σ\pm1\sigma regions calculated by Monte Carlo sampling of the individual posterior distributions. The higher redshift samples show a deviation from the log-normal shape seen at lower redshifts (e.g. [9] (Lee et al. 2007; Ly et al. 2010; Stark et al. 2013; Schenker et al. 2014; Begley et al. 2025)).

Extreme [O iii] + Hβ Emitters

The evolution in the fraction of so-called extreme emission-line galaxies, those galaxies displaying very high equivalent width rest-frame optical emission lines, has been investigated in a number of recent studies. Whilst rare in the local universe, JWST-based studies have shown that the fraction of EELGs was much higher at earlier times (e.g. Sun et al. (2023); Matthee et al. (2023); Topping et al. (2024); Boyett et al. (2024)), with their enhanced prevalence typically attributed to increasingly bursty star formation at high redshift.

Throughout our analysis, we define an EELG as an object with a rest-frame [O iii] + Hβ equivalent width of Wλ≥1000W_{\lambda} \ge1000 Å. Across our five redshift sub-samples, for galaxies in the range −21.6<MUV≤−19.75-21.6 < M_{\mathrm{UV}} \le-19.75, this results in EELG fractions of 0.14−0.02+0.020.14^{+0.02}_{-0.02}, 0.19−0.02+0.020.19^{+0.02}_{-0.02}, 0.31−0.03+0.030.31^{+0.03}_{-0.03}, 0.43−0.04+0.040.43^{+0.04}_{-0.04} and 0.50−0.10+0.100.50^{+0.10}_{-0.10}. These EELG fractions are shown as filled black circles in Figure 5. In Figure 5 we also show the EELG fractions for galaxies with −21.6<MUV≤−20.5-21.6 < M_{\mathrm{UV}} \le-20.5 as filled purple squares and galaxies with −20.5<MUV≤−19.75-20.5 < M_{\mathrm{UV}} \le-19.75 as filled orange squares. It can be seen that there is no indication that the evolution in the EELG fraction is a strong function of MUVM_{\mathrm{UV}}. This is reassuring, because it implies that within the MUVM_{\mathrm{UV}} range we have chosen to study, our redshift sub-samples are not strongly biased towards including only high equivalent width sources.

Figure 5. The fraction of extreme emission-line galaxies (EELGs) against redshift. The EELG fractions for the full −21.60<MUV≤−19.75-21.60 < M_{\mathrm{UV}} \leq-19.75 sample are shown as filled black circles. The EELG fractions for bright (−21.60<MUV≤−20.5-21.60 < M_{\mathrm{UV}} \leq-20.5) and faint (−20.5<MUV≤−19.75-20.5 < M_{\mathrm{UV}} \leq-19.75) sub-samples are shown as filled squares, coloured by their measure MUVM_{\mathrm{UV}}. The solid and dashed lines show quadratic fits to our full sample, and the bright and faint MUVM_{\mathrm{UV}} sub-samples, respectively. Also shown are results from Endsley et al. (2021), Boyett et al. [9], Endsley et al. (2023), Endsley et al. (2024), and Begley et al. (2025). The literature results cover a similar range of UV luminosities to our data, and have been corrected to a common EELG definition. The observed quadratic increase with redshift aligns with an increase in star-formation burstiness at high-redshift.

For comparison with our results, in Figure 5 we also show the EELG fractions measured by Endsley et al. (2021), Boyett et al. [9], Endsley et al. (2023), Endsley et al. (2024), and Begley et al. (2025), over a wide range of redshifts and at similar UV luminosities. For the studies that adopted lower EELG WλW_{\lambda} thresholds than employed here, or did not report EELG fractions (all except Begley et al. (2025)), we modelled their WλW_{\lambda} distributions as log-normal and calculated median, μ\mu, and standard deviation, σ\sigma, values. We then calculated the extreme-emitter fraction from this distribution, and repeated 1000 times. Taking the 50th and 16th/84th percentiles of the posterior distribution enabled us to obtain the median and 1σ1\sigma errors on the EELG fractions. For the Boyett et al. [9] data, we applied an additional correction factor of ×1.48\times1.48 to convert the [O iii]λ5007\lambda5007 measurements into [O iii] + Hβ measurements (Begley et al. 2025).

Combined with the existing literature results, our new set of selfconsistent measurements add to the consensus that, despite the scatter between the results of different studies, the fraction of EELGs is a strongly increasing function of redshift. The robustness of this evolution, along with the various quadratic fits shown in Figure 5, are discussed in § 5.1.

Stacked SEDs

In order to explore differences in the galaxy SEDs as a function of equivalent width, in Figure 6 we show stacked SEDs for our two highest redshift sub-samples (6.7≤z≤7.66.7 \leq z \leq7.6 and 8.0≤z≤8.58.0 \leq z \leq8.5). To create the stacked SEDs, each object’s bagpipes posterior SED distribution was sampled and the median flux at each wavelength taken to create a median spectrum. These spectra were then shifted to their rest-frame and normalised by the flux at λrest=4000\lambda_{\mathrm{rest}} = 4000 Å, after which the median was calculated to produce the final stacked SED. In Figure 6, stacked SEDs are presented for objects with Wλ≥1000W_{\lambda} \geq1000 Å, 300≤Wλ<1000300 \leq W_{\lambda} < 1000 Å and Wλ<300W_{\lambda} < 300 Å, along with the corresponding stacked photometry.

Figure 6. Stacked and normalised SEDs for the F410M sample at 6.7≤z≤7.66.7 \le z \le7.6 (left) and the F460M sample at 8.0≤z≤8.58.0 \le z \le8.5 (right). Stacked SEDs for galaxies with Wλ≥1000W_{\lambda} \ge1000 Å are shown in orange, with stacked SEDs for galaxies with 300≤Wλ<1000300 \le W_{\lambda} < 1000 Å and Wλ<300W_{\lambda} < 300 Å shown in blue and green, respectively. Stacked photometry for each sub-sample is overlaid as filled circles of the appropriate colour.

The stacked EELG SEDs, shown in orange, correspond to very young stellar populations, as indicated by the strong Balmer jump. The Balmer jump seen in the stacked bagpipes SED is not obvious in the stacked photometry due to the presence of strong line emission immediately long-ward of the jump. The strong nebular emission, both continuum and line, becomes less prominent in the lower equivalent width SED stacks, with the Wλ<300W_{\lambda} < 300 Å stack instead showing a clear Balmer break. Overall, this is consistent with the lower WλW_{\lambda} sources being systematically older and having lower specific star-formation rates.

This conclusion is supported by the results shown in Figure 7, where we present the bagpipes-derived median sSFR (10 Myr average) and mass-weighted age of the five redshift sub-samples as a function of [O III] + Hβ\beta equivalent width. Again, as one moves to lower WλW_{\lambda} the galaxies become systematically older and display lower sSFRs. We will return to examine the robustness of this conclusion in § 5.2.

Specific star-formation rate and median mass-weighted age for the galaxy sub-samples

Figure 7. The specific star-formation rate averaged over the last 10 Myr (top), and median mass-weighted age (bottom), for the F210M, F250M, F300M, F410M and F460M sub-samples (left-to-right), corresponding to z≃3.2z \simeq3.2, z≃4.0z \simeq4.0, z≃5.0z \simeq5.0, z≃7.2z \simeq7.2 and z≃8.3z \simeq8.3. The samples are divided-up by equivalent width from top-to-bottom in each panel and the quoted uncertainties are the standard errors.

W([Oiii] + H) - Stellar Mass Correlation

In Figure 8 we plot [O III]+Hβ\beta equivalent width against stellar mass for our full 2.95≤z≤8.52.95 \leq z \leq8.5 sample, with the data displayed as contours to aid visualisation. Overall there is a relatively strong correlation between the two parameters, with a Pearson correlation coefficient of rp=−0.61r_{\mathrm{p}} = -0.61 (pp-value of <0.0001< 0.0001). We do not fit a power-law to our sample because we do not expect it to be mass-complete, but our data set is consistent with, and bracketed by, the z≃2z \simeq2 and z≃8z \simeq8 power-law fits of Khostovan et al. (2016) and de Barros et al. (2019). It is clear from Figure 8 that the EELG population (grey shaded area) is confined to low stellar masses, with ≃70%\simeq70\% of EELGs having log⁡10(M⋆/M⊙)<8.0\log_{10}(M_{\star}/M_{\odot}) < 8.0.

Figure 8. WλW_{\lambda}([O III] + Hβ\beta) versus stellar mass, with contours showing a kernel density estimate visualisation of the overall data set. The power-law fits of de Barros et al. (2019) and Khostovan et al. (2016) are shown as dot-dash and dotted lines, respectively. We find that ≃70%\simeq70\% of EELGs (grey shaded region) have stellar masses of M⋆≲108M⊙M_{\star} \lesssim10^{8}M_{\odot}.

Discussion

The results presented in § 4 highlighted the redshift evolution of the [O III] + Hβ\beta equivalent width and the accompanying increase in the fraction of EELGs, along with evidence that the physical driver for this evolution is the increasing burstiness of star formation in the young Universe. In this section, we provide a more detailed discussion of each of these key issues.

Redshift evolution of [Oiii]+H equivalent width

We show the redshift evolution of the [O III]+Hβ\beta equivalent width (WλW_{\lambda}) distribution from z≃3z \simeq3 to z≃8z \simeq8 in Figure 4. Across the five redshift sub-samples, we see a clear monotonic increase in the typical WλW_{\lambda}, rising from ≃350\simeq350 Å at z≃3.2z \simeq3.2 to ≃1100\simeq1100 Å at z≃8.3z \simeq8.3.

For our F410M sub-sample at z≃7.2z \simeq7.2, we find a median equivalent width of μ=832−37+39\mu= 832^{+39}_{-37} Å (log⁡10(μ)=2.92±0.02\log_{10}(\mu) = 2.92 \pm0.02), in excellent agreement with Endsley et al. 2023 (μ=780−50+70\mu= 780^{+70}_{-50} Å, z≃7.3z \simeq7.3). Our result is also consistent with, although marginally higher than, the log⁡10(μ)=2.83−0.03+0.04\log_{10}(\mu) = 2.83^{+0.04}_{-0.03} obtained by Begley et al. (2025). The small offset can be explained by the additional SED constraints afforded by the MINERVA medium bands, as well as the SFH model adopted here allowing for larger bursts (see Section 3.2). We note that if we adopt the exact SFH prescription employed by Begley et al. (2025), our z≃7.2z \simeq7.2 result changes to log⁡10(μ)=2.78±0.02\log_{10}(\mu) = 2.78 \pm0.02.

The strong redshift evolution of the [O III]+Hβ\beta equivalent width is also clearly visible in Figure 5, where we show the fraction of sources classed as extreme emission-line galaxies (Wλ≳1000W_{\lambda} \gtrsim1000 Å; fEELGf_{\mathrm{EELG}}) for the five redshift sub-samples. We observe a strong increase in fEELGf_{\mathrm{EELG}}, rising from 0.14±0.020.14 \pm0.02 at z≃3.2z \simeq3.2 to 0.50±0.100.50 \pm0.10 at z≃8.3z \simeq8.3. We find that the evolution in fEELGf_{\mathrm{EELG}} is approximately quadratic in (1+z)(1+z), in agreement with past literature results (e.g. Endsley et al. 2021; Boyett et al. [9]; Endsley et al. 2023; Endsley et al. 2024; Begley et al. 2025). An exception to this trend is seen in Endsley et al. (2024), where the WλW_{\lambda} distribution median, and hence fEELGf_{\mathrm{EELG}}, actually decreases from z≃5.6z \simeq5.6 to z≃7.3z \simeq7.3.

In the selection of our redshift sub-samples, we implemented a SNR ≥3\geq3 threshold in the medium-band filter that sampled the [O III]+Hβ\beta complex. In principle, this could introduce an observational bias against the inclusion of the lowest WλW_{\lambda} sources. However, it can be seen from Figure 5 and Table 4 that, although there is a slight trend for the fainter MUVM_{\mathrm{UV}} sub-sample to display higher EELG fractions across our full redshift range, the difference is statistically marginal.

MUVM_{\mathrm{UV}} Subsamplef0EELGf^{\mathrm{EELG}}_0γ\gamma
−20.5<MUV≤−19.75-20.5 < M_{\mathrm{UV}} \le-19.750.0063±0.00100.0063 \pm0.00100.04±0.030.04 \pm0.03
−21.6<MUV≤−20.5-21.6 < M_{\mathrm{UV}} \le-20.50.0057±0.00140.0057 \pm0.00140.02±0.050.02 \pm0.05
Combined
−21.6<MUV≤−19.75-21.6 < M_{\mathrm{UV}} \le-19.750.0062±0.00100.0062 \pm0.00100.04±0.030.04 \pm0.03

Table 4. Parameters of our quadratic fit (fEELG=f0EELG(1+z)2+γf^{\mathrm{EELG}} = f^{\mathrm{EELG}}_0 (1 + z)^2 + \gamma) to the EELG fraction with redshift plotted in Figure 5. The parameters are listed for the two MUVM_{\mathrm{UV}} sub-samples and the combined sample.

Due to the previous lack of deep, longer-wavelength, JWST medium-band data, [O iii]+HβH\beta equivalent width estimates at z≳8z \gtrsim8 are scarce. Here, deep F460M imaging from MINERVA has enabled the most robust constraints to date of [O iii]+HβH\beta WλW_\lambda at these redshifts, where we infer log⁡10(μ)=3.05±0.03\log_{10}(\mu) = 3.05 \pm0.03 (≃1120−75+80\simeq1120^{+80}_{-75} Å) from a sample of N=49N = 49 galaxies at 8.0≤z≤8.58.0 \leq z \leq8.5. This is in good agreement with Greve et al. (2026), who obtain μ≃1260−260+330\mu\simeq1260^{+330}_{-260} Å at 9.4≤z<11.39.4 \leq z < 11.3, using MIRI/F560W imaging, although notably from a fainter sample (⟨MUV⟩≃−19.3\langle M_{\mathrm{UV}}\rangle\simeq-19.3) than our z≃8.3z \simeq8.3 sample (⟨MUV⟩≃−19.9\langle M_{\mathrm{UV}}\rangle\simeq-19.9). We also note that the z≃8.3z \simeq8.3 sample in this work is 16×16\times larger than the Greve et al. (2026) sample, further highlighting the vital statistical robustness afforded by large-scale medium band programmes such as MINERVA.

Evidence of Extremely Bursty Star-Formation

Our results provide multiple lines of evidence for increasingly bursty star-formation at z≃5−9z \simeq5 - 9. Firstly, the stacked SEDs shown in Figure 6 indicate that galaxies across the highest-redshift samples are dominated by young stellar populations, particularly amongst the most extreme emitters (Wλ≥1000W_\lambda\geq1000 Å). These sources exhibit both young mass-weighted ages and elevated sSFRs (Figure 7), consistent with intense episodes of recent star-formation following periods of reduced activity. Such rapidly varying SFHs are increasingly expected at early cosmic times, and are broadly consistent with both observational and simulation-based studies of high-redshift burstiness (e.g. (e.g. Faisst et al. 2019; Strait et al. 2023; Sun et al. 2023; Looser et al. 2024, 2025)).

Further evidence is provided by the evolution in the shape of the WλW_\lambda distributions shown in Figure 4. Galaxies observed during the peak of a star-formation burst are expected to exhibit strongly enhanced nebular emission and hence large [O iii] + HβH\beta equivalent widths. Conversely, a bursty population should also contain galaxies observed during lower-SFR phases with substantially weaker emission lines (see e.g. (see e.g. Looser et al. 2024; Faisst & Morishita 2024; Trussler et al. 2025). In this work, at z≃3.2z \simeq3.2 the inferred WλW_\lambda distribution is approximately symmetric, whereas progressively stronger asymmetry emerges towards higher redshift, becoming particularly pronounced by z≃8.3z \simeq8.3.

A similar departure from the approximately log-normal WλW_\lambda distributions seen at lower redshift was identified by Begley et al. (2025). While the typical equivalent width increases strongly with redshift, the increasingly extended tail towards low WλW_\lambda may trace galaxies observed during temporary downturns in their recent star-formation activity. Such a distribution is qualitatively expected if elevated-SFR bursts are short-lived relative to the intervening lower-SFR phases, and therefore provides independent evidence for increasingly variable SFHs towards early cosmic times. Physically, the connection between WλW_\lambda and burstiness is expected: nebular emission responds rapidly to changes in the instantaneous SFR, whereas the underlying stellar continuum traces star-formation over longer timescales. A rapid increase in SFR can therefore temporarily produce exceptionally large emission-line equivalent widths.

The stacked SED of the EELG subsample provides an additional signature of these young stellar populations, exhibiting a prominent Balmer jump at λrest≃3645\lambda_{\mathrm{rest}} \simeq3645 Å (Figure 6). This feature, which is mainly detected through SED modelling, arises from nebular free-bound emission and becomes particularly strong in young, nebular-dominated systems ((Cameron et al. 2024; Katz et al. 2025). A weaker (≲0.5×\lesssim0.5\times) discontinuity is also evident in the 300 Å ≤Wλ<1000\leq W_\lambda< 1000 Å subsample, consistent with somewhat older but still young stellar populations. In both cases, however, the Balmer jump is less apparent in the stacked photometry because of the strong contribution from emission lines red-ward of the discontinuity.

Although [O iii] + HβH\beta equivalent width provides a sensitive tracer of recent changes in star-formation activity, low WλW_\lambda does not uniquely imply a low instantaneous SFR. Recent studies have identified galaxies for which standard low-density direct-method assumptions imply extremely low metallicities, but where substantially higher electron densities alter the inferred gas conditions and metallicities ((Hayes et al. 2025; Hsiao et al. 2026; Berg et al. 2026; Arellano-Córdova et al. 2026), it should be noted that these studies focus on small numbers of extreme objects). At sufficiently high densities, collisional de-excitation can suppress [O iii] emission, producing weaker observed [O iii] + HβH\beta equivalent widths than would otherwise be expected. Some fraction of the low-WλW_\lambda population could therefore comprise actively star-forming, higher-metallicity systems in which the [O iii] emission is suppressed, rather than galaxies observed exclusively during post-burst phases or extremely metal-poor galaxies. Despite the data quality available to this study, it is clear that resolving this issue requires deep, medium resolution (R≃1000R \simeq1000), NIRSpec spectroscopy.

Associated Physical Properties

The results presented in Figure 8 show a clear anti-correlation between [O iii] +HβH\beta equivalent width and stellar mass. Although we do not fit a power-law relation to our sample due to mass-completeness concerns, our results are consistent with previous observational studies (e.g. (e.g. Reddy et al. 2018; Tang et al. 2019; Atek et al. 2022; Begley et al. 2025)) and the shape of our contours is in line with the (M⋆/M⊙)−0.48(M_\star/M_\odot)^{-0.48} scaling (dot-dash line) inferred by de Barros et al. (2019).

A correlation between [O iii] + HβH\beta equivalent width and stellar mass trend is expected if lower-mass galaxies experience increasingly variable star-formation histories, with burstiness becoming more prominent towards both lower stellar masses and higher redshifts (e.g. Looser et al. (2025)). In fact, ≃70%\simeq70\% of our EELG sample have stellar masses of M⋆≲108M⊙M_\star\lesssim10^8 M_\odot and ≃90%\simeq90\% have stellar masses of M⋆≲108.5M⊙M_\star\lesssim10^{8.5} M_\odot.

At fixed M⋆M_\star, however, we recover substantial scatter in WλW_\lambda, with several effects likely contributing to this dispersion. For example, in the absence of longer-wavelength constraints from MIRI, the stellar masses in the two highest-redshift bins remain relatively uncertain due to the potential impact of outshining (e.g. (e.g. Conroy & Gunn 2010). Some fraction of the scatter may also be intrinsic, with galaxies of similar stellar mass exhibiting substantially different equivalent widths, depending on whether they are observed during an elevated or suppressed phase of their recent SFH. Finally, covariance between the inferred stellar mass and emission-line strength within the SED fitting can further broaden the observed relation.

The combined PRIMER+MINERVA dataset also enables us to investigate the ionizing photon production efficiency, ξion\xi_{\mathrm{ion}}, across z≃3−8z \simeq3 - 8. We infer ξion\xi_{\rm ion} using the empirical relation from Begley et al. (2026), together with the conversion between [O iii] + Hβ\beta and [O iii]λ5007\lambda5007 equivalent widths described in § 4.2.

We recover median log⁡10(ξion /Hz−1erg)\log_{10}(\xi_{\rm ion} \, / \mathrm{Hz}^{-1}\mathrm{erg}) values of 25.14±0.0125.14 \pm0.01, 25.17±0.0125.17 \pm0.01, 25.23±0.0125.23 \pm0.01, 25.26±0.0125.26 \pm0.01 and 25.29±0.0325.29 \pm0.03 across the F210M, F250M, F300M, F410M, and F460M samples, respectively (Figure 9). These measurements imply only mild evolution with redshift, with dlog⁡10(ξion)/dlog⁡10(z)=0.031±0.003d\log_{10}(\xi_{\rm ion})/d\log_{10}(z) = 0.031 \pm0.003. This is shallower than the evolution reported by Simmonds et al. (2024a) and Stefanon et al. (2022), but broadly consistent with more recent measurements (e.g. Simmonds et al. 2024b; Begley et al. 2025). Indeed, Begley et al. (2025) showed that increasingly bursty SFHs towards high redshift can naturally flatten the apparent redshift evolution of ξion\xi_{\rm ion}. We emphasise, however, that ξion\xi_{\rm ion} depends on multiple correlated galaxy properties, and the weak evolution recovered here is therefore unlikely to reflect redshift alone (Begley et al. 2026).

Figure 9. The correlation between log⁡10(ξion)\log_{10}(\xi_{\rm ion}) and redshift, showing the median ionising photon production efficiency for each sub-sample. The best-fitting relation (black line) has the form: log⁡10(ξion)=(0.031±0.003)z+(25.053±0.013)\log_{10}(\xi_{\rm ion}) = (0.031 \pm0.003)z + (25.053 \pm0.013).

Conclusions

In this work, we explored the scientific gains enabled by combining JWST/NIRCam imaging from the PRIMER and MINERVA surveys in the UDS field. Together with the available HST/ACS imaging, the combined catalogue provides photometry across 19 broad-and medium-band filters, significantly improving the constraints provided by SED fitting. Using bagpipes, with non-parametric SFHs, we modelled a sample of N=1964N = 1964 galaxies with MUV≤−19.25M_{\rm UV} \leq-19.25 at z≃3−8.5z \simeq3 - 8.5 and inferred their [O iii] + Hβ\beta equivalent widths. Based on these measurements, we investigated the redshift evolution of [O iii] + Hβ\beta equivalent width and the relationship between equivalent width, age, stellar mass, and star-formation rate. Our main conclusions are as follows:

  1. Based on taking the median of 5 photometric redshift runs, we find that the combined PRIMER+MINERVA data set is capable of delivering highly accurate photometric redshift estimates. When compared to a robust sample of N=3145N = 3145 spectroscopic redshifts, drawn from ground-based surveys and the DJA, our final photometric redshifts have an outlier rate of η=4.71%\eta= 4.71\% and σ=0.015\sigma= 0.015, over the full redshift range and down to a magnitude limit of mF444W≤28m_{\mathrm{F444W}} \leq28. This represents a 42% reduction in outliers compared to photometric redshifts based on PRIMER data alone. When combined with a cut based on the standard deviation of the 5 photometric redshift runs, we showed that it is possible to reduce the outlier rate to η=2.69%\eta= 2.69\%, while culling <7%< 7\% of the sample.

  1. After exploring a variety of different star-formation history prescriptions with bagpipes, we determined that fitting the PRIMER+MINERVA photometry with BPASS SPS models and an adaptation of the Leja et al. (2019) non-parametric continuity SFH produced the most accurate recovery of the [O iii] + Hβ\beta equivalent width.

  1. As a final consistency check on our sample selection and equivalent width methodology, we presented new constraints on the bright-end of the [O iii] + Hβ\beta luminosity function at z≃7.3z \simeq7.3. Our results are in good agreement with Wold et al. (2025, 2026) and are consistent with those of Meyer et al. (2024) and Meyer et al. (2026). By combining our bright-end data with the fainter data of Wold et al. (2026), we fit a Schechter function with the following parameters: log⁡10(ϕ∗)=−4.3−0.5+0.3 Mpc−3 dex−1\log_{10}(\phi^*) = -4.3^{+0.3}_{-0.5}\,\mathrm{Mpc}^{-3}\,\mathrm{dex}^{-1}, log⁡10(L∗)=43.4−0.2+0.3 erg s−1 dex−1\log_{10}(L^*) = 43.4^{+0.3}_{-0.2}\,\mathrm{erg\,s}^{-1}\,\mathrm{dex}^{-1}, and α=−2.15−0.07+0.05\alpha= -2.15^{+0.05}_{-0.07}.

  1. Based on a consistent analysis of five redshift sub-samples spanning the redshift range z=3−8.5z = 3 - 8.5, we find that the median [O iii] + Hβ\beta equivalent width increases monotonically with redshift, rising from ≃350\simeq350 Å at z≃3z \simeq3 to ≃1100\simeq1100 Å at z≃8.3z \simeq8.3. These results are broadly consistent with, although more comprehensive than, previous measurements (e.g. Begley et al. 2025; de Barros et al. 2019; Endsley et al. 2023).

  1. In addition to the evolution in median equivalent width, we find that the equivalent width distribution becomes increasingly asymmetric with redshift. As a result, we see a roughly (1+z)2(1+z)^2 increase in the fraction of EELGs, rising from 0.14±0.020.14 \pm0.02 at z≃3.2z \simeq3.2 to 0.50±0.100.50 \pm0.10 at z≃8.3z \simeq8.3. This increase in the EELG fraction agrees well with previous constraints (e.g. Endsley et al. 2021, 2023; Begley et al. 2025; Greve et al. 2026).

  1. Stacked SEDs of the F460M and F410M samples were used to explore the differences between extreme emission-line galaxies and weaker emitters. This confirmed that EELGs are consistent with having very young stellar populations, as indicated by the presence of a strong Balmer jump. In contrast, those galaxies with low equivalent widths (i.e. Wλ<300W_{\lambda} < 300 Å), are consistent with older stellar populations, as indicated by the presence of a well developed Balmer break.

  1. The conclusion that the prevalence of EELGs at high redshift is connected with increasingly bursty star formation is further supported by considering the median mass-weighted age and specific star-formation rate of our galaxy sample over the full 10 Myr prior to observation. This confirms a consistent link between higher [O iii] + Hβ\beta equivalent width, high specific star-formation rate and young mass-weighted age. Over the full redshift range of our study, we find that EELGs are confined to a mass-weighted age of <23< 23 Myr and log⁡10(sSFR/Gyr)≥−7.1\log_{10}(\mathrm{sSFR}/\mathrm{Gyr}) \geq-7.1.

Acknowledgements

The authors would like to thank Johann Cohen-Tanugi and Olivier Ilbert for their help in improving the photometric redshift fitting through LePhare. IJBH acknowledges the support of the Science and Technology Facilities Council. RJM and TMS acknowledge support from a Science and Technology Facilities Council Grant (grant reference ST/Y001095/1). DJM and JSD acknowledge the support of the Royal Society through the award of a Royal Society Research Professorship to JSD. ACC, HL, ET, and SDS acknowledge support from a UKRI Frontier Research Grantee Grant (PI Carnall; grant reference EP/Y037065/1). FC acknowledges support from a UKRI Frontier Research Guarantee Grant (PI Cullen; grant reference EP/X021025/1). This paper is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The 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. This work made use of v2.2.1 of the Binary Population and Spectral Synthesis (BPASS) models as described in Eldridge et al. (2017) and Stanway & Eldridge (2018). This research has made use of the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/) supported from the Spanish MINECO through grant AYA2017-84089 (Rodrigo et al. 2012; Rodrigo & Solano 2020). This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

Data Availability

At the time of publication, all data relating to the images used in this work are available from the public Barbara A. Mikulski Archive for Space Telescopes (MAST).

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

References

  1. [9]Boyett K. N. K., Stark D. P., Bunker A. J., Tang M., Maseda M. V., 2022, MNRAS, 513, 4451

Paper details

vlnc:W_ypzjer9c
v1 · latestHistory (1)

published · Fresh — added to the Hub within the last 30 days

added ·

Usage

views · 0

downloads · 0

A view counts once per browser per day, a download each time the PDF is requested.

CC BY 4.0

Original source

  • Physics and Astronomy

Not recorded for this paper · View statement →

Looking for related papers…
Contents