The PHLEK Survey: A New Determination of the Primordial Helium Abundance

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 following:

• 1.  
  The temperature sensitive [O iii] λ4363 Å line, for a direct measurement of the oxygen abundance.
• 2.  
  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 Å.
• 3.  
  The NIR He i λ10830 Å line, whose emissivity is the most sensitive He i emission line to the density of the gas, relative to ${\rm{P}}\gamma \,\lambda$ 10940 Å.

In addition to these emission lines, we also detect in our spectra the [O ii] doublet at $\lambda \lambda$ 3727, 3729 Å, the [O iii] doublet at λλ4959, 5007 Å, the [N ii] doublet at λλ6548, 6584 Å, the [S ii] doublet at λλ6717, 6731 Å, and the Balmer series from Hα to at least H8.

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 near-infrared nights be allocated within a week of one another, so that our complementary observations for a given target would 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).

2.1.1. Optical Spectroscopy

2.1.2. Near-infrared Spectroscopy

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.⁴ We used a boxcar extraction technique to extract a single one-dimensional (1D) spectrum of each object.⁵ 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 1D 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 1D 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.

Zoom In Zoom Out Reset image size

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 required the systems to be: (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 calculated 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(λ)):

$$\begin{eqnarray*}&&\begin{array}{c}{\rm{S}}/{\rm{N}}([{\rm{O}}\,{\rm{III}}]\,\lambda 4363)\geqslant 5\\ {\rm{S}}/{\rm{N}}({\rm{He}}\,{\rm{I}}\,\lambda 5876)\geqslant 20\\ {\rm{S}}/{\rm{N}}({\rm{He}}\,{\rm{I}}\,\lambda 4471)\geqslant 3\\ {\rm{S}}/{\rm{N}}({\rm{He}}\,{\rm{I}}\,\lambda 6678)\geqslant 3\\ {\rm{S}}/{\rm{N}}({\rm{He}}\,{\rm{I}}\,\lambda 7065)\geqslant 3\end{array}\end{eqnarray*}$$

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 $({\rm{O}}/{\rm{H}})\times {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 $({\rm{O}}/{\rm{H}})\,\times {10}^{5}=4.82$. These values correspond to $12+{{\rm{log}}}_{10}({\rm{O}}/{\rm{H}})$ values of 8.12 and 7.68, respectively.

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 modeled using the Absorption LIne Software (ALIS); see Cooke et al. (2014) for a more detailed description of the software.⁶ ALIS simultaneously fits the emission line profile using a Gaussian model and the surrounding continuum using a 1 or 2D Legendre polynomial, and determines the best-fit parameters of the Gaussian and continuum model using a χ² minimization approach. Systems with high emission line fluxes, however, are not well-represented 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 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 in order 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:

$$\begin{eqnarray}&&\ F(\lambda )=\sum _{i}\,({F}_{i}-{h}_{i})d{\lambda }_{i},\end{eqnarray} \tag{ 1 }$$

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.⁷

The total reported error of the emission line fluxes is comprised of two terms added in quadrature: the measured error of the integrated emission line flux, and an assumed 2% 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; Izotov et al. 2007) 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 6 of Oke 1990).

To the PHLEK and SDSS samples, we also add the HeBCD sample of galaxies from Izotov & Thuan (2004) and Izotov et al. (2007), 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 only one in the NIR spectrum. 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 $({\rm{O}}/{\rm{H}})\times {10}^{5}\,=9.40$ or $12+{{\rm{log}}}_{10}({\rm{O}}/{\rm{H}})=$ 7.97. For these systems, we take the reported emission line flux ratios and equivalent widths but redetermine 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); these include the He i λ4026 flux as well as corrections to the original values found in Izotov et al. (2007).

The H i and He i emissivities, denoted by E(λ), provide a measure of the energy released per unit volume and time. The value of E(λ) is expressed in units of erg s⁻¹ cm⁻³ throughout.

3.1.1. Hydrogen Emissivity

3.1.2. Helium Emissivity

The collisional to recombination ratio, $\tfrac{C}{R}(\lambda )$, corrects for the amount of neutral hydrogen and helium atoms excited to higher-energy states due to collisions with electrons, as well as the emission detected as a result of the electrons subsequently cascading down to lower energy levels.

3.2.1. Hydrogen

Following the method of calculating the collisional to recombination correction factor in Aver et al. (2010), the $\tfrac{C}{R}(\lambda )$ ratio of an H i line is given by:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{C}{R}(\lambda ) & = & \displaystyle \frac{n({\rm{H}}\,{\rm{I}})\,(\sum _{i}{q}_{1\to i}{{BR}}_{i\to j})\,{{BR}}_{j\to n}{n}_{{\rm{e}}}}{n({\rm{H}}\,{\rm{II}})\,{\alpha }_{+\to j}{{BR}}_{j\to n}{n}_{{\rm{e}}}}\\ & = & \xi \times \displaystyle \frac{\sum _{i}{q}_{1\to i}{{BR}}_{i\to j}}{{\alpha }_{+\to j}}.\end{array}\end{eqnarray} \tag{ 10 }$$

Here, $n({\rm{H}}\,{\rm{I}})$ and $n({\rm{H}}\,{\rm{II}})$ are, respectively, the neutral and ionized hydrogen densities in units of cm⁻³. 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 ${\alpha }_{+\to j}$. The subsequent downward transition from $j\to n=2$ then gives rise to the Balmer wavelength of interest. For the Paschen series, the transition of interest becomes $j\to n=3$.

The numerator of Equation (10) expresses the contribution of emission stemming from collisional excitations. Here, ${q}_{1\to i}$ represents the rate coefficient of collisional excitation from the ground state n = 1 to a higher energy level i, in cm³ s⁻¹. The value of ${q}_{1\to i}$ depends on the effective collision strength of the transition, ϒ_1i, reported in Anderson et al. (2002, 2000) such that:

$$\begin{eqnarray}&&{q}_{1\to i}=4.004\times {10}^{-8}\,\sqrt{\displaystyle \frac{1}{{k}_{{\rm{B}}}T}}\,{\rm{\exp }}\left(\displaystyle \frac{-13.6{\rm{eV}}(1-\tfrac{1}{{i}^{2}})}{{k}_{{\rm{B}}}T}\right){{\rm{\Upsilon }}}_{1i}.\end{eqnarray} \tag{ 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^{\prime} l^{\prime} \to {nl}$, are expressed as branching ratios, ${{BR}}_{i\to j}$, and reported in Omidvar (1983). The $j\to 2$ transition of interest then occurs, with various path probabilities captured in the term ${{BR}}_{j\to 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({\rm{H}}\,{\rm{I}})$ and n_e.

Anderson et al. (2002, 2000) 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\to 2$), H8 ($8\to 2$), and Pγ ($6\to 3$):

$$\begin{eqnarray}&&\displaystyle \frac{C}{R}(\lambda )=\displaystyle \frac{C}{R}({\rm{H}}\gamma | {\rm{P}}\beta )\,{\rm{\exp }}\left(\displaystyle \frac{-13.6{\rm{eV}}(\tfrac{1}{{5}^{2}}-\tfrac{1}{{j}^{2}})}{{k}_{{\rm{B}}}T}\right).\end{eqnarray} \tag{ 12 }$$

For all other transitions n ≤ 5, the $n^{\prime} l^{\prime}$ orbitals we include are:

• 1.  
  Hα: $3s,\,3p,\,3d,\,4s,\,4p,\,4d,\,4f$;
• 2.  
  Hβ: $4s,\,4p,\,4d,\,4f,\,5s,\,5p,\,5d,\,5f,\,5g$;
• 3.  
  Hγ: $5s,\,5p,\,5d,\,5f,\,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 ${\alpha }_{+\to 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⁻¹. However, the value of ${\alpha }_{+\to j}$ must be proportional to all the emission that subsequently emanates from transitions out of energy level j:

$$\begin{eqnarray}&&{n}_{+}{n}_{{\rm{e}}}{\alpha }_{+\to j}h\nu =\sum _{k}E{\left(\lambda \right)}_{j\to k}.\end{eqnarray} \tag{ 13 }$$

Therefore, we can use emissivities of the subsequent transitions out of an energy level j as a proxy for ${\alpha }_{+\to 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\to 3$ transition or the $4\to 2$ transition. Thus, rather than using values of the recombination rates in the literature, we choose to substitute ${\alpha }_{+\to 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 $\xi ={10}^{-4}$ for illustration.

Zoom In Zoom Out Reset image size

3.2.2. Helium

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. 1999, 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 EWs 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.

3.3.1. 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}_{{\rm{H}}}({\rm{H}}8)$ 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.

$$\begin{eqnarray}\begin{array}{ll}{a}_{{\rm{H}}}({\rm{H}}\alpha )=0.942\qquad & {a}_{{\rm{H}}}({\rm{H}}\delta )=0.896\\ {a}_{{\rm{H}}}({\rm{H}}\beta )=1.000\qquad & {a}_{{\rm{H}}}({\rm{H}}8)=0.882\\ {a}_{{\rm{H}}}({\rm{H}}\gamma )=0.959\qquad & {a}_{{\rm{H}}}({\rm{P}}\gamma )=0.400\end{array}\end{eqnarray} \tag{ 14 }$$

3.3.2. 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}_{\mathrm{He}}(\lambda )$. 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}_{\mathrm{He}}({\rm{He}}\ {\rm{I}}\,\lambda 5015)$ given is determined using a linear fit to the coefficients of all other listed optical He i lines.

$$\begin{eqnarray}\begin{array}{rcl}{a}_{{\rm{He}}}({\rm{He}}\ {\rm{I}}\,\lambda 3889) & = & 1.400\quad {a}_{{\rm{He}}}({\rm{He}}\ {\rm{I}}\,\lambda 5876)=0.874\\ {a}_{{\rm{He}}}({\rm{He}}\ {\rm{I}}\,\lambda 4026) & = & 1.347\quad {a}_{{\rm{He}}}({\rm{He}}\ {\rm{I}}\,\lambda 6678)=0.525\\ {a}_{{\rm{He}}}({\rm{He}}\ {\rm{I}}\,\lambda 4471) & = & 1.000\quad {a}_{{\rm{He}}}({\rm{He}}\ {\rm{I}}\,\lambda 7065)=0.400\\ {a}_{{\rm{He}}}({\rm{He}}\ {\rm{I}}\,\lambda 5015) & = & 1.016\quad {a}_{{\rm{He}}}({\rm{He}}\ {\rm{I}}\,\lambda 10830)=0.800\end{array}\end{eqnarray} \tag{ 15 }$$

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 well-known 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({\rm{H}}\beta )$ 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 and 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({\rm{H}}\beta )$ 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 2001).

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 modeling 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⁻³. 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.

$$\begin{eqnarray}&&\ {f}_{\tau }(\lambda )=1+\displaystyle \frac{\tau }{2}[a+({b}_{0}+{b}_{1}{n}_{{\rm{e}}}+{b}_{2}{n}_{{\rm{e}}}^{2})\,{T}_{4}]\end{eqnarray} \tag{ 16 }$$

where T₄ = T_e/10,000 K.

To determine the best-fit parameters of each system via MCMC, our code yMCMC reads in a file containing the following 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γ):

$$\begin{eqnarray}&&\displaystyle \frac{F({\rm{He}}\ {\rm{I}}\,\lambda 10830)}{F(P\gamma )}=\displaystyle \frac{F({\rm{He}}\ {\rm{I}}\,\lambda 10830)}{F(H\beta )}/\displaystyle \frac{F(P\gamma )}{F(H\beta )},\end{eqnarray} \tag{ 17 }$$

where the right-hand side of the equation can be calculated using Equations (5) and (4), for the numerator and denominator, respectively.

Our model also 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 predicts individual H8 and He i λ3889 emission line flux ratios and sums the two for a blended flux:

$$\begin{eqnarray*}&&\displaystyle \frac{F({\rm{H8+He}}\,{\rm{I}}\,\lambda 3889)}{F(H\beta )}=\displaystyle \frac{F({\rm{H}}8)}{F({\rm{H}}\beta )}+\displaystyle \frac{F({\rm{He}}\ {\rm{I}}\,\lambda 3889)}{F({\rm{H}}\beta )}.\end{eqnarray*}$$

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 modeled 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({\rm{H}}\alpha )/{\rm{F}}({\rm{H}}\beta )$ 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({\rm{H}}\alpha )/{\rm{F}}({\rm{H}}\beta$). However, we find that, even without $F({\rm{H}}\alpha )/{\rm{F}}({\rm{H}}\beta$), 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\times {n}_{\mathrm{steps}}$ = 800 steps and dispose of all samples before the burn-in, leaving us with 10⁵ samples. To ensure that our MCMC chains have converged, we require the best recovered parameters from the two halves of the 10⁵ samples to agree to within a few percent. In this exercise, all best recovered ${y}^{+}$ values agree to within half a percent. We show an example contour plot and histogram of the recovered parameters of J0118+3512 in Figure 3.

Zoom In Zoom Out Reset image size

To identify the systems that are most suitable for determining the primordial helium abundance, we first require that all systems have a measured ${EW}({\rm{H}}\beta )\geqslant 50\,\mathring{\rm A}$. This ensures that systems have higher emission line flux-to-continuum level ratios, and 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; Izotov et al. 2007). 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 ${\tau }_{\mathrm{He}}\,\gt \,4$ and neutral to ionized hydrogen fractions $\xi \,\gt \,0.01$, if the 1σ lower bound on the recovered value of ξ does not encompass $\xi =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.

4.1.1. Sample 1

Our most stringent criterion 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.

Zoom In Zoom Out Reset image size

Table 4.  Best Recovered Parameters from MCMC Analysis

Galaxy	y+	Te	${{\rm{log}}}_{10}({n}_{{\rm{e}}}{/{\rm{cm}}}^{-3}$)	c(Hβ)	${a}_{{\rm{H}}}$	${a}_{\mathrm{He}}$	${\tau }_{\mathrm{He}}$	${\mathrm{log}}_{10}(\xi$)
		[K]			[Å]	[Å]
J0118+3512	${0.0763}_{-0.0073}^{+0.0076}$	${12700}_{-1400}^{+1500}$	${1.84}_{-0.82}^{+0.35}$	${0.153}_{-0.099}^{+0.132}$	${1.5}_{-1.0}^{+1.5}$	${0.39}_{-0.19}^{+0.20}$	${2.6}_{-1.2}^{+1.2}$	$-{3.3}_{-1.9}^{+1.8}$
J2030−1343	${0.0725}_{-0.0050}^{+0.0058}$	${12400}_{-1200}^{+1400}$	${1.50}_{-1.01}^{+0.89}$	${0.280}_{-0.061}^{+0.048}$	${0.53}_{-0.38}^{+0.70}$	${0.19}_{-0.12}^{+0.16}$	${0.55}_{-0.38}^{+0.65}$	$-{2.9}_{-2.1}^{+1.7}$
KJ29	${0.081}_{-0.015}^{+0.011}$	${11440}_{-940}^{+1320}$	${1.76}_{-1.01}^{+0.54}$	${0.096}_{-0.069}^{+0.106}$	${0.34}_{-0.24}^{+0.44}$	${0.80}_{-0.34}^{+0.28}$	${1.6}_{-1.1}^{+1.8}$	$-{3.1}_{-2.0}^{+1.9}$
spec-0301-51942-0531	${0.0915}_{-0.0064}^{+0.0059}$	${13700}_{-1500}^{+1500}$	${2.19}_{-1.42}^{+0.61}$	${0.179}_{-0.070}^{+0.069}$	${0.58}_{-0.43}^{+0.92}$	${0.22}_{-0.16}^{+0.24}$	${2.1}_{-1.1}^{+1.2}$	$-{1.60}_{-2.42}^{+0.71}$
spec-0364-52000-0187	${0.0863}_{-0.0029}^{+0.0045}$	${11790}_{-930}^{+1110}$	${1.25}_{-0.85}^{+0.86}$	${0.378}_{-0.047}^{+0.032}$	${0.45}_{-0.32}^{+0.49}$	${0.250}_{-0.090}^{+0.103}$	${0.82}_{-0.43}^{+0.45}$	$-{2.3}_{-2.6}^{+1.5}$
spec-0375-52140-0118	${0.0842}_{-0.0056}^{+0.0051}$	${14400}_{-1600}^{+1600}$	${1.89}_{-1.20}^{+0.76}$	${0.208}_{-0.039}^{+0.033}$	${0.53}_{-0.39}^{+0.75}$	${0.51}_{-0.25}^{+0.26}$	${1.84}_{-0.93}^{+0.90}$	$-{3.5}_{-1.6}^{+1.4}$
I Zw 18 SE1	${0.0763}_{-0.0028}^{+0.0031}$	${17900}_{-2000}^{+1900}$	${1.82}_{-0.15}^{+0.14}$	${0.016}_{-0.011}^{+0.018}$	${3.65}_{-0.59}^{+0.54}$	${0.24}_{-0.16}^{+0.22}$	${0.64}_{-0.42}^{+0.55}$	$-{4.71}_{-0.87}^{+0.92}$
SBS 0940+5442	${0.0804}_{-0.0023}^{+0.0033}$	${17100}_{-1400}^{+1400}$	${1.937}_{-0.095}^{+0.090}$	${0.053}_{-0.028}^{+0.025}$	${2.11}_{-0.93}^{+0.98}$	${0.39}_{-0.14}^{+0.15}$	${0.29}_{-0.21}^{+0.30}$	$-{3.75}_{-1.55}^{+0.94}$
Mrk 209	${0.0820}_{-0.0023}^{+0.0025}$	${17400}_{-1900}^{+1900}$	${1.85}_{-0.15}^{+0.14}$	${0.018}_{-0.012}^{+0.018}$	${1.93}_{-0.84}^{+0.79}$	${0.26}_{-0.12}^{+0.12}$	${1.26}_{-0.81}^{+1.02}$	$-{4.69}_{-0.90}^{+1.02}$

Note. The best recovered values of the eight parameters sampled with our MCMC analysis. These parameters describe a subset of galaxies from Sample 1, defined to be systems where the best recovered parameters can reproduce all the measured emission line flux ratios to within 2σ. This table lists three systems from our PHLEK sample (top three rows; see Section 2.1), the SDSS sample (middle three rows; see Section 2.3), and the HeBCD sample (bottom three rows; see Section 2.5).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

4.1.2. Sample 2

To calculate ionic abundances, we utilize the emission line analysis package PyNeb (Luridiana et al. 2015).⁸ To obtain a value and error on an ionic abundance, we calculate 10⁵ 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⁵ 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⁻³ (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}_{{\rm{e}}}\,\,{\mathrm{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.

4.2.1. Oxygen

The total oxygen abundance ${\rm{O}}/{\rm{H}}$ is the sum of the singly and doubly ionized ionic abundances:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{O}}}{{\rm{H}}}=\displaystyle \frac{{{\rm{O}}}^{+}}{{{\rm{H}}}^{+}}+\displaystyle \frac{{{\rm{O}}}^{++}}{{{\rm{H}}}^{+}}.\end{eqnarray} \tag{ 18 }$$

The values of ${{\rm{O}}}^{+}/{{\rm{H}}}^{+}$ and ${{\rm{O}}}^{++}/{{\rm{H}}}^{+}$ (hereafter abbreviated as O⁺ and O⁺⁺) of each galaxy depend on the measured emission line flux ratios relative to ${\rm{H}}\beta$, 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₃, where we calculate values of t₃ using the temperature sensitive [O iii] λ4363 line. Thus, the value of t₃ 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⁵ values of the O⁺⁺ abundance, each time adopting an electron density value as sampled in the 10⁵ 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₃ 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₂. A direct measure of t₂ 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₂ from t₃ following the relation from Pagel et al. (1992), which is based on the photoionization model grids by Stasińska (1990):

$$\begin{eqnarray}&&{t}_{2}=20,000\,{\rm{K}}/\left(\displaystyle \frac{10,000\,{\rm{K}}}{{{\rm{t}}}_{3}}+0.8\right).\end{eqnarray} \tag{ 19 }$$

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 are calculated by taking the mean and standard deviation of the 10⁵ 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.

Table 5.  Ionic and Total Abundances of Oxygen and Helium

Galaxy	O+/H+	O++/H+	$({\rm{O}}/{\rm{H}})$	${y}^{+}$	${y}^{++}$	y
	(×105)	(×105)	(×105)
J0118+3512	${1.44}_{-0.27}^{+0.38}$	${5.1}_{-1.3}^{+2.2}$	${6.5}_{-1.3}^{+2.2}$	${0.0764}_{-0.0075}^{+0.0076}$	${0.0028}_{-0.0002}^{+0.0002}$	${0.0792}_{-0.0075}^{+0.0076}$
J2030−1343	${2.28}_{-0.41}^{+0.55}$	${7.9}_{-2.0}^{+3.1}$	${10.2}_{-2.1}^{+3.1}$	${0.0725}_{-0.0050}^{+0.0058}$	${0.0018}_{-0.0001}^{+0.0001}$	${0.0743}_{-0.0050}^{+0.0058}$
KJ29	${2.10}_{-0.41}^{+0.43}$	${5.5}_{-1.5}^{+1.7}$	${7.6}_{-1.6}^{+1.8}$	${0.081}_{-0.015}^{+0.011}$		${0.081}_{-0.015}^{+0.011}$
spec-0301-51942-0531	${2.34}_{-0.33}^{+0.51}$	${7.1}_{-1.7}^{+2.6}$	${9.5}_{-1.7}^{+2.7}$	${0.0915}_{-0.0064}^{+0.0059}$	${0.0014}_{-0.0003}^{+0.0003}$	${0.0929}_{-0.0064}^{+0.0059}$
spec-0364-52000-0187	${2.73}_{-0.41}^{+0.51}$	${12.4}_{-2.9}^{+3.7}$	${15.1}_{-2.9}^{+3.7}$	${0.0863}_{-0.0029}^{+0.0045}$	${0.0007}_{-0.0001}^{+0.0001}$	${0.0870}_{-0.0030}^{+0.0045}$
spec-0375-52140-0118	${1.24}_{-0.19}^{+0.29}$	${7.7}_{-1.8}^{+2.9}$	${8.9}_{-1.8}^{+3.0}$	${0.0842}_{-0.0056}^{+0.0051}$	${0.0008}_{-0.0002}^{+0.0002}$	${0.0850}_{-0.0056}^{+0.0051}$
I Zw18 SE1	${0.465}_{-0.055}^{+0.083}$	${1.31}_{-0.27}^{+0.39}$	${1.78}_{-0.28}^{+0.40}$	${0.0763}_{-0.0028}^{+0.0031}$	${0.0008}_{-0.0002}^{+0.0002}$	${0.0772}_{-0.0028}^{+0.0031}$
SBS 0940+5442	${0.436}_{-0.045}^{+0.058}$	${3.37}_{-0.53}^{+0.71}$	${3.81}_{-0.53}^{+0.71}$	${0.0804}_{-0.0023}^{+0.0033}$	${0.0005}_{-0.0001}^{+0.0001}$	${0.0810}_{-0.0023}^{+0.0033}$
Mrk 209	${0.679}_{-0.087}^{+0.123}$	${4.38}_{-0.93}^{+1.30}$	${5.06}_{-0.94}^{+1.30}$	${0.0820}_{-0.0023}^{+0.0025}$	${0.0011}_{-0.0000}^{+0.0000}$	${0.0831}_{-0.0023}^{+0.0025}$

Note. The singly and doubly ionized oxygen abundances, total oxygen abundance, singly and doubly ionized helium abundances, and total helium abundance for a subset of galaxies from Sample 1. This table lists three systems from our PHLEK sample (top three rows; see Section 2.1), the SDSS sample (middle three rows; see Section 2.3), and the HeBCD sample (bottom three rows; see Section 2.5). The online version of this table also contains the remaining galaxies in Sample 1 as well as Sample 2.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

4.2.2. 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)). We recover y⁺ 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 nondetection 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⁵ 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 % 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

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{O}^{+}}{{S}^{+}}\times \displaystyle \frac{{S}^{++}}{{O}^{++}},\end{eqnarray} \tag{ 20 }$$

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.⁹ The ${S}^{++}$ abundance is extremely sensitive to temperature (Garnett 1992). Therefore, rather than assuming the value of ${t}_{{S}^{++}}$ to be t₃ or t₂, it is necessary to estimate the temperature of the ${S}^{++}$ zone following the relation from Garnett (1992):

$$\begin{eqnarray*}&&{t}_{{S}^{++}}=0.83{t}_{3}+0.17.\end{eqnarray*}$$

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 (four 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.

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 ⁴He, but this contribution is relatively minor compared to the quantity of ⁴He produced during BBN. Therefore, the post-BBN contribution to the ⁴He abundance can be modeled 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 sulfur abundance. Their work suggests that, as far as chemical evolution is concerned, sulfur may trace helium better than oxygen. However, in our work, we do not have access to the emission lines required to measure the sulfur abundance, and we therefore use the ${\rm{O}}/{\rm{H}}$ abundance in what follows.

Our determination of the primordial helium abundance, y_P, is based on a fit to the measured ${\rm{O}}/{\rm{H}}$ and $\mathrm{He}/{\rm{H}}\,\equiv \,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

$$\begin{eqnarray*}&&Y=\displaystyle \frac{4y(1-Z)}{1+4y}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&Z=c\times ({\rm{O}}/{\rm{H}}).\end{eqnarray*}$$

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:

• 1.  
  Observations of the helium abundance are intrinsically measuring a number abundance ratio.
• 2.  
  Calculations of the helium abundance are computed as a ratio of volume densities (${n}_{{}^{4}\mathrm{He}}/{n}_{{}^{1}{\rm{H}}}$), and later converted to a mass fraction to match the observationally reported mass fractions (see, e.g., Pitrou et al. 2018).
• 3.  
  The primordial helium mass fraction (Y_P) is not actually the fraction of mass in the form of ⁴He. It is defined as the ratio of volume densities ${Y}_{{\rm{P}}}=4\,n{(}^{4}\mathrm{He})/{n}_{{\rm{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.
• 4.  
  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 take c = 20 for Z = 20  × (O/H), Izotov et al. (2007) adopts c = 18.2, and Izotov et al. (2013) allow for a value of c that linearly scales with the metallicity, c = 8.64 × 12 + log₁₀(O/H)−47.44.

For these reasons, we have chosen to quote our primordial helium abundance in the form we most directly measure and the one most appropriate to compare to theoretical values—the primordial helium number abundance ratio, ${y}_{{\rm{P}}}$. However, a comparison of our measured y_P to previously reported values of Y_P can be simply calculated with the following equation:

$$\begin{eqnarray}&&\ {Y}_{{\rm{P}}}=\displaystyle \frac{4{y}_{{\rm{P}}}}{1+4{y}_{{\rm{P}}}}.\end{eqnarray} \tag{ 21 }$$

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:

$$\begin{eqnarray}&&{\rm{log}}({ \mathcal L })=-\displaystyle \frac{1}{2}\sum _{n}\left[\displaystyle \frac{{\left({y}_{n}-{{mx}}_{n}-b\right)}^{2}}{{\sigma }_{{y}_{n}}^{2}+{\sigma }_{\mathrm{intr}}^{2}}-{\rm{log}}({\sigma }_{{y}_{n}}^{2}+{\sigma }_{\mathrm{intr}}^{2})\right].\end{eqnarray} \tag{ 22 }$$

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 ${\rm{O}}/{\rm{H}}$ values, the slope $m\,\equiv$ d $y/{\rm{d}}({\rm{O}}/{\rm{H}})$, and the intercept $b\equiv {y}_{{\rm{P}}}$. We capture the error on our calculated ${\rm{O}}/{\rm{H}}$ abundances by drawing new values of ${\rm{O}}/{\rm{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 ${\sigma }_{{y}_{n}}$. We also introduce the term ${\sigma }_{\mathrm{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 Cooke et al. (2018).

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:

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{c}\begin{array}{c}0\leqslant \displaystyle \frac{{\rm{d}}y}{d\left({\rm{O}}/{\rm{H}}\right)}\leqslant 100\\ 0.06\leqslant {y}_{{\rm{P}}}\leqslant 0.10\\ 0\leqslant {\sigma }_{{\rm{intr}}}\leqslant \mathrm{0.01.}\end{array}\end{array}\end{array}\end{eqnarray*}$$

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, ${\sigma }_{{y}_{n}}$. We find a mean of $\langle \,{\sigma }_{{y}_{{\rm{n}}}}\,\rangle$ = 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 ${\sigma }_{{y}_{{\rm{n}}}}$. 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.

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{c}{y}_{{\rm{P}}}={0.0805}_{-0.0017}^{+0.0017}\\ \displaystyle \frac{{\rm{d}}y}{d({\rm{O}}/{\rm{H}})}={54}_{-16}^{+16}\\ {\sigma }_{{\rm{intr}}}\,\leqslant \,0.0019\,(2\sigma \,{\rm{CL}}).\end{array}\end{array}\end{eqnarray*}$$

This model has a ${\chi }^{2}/\mathrm{dof}=0.77$. For Sample 2, we recover:

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{c}{y}_{{\rm{P}}}={0.0813}_{-0.0013}^{+0.0013}\\ \displaystyle \frac{{\rm{d}}y}{d({\rm{O}}/{\rm{H}})}={40}_{-10}^{+11}\\ {\sigma }_{{\rm{intr}}}={0.0017}_{-0.0005}^{+0.0005},\end{array}\end{array}\end{eqnarray*}$$

with ${\chi }^{2}/\mathrm{dof}=0.82$. While the linear fits to Sample 1 and Sample 2 are comparable, we adopt ${y}_{{\rm{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 is nonzero 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).

Zoom In Zoom Out Reset image size

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 $\eta =(5.931\pm 0.051)\times {10}^{-10}$, which is based on the observationally measured abundance of primordial deuterium, ${({\rm{D}}/{\rm{H}})}_{{\rm{P}}}\,=(2.527\pm 0.030)\times {10}^{-5}$ (Cooke et al. 2018).

Table 6.  Primordial Helium Abundance Results Reported in the Literature

${y}_{{\rm{P}}}$	Observation/Method	Number of Systems	Citation
0.0856 ± 0.0010	H ii region	28	Izotov et al. (2014)
0.0811 ± 0.0018	H ii region	15	Aver et al. (2015)
0.0809 ± 0.0013	H ii region	5	Peimbert et al. (2016)
0.0802 ± 0.0022	H ii region	18	Fernández et al. (2019)
0.0812 ± 0.0011	H ii region in NGC 346	1	Valerdi et al. (2019)
${0.0805}_{-0.0017}^{+0.0017}$	H ii region	54	This work
0.0793 ± 0.011(2σ)	CMB	⋯	Planck Collaboration et al. (2018)
${0.085}_{-0.011}^{+0.015}$	Absorption line system	1	Cooke & Fumagalli (2018)
0.0820 ± 0.000074	SBBN calculation	⋯	Cyburt et al. (2016)
0.0820 ± 0.000075	SBBN calculation	⋯	Pitrou et al. (2018)

Note. A summary of primordial helium abundance results reported in recent literature, the method by which the values are measured or calculated, and their reference. The Planck measurement is the TT, TE, EE + low E value from Equation (80)(a) of Planck Collaboration et al. (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.

Download table as:  ASCIITypeset image

The primordial helium abundance that we report here is in 2.6σ disagreement with the Izotov et al. (2014) result, ${y}_{{\rm{P}}}=0.0856\pm 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 they suggest several possible reasons for the disagreement. Given that our model closely follows that of Aver et al. (2010, 2012, 2013, 2015), we expect many of their reasons to be equally 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, τ_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 modeling. 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₂ and ${t}_{{S}^{++}}$ from t₃. 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 uses this information to decide if a system qualifies. We refer readers to Izotov et al. (2006, 2013, 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 methods of analysis. 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.

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, ${{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}$, and the effective number of neutrino species, ${N}_{\mathrm{eff}}$. By adopting a measurement of ${{\rm{\Omega }}}_{{\rm{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. 2002, 2016; Pitrou 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, as well as of its assumptions.

Constraining the values of ${{\rm{\Omega }}}_{{\rm{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 Cooke et al. (2018),

$$\begin{eqnarray*}&&\begin{array}{c}{Y}_{{\rm{P}}}={0.2436}_{-0.0040}^{+0.0039}\\ {\left({\rm{D}}/{\rm{H}}\right)}_{{\rm{P}}}\times {10}^{5}=2.527\pm 0.030,\end{array}\end{eqnarray*}$$

and use calculations of BBN to infer the values of ${{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}$ and N_eff that best fit these abundances.

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 values of ${{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}$, ${N}_{\nu }$, and τ_n. We restate the formula for predicting Y_P here as an example (see their Equation (145), the surrounding text, and Table 6 of their paper for the values of the C_pqr coefficients that are referenced here):

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{\Delta }}{Y}_{{\rm{P}}}}{{Y}_{{\rm{P}}}}=\sum _{{pqr}}{C}_{{pqr}}{\left(\displaystyle \frac{{\rm{\Delta }}{{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}}{{{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}}\right)}^{p}{\left(\displaystyle \frac{{\rm{\Delta }}{N}_{\nu }}{{N}_{\nu }}\right)}^{q}{\left(\displaystyle \frac{{\rm{\Delta }}{\tau }_{{\rm{n}}}}{{\tau }_{{\rm{n}}}}\right)}^{r}.\end{eqnarray*}$$

We use the latest measurement of the mean neutron lifetime τ_n = 877.7 ± 0.7 (Pattie et al. 2018) to solve for ${{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}$ and N_ν. Furthermore, we use the scaling ${N}_{\mathrm{eff}}={N}_{\nu }\times 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:

$$\begin{eqnarray*}&&\begin{array}{c}0.0185\leqslant {{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}\leqslant 0.0267\\ 1.5\leqslant {N}_{\mathrm{eff}}\leqslant 4.5.\end{array}\end{eqnarray*}$$

With each step, a model set of primordial ${\rm{D}}/{\rm{H}}$ and ${Y}_{{\rm{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\times {n}_{\mathrm{steps}}=$ 1200 steps. Given the observed primordial element abundances, we report the following bounds on the effective number of neutrino species and the baryon density:

$$\begin{eqnarray*}&&\begin{array}{c}{N}_{\mathrm{eff}}={2.85}_{-0.25}^{+0.28}\\ {{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}={0.0215}_{-0.0005}^{+0.0005}.\end{array}\end{eqnarray*}$$

The result of our MCMC calculation is shown in Figure 7, together with the Planck bounds on these parameters.¹⁰ The best-fit value of N_eff is consistent with the value inferred by Planck (${N}_{\mathrm{eff}}={2.92}_{-0.18}^{+0.19};$ shown by the red contours in Figure 7) and the Standard Model value of N_eff = 3.046.

Zoom In Zoom Out Reset image size

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, to push on the lowest-metallicity regime; and (2) to obtain high-quality optical and NIR spectroscopy of a subset of the new, metal-poor galaxies, with priority on systems with metallicities determined to be $12+{{\rm{log}}}_{10}({\rm{O}}/{\rm{H}})\leqslant \,7.65$ based on strong-line calibration methods. The purpose of these specific goals was to better populate and constrain the metal-poor 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 subpercent-level accuracy.

5.2.1. Qualification Rates

5.2.2. Toward a Subpercent-level Measurement of ${y}_{{\rm{P}}}$