{'fraction_non_alphanumeric': 0.03496597722832311, 'fraction_numerical': 0.02718453142895641, 'mean_word_length': 4.455375191424196, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The PHLEK Survey: A New Determination of the Primordial Helium Abundance

We present Keck NIRSPEC and Keck NIRES spectroscopy of sixteen metal-poor galaxies that have pre-existing optical observations. The near-infrared (NIR) spectroscopy specifically targets the HeI 10830A emission line, due to its sensitivity to the physical conditions of the gas in HII regions. We use these NIR observations, combined with optical spectroscopy, to determine the helium abundance of sixteen galaxies across a metallicity range 12 + log (O/H) = 7.13 - 8.00. This data set is combined with two other samples where metallicity and helium abundance measurements can be secured: star-forming galaxies selected from the Sloan Digital Sky Survey spectroscopic database and existing low-metallicity systems in the literature. We calculate a linear fit to these measurements, accounting for intrinsic scatter, and report a new determination of the primordial helium number abundance, y_P = 0.0805+/-0.0017, which corresponds to a primordial helium mass fraction Y_P = 0.2436+/-0.0040. Using our determination of the primordial helium abundance in combination with the latest primordial deuterium measurement, (D/H)_P x 10^5 =2.527+/-0.030, we place a bound on the baryon density Omega_bh^2 = 0.0215+/-0.0005 and the effective number of neutrino species Neff = 2.85+/-0.28. These values are in 1.3-sigma agreement with those deduced from the Planck satellite observations of the temperature fluctuations imprinted on the Cosmic Microwave Background.


INTRODUCTION
The abundances of the light elements that were produced during Big Bang Nucleosynthesis (BBN) chiefly depend on: (1) the ratio of the baryon density to photon density, η 10 ≡ 10 10 (n B /n γ ) and (2) the expansion rate of the Universe (Hoyle & Tayler 1964;Peebles 1966). Baryonic matter in the Universe just prior to the onset of BBN mostly consisted of free neutrons and protons, which rapidly fused to form deuterium, and subsequently, other light elements. The freeze-out abundances of deuterium and the isotopes of helium and lithium, depend on a competition between the expansion rate of the Universe and the nuclear and weak interaction rates that govern the synthesis of the light elements (see the recent BBN reviews by Steigman 2007Steigman , 2012Cyburt et al. 2016;Pitrou et al. 2018).
The universal baryon density, Ω b h 2 η 10 /273.9 (Steigman 2006) is determined to ∼ 1 per cent precision via the temperature fluctuations of the Cosmic Microwave Background (CMB). The most recent determination of the baryon density inferred from the CMB is Ω b h 2 = 0.02236 ± 0.00016 (68 per cent confidence limits (CL) of the TT+TE,EE+lowE parameter estimation; see Table 2, Column 4 of Planck Collaboration et al. 2018). The expansion rate of the Universe is determined by the total energy density of the Universe. At the time of BBN, the total energy density was dominated by massless and relativistic particles, including photons, electrons, and the three Standard Model neutrinos (Steigman 2012;Mathews et al. 2017). The total radiation energy density is parameterized by the effective number of neutrino species, N eff = 3.046 + ∆N ν (equivalent to N ν = 3 + ∆N ν ). For the Standard Model of particle physics and cosmology, ∆N ν = 0. In the framework of the Standard Model in combination with the Planck measurement of Ω b h 2 , a mean neutron lifetime τ n , and cross-sections for the relevant reaction rates, the primordial element yields can be predicted to a precision of less than two per cent (Pitrou et al. 2018).
Similarly, observational measurements of the light element abundances in near-pristine environments provide an opportunity to infer the constituents of the early Universe. These observational measures of the primordial abundances offer an important test of standard Big Bang nucleosynthesis (SBBN); deviations from the SBBN light element abundances would indicate new physics in the early Universe. For example, if ∆N ν = 0, there may be a previously unrecognized particle that changes the total energy density of the Universe and thus the expansion rate of the early Universe (e.g., Di Valentino et al. 2013). To assess this possibility, reliable and precise observational measurements of the light element abundances must be made in order to firmly conclude the existence of physics beyond the Standard Model.
The light element nuclides deuterium D/H, helium-3 ( 3 He), helium-4 ( 4 He), and lithium-7 ( 7 Li) are made in as-trophysically measurable quantities, and have therefore been the targets of historic and current primordial abundance measurements. While all the primordial abundances depend on both the baryon density and the expansion rate of the Universe at the time of BBN, (D/H) and 7 Li are most sensitive to the baryon abundance whereas 4 He is primarily sensitive to the expansion rate of the Universe (see Figure 7 of Cyburt et al. 2016). 3 He is less sensitive to both the baryon density and the expansion rate than its peer primordial elements but provides orthogonal contours to (D/H) in the ∆N ν − Ω b h 2 plane (Cooke 2015). A 3 He abundance has been observed and measured in H II regions and planetary nebulae in the Milky Way, but these measures likely do not reflect the primordial 3 He composition, due to contamination by the complicated post-BBN production of 3 He (Olive et al. 1995;Vangioni-Flam et al. 2003).The primordial abundance of 7 Li can be inferred from the atmospheres of the most metal-poor dwarf stars in our Galaxy. The latest determinations (Aoki et al. 2009;Meléndez et al. 2010;Sbordone et al. 2010;Spite et al. 2015) are, however, in significant (∼ 6σ) disagreement with the SBBN value (Cyburt et al. 2008;Fields 2011), and has been famously dubbed the "lithium problem".
The primordial D/H ratio, (D/H) P , offers a sensitive probe of the baryon density and has a simple post-BBN chemical evolution. There are no pathways that net produce deuterium, so its abundance should decrease monotonically with increasing metallicity. Currently, the best environments to measure the primordial D/H ratio are high-redshift, nearpristine quasar absorption systems, where the current determination is at the 1 per cent level, (D/H) P = (2.527 ± 0.030) × 10 −5 , in agreement with SBBN .
The mass fraction of 4 He offers a sensitive test of physics beyond the Standard Model (Yang et al. 1979(Yang et al. , 1984Olive et al. 1981) due to its strong dependence on the effective number of neutrino species. Attempts to measure the primordial 4 He abundance, commonly denoted in the literature by the helium mass fraction, Y P , have most commonly utilized emission line observations of H II regions in low-metallicity dwarf galaxies, defined to have gas phase oxygen abundances less than a tenth solar metallicity, 12 + log 10 (O/H) ≤ 7.69. This method has shown the most promise to reach a ∼ 1 per cent inference on the helium abundance. Searle & Sargent (1972) presented an abundance analysis of the extragalactic H II regions I Zwicky 18 (I Zw18) and II Zwicky 40 and first suggested that metal-poor systems such as these would be crucial to pin down the primordial helium abundance. Finding new, metal-poor H II regions has historically been difficult, however. While all-sky surveys such as the Sloan Digital Sky Survey (SDSS) have provided a means to identify new, low-metallicity systems Izotov et al. 2013;Guseva et al. 2017), the number of metal-poor systems expected from the luminosity function greatly outnumbers the number of known metal-poor systems (Sánchez Almeida et al. 2017). It has been suggested that the most metal-poor systems tend to elude spectroscopic surveys, possibly due to their intrinsically low surface brightnesses as predicted by the luminosity-metallicity relation (James et al. 2015). Consistent with this line of reasoning, with the exception of the extremely metal-poor but more luminous systems such as I Zw 18 (Searle & Sargent 1972;Skillman & Kennicutt 1993) and SBS 0335-052 (Izotov et al. 1990), discoveries of new systems that push on the lowest-metallicity regime have been rare. Yet systems similar to these, i.e., at the hundredth solar metallicity level, are necessary for a precise extrapolation to the primordial helium value. Some recent exceptions include Leo P (Giovanelli et al. 2013;Skillman et al. 2013), AGC 198691 (Hirschauer et al. 2016, which were both initially found as H I gas rich regions in the Arecibo Legacy Fast ALFA Survey (Giovanelli et al. 2005), the Little Cub (Hsyu et al. 2017), J0811+4730 (Izotov et al. 2018), and HSC J1631+4426 (Kojima et al. 2019). Many of the latest efforts to significantly boost the number of low-metallicity H II regions have focused on using photometry to identify candidate systems, followed by spectroscopic confirmation combined with a direct measurement of the metallicity of the system. This method has yielded successful results, with 20 -60% of observed systems in these dedicated searches falling in the low-metallicity regime (James et al. 2015(James et al. , 2017Yang et al. 2017;Hsyu et al. 2018;Senchyna & Stark 2019).
Extracting a measure of the helium abundance of these near-pristine galaxies has its challenges. H II region modelling is believed to suffer from systematic uncertainties (for an incomplete list, see  and degeneracies among the model parameters, particularly between the electron density and temperature. This can lead to biases in the determination of the helium abundance (see Figure 3 of Aver et al. 2015). To help alleviate these biases, Izotov et al. (2014) included the near-infrared (NIR) He I λ10830Å line in their helium abundance analysis. The He I λ10830 line is very sensitive to the electron density, and helps to break the temperature-density degeneracy. Aver et al. (2015) confirmed the importance of He I λ10830 as an excellent density diagnostic -the addition of the He I λ10830 line to their analysis of 11 systems reduced the 1σ confidence interval on the electron density by 60 per cent. This reduction of the error on the electron density led to a reduction of the error on the helium abundance of each H II region ranging from 10-80 per cent.
However, these two works, which have systems in common in their analyses, report primordial helium abundances in mutual disagreement with one another. Izotov et al. (2014) reports Y P = 0.2551 ± 0.0022, which is higher than the SBBN predicted value, while Aver et al. (2015) finds Y P = 0.2449 ± 0.0040, consistent with the SBBN value of Y P = 0.24709 ± 0.00017 (Pitrou et al. 2018). Several other groups have recently reported competitive measurements of the primordial helium abundance in good agreement with the Aver et al. (2015) result, using a range of techniques. For example, Fernández et al. (2018) use sulphur (S) instead of oxygen (O) as a metallicity tracer, and find that the scatter in the Y P vs S/H plane is reduced compared with Y P vs O/H. These authors later employ probabilistic programming methods and find good agreement with their previous work (Fernández et al. 2019). Other groups have instead focused on modeling a small number of well-selected H II regions to infer the primordial value (Peimbert et al. 2016;Valerdi et al. 2019). It is perhaps promising that the different data sets used and the different modelling approaches employed yields mostly consistent results (with the exception of the value reported by Izotov et al. 2014). However, it is still necessary to take caution of confirmation bias (see e.g. Figure 8 of Steigman 2012), and understand why models are currently unable to simultaneously reproduce all of the observed H I and He I emission lines of some H II regions.
Motivated by the dearth of metal-poor systems that push on the lowest-metallicity regime and the need for more highquality, complementary optical and NIR spectra of external galaxies, we conducted a dedicated survey to identify new, metal-poor systems via SDSS photometry (Hsyu et al. 2018). Our follow-up spectroscopic survey of 94 objects found almost half of them to be in the low-metallicity regime, and our findings included one of the lowest-metallicity systems currently known, the Little Cub (Hsyu et al. 2017). After initial metallicity estimates, we obtained spectroscopy of a subset of the most promising systems, with a focus on obtaining high signal-to-noise (S/N) optical and NIR spectra. In this paper, we use this new sample, along with some previous systems in the literature, to report a new determination of the primordial helium abundance.
In Section 2, we describe the details of the full sample of galaxies that we use in this paper. This includes our own sample of new complementary optical and NIR data, for which we also include details of the observations, data reduction, and integrated emission line flux measurements. We supplement our data set with galaxies from the SDSS spectroscopic database and the HeBCD sample from Izotov & Thuan (2004);. The components of our model and the subsequent MCMC analysis used to solve for the best fit parameters of our H II regions are described in Section 3. In Section 4, we assess the potential systematics and select the most reliable set of H II regions to use in our determination of the primordial helium abundance. We discuss the implications of our work and consider future improvements to primordial helium work, both in observations of new systems and in model enhancements, in Section 5. Finally, we summarize our main conclusions in Section 6.

DATA COMPILATION AND PREPARATION
A well constrained measurement of the primordial helium abundance requires accurate measurements of the oxygen and helium abundance from a sizeable sample of galaxies that span a range of metallicities. In this section, we describe the observations of our galaxy sample, which populates the lowest-metallicity regime. Throughout this paper, we refer to our galaxy sample as the Primordial Helium Legacy Experiment with Keck (PHLEK) sample. We supplement our PHLEK sample with existing spectra from SDSS and the Izotov & Thuan (2004);  HeBCD data set. This combined sample provides a set of measurements that cover a broad range of metallicity. We note that the three data sets that make up our final, full sample of galaxies are thus likely heterogeneous data sets, and the degree of our involvement in processing each sample (e.g., converting the twodimensional, raw data into integrated emission line fluxes, varies).

Keck Observations
The primary goal of our observational program is to increase the sample size of very metal-poor galaxies where reliable oxygen and helium abundances can be determined. To this end, we acquired optical and near-infrared spectra of metal-poor H II regions in nearby dwarf galaxies using Keck Observatory, requiring that the spectra have confident detections of: • the temperature sensitive [O III] λ4363Å line for a direct measurement of the oxygen abundance • at least five optical He I emission lines to reliably determine the physical state of the H II regions, including: He I λ3889Å , λ4026Å , λ4471Å , λ5015Å , λ5876Å , λ6678Å , and λ7065Å • the NIR He I λ10830Å line, whose emissivity is the most sensitive He I emission line to the density of the gas, relative to Pγ λ10940Å In addition to these emission lines, we also detect in our spectra the To ensure that we observe the same region of each system either on multiple nights or on different instruments, we acquire each target by first centering on a bright nearby star, then applying an offset to the target based on SDSS astrometry. Additionally, we requested that our optical and nearinfrared nights be allocated within a week of one another such that our complementary observations for a given target be at similar airmass and parallactic angle. For the observations, we matched the slit widths of different instruments as best as possible. Spectroscopic observations of our metal-poor galaxy sample took place during semesters 2015B, 2016A, and 2018A (program IDs: U052LA/U052NI, U091LA/U091NS, U172).

Optical Spectroscopy
Optical spectroscopic observations of 32 metal-poor systems were made using the Low Resolution Imaging Spectrometer (LRIS) with the atmospheric dispersion corrector (ADC) at the W.M. Keck Observatory. LRIS has separate blue and red channels. On the blue side, our setup utilized the 600/4000 grism, which has an unbinned dispersion of 0.63 Å pix −1 . On the red side, we used the 600/7500 grating, which has an unbinned dispersion of 0.8 Å pix −1 . Using this instrument setup, the D560 dichroic, and a long slit, the full wavelength coverage achieved is ∼3200-8600 Å, with the separate blue and red channels covering ∼3200-5600 Å and ∼5400-8600 Å, respectively. We use 2×2 binning during readout. The blue and red channels have nominal FWHM resolutions of 2.6 Å and 3.1 Å for our adopted 0.70 slit. While the separate blue and red arms overlap in wavelength coverage, we take caution about the accuracy of the measurements here, as data near the region of overlap is compromised by the dichroic.
Our spectra were taken with a 175 × 0.70 slit, oriented at the parallatic angle. Our total exposure times range from 3 × 1200 s to 3 × 1800 s. We obtained bias frames, arc frames, and dome flats at the beginning of the night. For wavelength calibration on the blue side, we observed Hg, Cd, and Zn arc lamps; on the red side, we observed Ne, Ar, and Kr arc lamps. Photometric standard stars G191B2B, BD+284211, Feige 34, Feige 66, Feige 110, and/or HZ44 were observed at the start and end of each night for flux calibration. Excluding five previously unreported systems which are presented here, our observed and derived physical properties of the galaxies based on Keck+LRIS spectra are reported in Hsyu et al. (2018).

Near-Infrared Spectroscopy
We acquired complementary NIR observations for 16 of our 32 galaxies with optical spectroscopy. NIR observations were made using NIRSPEC in semesters 2015B and 2016A and the Near-Infrared Echellette Spectrometer (NIRES) in 2018A. Our NIRSPEC observations were done in low resolution mode using the NIRSPEC-1 filter, which offers a wavelength coverage of ∼9470-12100Å. NIRES covers wavelengths ∼9400-24500Å across five orders, with a gap between 18500-18800Å, though this wavelength gap does not affect our observation goals.
Our NIRSPEC observations were made using the 42× 0.72 slit to best match the slit width of our LRIS observations. The NIRES slit is fixed at 18× 0.55 . We observed all targets with the slit oriented at the parallactic angle. All NIR observations were made using an ABBA nod pattern for exposure times of 8 × 250 s to 8 × 360 s each. We obtained dome flats at the beginning of each night. An A0V calibration star near each of our science targets was observed following each observation for flux calibration.

Data Reduction
For optical LRIS observations, the two-dimensional raw images were individually bias subtracted, flat-field corrected, cleaned for cosmic rays, sky-subtracted, extracted, wavelength calibrated, and flux calibrated, all using PYPEIT (previously PYPIT), a Python based spectroscopic data reduction package. 1 We used a boxcar extraction technique to extract a single one-dimensional (1-D) spectrum of each object. 2 Multiple observations of the same target were coadded by weighting each exposure by the inverse variance at each pixel.
For our NIR data, PYPEIT combines a single set of ABBA observations during the reduction as A+A − (B+B), yielding an extracted 1-D spectrum at nod location A. Similarly, the frames are combined as B+B − (A+A) for a spectrum at nod location B. PYPEIT first flat-fields the individual frames, then combines and subtracts relevant frames, which removes the bias level and performs a first order sky subtraction. PYPEIT wavelength calibrates using the OH sky lines. Flux calibrations for NIR observations are performed separately from the automated reduction routine using the pypeit_flux_spec script. Our NIR observations of each target were acquired in two sets of ABBA observations, such that the final coadded spectrum consists of four 1-D extracted spectra, comprised of two spectra of A+A − (B+B) and two spectra of B+B − (A+A). We show an example of our reduced and coadded NIR spectra in Figure 1.

SDSS Sample
In addition to our new sample of metal-poor systems observed at Keck, we also use the SDSS spectroscopic database to identify additional emission-line galaxies that can be included in our primordial helium work. The SDSS sample complements our PHLEK sample described in Section 2.1 by providing a sample of higher metallicity galaxies. It also offers the potential to significantly increase the number of systems available for helium abundance analyses.
To take advantage of this database, we queried the SDSS specObj database for systems that are suitable to our analysis. Our query requires that the systems are: (1) classified as starburst galaxies and (2) within a redshift range of 0.02 < z < 0.15, such that the [O II] doublet and He I λ7065 lines, necessary for a metallicity and helium abundance, fall on the detector. Our SQL query can be found in Appendix A.
For the resulting galaxies, we calculate the emission line fluxes using the method described in Section 2.4 and filtered the systems to keep those with confident detections of: (1) the temperature sensitive [O III] λ4363Å line for a direct metallicity, and (2) multiple He I lines, to measure the helium abundance. We impose these criteria using the following S/N cuts, where S/N is defined to be the measured F(λ)/σ(F(λ)): Of these He I lines, the He I λ5876 line is typically the most significantly detected. We therefore require the strongest S/N condition on this line to ensure a confident detection of the weaker He I lines.
These steps filtered the SDSS spectroscopic database down to 1053 candidate systems to be included in our analysis. For reference, the peak of the metallicity distribution of this SDSS sample is (O/H) × 10 5 = 13.24, whereas the peak of the metallicity distribution of our PHLEK galaxies, including the systems presented in Hsyu et al. (2018) and here, is (O/H) × 10 5 = 4.82. These values correspond to 12 + log 10 (O/H) values of 8.12 and 7.68, respectively.

Emission Line Flux Measurements
For the Keck and SDSS samples, we calculate the integrated emission line fluxes by summing the total flux above the continuum level at each emission line, where the continuum level and its error are modelled using the Absorption LIne Software (ALIS, see Cooke et al. 2014 for a more detailed description of the software). 3 ALIS simultaneously fits the emission line profile using a Gaussian model and the surrounding continuum using a 1-or 2-D Legendre polynomial and determines the best fit parameters of the Gaussian and continuum model using a χ 2 minimization approach. Systems with high emission line fluxes, however, are not wellrepresented by a single Gaussian. We therefore adopt the continuum model and its associated error from the ALIS output, and use this to inform our calculation of the total flux  Table 2, as collected using NIRSPEC at Keck Observatory. The error spectra are shown in red. Only a small window of NIRSPEC's entire ∼9470-12100Å wavelength range is shown in these panels to best highlight the relevant emission lines of interest, He I λ10830 and Pγ λ10940, which are marked in the left panel.
above the continuum level. The width of the emission line included in the integrated flux is set to be ±5 pixels around the closest pixel to the redshifted central wavelength of the emission line. Two exceptions are the [O II] doublet, which has a width of ±7 pixels to encompass the full width of the blended doublet, and He I λ5015 where we take only 3 pixels (∼1.9Å) to the left of the central wavelength to avoid contamination from the [O III] λ5007 line (we still use 5 pixels to the right). We map the pixels to an array of change in wavelength at each pixel, dλ i , and determine the integrated flux: where F i is the flux and h i is the continuum level. The integrated flux measurements of our optical Keck spectra are published in Hsyu et al. (2018), except for five new systems, which are listed in Table 1. Our Keck NIR observations are described in Table 2. The measured emission line flux ratios of our systems, along with the 1053 systems derived from the SDSS galaxy sample that satisfy our S/N criteria, are also available on GitHub as MCMC input files as part of our primordial helium code, YMCMC. 4 The total reported error of the emission line fluxes comprises of two terms added in quadrature: the measured error of the integrated emission line flux and an assumed 2 percent relative flux uncertainty to account for the error of the flux calibration. The latter follows a common procedure in primordial helium work (Skillman et al. 1994;) and is taken from Oke (1990), which quantified the absolute flux uncertainties on a set of 25 standard stars now recognized as the Hubble Space Telescope spectrophotometric standards. Oke (1990) found these standard stars to be reliable to about 1-2% across the optical wavelength regime (see Table VI

HeBCD Sample
To the PHLEK and SDSS samples, we also add the HeBCD sample of galaxies from Izotov & Thuan (2004); , a fraction of which have follow-up NIR observations reported in Izotov et al. (2014). Their sample consists of 93 total systems, 21 of which have unique optical plus NIR spectroscopy, i.e., we do not consider systems with optical spectra reported for multiple regions but one NIR spectra. This is to ensure that the optical and NIR emission line fluxes originate from observations of the same part of a singular H II region. The HeBCD data set have metallicities that overlap with both our PHLEK sample and the SDSS sample, with a median metallicity of (O/H) × 10 5 = 9.40 or 12 + log 10 (O/H) = 7.97. For these systems, we take the reported emission line flux ratios and equivalent widths but re-determine their best-fit parameters, including the helium abundance, using our model, as described below in Section 3. Updated optical data of the HeBCD sample were obtained from E. Aver (2018;private communication) and include the He I λ4026 flux and corrections to the original values found in .

MODEL OVERVIEW
Most of the hydrogen and helium in an H II region is in an ionized state. Thus, the number abundance ratio of helium to hydrogen, y, of an H II region is given by the sum of the abundance ratios of singly and doubly ionized helium: The y + and y ++ abundances depend on the intrinsic helium to hydrogen ratio of the H II region, along with the detailed physical state of the ionized gas and the surrounding stellar population. Since the observed He I and H I relative line ratios depend on these physical parameters, we can measure the He I and H I line ratios to pin down the physical conditions of the ionized gas. Our analysis follows a similar approach NOTE-Optical emission line fluxes of systems observed using LRIS and previously unreported in Hsyu et al. (2018). The reported values are integrated flux measurements given relative to the Hβ flux, which is also quoted for reference. These fluxes are uncorrected for reddening, since reddening is a parameter we later solve for in our MCMC analysis.
to that described first by Aver et al. (2011) and subsequently by Aver et al. (2012Aver et al. ( , 2013Aver et al. ( , 2015. Our code YMCMC solves for the best fit parameters that reproduce the measured emission line ratios of our sample of galaxies described in Section 2. YMCMC closely follows the model and methods mentioned in the above works using a Python implementation of a Markov Chain Monte Carlo (MCMC) sampler, EMCEE (Foreman-Mackey et al. 2013), to survey an 8-dimensional parameter space: • the ionized helium abundance, y + • the electron temperature, T e [K] • the electron density, n e [cm −3 ] • the reddening parameter, c(Hβ) • the underlying hydrogen stellar absorption, a H [Å], normalized to the amount of absorption at Hβ • the underlying helium stellar absorption, a He [Å], normalized to the amount of absorption at He I λ4026 • the helium optical depth parameter, τ He , normalized to the value at He I λ3889 • the ratio of neutral to singly ionized hydrogen density, ξ ≡ n(H I)/n(H II) At each step of the MCMC chain, our model predicts the He I and H I emission line fluxes as a ratio relative to Hβ and calculates the log-likelihood function of the model: where σ(λ) is the uncertainty of the flux ratio of each emission line. The subscripts p and m represent the predicted and measured flux ratios, respectively. The predicted flux ratio of the hydrogen emission lines is given by: Here, E(λ) is the emissivity of an emission line at wavelength λ, EW (λ) is the measured equivalent width (EW ) of the emission line, C R (λ) is the collisional to recombination correction factor, and f (λ) is the reddening law. These individual components are discussed in further detail below. For helium emission lines, the predicted flux ratio is similarly given by: where f τ (λ) is the optical depth function. As shown in Equations 4 and 5, our model for predicting flux ratios depends on the measured quantity EW (Hβ), which has a corresponding uncertainty. To account for this uncertainty, at each step of the MCMC chain, we draw a new value for EW (Hβ) from a Gaussian distribution with a width equal to the measured uncertainty. This is the same approach adopted by Aver et al. (2011). Additionally, we perturb EW (Hα) and EW (Pγ) for our PHLEK sample and for systems with NIR data, respectively. In these two cases, we require EW (Hα) and EW (Pγ) to predict the theoretical F(Pγ) / F(Hβ) and F(Hα) / F(Hβ) ratios, which we use to match our predicted model fluxes to the format of our measured input fluxes (see Section 3.6 for details).
We further note that the equivalent width and the measured flux are not independent of one another. However, a conserved quantity is the height of the continuum around each emission line, h(λ). To ensure that the equivalent widths used in Equations 4 and 5 scale appropriately with the predicted fluxes, we introduce the following relation: which allows us to rewrite Equations 4 and 5, removing EW (λ) entirely, as follows: With Equations 7 and 8, YMCMC generates the model flux ratios given a set of parameters drawn from the MCMC. Motivated by physically meaningful limits, we impose the following uniform priors on the following parameters: The upper limit placed on log 10 (ξ) here is unrealistic for an H II region, as this upper bound would imply that only 55 per cent of the gas in the H II region is ionized. We allow our MCMC to explore this regime, but disqualify systems that have best recovered solutions that are unreasonable for H II regions (see Section 4.1).
To ensure that the electron temperature parameter explored by our MCMC stays within reasonable limits for the system, we include the following weak prior on T e : where we take σ to be 0.2T m and T m is the direct measurement of the electron temperature based on the [O III] λ4363 / (λ4959 + λ5007) ratio. This weak prior was also implemented by Aver et al. (2011), who demonstrated with synthetic data that the above prior improves the recovery of the input model parameters and removes local minima near the edges of the likelihood distributions. We also require that the electron temperature is within the range 10, 000 K ≤ T e ≤ 22, 000 K. In the following sections, we describe in detail the implementation of each term in Equations 7 and 8.

Emissivity
The H I and He I emissivities, denoted by E(λ), provide a measure of the energy released per unit volume and time. E(λ) is expressed in units of erg s −1 cm −3 throughout.

Hydrogen Emissivity
Our model determines the emissivity of an H I line at a given temperature and density following the hydrogen emissivity calculations made by P. Storey (2018; private communication) assuming Case B recombination. The Storey 2018 emissivities extend the Storey & Sochi (2015) hydrogen emissitivies down to the lowest density regime explored by our model, log 10 (n e /cm −3 ) = 0, and are available up to log 10 (n e /cm −3 ) = 5 at log 10 (n e /cm −3 ) = 1 intervals. The emissivities are calculated over the temperature range T e = 5,000 -25,000 K, at 1,000 K intervals. In our model, we interpolate linearly within this temperature and density grid using SCIPY's RectBivariateSpline().
The implementation of H I emissivities in our model assumes no error in the emissivity value. As an estimate of the uncertainty on these emissivities, we compare E(Hβ) from Storey 2018 with the parameterization of E(Hβ) by R. L. Porter (given in Eq. 3.1 of Aver et al. 2010; we note that this parameterization is independent of the electron density). Within the temperature and density ranges of interest, the emissivities differ by 0.10-0.55 per cent. At a fixed temperature, the difference in emissivities increases with increasing electron density. The ratio of the Hα, Hγ, and Hδ to Hβ emissivities from Storey 2018 differ by 0.10-0.20 percent compared to the parameterizations in Aver et al. (2010). That is, the extended Storey 2018 H I emissivities are not expected to significantly change our model and therefore the resulting best fit MCMC parameters. Rather, they make the H I emissivity grid more self-consistent, as it no longer relies on extrapolations to the lowest density regime.

Helium Emissivity
The He I line emissivity at a given temperature and density is determined in a similar manner to the H I emissivities. Our model adopts the He I emissivities introduced in Aver et al. (2013), which project the Porter et al. (2012Porter et al. ( , 2013 He I emissivities onto a finer grid. The Porter et al. (2012Porter et al. ( , 2013 emissivities assume Case B recombi-nation and are calculated for a grid of temperatures ranging from T e = 10,000 K -25,000 K and densities from log 10 (n e /cm −3 ) = 1 -5. We linearly interpolate the He I emissivities within this temperature and density grid using SCIPY's RectBivariateSpline() interpolator.
3.2. Collisional to Recombination Ratio, C R (λ) The collisional to recombination ratio, C R (λ), corrects for the amount of neutral hydrogen and helium atoms excited to higher energy states due to collisions with electrons, and the emission detected as a result of the electrons subsequently cascading down to lower energy levels.

Hydrogen
Following the method of calculating the collisional to recombination correction factor in Aver et al. (2010), the C R (λ) ratio of an H I line is given by: Here, n(H I) and n(H II) are, respectively, the neutral and ionized hydrogen densities in units of cm −3 . The ratio of these densities is defined as ξ and solved for as one of the free parameters in the MCMC.
The subscripts used in the numerator of Equation 10 represent energy level transitions to: the energy level i, which is above or equal to the transition level of interest, j (i.e., i ≥ j). In the denominator of Equation 10, the effective recombination rate from an ionized energy level above energy level j is given by α +→ j . The subsequent downward transition from j → n = 2 then gives rise to the Balmer wavelength of interest. For the Paschen series, the transition of interest becomes j → n = 3.
The numerator of Equation 10 expresses the contribution of emission stemming from collisional excitations. q 1→i represents the rate coefficient of collisional excitation from the ground state n = 1 to a higher energy level i, in cm 3 s −1 . The value of q 1→i depends on the effective collision strength of the transition, Υ 1i , reported in Anderson et al. (2000Anderson et al. ( , 2002 such that: (11) Once collisionally excited to the higher energy level i, the electron can cascade downward via various paths leading to energy level j. The probability of the different transition paths, n l → nl, are expressed as branching ratios, BR i→ j , and reported in Omidvar (1983). The j → 2 transition of interest then occurs, with various path probabilities captured in the term BR j→2 . Finally, the collisional contribution to emission depends both on the density of neutral hydrogen atoms and the density of electrons in the gas available for collisions, n(H I) and n e . Anderson et al. (2000Anderson et al. ( , 2002 only report collision strengths up to principle quantum number n = 5. While the contribution of collisional emission for transitions n > 5 is expected to be small, we apply scaling factors to quantify the collisional contributions of emission lines emanating from transitions n > 5, specifically Hδ (6 → 2), H8 (8 → 2), and Pγ (6 → 3): For all other transitions n ≤ 5, the n l orbitals we include are: • Hα: 3s, 3p, 3d, 4s, 4p, 4d, 4 f • Hβ: 4s, 4p, 4d, 4 f , 5s, 5p, 5d, 5 f , 5g • Hγ: 5s, 5p, 5d, 5 f , 5g The denominator of Equation 10 represents the amount of emission from the recombination of free electrons with ionized hydrogen atoms and the subsequent cascade down to less excited states. This is expressed by α +→ j , the rate that an ionized hydrogen atom recombines and transitions from higher energy levels, represented by +, down to the j energy level, in units of cm 3 s −1 . However, the value of α +→ j must be proportional to all the emission that subsequently emanates from transitions out of energy level j: And therefore, we can use emissivities of the subsequent transitions out of an energy level j as a proxy for α +→ j . For example, any recombination that brings an electron to energy level j = 4 will then transition out of j = 4 via either the 4 → 3 transition or the 4 → 2 transition. Thus, rather than using values of the recombination rates in the literature, we choose to substitute α +→ j with our latest hydrogen emissivities following Equation 13. This allows us to take advantage of the more refined temperature and density grid for which we have emissivity values. The functional form of the collisional corrections for our H I lines of interest, over a range of temperatures, are shown in Figure 2. In this figure, we use a neutral to ionized hydrogen density ratio of ξ = 10 −4 for illustration. The corrections for Hδ and H8 are scaled from the correction for Hγ following Equation 12, and similarly, the correction for Pγ is scaled from Pβ.

Helium
The C R (λ) correction for He I is folded in to the CLOUDY modelling done by Porter et al. (2012Porter et al. ( , 2013 in their latest emissivity work. The correctional factors are therefore included in our implementation of the interpolated He I emissivities described in Section 3.1.2. We refer readers to Section 3 of Aver et al. (2013) for a more detailed description of the collisional contribution included in the Porter et al. (2012Porter et al. ( , 2013 emissivities. We note that since the Porter et al. emissivities (E(λ)/E(Hβ) in Equation 8) include the collisional correction, we set the value of C R (λ) = 0 in the version of Equation 8 implemented in our code YMCMC.

Underlying Absorption
The observed H I and He I emission line fluxes are compromised by underlying stellar absorption from the atmospheres of the stars in the H II region. Failing to correct for the missing emission can lead to underestimating the total integrated flux of the H I and He I lines. The amount of underlying absorption depends on the particular stellar population in the galaxy. However, information about the specific stellar population, its age, and its metallicity, along with the possibility of multiple stellar populations, is difficult to extract from long-slit spectroscopy of the H II region. While older works assumed a constant EW of underlying stellar absorption at all H I and He I lines (Olive & Skillman 2001), it is now recognized that these values are wavelength dependent; assuming a constant amount of underlying absorption across the spectrum biases the derived value of the primordial helium abundance. Various works have estimated the average amount of stellar absorption expected at each H I and He I line based on synthetic spectra (González Delgado et al. 1999Delgado et al. , 2005. Our wavelength dependent underlying absorption corrections are given as coefficients normalized to the amount of absorption present at Hβ for H I lines and He I λ4471 for He I lines. Our model incorporates the coefficient values introduced by Aver et al. (2010) and repeated in Equations 4.2 and 4.3 of Aver et al. (2015) to include the stellar absorption at the NIR He I λ10830 and Pγ lines. The Aver et al. (2010) values represent the relative EW s of underlying absorption suitable over a range of stellar ages as calculated from a suite of stellar population models. These coefficients are summarized in the following subsections.

Hydrogen
The coefficients of underlying hydrogen stellar absorption are given below, normalized to the amount of absorption present at Hβ, referenced as the variable a H (λ), and in units of EW (Å). The value at a H (H8) is extrapolated from a linear fit to the wavelength and coefficients from Hβ to Hδ. We exclude Hα from the fit due to the decreasing nature of the underlying absorption at redder wavelengths.

Helium
We apply a correction to the optical and NIR He I emission lines to account for underlying stellar absorption. The following values are given in EW (Å) and normalized to the amount of underlying helium absorption at He I λ4471, denoted by the general variable a He (λ). That is, the amount of stellar absorption at a given He I line is the correctional coefficient at its wavelength, multiplied by a He (λ). The value of a He (He I λ5015) given is determined using a linear fit to the coefficients of all other listed optical He I lines.

Reddening Correction
Our observed emission line fluxes are expected to suffer from reddening due to dust along the line-of-sight. The theoretical emissivities of the H I recombination lines are wellknown and relatively insensitive to the temperature and density of the gas, and are therefore well-suited to the determination of the amount of reddening present in the observed spectrum. To correct for this effect, we include a logarithmic correction factor c(Hβ) in Equations 7 and 8. When combined with a reddening law, f (λ), the amount of extinction as a function of wavelength can be inferred. In our work, we assume the reddening law presented in Equations 2 and 3 of Cardelli et al. (1989). Using the formulation given by Cardelli et al. (1989), we generate a list of f (λ) values for a wavelength grid of 1000 values between 3100 -13000Å. We then linearly interpolate this functional form at the observed wavelengths of the H I and He I emission lines.
The best fit value of c(Hβ) includes reddening within our own Milky Way and in the observed system. For our sample of galaxies, however, the reddening correction is expected to be small because our candidate systems were selected to be away from the disk of the Milky Way and are expected to be of lower metallicity, where the effects of dust are less important. We note that it is typical to assume no error in the assumed reddening law (Olive & Skillman 2000).

Optical Depth Function
The optical depth function is a correction term that accounts for photons that are emitted but subsequently reabsorbed or scattered out of our line of sight. Accordingly, the correction depends on optical depth, the temperature, and the density of the gas. We use a set of optical depth corrections that are suited to the modelling of low metallicity H II regions (Benjamin et al. 2002). These assume Case B recombination, a spherically symmetric H II region with no systemic expansion or velocity gradients, and are valid for a temperature and density range of T e = 12,000 -20,000 K and n e = 1 -300 cm −3 . The coefficients of the fits to the optical depth correction are presented in Table 4 of Benjamin et al. (2002) and can be found listed in Equation A3 in the Appendix of Olive & Skillman (2004). The formulation of He I λ10830 is not included in the original work, but we apply the formula given by Equation 2.2 of Aver et al. (2015). For completeness, we give the functional form of the fits below in Equation 16, and the coefficients of individual He I lines are given in Table 3.

MCMC Details
To determine the best fit parameters of each system via MCMC, our code YMCMC reads in a file containing the fol-lowing four columns of measured values for a suite of emission lines: (1) the flux ratio, (2) the flux ratio uncertainty, (3) the equivalent width of the line in units of Å, and (4) the uncertainty of the equivalent width of the line (Å). The flux ratios and their corresponding errors are given relative to Hβ for all optical emission lines, while Pγ is used for the NIR He I λ10830 line. Since the input NIR flux ratio is not given relative to Hβ, our model separately calculates the predicted flux of He I λ10830 and Pγ relative to Hβ, and combines these two predicted values to match the input format, F(He I λ10830) / F(Pγ): where the right hand side of the equation can be calculated using Equations 5 and 4, for the numerator and denominator respectively. Additionally, our model predicts a total flux ratio of the blended H8+He I λ3889 lines, which differs from the deblending technique employed by Aver et al. (2010;see their Equation 4.1). To do this, our model individually predicts an H8 and He I λ3889 emission line flux ratio and sums the two for a blended flux: Finally, we note that our observations of the PHLEK sample used a dichroic at 5600 Å (see Section 2.1.1 for details), which means that the Hα and Hβ emission lines are detected on separate blue and red arms of LRIS. For these systems, we adapt the MCMC code to model the flux ratios relative to Hα for all optical emission lines that are detected on the red side of LRIS, in a manner equivalent to that of Equation 17 for the NIR emission lines. As standard, every emission line on the blue side of LRIS is modelled relative to Hβ. Because of the dichroic, we lose the Hα / Hβ Balmer line ratio in our analysis, and this has a minor impact on our ability to solve for the parameters. To test how the loss of F(Hα) / F(Hβ) affects our results, we generated a synthetic spectrum with emission line fluxes mirroring the format of our LRIS observations, and tested our MCMC's ability to recover the input model parameters. We show the results of this test in Appendix B. As expected, we do not constrain parameters that depend on the Balmer lines as tightly -for example, the 1σ errors on c(Hβ) double when we lose information on F(Hα) / F(Hβ). However, we find that even without F(Hα) / F(Hβ), our recovered parameters are within 1σ of the input parameters, and the errors on y + increase by a factor of just 1.06−1.20.
Our MCMC analysis uses 500 walkers and 1000 steps to determine the best fit model parameters of each system. We take our burn-in to be a conservative 0.8× n steps = 800 steps and dispose of all samples before the burn-in, leaving us with  . This functional form has been developed specifically for helium abundance measurements of H II regions and are valid only in the temperature and density range of Te = 12,000 -20,000 K and ne = 1 -300 cm −3 . There are no optical depth corrections for the singlet lines He I λ5015 and He I λ6678, i.e., fτ (λ) = 1.
10 5 samples. To ensure that our MCMC chains have converged, we require the best recovered parameters from the two halves of the 10 5 samples to agree to within a few percent. In this exercise, all best recovered y + values agree to within half a per cent. We show an example contour plot and histogram of the recovered parameters of J0118+3512 in Figure 3.

THE PRIMORDIAL HELIUM ABUNDANCE
In the following sections, we describe the sample definition, the calculation of the metal abundances, and our determination of the primordial helium abundance.

Qualification
To identify the systems that are most suitable for determining the primordial helium abundance, we first require that all systems have a measured EW (Hβ) ≥ 50 . This ensures that systems have higher emission line flux-to-continuum level ratios, thus weak emission lines are less affected by underlying absorption. In particular, this minimizes the effect of underlying He I absorption and ensures that the measured He I emission line ratios do not under predict the true helium abundance (Izotov & Thuan 2004;). This cut eliminates 3 systems from our PHLEK sample, 463 systems from the SDSS sample, and 4 from the HeBCD sample.
We also require that the recovered best fit parameters from the MCMC analysis are physical. Specifically, we remove all systems with recovered optical depths τ He > 4 and neutral to ionized hydrogen fractions ξ > 0.01, if the 1σ lower bound on the recovered value of ξ does not encompass ξ = 0.001. These specific values follow the work most recently highlighted by Aver et al. (2015), although we have opted to completely eliminate systems with recovered parameters in the regimes stated above, while Aver et al. (2015) consider some of these systems as part of their flagged data set. We then assert that the MCMC analysis recovers parameters that are able to successfully reproduce all measured emission line ratios, according to the criteria listed in Sections 4.1.1 and 4.1.2 below.

Sample 1
Our most stringent criteria is that all of the measured emission line ratios are reproduced to within 2σ, given the best recovered parameters from the MCMC. This qualification criteria is demonstrated for the galaxy J0118+3512 in Figure 4. We label the systems that qualify via these conditions "Sample 1". Sample 1 contains 3 galaxies from the PHLEK sample, 38 galaxies from the SDSS sample, and 13 galaxies from the HeBCD sample, resulting in a total of 54 systems. These systems and their best-recovered parameters are listed in part in Table 4 and are available in full online. The full MCMC chains for Sample 1 are available on GitHub as part of the HCPB20 branch of our primordial helium code, YMCMC.

Sample 2
We also consider a more lenient qualification criteria, which requires that all emission line ratios are reproduced to within 2σ, with the exception of one emission line ratio, which must be reproduced to within 3σ. We call this "Sample 2".
Sample 2 consists of all of the systems in Sample 1, plus an additional 4 galaxies from our PHLEK sample, 48 galaxies from the SDSS sample, and 13 galaxies from the HeBCD sample. Thus, Sample 2 contains a total of 119 galaxies. These systems and their best recovered parameters from the MCMC are available online as part of Table 4. The full y + = 0.0763 +0.0076  MCMC chains for Sample 2 are available on GitHub as part of the HCPB20 branch of YMCMC.

Abundance Measurements
To calculate ionic abundances, we utilize the emission line analysis package PYNEB (Luridiana et al. 2015). 5 To obtain  Hβ and Hα do not appear in this figure since our LRIS blue and red side emission lines are measured relative to those two emission lines, and therefore those emission lines do not carry any information. Objects that qualify in Sample 1 require that the emission line ratios reproduced by our model are within 2σ of the measured value. Note that the distributions shown by the histograms reflect the measurement errors. a value and error on an ionic abundance, we calculate 10 5 Monte Carlo realizations of each abundance by perturbing the measured flux ratios by their errors. For each realization, we use relevant parameters derived from our MCMC samples, namely the electron density, n e , and reddening parameter, c(Hβ). Our reported abundances and their errors are the mean and standard deviation of the 10 5 Monte Carlo realizations.
H II regions are expected to be in the low-density regime, where density diagnostics observed at optical wavelengths, such as the [S II] λλ6717, 6731 doublet, are not very sensitive to n e (see Figure 5.3 of Osterbrock 1989). As such, the n e value recovered by the MCMC analysis is loosely constrained when our observations do not include lines that are strongly sensitive to n e in the low density regime. Previous works in the literature usually assume n e = 100 for ionic abundance calculations instead of the measured electron density, an assumption that is within the 1σ bounds of their measured values. This choice is also within the range of densities expected of H II regions, n e = 100 -10, 000 cm −3 (Osterbrock 1989). However, n e can be pinned down when the density-sensitive He I λ10830 line is included in the analysis (see our distribution and best recovered value of log(n e /cm −3 ) in Figure 3 as an example, as well as the results of our trial MCMC runs on mock data including the He I λ10830 line in Appendix B). We therefore adopt the n e values sampled by our MCMC as input to PYNEB for the determination of the ionic abundances.

Oxygen
The total oxygen abundance O/H is the sum of the singly and doubly ionized ionic abundances: The values of O + /H + and O ++ /H + (hereafter abbreviated as O + and O ++ ) of each galaxy depend on the measured emission line flux ratios relative to Hβ, the electron temperature, and the electron density. We adopt a two-zone approximation of an H II region, with two distinct electron temperatures characterizing the high-and low-ionization zones. The O ++ abundance is calculated using the [O III] λλ4959, 5007 flux ratios in combination with the high ionization zone temperature, t 3 , where we calculate values of t 3 using the temperature sensitive [O III] λ4363 line. Thus, the value of t 3 differs from the electron temperature parameter in our MCMC model, T e , but the difference is expected to be small.
As mentioned previously, we calculate 10 5 values of the O ++ abundance, each time adopting an electron density value as sampled in the 10 5 density realizations of the MCMC chain. The measured [O III] flux ratios are perturbed each time by drawing a new value from a Gaussian distribution with a mean of the measured flux value and standard deviation of its measurement error. We also calculate a new value of t 3 at each step in the MCMC using the perturbed [O III] flux ratios.
The O + abundance is calculated using the [O II] λλ3727, 3729 doublet and the low-ionization zone temperature, t 2 . A direct measure of t 2 requires a detection of the [O II] λλ7320, 7330Å lines or the [N II] λ5755Å line (used in conjunction with the [N II] λλ6548, 6584 doublet). Since we do not detect these lines, we infer t 2 from t 3 following the relation from Pagel et al. (1992), which is based on the photoionization model grids by Stasińska (1990): t 2 = 20, 000 K 10, 000 K t 3 + 0.8 The total oxygen abundance of each system is calculated by summing the singly and doubly ionized oxygen abundances (i.e., Equation 18). The final reported oxygen abundance and its corresponding error, is calculated by taking the mean and standard deviation of the 10 5 Monte Carlo realizations. All abundance calculations are made using PYNEB's getIonAbundance() method. We report the ionic and total oxygen abundances of a subset of our systems in Table  5 and in full online.

Helium
The total helium abundance, y, is the sum of the abundances of singly ionized helium y + and doubly ionized helium y ++ (see Equation 2). y + is recovered as a parameter of the MCMC analysis, and the presence of y ++ in an H II region can be inferred via emission at He II λ4686Å (Pagel et al. 1992;Skillman et al. 2013). Therefore, if the He II λ4686 line is detected, we calculate and include the y ++ abundance in the total helium abundance. A non-detection of He II λ4686 in the spectrum is assumed to indicate a negligible y ++ abundance.
As with calculating the oxygen abundance, we assume an electron density n e as recovered by the MCMC. However, for the helium abundances, we also assume the electron temperature T e from the MCMC chains. We make 10 5 realizations of the y ++ abundance by perturbing the measured He II flux ratios by the error in its measurement. While we expect doubly ionized helium to occupy a region of higher temperatures than T e (i.e., the temperature at which singly ionized helium is found), this assumption has a negligible effect on the total helium abundance, since the y ++ abundance typically contributes a ∼ 1 per cent correction to the overall helium abundance.
Furthermore, some of the helium in H II regions may be in the neutral state; thus, the total helium abundance may require a correctional factor for undetected neutral helium. To assess whether a neutral helium component is present, we follow the use of the radiation softness parameter, η (Vilchez & Pagel 1988), defined as: to estimate the hardness of the ionizing radiation. The S ++ abundance depends on the temperature of the gas, t S ++ . A direct measure of t S ++ requires the detection of [S III] emission at λ6312Å, λ9069Å, and λ9532Å, the latter two of which fall outside the wavelength coverage of our instrument setup. 6 The S ++ abundance is extremely sensitive to temperature (Garnett 1992), therefore, rather than assuming the value of t S ++ to be t 3 or t 2 , it is necessary to estimate the temperature of the S ++ zone following the relation from Garnett (1992): We assume t S + = t S ++ , following the expected ionization structure in a two-zone photoionization model (see e.g., Figure 2 of Garnett 1992). Adopting this temperature, the S ++ and S + abundances can be calculated with the [S III] λ6312 and [S II] λλ6717, 6731 emission line fluxes. Based on photoionization models, Pagel et al. (1992) concluded η to be suitable for determining whether a correctional factor is necessary for undetected neutral helium; if log(η) < 0.9, the neutral helium abundance can be assumed to be negligible (see Figure 6 of Pagel et al. 1992). We choose to exclude the systems from our sample that were found to have a non-negligible neutral helium abundance following this metric (4 total systems, all of which are from the SDSS sample), due to the additional uncertainties introduced when assuming a correctional factor.
The ionic and total helium abundances of our systems are partially listed in Table 5 and are available in full online.

Extrapolation to y P
The standard approach for determining the primordial helium abundance is to perform a linear regression to a set of measured oxygen and helium abundances. This technique was initially proposed by Peimbert & Torres-Peimbert (1974, 1976 and is still used by the most recent primordial helium abundance investigations. The analysis follows the expectation from BBN calculations that most of the helium in the Universe is produced during BBN, while essentially no oxygen is produced. Through the chemical evolution of stars, there is a net production of 4 He, but this contribution is relatively minor compared to the quantity of 4 He produced during BBN. Therefore, the post-BBN contribution to the 4 He abundance can be modelled as a small (linear) deviation from the BBN value that increases with increasing metallicity. We also note that Fernández et al. (2018) have recently proposed that a tighter relation exists between the helium abundance and the sulphur abundance. Their work suggests that, as far as chemical evolution is concerned, sulphur may trace helium better than oxygen. However, in our work, we do not have access to the emission lines required to measure the sulphur abundance, and we therefore use the O/H abundance in what follows.
Our determination of the primordial helium abundance, y P , is based on a fit to the measured O/H and He/H ≡ y number abundance ratios of the galaxies that qualify for Sample 1 and 2. We note that our choice to use the helium number abundance ratio differs from the typical format historically found in the literature, where the primordial helium abundance is expressed as the primordial helium mass fraction, Y P . For reference, the helium mass fraction, Y , can be converted from y using: Here, Z is the metallicity dependent heavy element mass fraction, which is linearly proportional to the constant c, which depends on chemical evolution (see directly below for a further discussion of this constant).
We have decided to use the helium number abundance y instead of the helium mass fraction Y , as done historically, for the following reasons: • Observations of the helium abundance are intrinsically measuring a number abundance ratio.
• Calculations of the helium abundance are computed as a ratio of volume densities (n4 He /n1 H ), and later con- verted to a mass fraction to match the observationally reported mass fractions (see e.g., Pitrou et al. 2018).
• The primordial helium mass fraction (Y P ) is not actually the fraction of mass in the form of 4 He. It is defined as the ratio of volume densities Y P = 4 n( 4 He)/n b , where n b is the baryon density. Thus, the term "mass fraction" is a misnomer that should probably be avoided as we enter the era of precision cosmology.
• Our choice eliminates the dependence on Z, whose value has varied across primordial helium works. For reference, Pagel et al. (1992) and Aver et al. (2015) both For these reasons, we have therefore chosen to quote our primordial helium abundance in the form we most directly measure and most appropriate to compare to theoretical valuesthe primordial helium number abundance ratio, y P . However, a comparison of our measured y P to previously reported values of Y P can be simply calculated with the following equation: Y P = 4y P 1 + 4y P Our linear fits to the two galaxy samples described in Sections 4.1.1 (Sample 1) and 4.1.2 (Sample 2) are optimized using EMCEE given the likelihood function of our linear model: Here, the summation is over all individual galaxies in each sample. Our linear model is given by mx n + b, where the x n are our measured O/H values, the slope m ≡ dy/d(O/H), and the intercept b ≡ y P . We capture the error on our calculated O/H abundances by drawing new values of O/H from a Gaussian with a mean of the calculated values and standard deviation of the calculated errors during each step of the MCMC procedure. The total measured error of y n is captured by the term σ yn . We also introduce the term σ intr to our likelihood function to quantify the intrinsic scatter of our sample of y measurements to account for unknown systematic uncertainties that are introduced by our model, following the method presented in Section 4.3 of .
To solve for the parameters that best describe our linear model and the intrinsic scatter, we use 1000 walkers each taking 1000 steps in the MCMC. We set the following uniform priors on the model parameters: The range of allowed y P values matches the range of y + values of our model described in Section 3. Similarly, we allow a generous range of possible dy/d(O/H) values. The range in σ intr is chosen to be comparable to the measurement error of the y values, σ yn . We find a mean of σ yn = 0.005, and allow for the range of σ intr to be twice that value, although it is desirable for this parameter to be less than σ yn . We conservatively use a burn-in of 800 steps. The distribution of the explored parameter space and the best recovered parameters using Sample 1 and Sample 2 are shown in Figure 5. In Figure 6, we plot Samples 1 and 2, along with their best fit linear models and extrapolations to y P . The optimal parameter values recovered from the MCMC for Sample 1 are: with χ 2 /dof = 0.82. While the linear fits to Sample 1 and Sample 2 are comparable, we adopt y P from Sample 1 as our reported value and in all further analyses. This choice is motivated by our model being able to more confidently reproduce all of the observed emission line fluxes of Sample 1. This increases our confidence in the recovered parameters, including the primordial helium abundance. This confidence is also reflected in the recovered value of σ intr , which is consistent with zero for Sample 1 but non-zero for Sample 2. The recovered y P values from Sample 1 and Sample 2 are within 1σ of each other, but we note that Sample 1 and Sample 2 are not independent of one another (i.e., Sample 1 is a subset of Sample 2).

Comparison to Existing Measurements of Y P
We now compare our result to existing primordial helium abundance measurements that are reported in the literature. To allow for a comparison of the primordial helium number abundance ratio, y P , we convert all literature measurements of Y P to y P using Equation 21. The literature results are summarized in Table 6. Our result agrees with measurements derived from emission line observations of H II regions in nearby galaxies (Aver et al. 2015;Peimbert et al. 2016;Fernández et al. 2019;Valerdi et al. 2019), absorption line observations of a near-pristine gas cloud along the line-of-sight to a background quasar (Cooke & Fumagalli 2018), the primordial helium abundance derived from the damping tail of the CMB recorded by the Planck satellite (Planck Collaboration et al. 2016), and SBBN calculations of the primordial abundances (Cyburt et al. 2016;Pitrou et al. 2018) that assume a baryon-to-photon ratio η = (5.931 ± 0.051) × 10 −10 , which is based on the observationally measured abundance of primordial deuterium, (D/H) P = (2.527 ± 0.030) × 10 −5 .
The primordial helium abundance that we report here is in 2.6σ disagreement with the Izotov et al. (2014) result, y P = 0.0856 ± 0.0010. We note that the HeBCD sample included in this work was compiled by Izotov et al. (2014) and was subsequently the sample analyzed by Aver et al. (2015). Aver et al. (2015) also find a discrepancy with the Izotov et al. results (2.2σ), and suggest several possible reasons for the disagreement. Given that our model closely follows that of Aver et al. (2010Aver et al. ( , 2012Aver et al. ( , 2013Aver et al. ( , 2015, we expect many of these reasons to equally be relevant in the comparison between our result and that of Izotov et al. (2014). For example, Izotov et al. first use the observed Balmer line ratios to solve for the amount of reddening and underlying hydrogen stellar absorption present, while assuming that the underlying absorption is the same for all hydrogen lines. After correcting observations for reddening, they then use Monte Carlo to find the best fit value of y + , given T e , n e , and τ He and the observed He I lines. This differs from the MCMC method adopted in this work, which solves for all parameters using all observed emission lines simultaneously. Within their model, Izotov et al. implement a correction for hydrogen emission resulting from collisional excitation based on CLOUDY photoionization modelling. There are also slight differences in the assumed coefficients for underlying stellar absorption between our models, the incorporation and scaling of the NIR lines to Hβ, and the calculation of t 2 and t S ++ from t 3 . Finally, Izotov et al. apply a cut on their sample prior to solving for the best fit parameters, based on properties such as their measured EW (Hβ) and ionization parameter. Our approach, on the other hand, solves for the best fit parameters of every galaxy and subsequently use this information to decide if a system qualifies. We refer readers to Izotov et al. (2006Izotov et al. ( , 2013Izotov et al. ( , 2014 for details of their model and approach.
It is reassuring that our extrapolation to y P is in agreement with numerous existing values reported in the literature. Most of these are based on distinct samples, a variety of sample sizes, and adopt different analysis methods. This does not, however, rule out the need for improvements in future primordial helium research; we discuss current model deficiencies that could warrant additional enhancements in Section 5.2.

DISCUSSION
In this section, we use our determination of y P to place a limit on physics beyond the Standard Model and discuss future improvements that can be made to push measurements of y P to sub-percent level accuracy.  (2018) and is BBN-independent. All values are quoted with 1σ confidence limit, except the CMB value, which is quoted with 2σ confidence limit, as indicated.

Implications for the Standard Model -BBN bounds on Ω b h 2 and N eff
Physics beyond the Standard Model at the time of BBN can be identified by comparing observational measurements of the primordial abundances with the SBBN predicted values. The primordial element abundances produced during BBN are captured primarily by two parameters: the baryon density, Ω b h 2 , and the effective number of neutrino species, N eff . By adopting a measurement of Ω b h 2 from the CMB (Planck Collaboration et al. 2018) and assuming N eff = 3.046 (i.e. the Standard Model value; Cyburt et al. 2002Cyburt et al. , 2016Pitrou et al. 2018), BBN is a parameter free theory. Note that the SBBN predicted abundances are still subject to other uncertainties, such as the mean neutron lifetime τ n and nuclear reaction rates, but these values are measured in laboratories or inferred using ab initio calculations. Primordial abundances deduced from observations of astrophysical regions thus provide a valuable test of the Standard Model of particle physics and cosmology and its assumptions.
Constraining the values of Ω b h 2 and N eff using observations requires using two or more measurements of the primordial abundances. For this exercise, we take our measurement of the primordial helium abundance in conjunction with the latest primordial deuterium abundance reported by : In what follows, we use the detailed primordial abundance calculations reported by Pitrou et al. (2018). These authors provide formulae for calculating the primordial abundances, given a value of Ω b h 2 , N ν , and τ n . We restate the formula for predicting Y P here as an example (see their Equation 145, the surrounding text, and Table VI of their paper for the values of the C pqr coefficients that are referenced here): We use the latest measurement of the mean neutron lifetime τ n = 877.7 ± 0.7 (Pattie et al. 2018) to solve for Ω b h 2 and N ν . Furthermore, we use the scaling N eff = N ν × 3.046/3 (see Pitrou et al. 2018). This choice of scale is commonly used, and allows us to fairly compare the BBN results to the CMB. We use EMCEE with 100 walkers taking 1500 steps each and sample the parameter space: With each step, a model set of primordial D/H and Y P abundances are predicted. The optimal parameters are solved for assuming a Gaussian likelihood function. We take the burn in to be at 0.8× n steps = 1200 steps. Given the observed primordial element abundances, we report the following bounds on the effective number of neutrino species and   The result of our MCMC calculation is shown in Figure 7, together with the Planck bounds on these parameters. 7 The best fit value of N eff is consistent with the value inferred by Planck (N eff = 2.92 +0.19 −0.18 ; shown by the red contours in Figure  7) and the Standard Model value of N eff = 3.046.

Future Improvements
The goal of the spectroscopic survey reported by Hsyu et al. (2018) was twofold: (1) to increase the number of known systems in the low-metallicity regime and in particular, push on the lowest-metallicity regime, and (2) to obtain high-quality optical and NIR spectroscopy of a subset of 7 We use the Planck Release 3 data with the prefix 'plikHM_TTTEEE_lowl_lowE'. the new, metal-poor galaxies, with priority on systems with metallicities determined to be 12 + log 10 (O/H) ≤ 7.65 based on strong-line calibration methods. The purpose of the specific goals was to better populate and constrain the metalpoor end in the extrapolation to a primordial helium abundance. In the following text, we discuss the current limitations of the PHLEK survey and future improvements that need to be explored to push y P to sub-percent level accuracy.

Qualification Rates
One of the main obstacles of measuring the primordial helium abundance is the difficulty of accurately modelling a large fraction of emission line observations. To give an overview of the current status of modelling the H II region emission lines, in Table 7 we have compiled the qualification rates of the three survey samples considered in this paper. Our results show that the qualification rates of H II regions that only have optical data are consistent among the PHLEK, SDSS, and HeBCD data sets -about 10% of systems make it into Sample 1. The meager number of currently known, nearpristine systems that push on the lowest-metallicity regime hinders our ability to constrain the slope (and thus intercept) of the linear extrapolation to the primordial value. The effect of a dearth of the most metal-poor systems is multiplied when these systems are unsuccessfully modelled and consequently excluded from primordial helium analyses after quality screening.
However, Table 7 shows that systems with complementary optical and NIR data are more successfully modelled. Our PHLEK sample sees an increase from 10% to 13% of systems qualifying for Sample 1 and more noticeably, the HeBCD sample success rate increases to 29% when NIR data are included. An aspect that contributes to the small fraction of systems that can be well-modelled lies in the difficulty of confidently detecting the weak optical He I lines necessary for accurately determining the physical conditions of the H II region, including the helium abundance. The addition of the NIR He I λ10830 line to primordial helium work saw an appreciable reduction in the errors on the recovered helium abundances and electron densities (Izotov et al. 2014;Aver et al. 2015). The value of the He I λ10830 line is likely the reason behind the higher success rates we see in Table 7 for systems with optical and NIR spectroscopy. The sensitivity of the He I λ10830 line emissivity to the electron density eliminates degeneracies between the electron temperature and electron density when modelling systems. Although the He I λ10830 line is the brightest emission line detected in our NIR observations, it still sometimes eluded detection completely. From our experience, it is useful to target systems with F(Hβ) 10 −15 erg s −1 cm −2 for complementary NIR spectroscopy. This assumes existing optical data, but the criteria increases the chance of acquiring high S/N NIR data.
Even with our systems that satisfy this criteria, however, we recover a lower success rate in modelling the PHLEK sample compared to the HeBCD sample. Table  7 shows that a total of 27% of our systems with optical plus NIR spectroscopy qualify in Sample 2 (which includes the systems that qualify in Sample 1), compared to 67% for the HeBCD sample. We presume the difference comes from these two data sets consisting of different types of galaxies. The correlation between F(Hβ) and a NIR detection mentioned above is not unlike the criteria imposed by Izotov & Thuan as part of their HeBCD sample selection. Specifically, the construction of the HeBCD data set was based on existing observations and used a selection criteria of high EW (Hβ), quoted to be generally EW (Hβ) ≥ 200 Å, and have metallicities ranging from 12 + log 10 (O/H) = 7.00 − 8.21 (Izotov & Thuan (2004); although we note that only 35 of the 93 HeBCD sample satisfies the EW (Hβ) condition, and 22 of the 93 fall in the lowmetallicity regime). Subsequent analysis of the HeBCD data set by  was combined with SDSS DR5 spectroscopy, with the requirement that only SDSS galaxies with EW (Hβ) ≥ 50 Å, and F(Hβ) 10 −14 erg s −1 cm −2 are included in the analysis. While we also impose a EW (Hβ) ≥ 50 Å criteria (see Section 4.1), a comparison of the F(He I λ10830) / F(Pγ) detection levels shows that all HeBCD systems have He I λ10830Å to Pγ ratios detected with S/N ≥ 50, whereas the PHLEK sample have S/N ≤ 50, regardless of the measured EW (Hβ) and the status of the modelling success rate.
Thus, it is evident that a sample selection based on measured high EW (Hβ), as with the HeBCD data set, versus one based on low estimated metallicities via strong-line calibrations, as with the PHLEK sample, yield different data sets. The former yielded systems with higher significance detections of the necessary NIR He I λ10830 emission line and a higher modelling success rate. Meanwhile, the latter successfully populates the lowest-metallicity end of the galaxy sample, but currently face more significant limitations in accurately modelling their physical conditions and characteristics, even with He I λ10830. Until we identify and improve current shortcomings, which may lie in data processing or in model simplicities and deficiencies, we are not equipped to equally include all types of metal-poor systems to deduce the primordial helium abundance. Now that a substantial sample of metal poor star-forming galaxies are known, we suggest that a more detailed analysis of individual systems may allow us to better model the complicated physics of H II regions. This in turn may allow us to construct a model with improved capabilities of recovering the helium abundance in a wider variety of star-forming galaxies.

Towards a Sub-percent Measurement of yP
Other potential obstacles faced by primordial helium analyses come post-data collection, a facet of which is in data processing. Currently, fluxing our emission line spectra using observations of spectrophotometric standards introduces uncertainties in relative flux measurements between 1-2%, shown by Oke (1990). The weakest He I lines and their measured flux ratios are therefore easily susceptible to errors introduced during flux calibration. To push observational primordial helium measurements to the sub-percent level will require flux calibrations to mirror this precision. Izotov & Thuan (2004);  take precaution to derive sensitivity curves using only hot white dwarf standard stars that show relatively weak absorption features, such as Feige 34, Feige 110, and HZ 44. Such stars allow for sensitivity curves that are accurate to 1% over the optical wavelength range (Oke 1990). Following this necessity of subpercent flux calibrations, we propose the use of the nearperfect blackbody stars from Suzuki & Fukugita (2018) to flux calibrate observations of metal-poor H II regions. The spectra of these stars, thought to be white dwarfs, are nearly featureless, and their blackbody nature ranges from the ultraviolet to infrared. These stars offer the potential to improve the precision of flux calibrations even further and can bring primordial helium abundances closer to the sub-percent level.
Additionally, the model we assume in this work to describe our emission line observations is subject to deficiencies. As part of our analysis, we investigated obvious shortcomings in our model, such as the inability to model a specific emission line or the inability to model systems when their parameters fall in a particular regime. Reassuringly, we found no obvious parameters or combination of parameters that perform poorly for non-qualifying systems. However, since the model is unable to reproduce a high fraction of the initial galaxy sample, we conclude that some aspects of H II region modelling are currently unaccounted for.
Helium abundance measurements have historically been derived from longslit observations. These observations are assumed to be representative of the entire H II region. In reality, the integrated light that enters the slit likely samples multiple radii of an H II region and can also be the result of multiple, overlapping H II regions. Such simplistic assumptions likely affect our ability to fit the observed data with a single set of parameters. Possible model enhancements include dropping the simple two-zone photoionization model characterized by two temperatures, and introducing a temperature structure to our model. However, we note that we do not anticipate the temperature structure within an H II region to vary much beyond the limits we place on our temperature prior (selected to be σ = 0.2T m , where T m is the direct measurement of the electron temperature from the [O III] lines). Similarly, the density of the H II region is likely a function SDSS · · · · · · · · · 38 (7%) 85 (15%) 578 HeBCD 6 (29%) 14 (67%) 21 7 (10%) 12 (17%) 69 NOTE-The number (and percentage) of systems from our PHLEK sample, the SDSS sample, and the HeBCD sample that qualify for Sample 1 and Sample 2, out of the total number of systems available in each sample. The statistics are separated by systems for which optical and NIR spectroscopy are available and systems for which only optical spectroscopy exists. We remind readers that Sample 1 is included in Sample 2, and the two are described in Sections 4.1.1 and 4.1.2 respectively. The HeBCD optical+NIR sample has the highest rate of satisfying our criteria for Sample 1 and Sample 2, likely due to the higher S/N NIR data that exist for the HeBCD sample. The qualification rate for systems with only optical data are comparable across all data sets and lower than the rate for systems with complementary NIR data in the same data set when the comparison is available. These rates illustrate the difficulty of modelling systems well without complementary NIR spectroscopy and more importantly, the need for high-quality NIR data.
of distance from the central star, and introducing a density structure may improve our model as well.

SUMMARY AND CONCLUSION
We present a sample of NIR observations of several metalpoor galaxies reported by Hsyu et al. (2018). Using this sample along with galaxies from the SDSS spectroscopic database and existing metal-poor galaxies in the literature, we report a new determination of the primordial helium abundance. We summarize the main results of our analysis as follows: 1. We obtain near-infrared (NIR) spectra of a sample of sixteen galaxies to complement optical spectroscopy presented by Hsyu et al. (2018). The NIR observations are taken using NIRSPEC or NIRES at Keck Observatory and are designed to obtain a measurement of the He I λ10830Å to Pγ flux ratio. We supplement this sample with 1053 starburst galaxies in the SDSS spectroscopic database, selected based on their star-forming nature and sufficiently high S/N data on a suite of optical He I lines, and 93 systems from the  HeBCD sample, a subset of which include follow-up NIR observations reported by Izotov et al. (2014).
2. We outline our Python-based code YMCMC, which uses a Markov Chain Monte Carlo (MCMC) approach to find parameters that best describe the observed emission line flux ratios of each galaxy. The parameters we solve for in the analysis are: the singly ionized helium abundance, the electron temperature and density, the reddening parameter, the underlying H I and He I stel-lar absorption, the helium optical depth parameter, and the ratio of neutral to ionized hydrogen densities. Our method is largely based on the approach developed by Aver et al. (2011Aver et al. ( , 2012Aver et al. ( , 2013Aver et al. ( , 2015. Our implementation of the techniques include new H I emissivities that extend to the lowest density regime, n e = 1 cm −3 , a different treatment of the blended H8+He I λ3889 lines, and a different method of correcting for H I emission stemming from collisional excitation. We also employ an alternative qualification approach, which requires that we statistically reproduce all emission lines used in the analysis. 3. Using YMCMC, we solve for the best fit parameters that describe our sample of galaxies. We construct two qualifying samples, named "Sample 1" and "Sample 2". Sample 1 contains all galaxies whose H I and He I emission lines are reproduced by our model to within 2σ. Sample 2 is defined such that all except one observed emission line ratios are reproduced to within 2σ (however, the one emission line that fails this 2σ limit must be reproduced to within 3σ).
4. We calculate ionic abundances of O + , O ++ , and y ++ , and combine these with the y + values recovered from the MCMC analysis to calculate total number abundance ratios, O/H and y. We fit a linear model to the O/H versus y abundances of Sample 1 and Sample 2, and extrapolate to zero metallicity to infer the primordial helium abundance, y P . Our linear model allows for the presence of an intrinsic scatter of the measurements due to systematic uncertainties that may currently be unaccounted for. We find that Sample 1 contains no evidence of intrinsic scatter, while Sample 2 contains some intrinsic scatter. Both samples yield primordial helium abundances that are in mutual agreement with one another. However, we adopt the y P determination based on Sample 1 due to our increased confidence in the model. We report a primordial helium number abundance ratio y P = 0.0805 +0.0017 −0.0017 , which corresponds to a primordial helium mass fraction Y P = 0.2436 +0.0039 −0.0040 . 5. Combining our determination of y P with (D/H) P from , we find Ω b h 2 = 0.0215 +0.0005 −0.0005 and N eff = 2.85 +0.28 −0.25 . This value of N eff is within 1σ agreement with the Standard Model value of N eff = 3.046. Our value of Ω b h 2 is in 1.3σ agreement with the value measured by Planck.
Observational measurements of the primordial light element abundances provide a unique window to study the conditions of the early Universe and offer the potential to identify non-standard physics at the time of BBN. The latest (D/H) P determination has reached the percent level, comparable to the precision achieved by the latest CMB constraints. y P determinations are reaching similar precision; in particular, the recent addition of the NIR He I λ10830 line to helium abundance analyses has led to an improvement in the helium abundance measurements of individual galaxies.
From a theoretical perspective, the primordial helium abundance can be reliably calculated; a high precision observational determination of the primordial helium abundance will therefore provide the most sensitive test of the Standard Model. As we move towards this era of high precision cosmology, we advocate that it will become necessary to understand the nuances of our current limitations of H II region modelling and push data processing techniques to higher accuracy. These can potentially be done by conducting detailed observations of individual systems and improving flux calibration using near-featureless blackbody stars as standards. Finally, we strongly suggest that observational primordial helium works, and also BBN calculations, shift towards re-porting helium number abundances, as opposed to a helium mass fraction commonly adopted in literature, for the most direct comparison to theoretical works. The following text shows the CasJobs query submitted to retrieve the full sample of galaxies in the SDSS spectroscopic database potentially suited to be included in our primordial helium abundance determination. The query requires that these systems be star-forming galaxies within a redshift range of 0.02 < z < 0.15. This guarantees that the [O II] doublet and He I λ7065 lines, necessary for oxygen and helium abundance measurements, are detected.
SELECT e.specObjID, e.ra, e.dec, e.z, e.zErr, e.subClass, dbo.fGetUrlFitsSpectrum(e.specObjID) as urlfits into mydb.specobj_starburst FROM specObj AS e WHERE e.z < 0.15 AND e.z > 0.02 AND e.subClass = 'STARBURST' B. MOCK DATA AND MCMC RECOVERY Given a set of input parameter values and EW s, we generate mock data (i.e., the would-be observed flux ratios) to test how well YMCMC can recover the input parameters. In all these test runs, we adopt a weak temperature prior equal to the input temperature of T e = 18,000 K. The results are shown in the following table. NOTE-MCMC results on the recovery of parameters, given mock data. The corresponding contours and histograms showing the best recovered model parameters on these trial runs are shown in Figures 8, 9, 10, and 11. Optical+NIR, All Figure 8. Contours (off-diagonal panels) and histograms (diagonal panels) showing the best recovered model parameters on mock data including optical and near-infrared data, with the F(Hα)/F(Hβ) ratio, i.e., the first column of recovered values in Table 8. The contours show the 1σ, 2σ, and 3σ levels. The solid green line in the histograms show the best recovered parameter value, and the dotted green lines show the ±1σ values. In the panel showing the results for the log10(ξ) parameter, the solid vertical line represents a 2σ upper limit. Optical+NIR, no F (Hα)/F(Hβ) Figure 9. Same as Figure 8, but on mock data including optical and near-infrared data, without the F(Hα)/F(Hβ) ratio, i.e., the second column of recovered values in Table 8. Optical Only, All Figure 10. Same as Figure 8, but on mock data including only optical data, with the F(Hα)/F(Hβ) ratio, i.e., the third column of recovered values in Table 8. Optical Only, no F (Hα)/F(Hβ) Figure 11. Same as Figure 8, but on mock data including only optical data, without the F(Hα)/F(Hβ) ratio, i.e., the fourth column of recovered values in Table 8.