WorksPhysics and Astronomy
A combined PRIMER+MINERVA study of the evolution of \mathbf{[O\mbox{III}]+Hβ} equivalent width from to
Abstract
We present an analysis of the evolution of equivalent width () over the redshift range to , 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 at , a reduction compared with PRIMER imaging alone. We also demonstrate that the combined data set facilitates the accurate recovery of equivalent widths in redshift bins centred on ; corresponding to the F210M, F250M, F300M, F410M and F460M filters. Based on a final sample of galaxies with , we first provide new constraints on the bright end of the luminosity function (, ), before proceeding to investigate the redshift evolution of equivalent width. Using SED modelling based on non-parametric continuity star-formation histories, we find that the median equivalent width monotonically increases with redshift, rising from at to at . In tandem, the distribution of is observed to become increasingly asymmetric with redshift, producing a fraction of extreme emission-line galaxies (EELGs, Å) that reaches by . 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] Å reflects the production rate of massive O- and B-type stars on a timescale of 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 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 (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 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 (), the equivalent width () of [O iii] is known to be highly correlated with ionizing photon production efficiency (; Tang et al. 2019). Recent NIRSpec studies have found a similar − relation at , together with a strong correlation between [O iii] and , 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 for large, statistically complete samples. Recent photometric studies have demonstrated that extreme emission-line galaxies (EELGs: sources with [O iii] + H Å) were much more common in the early Universe, and that the typical [O iii] + H increases from Å at (Boyett et al. [9]) to Å at (e.g. Begley et al. 2025; Endsley et al. 2025).
Until recently, photometric studies of [O iii] + H 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 –.
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 equivalent width distribution and the fraction of EELGs, within a –complete galaxy sample spanning the redshift range . 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 luminosity function at . 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 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 [O iii] + H luminosity function. We present our [O iii] + H 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: , , , 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 arcmin 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 arcmin, while an area of arcmin 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 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 and -diameter circular apertures, and corrected for galactic extinction using the Gordon et al. (2023) extinction curve and , 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 -diameter circular apertures (corrected to total) are presented in Table 1.
| Filter | Filter | ||
| F435W† | 27.39 | F140M | 27.10 |
| F606W† | 28.08 | F162M | 27.16 |
| F814W† | 28.09 | F182M | 27.43 |
| F090W | 27.71 | F210M | 27.25 |
| F115W | 27.75 | F250M | 27.33 |
| F150W | 27.95 | F300M | 27.69 |
| F200W | 28.12 | F360M | 27.69 |
| F277W | 28.38 | F460M | 26.71 |
| F356W | 28.37 | ||
| F410M | 27.72 | ||
| F444W | 27.83 |
Table 1. The global depths in all 19 filters. The depths are calculated in a -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 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 and -diameter aperture photometry. For the vast majority of sources (87% of the catalogue), the final photometric redshifts were based on the higher SNR -diameter photometry. However, for bright objects () we adopted the photometric redshifts derived from the -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 high-quality spectroscopic redshifts within the 234 arcmin 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: , the catastrophic outlier fraction, (where ), the standard deviation, and the bias, . The full set of statistics are reported in Table 2 for those objects that have SNR 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 and an outlier fraction of 8.1%. However, when the medium-band data from MINERVA are included, the overall accuracy improves to 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 | |||||
| Photometry | N | ||||
| PRIMER | 2599 | 8.12% | 0.033 | 0.008 | |
| PRIMER+MINERVA | 3145 | 4.71% | 0.015 | 0.001 | |
| PRIMER+MINERVA | 2972 | 2.69% | 0.014 | 0.001 | |
| (b) Ground | |||||
| Photometry | N | ||||
| PRIMER | 654 | 2.91% | 0.023 | 0.004 | |
| PRIMER+MINERVA | 654 | 0.92% | 0.012 | -0.002 | |
| PRIMER+MINERVA | 650 | 0.92% | 0.012 | -0.001 | |
| (c) DJA | |||||
| Photometry | N | ||||
| PRIMER | 1945 | 9.87% | 0.039 | 0.010 | |
| PRIMER+MINERVA | 2491 | 5.70% | 0.016 | 0.002 | |
| PRIMER+MINERVA | 2322 | 3.19% | 0.015 | 0.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 () 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 (). 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 , resulting in an outlier rate of only , down to a magnitude limit of . This criterion is therefore capable of reducing the outlier fraction by a further 43%, while excluding 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 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 (), 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 () when the MINERVA photometry is added to the PRIMER photometry. We see a further 43% reduction in the outlier rate when including the 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 cut). Moreover, isolating only objects from the RUBIES survey (PID-4233; PIs: de Graaff & Brammer), we measure ; compared to the obtained by de Graaff et al. (2025) using only PRIMER data.
[Oiii] + H 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 equivalent width out to , while the F410M PRIMER imaging provides an opportunity to benchmark our results against previous studies at (e.g. Meyer et al. 2024; Begley et al. 2025).
To determine the [O iii] + H equivalent widths () 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- 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 equivalent widths for 115 galaxies drawn from our full sample that have publicly available NIRSpec PRISM spectra on the DJA. These sources have SNR in the F356W filter and , 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 (). To ensure accurate spectroscopic measurements of the [O iii] + H equivalent widths, these spectra were chosen to have SNR over the [O iii] + H complex, together with significant detections of the adjacent continuum: median SNR/pix 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 prior scale from 0.3 to 1, and decreased the prior d.f. from 2 to 1.5. These hyper-parameters dictate and , respectively, in the Student’s -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.
| Parameter | Range | Prior | Hyper-parameters |
| Redshift | Gaussian | ||
| Total stellar mass formed | |||
| Stellar Metallicity | Uniform | ||
| SFR Change () | Student-t | prior scale = 1, prior d.f. = 1.5 | |
| V-band attenuation | Uniform | ||
| Nebular ionisation parameter | Uniform |
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 equivalent widths. We take the 50th percentile of the resulting distribution to be our measurement, and the 16th and 84th percentiles as our 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 Å. To account for the minor discrepancy between our photometric measurements and their spectroscopic counterparts, we fit the data with a linear relationship, finding . All future measurements are corrected using this relation.
!Photometric versus spectroscopic measurements of rest-frame [O III] + Hβ equivalent width bbox=0.07,0.07,0.47,0.53
Figure 2. Photometric versus spectroscopic measurements of rest-frame [O III] + equivalent width for our sub-sample of 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 against the spectroscopic . We find a relationship of .
Sample Selection
Our initial sample of galaxies comprised objects with SNR in the F356W filter and a photometric redshift in the ranges , , , , or , corresponding to the [O III] + 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 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] + equivalent widths in the redshift range . This is because the F360M filter sits within the broad F356W filter, meaning that our parent photometry catalogue is biased towards high [O III] + equivalent width objects.
To exclude potential active galactic nuclei (AGN) we removed sources detected in X-rays within 0.5″ (using Kovcevski et al. 2018), and 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 objects, over the redshift range . However, in order to mitigate the evolving limit, we restrict our final sample to those objects with , 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 objects, with in the bin, in the bin, in the bin, in the bin and in the bin.
[Oiii] + H Luminosity Function
As a final consistency check on our sample selection and SED-fitting methodology, we computed the [O III] + luminosity function (LF) at by combining our F410M and F460M sub-samples.
LF estimate is based on the 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 , together with recent literature results at 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 LFs by Wold et al. (2026) is also shown as a pink dotted line. Where necessary, we have corrected [O iii]5007 luminosities to [O iii]+H luminosities using a factor of 1.48 (Begley et al. 2025). Our new LF estimate is plotted at , which corresponds to our conservative completeness limit, based on the 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] + luminosity function at redshift (filled blue circles), based on combining our and sub-samples. We also show two LF estimates at from Meyer et al. (2024) and Meyer et al. (2026), an LF estimate at by Wold et al. (2025) and the two highest luminosity data points from the LF estimate of Korber et al. (2026). Where necessary we have corrected from [O III]5007 to [O III] + luminosities. The dashed line shows the Schechter function fit at 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 LF of Wold et al. (2026) and also follows cleanly on from the 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 , 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 . Our best-fitting Schechter function is shown as the solid black line in Figure 3 and has the following parameters: , , and .
It can be seen that our new Schechter function is 0.5 dex higher than the Wold et al. (2026) fit at . 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 galaxies. In this section we investigate the evolving [O iii] + H equivalent width distributions and underlying physical properties of our galaxy samples at , , , , and .
As [O iii] + H 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 equivalent widths for all redshift samples in Figure 4. For each object in the samples, a random 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 , from which the 16, 50, and 84 percentiles were taken. The median values of for the F210M, F250M, F300M, F410M, and F460M samples are listed in the Figure 4 legend. It is clear that the median value of is monotonically increasing with redshift. Moreover, it can also be seen that the 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 equivalent width distributions for the F210M (purple, ), F250M (red, ), F300M (green, ), F410M (orange, ) and F460M (blue, ) sub-samples. The shaded regions surrounding the lines are the 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 Å. Across our five redshift sub-samples, for galaxies in the range , this results in EELG fractions of , , , and . These EELG fractions are shown as filled black circles in Figure 5. In Figure 5 we also show the EELG fractions for galaxies with as filled purple squares and galaxies with 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 . This is reassuring, because it implies that within the 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 sample are shown as filled black circles. The EELG fractions for bright () and faint () sub-samples are shown as filled squares, coloured by their measure . The solid and dashed lines show quadratic fits to our full sample, and the bright and faint 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 thresholds than employed here, or did not report EELG fractions (all except Begley et al. (2025)), we modelled their distributions as log-normal and calculated median, , and standard deviation, , 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 errors on the EELG fractions. For the Boyett et al. [9] data, we applied an additional correction factor of to convert the [O iii] 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 ( and ). 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 Å, after which the median was calculated to produce the final stacked SED. In Figure 6, stacked SEDs are presented for objects with Å, Å and Å, along with the corresponding stacked photometry.
Figure 6. Stacked and normalised SEDs for the F410M sample at (left) and the F460M sample at (right). Stacked SEDs for galaxies with Å are shown in orange, with stacked SEDs for galaxies with Å and Å 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 Å stack instead showing a clear Balmer break. Overall, this is consistent with the lower 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 equivalent width. Again, as one moves to lower the galaxies become systematically older and display lower sSFRs. We will return to examine the robustness of this conclusion in § 5.2.

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 , , , and . 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 equivalent width against stellar mass for our full 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 (-value of ). 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 and 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 of EELGs having .
Figure 8. ([O III] + H) 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 of EELGs (grey shaded region) have stellar masses of .
Discussion
The results presented in § 4 highlighted the redshift evolution of the [O III] + H 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 equivalent width () distribution from to in Figure 4. Across the five redshift sub-samples, we see a clear monotonic increase in the typical , rising from Å at to Å at .
For our F410M sub-sample at , we find a median equivalent width of Å (), in excellent agreement with Endsley et al. 2023 ( Å, ). Our result is also consistent with, although marginally higher than, the 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 result changes to .
The strong redshift evolution of the [O III]+H equivalent width is also clearly visible in Figure 5, where we show the fraction of sources classed as extreme emission-line galaxies ( Å; ) for the five redshift sub-samples. We observe a strong increase in , rising from at to at . We find that the evolution in is approximately quadratic in , 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 distribution median, and hence , actually decreases from to .
In the selection of our redshift sub-samples, we implemented a SNR threshold in the medium-band filter that sampled the [O III]+H complex. In principle, this could introduce an observational bias against the inclusion of the lowest sources. However, it can be seen from Figure 5 and Table 4 that, although there is a slight trend for the fainter sub-sample to display higher EELG fractions across our full redshift range, the difference is statistically marginal.
| Subsample | ||
| Combined | ||
Table 4. Parameters of our quadratic fit () to the EELG fraction with redshift plotted in Figure 5. The parameters are listed for the two sub-samples and the combined sample.
Due to the previous lack of deep, longer-wavelength, JWST medium-band data, [O iii]+ equivalent width estimates at are scarce. Here, deep F460M imaging from MINERVA has enabled the most robust constraints to date of [O iii]+ at these redshifts, where we infer ( Å) from a sample of galaxies at . This is in good agreement with Greve et al. (2026), who obtain Å at , using MIRI/F560W imaging, although notably from a fainter sample () than our sample (). We also note that the sample in this work is 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 . 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 ( Å). 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 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] + 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 the inferred distribution is approximately symmetric, whereas progressively stronger asymmetry emerges towards higher redshift, becoming particularly pronounced by .
A similar departure from the approximately log-normal 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 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 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 Å (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 () discontinuity is also evident in the 300 Å Å 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] + equivalent width provides a sensitive tracer of recent changes in star-formation activity, low 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] + equivalent widths than would otherwise be expected. Some fraction of the low- 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 (), NIRSpec spectroscopy.
Associated Physical Properties
The results presented in Figure 8 show a clear anti-correlation between [O iii] + 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 scaling (dot-dash line) inferred by de Barros et al. (2019).
A correlation between [O iii] + 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, of our EELG sample have stellar masses of and have stellar masses of .
At fixed , however, we recover substantial scatter in , 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, , across . We infer using the empirical relation from Begley et al. (2026), together with the conversion between [O iii] + H and [O iii] equivalent widths described in § 4.2.
We recover median values of , , , and across the F210M, F250M, F300M, F410M, and F460M samples, respectively (Figure 9). These measurements imply only mild evolution with redshift, with . 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 . We emphasise, however, that 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 and redshift, showing the median ionising photon production efficiency for each sub-sample. The best-fitting relation (black line) has the form: .
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 galaxies with at and inferred their [O iii] + H equivalent widths. Based on these measurements, we investigated the redshift evolution of [O iii] + H equivalent width and the relationship between equivalent width, age, stellar mass, and star-formation rate. Our main conclusions are as follows:
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 spectroscopic redshifts, drawn from ground-based surveys and the DJA, our final photometric redshifts have an outlier rate of and , over the full redshift range and down to a magnitude limit of . 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 , while culling of the sample.
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 equivalent width.
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 luminosity function at . 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: , , and .
Based on a consistent analysis of five redshift sub-samples spanning the redshift range , we find that the median [O iii] + H equivalent width increases monotonically with redshift, rising from Å at to Å at . 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).
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 increase in the fraction of EELGs, rising from at to at . 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).
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. Å), are consistent with older stellar populations, as indicated by the presence of a well developed Balmer break.
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 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 Myr and .
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
- [9]Boyett K. N. K., Stark D. P., Bunker A. J., Tang M., Maseda M. V., 2022, MNRAS, 513, 4451





