PRECISION MEASURES OF THE PRIMORDIAL ABUNDANCE OF DEUTERIUM^*

Most newly discovered DLAs are now identified using the low-resolution spectra provided by the Sloan Digital Sky Survey (SDSS). Recent searches for DLAs have yielded several thousand systems (Prochaska & Wolfe 2009; Noterdaeme et al. 2012). At the modest resolution (R ∼ 2000) and signal-to-noise ratio (S/N ∼ 20) of the SDSS spectra, the metal absorption lines of the most metal-poor DLAs are unresolved and undetected. Follow-up, high spectral resolution observations are thus required to pin down their chemical abundances (see Cooke et al. 2011).

Penprase et al. (2010) were the first to recognize the very metal-poor DLA at redshift z_abs ≃ 3.0674 toward the z_em ≃ 3.173 quasar SDSS J1358+6522. On the basis of their intermediate-resolution observations (60 km s⁻¹ FWHM) with the Keck Echellette Spectrograph and Imager, these authors concluded that this DLA has a metallicity of [Fe/H] ⩽−3.02 and is thus among the most metal-poor known (Cooke et al. 2011). We re-observed this QSO with the Keck High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) using the red cross-disperser to perform an accurate chemical abundance analysis (Cooke et al. 2012). The HIRES data pinned down the metallicity of the DLA ([Fe/H] =−2.84 ± 0.03) and provided the first indication for the presence of a handful of resolved D i absorption lines.

Taking advantage of the high efficiency of the HIRES UV cross-disperser at blue wavelengths, we then performed dedicated HIRES observations with the goal of obtaining a precise measure of the deuterium abundance. Our observations included a total of 24, 300 s on 2012 May 11, (divided equally into nine exposures) and an additional 29, 750 s on 2013 May 3 and 5. An identical instrument setup was employed for both runs; we used the C1 decker (7'' × 0861, well matched to the seeing conditions of 06–08), which delivers a nominal spectral resolution of R ∼ 48, 000 (≡ 6.3 km s⁻¹ FWHM) for a uniformly illuminated slit. Our setup covers the wavelength range 3650–6305 Å, with small gaps near 4550 Å and 5550 Å due to the gaps between the HIRES detector chips. The exposures were binned on-chip 2 × 2.

The data were reduced following the standard procedures of bias subtraction, flat-fielding, order definition, and extraction. The data were mapped onto a vacuum, heliocentric wavelength scale from observations of a ThAr lamp that bracketed the science exposures. To test the quality of the data reduction, we reduced the data using two pipelines: makee, which is maintained by T. Barlow,⁸ and HIRedux, which is maintained by J. X. Prochaska.⁹ We found that the skyline subtraction performed by HIRedux was superior to that of makee and substantially improved the S/N of the spectra at blue wavelengths. This was very important for the removal of telluric features near the redshifted D i and H i transitions of the DLA. All data analyzed and presented herein for J1358+6522 were reduced with HIRedux. The individual orders were combined with UVES_popler, which is maintained by M. T. Murphy.¹⁰ As a final step, we flux-calibrated the echelle data by comparison to the SDSS data. The combined S/N of the final spectrum per 2.5 km s⁻¹ pixel is ∼30 near 4000 Å, ∼45 at 5000 Å, and ∼20 near 6000 Å. The reduced, flux-calibrated spectrum is presented in Figure 1, where the red tick marks indicate the wavelengths of the redshifted H i Lyman series transitions.

Zoom In Zoom Out Reset image size

2.2.1. Selection Criteria

2.2.2. HS0105+1619, z_abs = 2.53651

2.2.3. Q0913+072, z_abs = 2.61829

2.2.4. J1419+0829, z_abs = 3.04973

2.2.5. J1558−0031, z_abs = 2.70242

To account for the dominant systematics that could potentially affect the measurement of D/H, we have developed a new suite of software routines written in the python software environment. Our Absorption LIne Software (alis; described in more detail by R. J. Cooke, in preparation) uses a modified version of the mpfit package (Markwardt 2009). mpfit employs a Levenberg–Marquardt least-squares minimization algorithm to derive the model parameters that best fit the data (i.e., the parameters that minimize the difference between the data and the model, weighted by the error in the data). Unlike the software packages that have previously been used to measure D/H, alis has the advantage of being able to fit an arbitrary emission profile for the quasar while simultaneously fitting to the absorption from the DLA; any uncertainty in the continuum is therefore automatically folded into the final uncertainty in the D/H ratio. For all transitions, we use the atomic data compiled by Morton (2003).

We first highlight one of the main strengths of our analysis that, to our knowledge, has not been implemented in previous D/H studies. In order to remove human bias from the results, we have adopted a blind analysis strategy, whereby the only information revealed during the profile-fitting process is the best-fitting χ² value, and plots of the corresponding model fits to the data; all of the parameter fitting-results are entirely hidden from view. Only after the analysis converged on the best-fitting model were the full results of the model fit uncovered, and these are the final values quoted here. No further changes were made to the model fit nor data reduction after the results were revealed.

To allow for an efficient computation, a small wavelength interval (typically ±300 km s⁻¹) around each absorption line of interest was extracted. For each of these intervals, we manually identified the regions to be used in the fitting process by selecting pixels that we attribute to either the quasar continuum or the DLA's absorption lines. In some cases, it was also necessary to model unrelated absorption lines in the fitting procedure if there was significant line blending.

We have also attempted to account for the main systematics and limitations that might affect the reported D/H value. The dominant systematics likely include the instrumental broadening function, the cloud model, and the placement of the quasar continuum and zero level. We marginalize over the first of these systematics by assuming that the instrumental profile is a Gaussian,¹² where the Gaussian FWHM was allowed to be a free model parameter during the χ² minimization process.

Any uncertainty or bias in the choice of the cloud model is overcome by our selection criteria listed in Section 2.2.1. Specifically, DLAs have two advantages over other absorption-line systems for the determination of the D/H ratio, both stemming from their high H i column densities: (1) the H i Lyα line exhibits well-defined Lorentzian damping wings¹³ from which the total H i column density can be deduced precisely and independently of the cloud model, and (2) many D i transitions in the Lyman series are normally detectable. With their widely varying oscillator strengths, these multiple D i lines can pin down the D i column density very precisely; the availability of unsaturated, high-order lines ensures that the value of D/H thus deduced does not depend on the detailed kinematic structure of the gas (i.e., the cloud model). Neither of these two advantages applies to Lyman limit systems whose H i column density can be three orders of magnitude lower than that of DLAs.

Finally, if the background level is not accurately determined, a small offset can systematically affect the derived D/H value. We have overcome this concern by allowing the zero level to be fit during the minimization process (assuming the offset to be the same at all wavelengths).

Arguably, the above systematic effects are sub-dominant compared to the choice of the continuum level. One of the improvements over previous analyses is that we have not specified a priori the level of the quasar continuum. Instead, we have marginalized over the continuum level during the fitting process by fitting Legendre polynomials to the flux-calibrated quasar continuum. The order of the polynomial was chosen to be sufficiently large to allow flexibility in arbitrarily defining the continuum level, while ensuring that the continuum was not being overfit to the noise (examples are shown in the top panels of Figures 2, 8, 9, and 10).

Zoom In Zoom Out Reset image size

We emphasize that by including the above systematics into the fitting process, the uncertainties in the continuum, zero level, and instrumental broadening are incorporated into the final uncertainty in the derived D/H value. To summarize, we have marginalized over the continuum, zero level, and instrumental FWHM, while simultaneously fitting for all of the DLA's unblended metal lines, the DLA's Lyα transition, and the entire H i and D i Lyman series. Our analysis therefore includes all of the available information that can help pin down the D/H value, and it accounts for the systematics that are currently thought to contribute the most to the overall error budget. We finally point out that every H i component used to fit the DLA was forced to have the same D i/H i ratio. We have therefore fit directly to the D i/H i value, rather than deriving separately the column densities of D i and H i. All absorption lines were modeled as Voigt profiles where the line broadening has both thermal and turbulent contributions. All of the DLAs in our survey were first modeled with a single absorption component. If the metal lines exhibited an asymmetric profile, an additional component was included to improve the quality of the fit.

To confirm that the model parameters had converged to their best-fitting values, we performed a convergence test; once the difference in the χ² between successive iterations was <0.01, the parameter values were stored and the χ² minimization recommenced with a tolerance of 10⁻³. Once a successive iteration reduced the χ² by <10⁻³, the minimization was terminated and the parameter values from the two convergence criteria were compared. If all parameter values differed by <0.2σ (i.e., 20% of the parameter error), then the model fit has converged.

As a final step, we repeated the χ² minimization process 20 times, perturbing the starting parameters of each run by the covariance matrix. This exercise ensures that our choice of starting parameters does not influence the final result. We found that the choice of starting parameters has a negligible contribution to the error on D i/H i (typically 0.002 dex), but can introduce a small bias (again, typically 0.002 dex). We have accounted for this small bias in all of the results quoted herein.

Most of the narrow, low-ionization metal lines of the DLA toward J1358+6522 consist of a single component at z_abs = 3.067259. A second weaker component, blueshifted by 17.4 km s⁻¹ (z_abs = 3.06702), contributes to Si iii λ1206.5 and to the strongest C ii and Si ii lines. Evidently, this weaker absorption arises in nearby ionized gas.

In fitting the absorption lines, we tied the redshift, turbulent Doppler parameter, and kinetic temperature of the gas to be the same for the metal, D i, and H i absorption lines. We allowed all of the cloud model parameters to vary, while simultaneously fitting for the continuum near every absorption line. Relevant parameters of the best-fitting cloud model so determined are collected in Table 1. Figures 2, 3, and 4 compare the data and model fits for, respectively, the damped Lyα line, the full Lyman series, and selected metal lines. [Since the metal lines analyzed here are the same as those shown in Figure 1 of Cooke et al. (2012), albeit now with a higher S/N, we only present a small selection of them in Figure 4 to avoid repetition]. The best-fitting chi-squared value for this fit is also provided for completeness.¹⁴

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Returning to Table 1, it can be seen that we found it necessary to separate the main absorption into two separate components, labeled 1a and 1b in the table. A statistically acceptable fit¹⁵ to the metal, D i and H i lines could not be achieved with a single absorbing cloud in which the turbulent broadening is the same for all species and the thermal broadening is proportional to the square root of the ion mass (i.e., $b_{\rm th}^2 = 2KT/m$, where K is the Boltzmann constant). The main absorption component of this DLA appears to consist of two "clouds" with very similar redshifts, temperatures, and H i column densities, but with significantly different turbulence parameters (see Table 1). The turbulent broadening for component 1a is bounded by the metal lines, whereas the thermal broadening is bounded by the relatively narrow H i line profiles. This combination of turbulent and thermal broadening is unable to reproduce the observed widths of the strongest D i lines, which require an additional component with a larger contribution of turbulent broadening. Surprisingly, metal absorption is only seen in the low-turbulence cloud (component 1a); component 1b has less than 1/100 of the metallicity of component 1a. We propose three possible interpretations for this unusual finding: (1) there exists an essentially pristine cloud of gas near the very metal-poor DLA; (2) the metals see a different potential to D i and H i (i.e., the metals have not been entirely mixed into the H i); (3) the assumption of a constant kinetic temperature for all elements through the entire sightline is incorrect.

While each of these explanations is interesting in its own right, we stress here that the choice of the cloud model does not affect the derived value of D/H in this system, for the reasons discussed in Section 3.2. Provided that there exists a series of optically thin D i lines, and H i Lyα exhibits damping wings (as indeed is the case for the DLA toward J1358+6522—see Figures 2 and 3), the derived total D i and H i column densities will be independent of the cloud model. We therefore only require that the model provides a good fit to the data.

The final, best-fitting value of the deuterium abundance for the DLA toward J1358+6522 is (expressed as a logarithmic and linear quantity)

$$\begin{equation} \log \,({\rm D}\, {\scriptsize I}/{\rm H}\,{\scriptsize I}) = -4.588 \pm 0.012, \end{equation} \tag{ 1 }$$

$$\begin{equation} 10^5\,{\rm D}\,{\scriptsize I}/{\rm H}\,{\scriptsize I} = 2.58 \pm 0.07, \end{equation} \tag{ 2 }$$

where the quoted error term includes the random (observational) error in addition to the systematic uncertainties that were marginalized over during the fitting process.

Using the most up-to-date calculations of the network of nuclear reactions involved in BBN, the primordial abundance of deuterium is related to the cosmic density of baryons (in units of the critical density), Ω_{b, 0}, via the following relations (Steigman 2012; G. Steigman 2013, private communication):

$$\begin{equation} {\rm (D\,/\,H)}_{\rm p} = 2.55\times 10^{-5}\,(6/\eta _{\rm D})^{1.6}\times (1\pm 0.03), \end{equation} \tag{ 5 }$$

$$\begin{equation} \eta _{\rm D} = \eta _{\rm 10} - 6(S-1) + 5\xi /4, \end{equation} \tag{ 6 }$$

where η₁₀ = 273.9 Ω_{b, 0} h², S = [1 + 7(N_eff − 3.046)/43]^1/2 is the expansion factor, and ξ is the neutrino degeneracy parameter (related to the lepton asymmetry by Equation (14) from Steigman 2012). The rightmost term in Equation (5) represents the current 3% uncertainty in the conversion of (D/H)_p to η_D due to the uncertainties in the relevant nuclear reaction rates (see Section 4.2). For the standard model, N_eff ≃ 3.046 and ξ = 0. In this case, the Precision Sample of D/H measurements implies a cosmic density of baryons:

$$\begin{eqnarray} 100\,\Omega _{\rm b,0}\,h^2 ({\rm BBN}) &=& 2.202\pm 0.020\; ({\rm random}) \nonumber \\ && \pm 0.041\; ({\rm systematic}), \end{eqnarray} \tag{ 7 }$$

where we have decoupled the error terms from our measurement (i.e., the random error term) and the systematic uncertainty in converting the D abundance into the baryon density parameter.

As can be seen from Figure 5, this value of Ω_{b, 0} h² is in excellent agreement with that derived from the analysis of the CMB temperature fluctuations measured by the Planck satellite (Planck Collaboration 2013):

$$\begin{equation} 100 \, \Omega _{\rm b,0} \, h^2 {\rm (CMB)} =2.205 \pm 0.028. \end{equation} \tag{ 8 }$$

In the era of high-precision cosmology, we feel that it is important to highlight the main limitations affecting the use of (D/H)_p in the estimation of cosmological parameters. As can be seen from Equation (7), the main source of error is in the conversion of (D/H)_p to the baryon density parameter (η_D, and hence Ω_{b, 0} h²). In large part, this systematic uncertainty is due to the relative paucity of experimental measures for several nuclear cross sections that are important in the network of BBN reactions, particularly deuteron–deuteron reactions and the d(p, γ)³He reaction rate at the relevant energies (Fiorentini et al. 1998; Nollett & Burles 2000; Cyburt 2004; Serpico et al. 2004). Since these studies, estimates for the deuteron–deuteron reaction cross sections (Leonard et al. 2006) have improved and their contribution to the error budget has been reduced. The main lingering concern involves the reaction rate d(p, γ)³He, for which only a single reliable data set of the S-factor is currently available in the relevant energy range (Ma et al. 1997).¹⁶ The concern for d(p, γ)³He is made worse by the fact that the theoretical and experimental values of the S-factor do not agree. This paucity of data, in addition to the difficulties of obtaining an accurate and precise measure of d(p, γ)³He at BBN energies, led Nollett & Holder (2011) to adopt a theoretical curve for the S-factor between 50 and 500 keV. This has resulted in the improved conversion given in Equation (5).

It is worth pointing out here that a further reduction by a modest factor of two in the uncertainty of the conversion from (D/H)_p to η_D would make the precision of Ω_{b, 0} h²(BBN) from the DLA Precision Sample comparable to that of Ω_{b, 0} h²(CMB) achieved by the Planck mission (see Equations (7) and (8)).

It has long been known that the mass fraction of ⁴He, Y_P, is potentially a very sensitive probe of additional light degrees of freedom (Steigman et al. 1977). Unfortunately, systematic uncertainties in the determination of Y_P have limited its use as an effective probe of N_eff (see Figure 8 of Steigman 2012). However, Nollett & Holder (2011) (see also Cyburt 2004) have recently highlighted the potential of using precise measurements of the primordial deuterium abundance in conjunction with observations of the CMB to place a strong, independent bound on N_eff.

In the left panel of Figure 6, we show the current 1σ and 2σ confidence contours for N_eff and Ω_{b, 0} h² derived from the Planck+WP+highL CMB analysis¹⁷ (green) and from the BBN-derived (D/H)_p reported here (blue). The combined confidence bounds are overlaid as red contours. In what follows, it is instructive to recall that the CMB-only bounds are (Planck Collaboration 2013)

$$\begin{equation} 100\,\Omega _{\rm b,0}\,h^2 = 2.23\pm 0.04, \end{equation} \tag{ 9 }$$

$$\begin{equation} N_{\rm eff} = 3.36\pm 0.34 . \end{equation} \tag{ 10 }$$

(Note that solving simultaneously for Ω_{b, 0} h² and N_eff leads to a slightly different best-fitting value of Ω_{b, 0} h² than that obtained for the standard model; see Equations (8) and (9)). For comparison, from the joint BBN+CMB analysis we deduce

$$\begin{equation} 100\,\Omega _{\rm b,0}\,h^2 = 2.23\pm 0.04, \end{equation} \tag{ 11 }$$

$$\begin{equation} N_{\rm eff} = 3.28 \pm 0.28 . \end{equation} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

Thus, combining (D/H)_p with the CMB does not significantly change the uncertainty in Ω_{b, 0} h², but does reduce the error on N_eff by ∼20%. The uncertainty on N_eff could be reduced further by an improvement in the cross section of the d(p, γ)³He (see right panel of Figure 6, and Section 4.2). Based on the current bound on N_eff from CMB+(D/H)_p, we can nevertheless rule out the existence of an additional (sterile) neutrino (i.e., N_eff = 4.046) at 99.3% confidence (i.e., ∼2.7σ), provided that N_eff and η₁₀ remained unchanged between BBN and recombination. However, as noted recently by Steigman (2013), if the CMB photons are heated after the neutrinos have decoupled [for example, by a weakly interacting massive particle that annihilates to photons], N_eff will be less than 3.046 for three standard model neutrinos; a sterile neutrino can in principle exist even when N_eff <4.046.

Looking to the future, Y_P has contours that are almost orthogonal to those of the CMB and (D/H)_p (see, e.g., Steigman 2007). Thus, measures of Y_P that are not limited by systematic uncertainties could potentially provide a very strong bound, when combined with (D/H)_p, on the number of equivalent neutrinos during the epoch of BBN, independently of CMB observations.

Following improvements in the He i emissivity calculations (Porter et al. 2012, 2013), there have been two recent reassessments of Y_P, by Izotov et al. (2013) and by Aver et al. (2013), respectively. Izotov et al. (2013) find Y_P = 0.253 ± 0.003; this value includes a small correction of −0.001 to reflect the fact that their Cloudy modeling overpredicts Y_P by this amount at the lowest metallicities—see their Figure 7(a). Aver et al. (2013) find Y_P = 0.254 ± 0.004 from the average of all the low-metallicity H ii regions in their sample. Thus, the values deduced by these two teams are in good mutual agreement; in the analysis that follows we adopt Y_P = 0.253 ± 0.003 from Izotov et al. (2013). Using the following conversion relation for Y_P (Steigman 2012; G. Steigman 2013, private communication):

$$\begin{equation} {Y}_{\rm P} = 0.2469 \pm 0.0006 + 0.0016\,(\eta _{\rm He} - 6), \end{equation} \tag{ 13 }$$

$$\begin{equation} \eta _{\rm He} = \eta _{\rm 10} + 100(S-1) - 575\xi /4, \end{equation} \tag{ 14 }$$

combined with the Izotov et al. (2013) measure of Y_P, we derive the following BBN-only bound on the baryon density and the effective number of neutrino species:

$$\begin{equation} 100\,\Omega _{\rm b,0}\,h^2 = 2.28\pm 0.05, \end{equation} \tag{ 15 }$$

$$\begin{equation} N_{\rm eff} = 3.50 \pm 0.20 . \end{equation} \tag{ 16 }$$

The corresponding contours are shown in Figure 7. Thus, it appears that even with the most recent reappraisals of the primordial abundance of ⁴He by Izotov et al. (2013) and Aver et al. (2013), there is better agreement (within the standard model) between (D/H)_p and the CMB than between (D/H)_p and Y_P.

Zoom In Zoom Out Reset image size

In the past, the primordial deuterium abundance has been commonly used as a tool for measuring the present-day universal density of baryons (see, e.g., Steigman 2007), and more recently as a probe of the effective number of neutrino families (Cyburt 2004; Nollett & Holder 2011; Pettini & Cooke 2012; see also Section 5.1). Here, we demonstrate that precise measures of the primordial deuterium abundance (in combination with the CMB) can also be used to estimate the neutrino degeneracy parameter, ξ, which is related to the lepton asymmetry by Equation (14) from Steigman (2012).

Steigman (2012) recently suggested that combined estimates for (D/H)_p, Y_P, and a measure of N_eff from the CMB can provide interesting limits on the neutrino degeneracy parameter (ξ ⩽ 0.079, 2σ; see also Serpico & Raffelt 2005; Popa & Vasile 2008; Simha & Steigman 2008). By combining (D/H)_p and Y_P, this approach effectively removes the dependence on Ω_{b, 0} h². Using the conversion relations for (D/H)_p and Y_P (Equations (5)–(6) and 13–14) and the current best determination of Y_P (0.253 ± 0.003; Izotov et al. 2013), in addition to the Planck+WP+highL¹⁸ constraint on N_eff and the precise determination of (D/H)_p reported here, we derive a 2σ upper limit on the neutrino degeneracy parameter, |ξ| ⩽ 0.064, based on the approach by Steigman (2012).

We propose that an equally powerful technique for estimating ξ does not involve removing the dependence on Ω_{b, 0} h² by combining (D/H)_p and Y_P, as in Steigman (2012). Instead, one can obtain a measure of both Ω_{b, 0} h² and N_eff from the CMB and use either (D/H)_p or Y_P to obtain two separate measures of ξ. This has the clear advantage of decoupling (D/H)_p and Y_P; any systematic biases in either of these two values could potentially bias the measure of ξ. Separating (D/H)_p and Y_P also allows one to check that the two estimates agree with one another.

Our calculation involved a Monte Carlo technique, whereby we generated random values from the Gaussian-distributed primordial D/H abundance measurements, while simultaneously drawing random values from the (correlated) distribution between Ω_{b, 0} h² and N_eff from the Planck+WP+highL CMB data (Planck Collaboration 2013).¹⁹ Using Equation (19) from Steigman (2012, equivalent to Equation (6) here), we find ξ_D = +0.05 ± 0.13 for (D/H)_p, leading to a 2σ upper limit of |ξ_D| ⩽ 0.31.

With the technique outlined above, we have also computed the neutrino degeneracy parameter from the current observational bound on Y_P. For this calculation, we have used the MCMC chains from the Planck+WP+highL CMB base cosmology with N_eff and Y_P added as free parameters. In this case, the CMB distribution was weighted by the observational bound on Y_P (Y_P = 0.253 ± 0.003; Izotov et al. 2013). Using Equations (19)–(20) from Steigman (2012, equivalent to Equations (6) and (14) here), we find ξ_D = +0.04 ± 0.15 for (D/H)_p and ξ_He = −0.010 ± 0.027 for Y_P. These values translate into corresponding 2σ upper limits |ξ_D| ⩽ 0.34 and |ξ_He| ⩽ 0.064. Combining these two constraints then gives ξ = −0.008 ± 0.027, or |ξ| ⩽ 0.062 (2σ).

Alternatively, if we assume that the effective number of neutrino species is consistent with three standard model neutrinos (i.e., N_eff ≃ 3.046), we obtain the following BBN-only bound on the neutrino degeneracy parameter by combining (D/H)_p and Y_P, ξ = −0.026 ± 0.015, or |ξ| ⩽ 0.056 (2σ). We therefore conclude that all current estimates of the neutrino degeneracy parameter, and hence the lepton asymmetry, are consistent with the standard model value, ξ = 0.

From the above calculations, it is clear that Y_P is the more sensitive probe of any lepton asymmetry, as is already appreciated. However, (D/H)_p offers an additional bound on ξ that is complementary to that obtained from Y_P. We note further that, if the uncertainty in the conversion of (D/H)_p to η_D could be reduced by a factor of two, the bound on ξ_D would be reduced by 15%, corresponding to a 1σ uncertainty on ξ of ∼0.11 from the CMB and deuterium measurements alone.