A Bayesian direct method implementation to fit emission line spectra:   Application to the primordial He abundance determination

Vital Fern\'andez; Elena Terlevich; Angeles I. D\'iaz; Roberto; Terlevich

arXiv:1905.09215·astro-ph.GA·June 5, 2019

A Bayesian direct method implementation to fit emission line spectra: Application to the primordial He abundance determination

Vital Fern\'andez, Elena Terlevich, Angeles I. D\'iaz, Roberto, Terlevich

PDF

TL;DR

This paper introduces a Bayesian algorithm using probabilistic programming to fit emission line spectra, enabling precise primordial helium abundance measurement consistent with Big Bang predictions.

Contribution

The work develops a novel Bayesian method with a 14-dimensional model for emission spectra fitting, improving analysis speed and accuracy over traditional approaches.

Findings

01

Primordial helium abundance measured as 0.243±0.005

02

Method shows excellent convergence and speed with PyMC3

03

Results align with standard Big Bang cosmology

Abstract

This work presents a Bayesian algorithm to fit the recombination and collisionally excited line spectra of gas photoionized by clusters of young stars. The current model consists in fourteen dimensions: two electron temperatures, one electron density, the extinction coefficient, the optical depth on the $H e I$ recombination lines and nine ionic species. The results are in very good agreement with those previously published using the traditional methodology. The probabilistic programming library PyMC3 was chosen to explore the parameter space via a NUTs sampler. These machine learning tools provided excellent convergence quality and speed. The primordial helium abundance measured from a multivariable regression using oxygen, nitrogen and sulfur was $Y_{P, O - N - S} = 0.243 \pm 0.005$ in agreement with a standard Big Bang scenario.

Tables12

Table 1. Table 1: Priors and likelihood distributions in our model. The term X i + superscript 𝑋 limit-from 𝑖 X^{i+} includes all the ionic metal abundances: A r 2 + 𝐴 superscript 𝑟 limit-from 2 Ar^{2+} , A r 3 + 𝐴 superscript 𝑟 limit-from 3 Ar^{3+} , y + superscript 𝑦 y^{+} , y 2 + superscript 𝑦 limit-from 2 y^{2+} , O + superscript 𝑂 O^{+} , O 2 + superscript 𝑂 limit-from 2 O^{2+} , N + superscript 𝑁 N^{+} , S + superscript 𝑆 S^{+} and S 2 + superscript 𝑆 limit-from 2 S^{2+} . These abundances are define in 12 + l o g ( X i + ) 12 𝑙 𝑜 𝑔 superscript 𝑋 limit-from 𝑖 12+log\left(X^{i+}\right) scale.

Parameter	Prior distribution
$T_{l o w}$	$N o r m a l (μ = 15000 K, σ = 5000 K)$
$T_{h i g h}$	$N o r m a l (μ = 15000 K, σ = 5000 K)$
$n_{e}$	$N o r m a l (μ = 150 c m^{- 3}, σ = 50 c m^{- 3})$
$c (H β)$	$l o g N o r m a l (μ = 0, σ = 1)$
$X^{i +}$	$N o r m a l (μ = 5, σ = 5)$
$y^{+}$	$1 / 10 \cdot l o g N o r m a l (μ = 0, σ = 1)$
$y^{2 +}$	$1 / 1000 \cdot l o g N o r m a l (μ = 0, σ = 1)$
$τ$	$l o g N o r m a l (μ = 0, σ = 0.4)$
Parameter	Likelihood distribution
$\frac{F_{X^{i +}, λ}}{F_{H β}}$	$N o r m a l (μ = \frac{F_{X^{i +}, λ, o b s}}{F_{H β}}, σ = \frac{σ_{X^{i +}, λ, o b s}}{F_{H β}})$

Table 2. Table 2: Fitting results for the test cases with increasing number of parameters described in the text

Parameter	True value	Test 1	Test 2	Test 3	Test 4
$T_{l o w}$	15590	$15590 \pm 244$	$15750 \pm 252$	$16470 \pm 787$	$15640 \pm 282$
$n_{e}$	500	$491 \pm 39$	$492 \pm 39$	$494 \pm 39$	$492 \pm 36$
$S^{+}$	5.48	$5.479 \pm 0.014$	$5.459 \pm 0.028$	$5.402 \pm 0.064$	$5.475 \pm 0.019$
$S^{2 +}$	6.36	$6.360 \pm 0.015$	$6.332 \pm 0.036$	$6.258 \pm 0.084$	$6.354 \pm 0.022$
$O^{+}$	7.80	✗	$7.768 \pm 0.047$	$7.663 \pm 0.116$	$7.794 \pm 0.037$
$O^{2 +}$	8.05	✗	$8.037 \pm 0.022$	$8.027 \pm 0.026$	$8.048 \pm 0.017$
$A r^{2 +}$	5.72	✗	$5.696 \pm 0.031$	$5.634 \pm 0.071$	$5.715 \pm 0.019$
$A r^{3 +}$	5.06	✗	$5.049 \pm 0.019$	$5.043 \pm 0.021$	$5.058 \pm 0.016$
$N^{+}$	5.84	✗	$5.820 \pm 0.027$	$5.762 \pm 0.064$	$5.835 \pm 0.019$
$c (H β)$	0.100	✗	$0.137 \pm 0.047$	$0.21 \pm 0.09$	$0.11 \pm 0.018$
$T_{h i g h}$	16000	✗	✗	$16310 \pm 334$	$16020 \pm 193$
$y^{+}$	0.0850	✗	✗	✗	$0.0850 \pm 0.001$
$y^{2 +}$	0.00088	✗	✗	✗	$0.00088 \pm 0.00001$
$τ$	1.0	✗	✗	✗	$0.991 \pm 0.225$

Table 3. Table 3: Sample fitting results. The values displayed correspond to the mean and standard deviation for n e subscript 𝑛 𝑒 n_{e} , the low and high ionization region T e subscript 𝑇 𝑒 T_{e} and the extinction coefficient.

	$n_{e} [S I I]$	$T_{e} [S I I I]$	$T_{e} [O I I I]$
HII Galaxy	$(c m^{- 3})$	$(K)$	$(K)$	$c (H β)$
FTDTR-1	42.1 $\pm$ 20	14500 $\pm$ 740	-	0.233 $\pm$ 0.024
IZw18	31.1 $\pm$ 20	19000 $\pm$ 1300	-	0.0567 $\pm$ 0.019
MRK36-A1	30.4 $\pm$ 20	14900 $\pm$ 610	14700 $\pm$ 780	0.0743 $\pm$ 0.022
MRK36-A2	64.1 $\pm$ 27	13000 $\pm$ 290	15900 $\pm$ 320	0.0516 $\pm$ 0.019
MRK475	36.5 $\pm$ 19	14100 $\pm$ 390	14600 $\pm$ 320	0.0787 $\pm$ 0.022
FTDTR-2	39.5 $\pm$ 22	-	13900 $\pm$ 680	0.109 $\pm$ 0.025
IZw70	38.4 $\pm$ 20	11200 $\pm$ 500	-	0.0873 $\pm$ 0.023
MRK689	41.8 $\pm$ 20	10300 $\pm$ 390	-	0.266 $\pm$ 0.028
MRK67	70.5 $\pm$ 31	12900 $\pm$ 470	-	0.238 $\pm$ 0.027
FTDTR-3	45.7 $\pm$ 19	13000 $\pm$ 430	13000 $\pm$ 320	0.101 $\pm$ 0.023
SHOC022	37.0 $\pm$ 20	11600 $\pm$ 600	12300 $\pm$ 610	0.123 $\pm$ 0.026
FTDTR-4	45.0 $\pm$ 20	12400 $\pm$ 590	15100 $\pm$ 310	0.186 $\pm$ 0.024
FTDTR-5	146 $\pm$ 86	-	13600 $\pm$ 1600	0.122 $\pm$ 0.027
FTDTR-6	45.9 $\pm$ 20	12700 $\pm$ 500	13000 $\pm$ 380	0.181 $\pm$ 0.025
FTDTR-7	133 $\pm$ 18	13300 $\pm$ 350	13900 $\pm$ 210	0.289 $\pm$ 0.023
MRK627	44.2 $\pm$ 21	10500 $\pm$ 470	11400 $\pm$ 760	0.151 $\pm$ 0.027
PHL293B	113 $\pm$ 27	15600 $\pm$ 640	16600 $\pm$ 230	0.225 $\pm$ 0.027
FTDTR-8	49.7 $\pm$ 22	-	13100 $\pm$ 1500	0.151 $\pm$ 0.029
SHOC263	233 $\pm$ 60	-	11000 $\pm$ 340	0.0685 $\pm$ 0.021
FTDTR-9	145 $\pm$ 48	11800 $\pm$ 870	11800 $\pm$ 740	0.208 $\pm$ 0.028
FTDTR-10	582 $\pm$ 80	-	16500 $\pm$ 630	0.0850 $\pm$ 0.023

Table 4. Table 4: Helium, oxygen, and nitrogen ionic abundances.

HII Galaxy	$H e^{+} / H^{+}$	$H e^{2 +} / H^{+}$	$12 + l o g (O^{+} / H^{+})$	$12 + l o g (O^{2 +} / H^{+})$	$12 + l o g (N^{+} / H^{+})$
FTDTR-1	0.0975 $\pm$ 0.0018	0.00117 $\pm$ 0.00017	7.22 $\pm$ 0.11	7.72 $\pm$ 0.060	5.31 $\pm$ 0.056
IZw18	0.0752 $\pm$ 0.0018	0.000527 $\pm$ 0.00	6.59 $\pm$ 0.11	7.00 $\pm$ 0.065	4.79 $\pm$ 0.057
MRK36-A1	0.0767 $\pm$ 0.0018	0.00180 $\pm$ 0.00013	7.27 $\pm$ 0.084	7.72 $\pm$ 0.060	5.56 $\pm$ 0.048
MRK36-A2	0.0784 $\pm$ 0.0015	0.00124 $\pm$ 0.00	7.60 $\pm$ 0.052	7.66 $\pm$ 0.025	5.68 $\pm$ 0.039
MRK475	0.0821 $\pm$ 0.0018	0.00171 $\pm$ 0.00	7.28 $\pm$ 0.062	7.79 $\pm$ 0.028	5.68 $\pm$ 0.042
FTDTR-2	0.0750 $\pm$ 0.0022	0.000963 $\pm$ 0.00024	-	7.84 $\pm$ 0.061	5.53 $\pm$ 0.059
IZw70	0.0939 $\pm$ 0.0020	0.000570 $\pm$ 0.00014	7.99 $\pm$ 0.11	7.96 $\pm$ 0.066	6.37 $\pm$ 0.050
MRK689	0.0795 $\pm$ 0.0018	-	8.24 $\pm$ 0.10	8.19 $\pm$ 0.060	6.28 $\pm$ 0.047
MRK67	0.0836 $\pm$ 0.0020	0.000964 $\pm$ 0.00	7.63 $\pm$ 0.083	7.95 $\pm$ 0.048	5.84 $\pm$ 0.049
FTDTR-3	0.0833 $\pm$ 0.0019	-	7.33 $\pm$ 0.077	8.09 $\pm$ 0.032	5.59 $\pm$ 0.046
SHOC022	0.0833 $\pm$ 0.0020	0.00134 $\pm$ 0.00014	7.83 $\pm$ 0.13	7.94 $\pm$ 0.065	6.21 $\pm$ 0.057
FTDTR-4	0.0865 $\pm$ 0.0017	0.000744 $\pm$ 0.00	7.73 $\pm$ 0.11	7.82 $\pm$ 0.025	6.00 $\pm$ 0.050
FTDTR-5	0.0764 $\pm$ 0.0040	-	-	8.04 $\pm$ 0.10	-
FTDTR-6	0.0821 $\pm$ 0.0017	0.000884 $\pm$ 0.00	7.58 $\pm$ 0.094	7.96 $\pm$ 0.038	5.76 $\pm$ 0.051
FTDTR-7	0.0840 $\pm$ 0.0015	0.000407 $\pm$ 0.00	7.39 $\pm$ 0.062	8.03 $\pm$ 0.021	5.69 $\pm$ 0.029
MRK627	0.0869 $\pm$ 0.0022	-	8.16 $\pm$ 0.12	8.00 $\pm$ 0.094	6.52 $\pm$ 0.053
PHL293B	0.0669 $\pm$ 0.0018	0.00174 $\pm$ 0.00017	6.89 $\pm$ 0.086	7.62 $\pm$ 0.019	-
FTDTR-8	0.0737 $\pm$ 0.0035	-	7.30 $\pm$ 0.27	7.85 $\pm$ 0.040	5.32 $\pm$ 0.16
SHOC263	0.0874 $\pm$ 0.0042	-	8.04 $\pm$ 0.093	8.01 $\pm$ 0.047	6.59 $\pm$ 0.045
FTDTR-9	0.0876 $\pm$ 0.0032	0.00160 $\pm$ 0.00013	7.82 $\pm$ 0.18	8.01 $\pm$ 0.086	6.21 $\pm$ 0.080
FTDTR-10	0.0736 $\pm$ 0.0031	0.000976 $\pm$ 0.00014	6.98 $\pm$ 0.095	7.53 $\pm$ 0.046	5.84 $\pm$ 0.048

Table 5. Table 5: Sulphur and argon ionic abundances

HII Galaxy	$12 + l o g (S^{+} / H^{+})$	$12 + l o g (S^{2 +} / H^{+})$	$I C F (S^{3 +})$	$12 + l o g (A r^{2 +} / H^{+})$	$12 + l o g (A r^{3 +} / H^{+})$
FTDTR-1	5.15 $\pm$ 0.046	5.83 $\pm$ 0.046	1.93 $\pm$ 0.10	5.13 $\pm$ 0.052	5.14 $\pm$ 0.074
IZw18	4.70 $\pm$ 0.048	5.15 $\pm$ 0.047	1.19 $\pm$ 0.041	4.59 $\pm$ 0.057	3.85 $\pm$ 0.12
MRK36-A1	5.32 $\pm$ 0.036	5.98 $\pm$ 0.037	1.22 $\pm$ 0.045	5.38 $\pm$ 0.036	4.66 $\pm$ 0.098
MRK36-A2	5.43 $\pm$ 0.024	6.16 $\pm$ 0.026	1.17 $\pm$ 0.018	5.49 $\pm$ 0.024	4.64 $\pm$ 0.049
MRK475	5.38 $\pm$ 0.027	6.13 $\pm$ 0.030	1.15 $\pm$ 0.021	5.51 $\pm$ 0.028	4.57 $\pm$ 0.068
FTDTR-2	5.26 $\pm$ 0.047	5.98 $\pm$ 0.043	1.53 $\pm$ 0.098	5.31 $\pm$ 0.056	5.02 $\pm$ 0.11
IZw70	5.87 $\pm$ 0.048	6.36 $\pm$ 0.046	-	5.73 $\pm$ 0.048	-
MRK689	5.92 $\pm$ 0.045	6.55 $\pm$ 0.044	-	5.78 $\pm$ 0.045	-
MRK67	5.51 $\pm$ 0.037	6.17 $\pm$ 0.039	1.28 $\pm$ 0.045	5.60 $\pm$ 0.038	5.00 $\pm$ 0.095
FTDTR-3	5.21 $\pm$ 0.034	6.14 $\pm$ 0.035	1.93 $\pm$ 0.086	5.51 $\pm$ 0.034	5.48 $\pm$ 0.033
SHOC022	5.74 $\pm$ 0.054	6.40 $\pm$ 0.052	-	5.69 $\pm$ 0.055	-
FTDTR-4	5.61 $\pm$ 0.048	6.22 $\pm$ 0.046	1.15 $\pm$ 0.024	5.58 $\pm$ 0.048	4.67 $\pm$ 0.066
FTDTR-5	5.25 $\pm$ 0.11	6.01 $\pm$ 0.098	1.88 $\pm$ 0.35	5.37 $\pm$ 0.12	5.30 $\pm$ 0.17
FTDTR-6	5.41 $\pm$ 0.040	6.27 $\pm$ 0.040	1.23 $\pm$ 0.039	5.60 $\pm$ 0.040	4.88 $\pm$ 0.076
FTDTR-7	5.27 $\pm$ 0.028	6.08 $\pm$ 0.030	1.70 $\pm$ 0.050	5.44 $\pm$ 0.028	5.28 $\pm$ 0.023
MRK627	5.96 $\pm$ 0.051	6.53 $\pm$ 0.049	-	5.83 $\pm$ 0.053	-
PHL293B	4.96 $\pm$ 0.036	5.73 $\pm$ 0.039	1.60 $\pm$ 0.10	5.14 $\pm$ 0.037	4.91 $\pm$ 0.082
FTDTR-8	5.27 $\pm$ 0.11	6.00 $\pm$ 0.10	1.90 $\pm$ 0.24	5.41 $\pm$ 0.11	5.38 $\pm$ 0.078
SHOC263	5.96 $\pm$ 0.036	6.56 $\pm$ 0.043	-	5.81 $\pm$ 0.036	-
FTDTR-9	5.70 $\pm$ 0.076	6.47 $\pm$ 0.070	-	5.66 $\pm$ 0.081	-
FTDTR-10	5.39 $\pm$ 0.038	6.07 $\pm$ 0.034	-	5.26 $\pm$ 0.043	-

Table 6. Table 6: Element abundances and helium mass fractions using either oxygen Y O / H subscript 𝑌 𝑂 𝐻 Y_{\nicefrac{{O}}{{H}}} or sulphur Y S / H subscript 𝑌 𝑆 𝐻 Y_{\nicefrac{{S}}{{H}}} .

HII Galaxy	$H e / H$	$Y_{(O / H)}$	$Y_{(S / H)}$	$12 + l o g (O / H)$	$12 + l o g (N / H)$	$12 + l o g (S / H)$
FTDTR-1	0.0987 $\pm$ 0.0018	0.283 $\pm$ 0.0037	0.283 $\pm$ 0.0037	7.84 $\pm$ 0.070	5.93 $\pm$ 0.041	6.20 $\pm$ 0.056
IZw18	0.0757 $\pm$ 0.0018	0.232 $\pm$ 0.0042	0.232 $\pm$ 0.0042	7.15 $\pm$ 0.077	5.35 $\pm$ 0.038	5.36 $\pm$ 0.049
MRK36-A1	0.0785 $\pm$ 0.0018	0.239 $\pm$ 0.0041	0.239 $\pm$ 0.0041	7.85 $\pm$ 0.050	6.15 $\pm$ 0.063	6.15 $\pm$ 0.035
MRK36-A2	0.0796 $\pm$ 0.0015	0.241 $\pm$ 0.0035	0.242 $\pm$ 0.0035	7.93 $\pm$ 0.029	6.01 $\pm$ 0.035	6.30 $\pm$ 0.024
MRK475	0.0838 $\pm$ 0.0018	0.251 $\pm$ 0.0040	0.251 $\pm$ 0.0040	7.91 $\pm$ 0.027	6.31 $\pm$ 0.045	6.26 $\pm$ 0.028
FTDTR-2	0.0760 $\pm$ 0.0022	-	0.233 $\pm$ 0.0052	-	-	6.24 $\pm$ 0.055
IZw70	0.0944 $\pm$ 0.0021	0.273 $\pm$ 0.0044	0.274 $\pm$ 0.0043	8.27 $\pm$ 0.090	6.65 $\pm$ 0.032	6.48 $\pm$ 0.046
MRK689	0.0795 $\pm$ 0.0018	0.240 $\pm$ 0.0042	0.241 $\pm$ 0.0042	8.51 $\pm$ 0.082	6.56 $\pm$ 0.029	6.64 $\pm$ 0.044
MRK67	0.0846 $\pm$ 0.0020	0.252 $\pm$ 0.0043	0.253 $\pm$ 0.0044	8.12 $\pm$ 0.059	6.33 $\pm$ 0.037	6.36 $\pm$ 0.043
FTDTR-3	0.0833 $\pm$ 0.0019	0.249 $\pm$ 0.0042	0.250 $\pm$ 0.0042	8.16 $\pm$ 0.031	6.41 $\pm$ 0.053	6.47 $\pm$ 0.027
SHOC022	0.0846 $\pm$ 0.0020	0.252 $\pm$ 0.0044	0.253 $\pm$ 0.0044	8.19 $\pm$ 0.069	6.57 $\pm$ 0.043	6.49 $\pm$ 0.052
FTDTR-4	0.0873 $\pm$ 0.0017	0.258 $\pm$ 0.0038	0.259 $\pm$ 0.0038	8.08 $\pm$ 0.053	6.36 $\pm$ 0.024	6.37 $\pm$ 0.042
FTDTR-5	0.0764 $\pm$ 0.0040	-	0.234 $\pm$ 0.0093	-	-	6.35 $\pm$ 0.091
FTDTR-6	0.0830 $\pm$ 0.0017	0.248 $\pm$ 0.0039	0.249 $\pm$ 0.0039	8.12 $\pm$ 0.040	6.30 $\pm$ 0.051	6.41 $\pm$ 0.037
FTDTR-7	0.0844 $\pm$ 0.0015	0.252 $\pm$ 0.0033	0.252 $\pm$ 0.0034	8.12 $\pm$ 0.021	6.42 $\pm$ 0.031	6.37 $\pm$ 0.023
MRK627	0.0869 $\pm$ 0.0022	0.257 $\pm$ 0.0048	0.258 $\pm$ 0.0048	8.39 $\pm$ 0.085	6.75 $\pm$ 0.045	6.63 $\pm$ 0.049
PHL293B	0.0687 $\pm$ 0.0018	0.215 $\pm$ 0.0044	0.215 $\pm$ 0.0044	7.70 $\pm$ 0.022	-	6.00 $\pm$ 0.038
FTDTR-8	0.0737 $\pm$ 0.0035	0.227 $\pm$ 0.0084	0.228 $\pm$ 0.0084	7.97 $\pm$ 0.082	6.00 $\pm$ 0.16	6.35 $\pm$ 0.071
SHOC263	0.0874 $\pm$ 0.0042	0.258 $\pm$ 0.0092	0.259 $\pm$ 0.0092	8.32 $\pm$ 0.068	6.88 $\pm$ 0.042	6.65 $\pm$ 0.040
FTDTR-9	0.0892 $\pm$ 0.0032	0.262 $\pm$ 0.0070	0.263 $\pm$ 0.0070	8.23 $\pm$ 0.095	6.63 $\pm$ 0.069	6.54 $\pm$ 0.071
FTDTR-10	0.0746 $\pm$ 0.0031	0.229 $\pm$ 0.0072	0.230 $\pm$ 0.0073	7.64 $\pm$ 0.053	6.50 $\pm$ 0.057	6.15 $\pm$ 0.034

Table 7. Table 7: Comparison between the synthetic test case published in Aver et al. ( 2015 ) and one fitted with our algorithm. The metal emission parameters have been excluded from this comparison since they are not fitted by those authors methodology. Similarly, the absorption on the hydrogen and helium lines are not included since they are not fitted by our algorithm.

Parameter	True value	Aver et al (2015)	Current work
$y^{+}$	0.085	$0.0858 \pm 0.0027$	$0.0851 \pm 0.0009$
$y^{2 +}$	0.00088	✗	$0.00088 \pm 0.00002$
$n_{e} (H e)$	500.0	$473 \pm 67$	$487 \pm 34$
$τ (H e I)$	1.0	$0.78 \pm 0.31$	$0.95 \pm 0.18$
$T_{e} (H e)$	16000.0	$17320 \pm 1090$	$16311 \pm 824$
$c (H β)$	0.10	$0.09 \pm 0.03$	$0.106 \pm 0.015$
$ξ (\times 10^{4})$	1.0	$13_{+ 13}^{+ 19}$	✗

Table 8. Table 8: Primordial helium abundance determinations from all the linear regression combinations and comparison with the literature.

Element regression	Value	Number of objects
$Y_{P, O}$	$0.241 \pm 0.004$	16
$Y_{P, N}$	$0.244 \pm 0.003$	17
$Y_{P, S}$	$0.240 \pm 0.004$	19
$Y_{P, O - N - S}$	$0.243 \pm 0.005$	16
$Y_{P, O}^{1}$	$0.246 \pm 0.005$	18
$Y_{P, N}^{1}$	$0.251 \pm 0.005$	18
$Y_{P, S}^{1}$	$0.244 \pm 0.006$	21
$Y_{P, O - N - S}^{1}$	$0.245 \pm 0.007$	17
$Y_{P, O}^{2}$	$0.2446 \pm 0.0029$	5
$Y_{P, O}^{3}$	$0.2449 \pm 0.0040$	15
$Y_{P, O}^{4}$	$0.2551 \pm 0.0022$	28
$Y_{P, P l a n c k B B N}^{5}$	$0.24467 \pm 0.0002$	-

Table 9. Table 9: SHOC022 fitting results.

Parameter	F2018 value	Mean	Standard deviation	Number of points	Median	$16^{t h}$ percentile	$84^{t h}$ percentile	Difference $%$
$T_{l o w}$	11390	11590	597.2	12000	11590	11000	12190	1.703
$n_{e}$	49.46	37.04	20.36	12000	36.91	17.05	57.30	-34.01
$T_{h i g h}$	12170	12270	609.4	12000	12280	11660	12880	0.9137
$A r^{2 +}$	5.736	5.686	0.05484	12000	5.684	5.632	5.741	-0.9107
$N^{+}$	6.238	6.210	0.05657	12000	6.208	6.153	6.266	-0.4898
$O^{+}$	7.966	7.827	0.1275	12000	7.822	7.699	7.953	-1.836
$O^{2 +}$	7.972	7.937	0.06491	12000	7.933	7.874	8.001	-0.4868
$S^{+}$	5.767	5.744	0.05402	12000	5.742	5.690	5.797	-0.4333
$S^{2 +}$	6.419	6.399	0.05150	12000	6.398	6.347	6.451	-0.3336
$c (H β)$	0.1205	0.1229	0.02568	12000	0.1227	0.09705	0.1488	1.816
$τ$	None	1.075	0.4445	12000	0.9978	0.6639	1.479	None
$y^{+}$	0.08388	0.08328	0.001974	12000	0.08327	0.08131	0.08521	-0.7379
$y^{2 +}$	0.001362	0.001340	0.0001416	12000	0.001339	0.001200	0.001479	-1.764
$\frac{O I}{H I}$	8.270	8.192	0.06910	12000	8.188	8.124	8.261	-0.9948
$\frac{N I}{H I}$	6.552	6.575	0.04346	12000	6.571	6.533	6.617	0.2797
$\frac{S I}{H I}$	6.612	6.486	0.05178	12000	6.485	6.434	6.538	-1.954
$\frac{H e I}{H I}$	None	0.08462	0.001979	12000	0.08460	0.08263	0.08656	None
$Y_{O}$	0.2532	0.2521	0.004405	12000	0.2520	0.2477	0.2564	-0.4633
$Y_{S}$	0.2542	0.2528	0.004418	12000	0.2528	0.2484	0.2572	-0.5568

Table 10. Table 10: Atomic data references for the emission lines considered along with the parametrised relations.

Ion	Atomic data		Emissivity parametrisation
$H$	Storey & Hummer (1995)		$a + b \cdot l o g (T_{e}) + c \cdot l o g^{2} (T_{e})$
$H e$	Porter et al. (2013)		$(a + b \cdot n_{e}) \cdot l o g (t_{e}) - l o g (c + d \cdot n_{e})$
$H e^{+}$	Storey & Hummer (1995)		$a + b \cdot l o g (t_{e})$
Ion	Collision Strengths	Transition probabilities	Emissivity parametrisation
$O^{+}$	Pradhan et al. (2006); Tayal (2007)	Zeippen (1982); Wiese et al. (1996)	$a + \frac{b}{T_{e}} \cdot l o g (T_{e}) + c \cdot l o g (T_{e}) + l o g (1 + d \cdot n_{e})$
$S^{+}$	Tayal & Zatsarinny (2010)	Podobedova et al. (2009)
$O^{+ 2}$	Aggarwal & Keenan (2000)	Storey & Zeippen (2000); Wiese et al. (1996)	$a + b \cdot l o g (T_{e}) + c \cdot l o g (T_{e})$
$N^{+}$	Tayal (2011)	Wiese et al. (1996); Galavís et al. (1997)
$S^{+ 2}$	Hudson et al. (2012)	Podobedova et al. (2009)
$S^{+ 3}$	Tayal (2000)	Dufton et al. (1982); Johnson et al. (1986)
$A r^{+ 2}$	Galavis et al. (1995)	Kaufman & Sugar (1986); Galavis et al. (1995)
$A r^{+ 3}$	Ramsbottom & Bell (1997)	Mendoza & Zeippen (1982)

Table 11. Table 11: Emissivity parametrisation coefficients for the relations given in Table 10 . The Fitting precision for the HI, HeII and metal emissivities is better than 1 % percent 1 1\% for all the temperature and density ranges. In the HeI emissivities the precision decreases with up to a 5 % percent 5 5\% discrepancy in the n e subscript 𝑛 𝑒 n_{e} and T e subscript 𝑇 𝑒 T_{e} range of interest.

$L i n e$	a	b	c	d
$4341$ $H γ$	-0.846	0.230	-0.0251	-
$4363$ $[O I I I]$	4.97	-2.71	0.530	-
$4471$ $H e I$	0.0153	0.000383	2.20	0.000041
$4686$ $H e I I$	1.09	-0.0630	-	-
$4470$ $[A r I V]$	5.68	-1.21	0.736	-
$4959$ $[O I I I]$	5.23	-1.26	0.570	-
$5007$ $[O I I I]$	5.71	-1.26	0.570	-
$5876$ $H e I$	0.000	0.000791	0.828	0.000052
$6312$ $[S I I I]$	5.25	-1.64	0.712	-
$6548$ $[N I I]$	5.12	-0.906	0.545	-
$6563$ $H α$	1.59	-0.492	0.0522	-
$6583$ $[N I I]$	5.59	-0.906	0.545	-
$6678$ $H e I$	0.000637	0.000833	2.93	0.000109
$6717$ $[S I I]$	6.29	-0.924	0.423	0.00002
$6731$ $[S I I]$	6.14	-0.910	0.399	0.0001
$7136$ $[A r I I I]$	5.90	-0.820	0.506	-
$7319 + 7330$ $[O I I]$	4.84	-2.48	0.459	0.00021
$9069$ $[S I I I]$	5.47	-0.675	0.584	-
$9531$ $[S I I I]$	5.86	-0.675	0.584	-

Table 12. Table 12: Fitting results for additional test cases with different priors design described in the text.

Parameter	True value	Test 5	Test 6	Test 7	Test 8	Test 9
$T_{l o w}$	15590	$15640 \pm 283$	$15580 \pm 280$	$15640 \pm 281$	$15630 \pm 286$	$15630 \pm 282$
$n_{e}$	500	$493 \pm 36$	$493 \pm 36$	$489 \pm 38$	$498 \pm 38$	$497 \pm 38$
$S^{+}$	5.48	$5.474 \pm 0.020$	$5.478 \pm 0.019$	$5.474 \pm 0.020$	$5.476 \pm 0.020$	$5.476 \pm 0.020$
$S^{2 +}$	6.36	$6.354 \pm 0.023$	$6.359 \pm 0.022$	$6.354 \pm 0.022$	$6.355 \pm 0.023$	$6.355 \pm 0.022$
$O^{+}$	7.80	$7.793 \pm 0.037$	$7.802 \pm 0.037$	$7.795 \pm 0.037$	$7.794 \pm 0.037$	$7.794 \pm 0.037$
$O^{2 +}$	8.05	$8.048 \pm 0.017$	$8.047 \pm 0.017$	$8.049 \pm 0.017$	$8.049 \pm 0.017$	$8.049 \pm 0.017$
$A r^{2 +}$	5.72	$5.715 \pm 0.020$	$5.719 \pm 0.019$	$5.715 \pm 0.019$	$5.716 \pm 0.019$	$5.716 \pm 0.020$
$A r^{3 +}$	5.06	$5.058 \pm 0.016$	$5.057 \pm 0.016$	$5.059 \pm 0.016$	$5.059 \pm 0.015$	$5.059 \pm 0.015$
$N^{+}$	5.84	$5.835 \pm 0.020$	$5.839 \pm 0.019$	$5.836 \pm 0.019$	$5.836 \pm 0.019$	$5.836 \pm 0.019$
$c (H β)$	0.100	$0.106 \pm 0.018$	$0.104 \pm 0.018$	$0.106 \pm 0.018$	$0.105 \pm 0.018$	$0.105 \pm 0.018$
$T_{h i g h}$	16000	$16020 \pm 195$	$16040 \pm 198$	$16010 \pm 196$	$16020 \pm 191$	$16020 \pm 193$
$y^{+}$	0.0850	$0.085 \pm 0.001$	$0.085 \pm 0.001$	$0.085 \pm 0.001$	$0.085 \pm 0.001$	$0.085 \pm 0.001$
$y^{2 +}$	0.00088	$0.00088 \pm 0.00001$	$0.00088 \pm 0.00002$	$0.00088 \pm 0.00002$	$0.00088 \pm 0.00002$	$00088 \pm 00002$
$τ$	1.0	$0.987 \pm 0.228$	$0.990 \pm 0.228$	$1.030 \pm 0.281$	$0.972 \pm 0.230$	$0.976 \pm 0.228$

Equations10

\frac{F _{X^{i +}, λ}}{F _{H β}} = X^{i +} \frac{ϵ _{X^{i +}, λ} ( T _{e} , n _{e} )}{ϵ _{H β} ( T _{e} , n _{e} )} \cdot 1 0^{- c (H β) \cdot f_{λ}} \cdot k_{H e I, λ} (τ, T_{e}, n_{e})

\frac{F _{X^{i +}, λ}}{F _{H β}} = X^{i +} \frac{ϵ _{X^{i +}, λ} ( T _{e} , n _{e} )}{ϵ _{H β} ( T _{e} , n _{e} )} \cdot 1 0^{- c (H β) \cdot f_{λ}} \cdot k_{H e I, λ} (τ, T_{e}, n_{e})

P r (θ ∣ y) = \frac{P r ( y ∣ θ ) P r ( θ )}{P r ( y )}

P r (θ ∣ y) = \frac{P r ( y ∣ θ ) P r ( θ )}{P r ( y )}

t_{e} [O I I I] = 1.0807 t_{e} [S I I I] - 0.0846

t_{e} [O I I I] = 1.0807 t_{e} [S I I I] - 0.0846

Y_{P, S} = \frac{4 \frac{H e}{H} ( 1 - 20 \cdot \frac{O}{S} \frac{S}{H} )}{1 + 4 \frac{H e}{H}}

Y_{P, S} = \frac{4 \frac{H e}{H} ( 1 - 20 \cdot \frac{O}{S} \frac{S}{H} )}{1 + 4 \frac{H e}{H}}

Y_{P, O - N - S} = 0.243 \pm 0.005

Y_{P, O - N - S} = 0.243 \pm 0.005

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

\savesymbol

iint \savesymboliiint \savesymboliiiint

A Bayesian direct method implementation to fit emission line spectra: Application to the primordial He abundance determination

Vital Fernández1, Elena Terlevich1, Angeles I. Díaz2,3,4, Roberto Terlevich1,5

1Instituto Nacional de Astrofísica, Óptica y Electrónica, Luis E. Erro 1, 72840 Tonantzintla, Puebla, Mexico

2Departamento de Física Teórica, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

3Centro de Investigación Avanzada en Física Fundamental CIAFF-UAM

4Astro-UAM, UAM, Unidad Asociada CSIC

5 Institute of Astronomy, University of Cambridge, Madingley Rd., Cambridge CB3 0HA , UK E-mail: [email protected] (Vital Fernández)

(Accepted XXX. Received YYY; in original form ZZZ)

Abstract

This work presents a Bayesian algorithm to fit the recombination and collisionally excited line spectra of gas photoionized by clusters of young stars. The current model consists in fourteen dimensions: two electron temperatures, one electron density, the extinction coefficient, the optical depth on the $HeI$ recombination lines and nine ionic species. The results are in very good agreement with those previously published using the traditional methodology. The probabilistic programming library PyMC3 was chosen to explore the parameter space via a NUTs sampler. These machine learning tools provided excellent convergence quality and speed. The primordial helium abundance measured from a multivariable regression using oxygen, nitrogen and sulfur was $Y_{P,\,O-N-S}=0.243\pm 0.005$ in agreement with a standard Big Bang scenario.

keywords:

cosmology:primordial helium abundance – ISM:HII regions – sulfur abundance

††pubyear: 2015††pagerange: A Bayesian direct method implementation to fit emission line spectra: Application to the primordial He abundance determination–A

1 Introduction

The Standard Big Bang Nucleosynthesis (SBBN) is the model of choice for explaining the non-zero helium mass fraction observed in the universe for objects with almost zero metals. Inside its framework for a certain set of atomic parameters, the primordial nucleosynthesis by-products $(D,\,^{4}He,\,{}^{3}He$ and ${}^{7}Li)$ can be modelled as a function of the baryons-to-photons ratio, $\eta$ . These are commonly referred to as light elements and an empirical confirmation of their primordial abundances can provide very strong constraints on the early universe. An indirect measurement for these primordial abundances can be accomplished from the analysis of the Cosmic Microwave Background (CMB) since $\eta$ can be measured from its spectrum. Recently, the Planck Collaboration et al. (2018a) combined their CMB measurements with the algorithms by Pisanti et al. (2008) to predict the light element abundances under a $\Lambda CDM$ paradigm. Their prediction for the primordial helium mass fraction is $Y_{P}=0.24672_{-\left(0.00012\right)0.00061}^{+\left(0.00011\right)0.00061}$ . However, once the fitting is repeated allowing $Y_{P}$ and $N_{eff}$ (the number of neutrino species) to vary, the prediction uncertainty increases to $Y_{P}^{BBN}=0.247_{-0.036}^{+0.034}$ with $N_{eff}=2.89_{-0.57}^{+0.63}$ using a neutron lifetime without uncertainties of $\tau_{n}=880.2s$ .

The measurement of $Y_{P}$ involves the analysis of chemically unevolved gas reservoirs. In the recent work by Cooke & Fumagalli (2018), the helium abundance was measured from the absorption features of intergalactic gas observed against the light of a background quasar. Their $Y_{P}=0.250_{-0.025}^{+0.033}$ value is in agreement with the SBBN Planck Collaboration et al. (2018b) prediction. In the past, this technique has provided very accurate measurements for the primordial deuterium mass fraction (see Carswell et al., 1994) which may extend to the $Y_{P}$ estimation. Currently, however, state of the art measurements for $Y_{P}$ are obtained from the chemical analysis of the photo-ionised gas in extragalactic HII regions.

Over 40 years ago, Peimbert & Torres-Peimbert (1974) and Lequeux et al. (1979) published what is now known as the traditional method to determine $Y_{P}$ . These authors argued that in low metallicity regions the chemical enrichment $\left(\nicefrac{{dY}}{{dZ}}\right)$ is linear. They proposed to use oxygen as the tracer for the total metal content $\nicefrac{{dY}}{{dZ}}\approx\nicefrac{{dY}}{{dO}}$ . This is a very convenient choice for three reasons: firstly oxygen is the most common element after hydrogen and helium, accounting for almost half the metals mass fraction, $Z$ . Secondly, in the optical range oxygen has very intense lines, which provide an accurate measurement for its main ionic populations $\left(O^{+},\,O^{2+}\right)$ . Finally, oxygen is a heavy element produced by very massive stars in very short time scales. This guarantees that the gas reservoir is not largely contaminated by the current star-forming burst. These assumptions make it possible to apply a linear regression on the $Y$ vs $\nicefrac{{O}}{{H}}$ abundances relation for a given data sample. The coordinate $\nicefrac{{O}}{{H}}=0$ corresponds to $Y=Y_{P}$ : the helium abundance produced by the primordial nucleosynthesis.

Recent primordial helium abundance determinations using this technique have been published by Izotov et al. (2014) with $Y_{P}=0.2551\pm 0.0022$ , Aver et al. (2015) with $Y_{P}=0.2449\pm 0.004$ , Peimbert et al. (2016) with $Y_{P}=0.2446\pm 0.0029$ and Valerdi et al. (2019) with $Y_{P}=0.2451\pm 0.0026$ . The last works from Peimbert and collaborators include a computation of the neutron lifetime using the measured value $Y_{P}$ and a SBBN model. The estimated values: $\tau_{n}=870\pm 14$ and $\tau_{n}=873\pm 14$ are in good agreement with the value measured by Pattie et al. (2018) $\tau_{n}=877.7\pm 0.7_{stat}+0.4/-0.2_{sys}$ . Additionally, Fernández et al. (2018) (hereafter F2018) published the first $Y_{P}$ determination using sulfur as the metallicity tracer resulting in $Y_{P,S}=0.244\pm 0.006$ . A multivariable linear regression using oxygen, nitrogen and sulfur was proposed giving $Y_{P,\,O-N-S}=0.245\pm 0.007$ . It can be appreciated that the latter three groups measurements agree with a SBBN scenario, while the value from Izotov et al. (2014) is $2.2\sigma$ above the standard cosmology predictions. This discrepancy may be explained by the complexity in the chemical analysis. Peimbert & Peimbert (2011) summarised thirteen uncertainty sources, most of them systematic, which affect the helium abundance calculation. Many of these error sources also impact the metals abundance determination. These include the uncertainty on the atomic data, the temperature and density distributions or the accuracy on the ionisation correction factors.

The methodology described in F2018 follows a standard approach to derive the chemical composition of ionised gas. This procedure, recently reviewed by Peimbert et al. (2017) and Pérez-Montero (2017) is designed to isolate each phenomenon contributing to the observed emission. However, as both the astronomical data and models complexity keep increasing, this methodology becomes more cumbersome to apply. This is the reason why many researchers are embracing methodologies capable of exploring a multidimensional space. Examples of chemical analysis in the optical spectral range include the pioneer work by Olive & Skillman (2001) and Olive & Skillman (2004) which fitted simultaneously five parameters contributing to the helium lines emission. The aim of this new methodology is to provide a better quantification of the errors for the helium abundance determination. Later, this methodology was further enhanced by Aver et al. (2011); Aver et al. (2013); Aver et al. (2015), who included the hydrogen emission lines and increased the number of dimensions to eight. They also discussed the impact of the computational technique in order to explore the parameter space. Their chosen sampler was COSMOMC (see Li et al., 2014), a Fortran based program which applies a Markov Chain Monte Carlo (MCMC) process via a Metropolis-Hastings algorithm.

In this paper, we propose a methodology to fit the recombination and collisionally excited emission spectra. As it was done in F2018, the electron temperature and density are determined using the direct method (see Osterbrock & Ferland, 2005) which are anchored by appropriate collisionally excited line ratios. The main advantage of the new method is that all the model parameters contributing to the observed emission are solved simultaneously. This represents a space with up to fourteen dimensions. These parameters include nine ionised species of argon, helium, oxygen, nitrogen and sulfur in the optical and near infrared regions. To solve this system, we apply a HMC (Hamiltonian Monte Carlo) sampler based in machine learning algorithms. These tools provide a remarkable convergence quality and speed, which would have been impossible with standard MCMC algorithms. The data used is the one described in F2018.

The paper is organised as follows: Section 2 briefly reviews the F2018 sample and the treatment to account for the stellar and nebular continua. Section 3 describes how the direct method was adapted to a Bayesian paradigm and a technical description of the HMC sampler. Section 4 overviews a set of test cases to evaluate the accuracy and physical viability of the new technique. Section 5 presents the results and discusses merits and limitations of this technique. Finally, Section 6 summarises the main results.

2 The data

The data used in this paper are the 27 HII galaxies presented in Table 1 of F2018. In that paper we discussed the selection of the objects, the data reduction process and the derivation of their chemical abundances. The trends $Y-O/H$ , $Y-N/H$ and $Y-S/H$ and the obtained values for $Y_{P}$ are also discussed in that paper.

The measurement of the line fluxes has been described in detail in F2018. The process takes into account corrections for extinction and for underlying stellar and nebular continua. For the starforming ionised regions, the nebular continuum can reach comparable values to the underlying stellar continuum (see Reines et al., 2009).

The impact of this treatment of the spectra is illustrated in Fig. 1, where we plot the percentage increase in the line fluxes after the corrections have been applied for the most relevant recombination lines in our spectra. Even though the contribution for the strong hydrogen Balmer lines may be considered negligible, we observe about a 20 % increase in the HeI lines intensity once the continuum correction is taken into account and this percentage changes with wavelength.

3 Bayesian direct method.

The flux from an emission line can be expressed as a function of the transition emissivity relative to a Balmer line, e.g. $H\beta$ :

[TABLE]

where $\nicefrac{{\epsilon_{X^{i+},\,\lambda}}}{{\epsilon_{H\beta}}}$ is the relative emissivity at the transition wavelength $\lambda$ , for an ion with abundance $X^{i+}$ , at certain electron temperature $T_{e}\,\left(K\right)$ and electron density $n_{e}$ $\left(cm^{-3}\right)$ . The term $c\left(H\beta\right)$ is the logarithmic extinction at $H\beta$ for a reddening law $f_{\lambda}$ . The parameter $k_{HeI,\,\lambda}$ represents a fluorescence excitation correction for the $HeI$ lines calculated by Benjamin et al. (2002) as a function of the optical depth $\tau$ .

Traditionally, as already mentioned, the parameters in eq. 1 are solved individually by isolating the different processes contributing to the observed flux. Usually, the first step consists in computing $c\left(H\beta\right)$ from the difference between the theoretical and observed fluxes of the hydrogen Balmer series. This is a valid approach since most of the hydrogen is ionised. Additionally, their emissivity remains almost constant for the expected $T_{e}$ and $n_{e}$ values.

The next step in the traditional method consists in computing the ionic abundances from the observed emission lines. In the so called direct method, we calculate the emissivity coefficient by providing values for $T_{e}$ and $n_{e}$ obtained by measuring temperature-sensitive and density-sensitive line ratios. Finally, using these parameters back in eq. 1 the ionic abundance can be computed from each emission line.

As we discussed in F2018, however, there are problems inherent to this method. Firstly, it can underestimate the uncertainty in the model parameters since these are fitted individually. To properly propagate the error to the abundances, several iterations may be necessary, as well as complex networks of chained Monte Carlo algorithms. Secondly, eq. 1 represents the most basic interpretation. As we introduce more complexity in the physics of the problem, like e.g. $k_{\tau}$ in eq. 1, it becomes harder to solve the model parameters. Finally, with the large amount of data available today, plus its larger wavelength coverage, this methodology becomes harder to apply.

In this paper, we use a method that allows to fit all the parameters in eq. 1 simultaneously. These are two electron temperatures for the low and high ionisation regions, $T_{low}$ and $T_{high}$ , one electron density $n_{e}$ , which is assumed to remain constant within the star forming region, the extinction coefficient $c\left(H\beta\right)$ and nine ionic abundances $Ar^{2+}$ , $Ar^{3+}$ , $y^{+}$ , $y^{2+}$ , $O^{+}$ , $O^{2+}$ , $N^{+}$ , $S^{+}$ and $S^{2+}$ . Here we add the optical depth, $\tau$ , to account for the fluorescence excitation in the helium lines as tabulated by Benjamin et al. (2002). This brings the maximum number of dimensions to fourteen and in order to explore this relatively large parameter space, it becomes necessary to apply machine learning tools.

Currently, data science has become very popular due to the computational advances in three domains: big data, deep learning and probabilistic programming. In the first field, the user is interested in finding patterns in large data sources. In deep learning, the user defines the parameter space via neural networks which have several ‘depths’. Finally, probabilistic programming focuses in declaring models, whose parameters and outputs are probability functions. We applied the probabilistic programming package PyMC3 by Salvatier et al. (2016), which makes use of the deep learning library Theano by The Theano Development Team et al. (2016) to define the physical model. This package includes a NUTs (No-U-Turns) sampler (see Hoffman & Gelman, 2011) to explore the parameter space. This sampler follows a Hamiltonian Monte Carlo (HMC) paradigm, which drops the stochastic jumps from the Markov Chain by an informed sampling guided by the mathematical model derivatives. This algorithm implementation involves a more challenging programming experience than a MCMC sampler. Still, this HMC algorithm provides excellent convergence quality for large parameter spaces. Moreover, this sampler decreases the fitting computational time from several hours $\left(>6\,h\right)$ to a couple of minutes.

It may be inferred from the previous definition that any probabilistic programming implementation is intrinsically Bayesian. This paradigm is characterised by the application of the Bayesian theorem:

[TABLE]

where $Pr\left(\theta\right)$ is referred to as the prior. This term represents the probability distribution of the model parameters $\theta$ before the fitting. $Pr\left(y|\theta\right)$ is the likelihood of the model and it provides an evaluation on how the observational data $y$ adjusts to the theoretical model. $Pr\left(y\right)$ is referred to as the evidence. This parameter is actually the integral of the numerator in eq. 2 over the complete parameter space and it represents the probability that the observed data has been generated by the processes described by the model. Finally, $P\left(\theta|y\right)$ is the probability of the model parameters given the observations. In a Bayesian inference, the posterior provides the user with a credible region: given the observed data, there is a $95\%$ probability that the true value of $\theta$ falls within this credible region $CR_{\theta}$ . To successfully compute this diagnostic, however, it is essential to properly define the terms in the Bayesian theorem.

Recently, Tak et al. (2018) reviewed the quality of Bayesian models in astronomical literature. These authors emphasised how improper priors can result in output distributions not meeting the posterior propriety. An example of improper prior includes a uniform distribution covering the real space $\left(-\infty,\,\infty\right)$ . Moreover, in the particular case of a NUTs sampler a uniform prior can dramatically affect both the convergence quality and the simulation speed. To deal with this issue, in the Bayesian approach, one can define a joint distribution using proper priors to ensure the posterior propriety. On the one hand, we can rely on the scientific evidence to implement very informative (or constrained) prior distributions. On the other hand, in cases where little physical knowledge is available, it is acceptable to provide uninformative (or wide) prior distributions. In practice, this provides a uniform probability distribution for a region of interest. These two strategies are considered here and the input priors are displayed in Table 1. The following paragraphs describe how the traditional direct method applied in F2018 was adjusted to a Bayesian paradigm using proper priors:

•

We consider in this study two ionisation regions. These are characterised by two electron temperatures: $T_{low}$ and $T_{high}$ . In the low ionisation region the ionised species are $Ar^{2+}$ , $H^{+},$ $O^{+}$ , $N^{+}$ , $S^{+}$ and $S^{2+}$ . The high ionisation species are $Ar^{3+}$ , $He^{+}$ , $He^{2+}$ and $O^{2+}$ . The two temperatures share the same prior design: a Gaussian distribution with $\mu=15000\,K$ and $\sigma=5000\,K$ . This range provides a good coverage for the temperatures commonly encountered in HII galaxies (as it was shown in F2018). The temperature depends on the auroral lines observed: $\left[SIII\right]6312\text{\AA}$ for $T_{low}$ and $\left[OIII\right]4363\text{\AA}$ for $T_{high}$ . In the objects where only one auroral line was available, only the corresponding temperature priors is declared. The other temperature is calculated using the linear relation provided by Pérez-Montero (2017):

[TABLE]

where $t_{e}$ is in $10^{4}K$ units.

•

The electron density calculation from the $\left[SII\right]\lambda\lambda 6717$ , $6731\text{\AA}$ line ratio falls within the very low density regime for most objects $\left(<\,100\,cm^{-3}\right)$ . This is actually the regime at which this line ratio becomes less sensitive to density. This means that the $\left[SII\right]$ density cannot be computed to high accuracy. This is not, however, a big issue from the chemical analysis point of view as for most ionised species the emissivity also remains independent with the density at these regimes. This can be an issue from the mathematical point of view, though. If the density cannot be fitted from the available data, the resulting posterior distribution will actually be the prior one. In order to force the simulation to stay in a representative density region, we define the density prior via a Gaussian distribution with $\mu=n_{\left[SII\right]}$ and $\sigma=\sigma_{n_{e}\left[SII\right]}$ measured in F2018.

•

The reddening curve chosen for this work is the ’LMC average’ published by Gordon et al. (2003) with a $R_{V}=3.4$ appropriate for a star forming region. The gas extinction in these objects is usually very low (e.g. Terlevich et al., 1991). This parameter, however, cannot be negative. To account for this physical limit a log-normal distribution was chosen with $\sigma=1.0$ . This guarantees higher probability to the expected extinction coefficient range.

•

In F2018, all calculations involving the line emissivities were accomplished using PyNeb (see Luridiana et al., 2015). Third party libraries, however, are not easily imported into deep learning algorithms. Implementing a 2-dimensional interpolation on emissivity grids is beyond the scope of this project. The alternative chosen is a strategy commonly used in the traditional chemical analysis: Parametrised equations for the emissivity as a function of the electron temperature and density. For each emission line an emissivity grid was computed using PyNeb for the references in Table 10 in the appendix and which also includes the emissivity parametrisation for each ion. Table 11 displays the fitted coefficients for the $8,000$ - $25,000\,K$ and $1$ - $600\,cm^{-3}$ surface for each emission line. These fittings provide the line emissivity in log scale. This was done for two reasons: first, the emissivity grids are easier to fit in a log scale surface. Second, to improve the metal abundance sampling it was desirable to use the standard $12+log\left(X^{+}\right)$ notation. For the ionic abundances uninformative priors can be applied. For the metal species Gaussian distributions are considered with $\mu=5$ and $\sigma=5$ . The $y^{+}$ and $y^{2+}$ abundances are given in the linear scale and a log-normal distribution is used to model their priors with $\mu=0$ and $\sigma=1$ . The values drawn from these priors are scaled by the expected $y^{+}$ and $y^{2+}$ regime via coefficients $k_{y^{+}}=\nicefrac{{1}}{{10}}$ and $k_{y^{2+}}=\nicefrac{{1}}{{1000}}$ .

•

To account for fluorescence excitation on the helium lines, a new correction has been included on the analysis. The radiative transfer calculations from Benjamin et al. (2002) are included in eq. 1 as a function of $T_{high}$ and $n_{e}$ . A successful fitting for this parameter depends on the availability of $HeI$ lines which are sensitive to this effect. As it will be showed in the next section test cases, this parameter fits properly in synthetic observations for both prior probability distributions considered. However, our current observational data lacks emission lines which are strongly affected by this phenomenon. Therefore, it cannot be quantified. Consequently, we defined a prior with a log-normal distribution with $\mu=0$ and $\sigma=0.4$ . This choice is justified by the following arguments: 1) In Aver (2012), it was shown that more than half of the fifty successfully fitted spectra using the same optical depth model have an opacity below one and only three objects displayed opacities above four. Therefore, this distribution guarantees a proper prior, even if very informative, which is justified by known data. This is an acceptable practice for statistical inference in scientific models. 2) Using synthetic test cases, it was confirmed that this probability distribution does not result in degeneracies with other parameters, in particular $y^{+}$ . This was not the case of prior designs, which covered large $\tau$ values.

The last element to define in eq. 2 is the model likelihood. At this point, it is important to remember that in the Bayesian paradigm there is not uncertainty in the model data $y$ : any randomness in the observables is caused by the model parameters, which behave as a probability distribution. This interpretation has a physical foundation in models such as this one. For example: $T_{e}$ is generally interpreted as the result of a Maxwellian velocity distribution of the electrons. In practice, it is not common for spectra to show the same uncertainty everywhere along the wavelength range. For example, in the present case it is essential to assert the difference in uncertainty between the $\left[OIII\right]$ auroral and nebular lines. A valid approach to account for this uncertainty consists in establishing a normal distribution, whose standard deviation is weighted by the observational error, as the likelihood. This was the approach considered in this model (shown in Table 1) with a normal distribution likelihood with $\mu=F_{\lambda}$ and $\sigma=\sigma_{F_{\lambda}}$ for each emission line where $F_{\lambda}$ represents the flux relative to $H\beta$ including the error propagation in the ratio.

4 Model test cases

Possibly the most discouraging feature of Monte Carlo simulations is assessing the convergence quality. Fortunately, in empirical models such as this one, the convergence can be evaluated graphically. This process has two phases:

In the first one, a set of values is chosen for the model parameters from which the synthetic observables are generated. Afterwards, this synthetic spectrum is fed into the algorithm: The output fittings are then compared with the true values. In addition to a good accuracy it is expected that the credible region width is of the same order as the uncertainty found in the traditional methodology. Table 2 displays in column 1 a set of values for the model parameters, which were used to generate an array of synthetic emission lines using eq. 1. In the test cases below, a uniform $2\%$ uncertainty was adopted for the emission lines flux considered. The following paragraphs describe these test cases.

In the simplest example only four model parameters are fitted $\left(T_{low},\,n_{e},\,S^{+},S^{2+}\right)$ using five emission lines: $\left[SII\right]\lambda\lambda 6717$ , $6731\text{\AA}$ , $\left[SIII\right]\lambda 6312\text{\AA}$ and $\left[SIII\right]\lambda\lambda 9069,9531\text{\AA}$ . In a Monte Carlo simulation, such as this one, the user declares the number of iterations and the number of tuning steps, which are the number of initial steps not included in the final statistics. This is because they account for the time the algorithm requires to reach the optimum jump step size to the solution region. An additional good practice consists in repeating the simulation to confirm that the initial conditions do not alter the results. Fortunately, packages such as PyMC3 can parallelize the models easily by running one simulation per machine processor. Once all the fittings finish, the parameters measurement includes the traces from all the machine cores. This enables shorter simulations and a better resource management. All this information can be seen in Fig.2. The four plots on the left side represent the model parameter traces evolution with time. It can be appreciated that during the simulation steps, all the traces have a white noise-like appearance. This signals a simulation which has successfully converged to a parameter space coordinate (not necessarily the right one). This simulation configuration consists in 6000 simulation steps and 2000 tuning steps parallelised in a 1.8-GHz i7-5400U dual-core processor laptop. These simulation pairs are joined in Fig.2 where we can see that the traces are indistinguishable. These traces provide us with the statistics for the parameters measurement. The value written over the traces are the mean and standard deviation from the traces. These are the values in Table 2. The $16^{th}$ and $84^{th}$ percentiles from the posterior distribution are written on the ordinate axis. The median value is given in the histograms. It can be concluded from this initial test that the simulation fittings are in good agreement with the true values.

Fig. 3 displays the results for the test cases 2 and 3 similarly to Fig.2, but now all the metal emission lines are included as input. Consequently, all metals ionic abundances are included as well. In this simulation, we still only have one temperature prior $(T_{low})$ . In order to compute the physically correct temperature for the high ionization ions, the algorithm automatically includes the empirical linear model in eq. 3. Hence, at each iteration the temperature value drawn from this distribution is applied directly to compute the emissivity for the low ionization species emissivity. For the high ionization ones, however, this temperature is previously corrected according to eq. 3. It may be concluded from the results that all the abundances were properly fitted even for elements, such as the argon ions, where only one emission line was observed for each ion. The extinction coefficient $c(H\beta)$ is also considered in the analysis. Usually, the extinction cannot be determined using collisionally excited lines due to their sensitivity to temperature and density. In this case, however, we can see from Fig. 3 that $c\left(H\beta\right)$ was successfully measured using only collisionally excited lines. There are two reasons for that: first, two auroral lines, $\left[OIII\right]\lambda 4363\text{\AA}$ and $\left[SIII\right]\lambda 6312\text{\AA}$ , are anchoring the electron temperature. Second, this fitting includes emission lines from $\left[OIII\right]\lambda 4363\text{\AA}$ to $\left[SIII\right]\lambda 9531\text{\AA}$ . This wide wavelength range guarantees data points with varying sensitivity to the extinction, which improves the sampling of $c\left(H\beta\right)$ .

In the $3^{rd}$ test case, the simulation includes a prior $T_{high}$ for the high ionization region electron temperature. It should be emphasised that in the previous test cases the ionised species had the same assignment the low or high ionisation regions. The difference now is that $T_{high}$ is fitted along the rest of the model parameters instead of being calculated using eq. 3. It can be appreciated in the right hand side of Fig. 3 that in this simulation there is a greater uncertainty in all the parameter fittings. This is because the previous conditions are no longer met. To begin with only one auroral line is anchoring each electron temperature. Moreover, in this simulation only the $\left[OIII\right]$ and $\left[ArIV\right]$ emissions are located in the high ionisation region. These lines cover a narrow wavelength range to sample the parameter space which makes the $c\left(H\beta\right)$ sampling harder. Finally, some small divergence on the traces for some parameters during the simulation can be appreciated. We conclude that this test case does not have enough data to fit the model parameters.

In the $4^{th}$ test case, the simulation runs the complete model: three more dimensions $\left(y^{+},\,y^{2+},\tau\right)$ are included reaching a total of fourteen. As all the hydrogen and helium lines in the covered spectral region are included, the number of inputs has also increased. Moreover, the hydrogen recombination lines provide a firm anchor to the dust extinction since this parameter has the greatest impact in their intensity. Fig. 4 shows that in this case the simulation has properly converged. Additionally, it can be appreciated that all the parameter fittings are very close to the true values tabulated in Table 2.

An important practice in Bayesian models is to check the impact of the priors design in the sampling process and results. The results from additional test cases for the complete model can be found in Table 12. In the $5^{th}$ test case the $T_{high}$ temperature prior has a uniform distribution which covers the complete temperature range: from $8000$ to $22000K$ . In the $6^{th}$ test case both temperature priors have a Gaussian distribution but $T_{low}$ is centred at $10000K$ while $T_{high}$ is centred at $20000K$ . Both distributions have a standard deviation of $2500K$ which means that the true values at $15590K$ and $16000K$ have a relatively low probability. In the $7^{th}$ and $8^{th}$ test cases we are assigning uniform priors for the optical depth and the electron density respectively. In the former case, the uniform priors lower and upper limits are 0 and 10 while in the latter the limits are 1 and 1000 $cm^{-3}$ . Finally, in the $9^{th}$ test case the ionic helium abundance priors are defined by the relation: $y^{i+}=k_{y^{i+}}+Normal\left(\mu=0.0,\,\sigma=1.0\right)\cdot k_{y^{i+}}$ , where $k_{y^{i+}}=0.1$ and $k_{y^{i+}}=0.001$ for the $y^{+}$ and $y^{2+}$ abundances respectively. These priors represent a basic reparameterization: Distribution centred at zero with a certain offset. This prior design is not easy to read from a researcher point of view. Indeed, it may seem the obvious choice to assign a probability distribution which represents the physical parameter magnitude and behaviour as close as possible. Nevertheless, from the computational point of view, it is encouraged to use the same probability distributions centred at zero for all the model parameters and parameterizing the theoretical model instead. This is because the sampling process becomes more efficient. We can confirm in Table 12 that the prior design is not affecting the results and the simulations are consistently converging at the same coordinate. These test cases represent an ideal scenario, where all the necessary inputs are available. In real observations where some emission lines may not be available the prior design impact needs to be checked again.

An additional graphical tool to establish a Monte Carlo simulation convergence quality is a scatter plot matrix. This is shown in Fig. 5, in which the parameter traces from the complete test case are plotted against each other as surface distributions. The red lines mark the true value location, which are within the algorithm sampled region. It can be appreciated that there is some degeneracy between the metals ionic abundance and electron temperature. This degeneracy, however, is intrinsic for the physical model since for the temperature and abundance range these parameters are correlated. This degeneracy might be decreased by parameterizing the mathematical parameter space, thus making it easier to sample. This would also improve the speed of the model, which for the current implementation is over 100 steps per second. However, the current algorithm convergence quality and speed is good and the fittings presented display accuracy and precision in the measurement. In the second phase, the algorithm is tested on real observations and its convergence is compared to the synthetic ones. This is done in the following section using the data values from F2018.

5 Results and discussion

The initial implementations of this Bayesian algorithm in real spectra resulted in a slow sampling process or even failure to launch. This was caused by some emission lines which had a very narrow likelihood in the inference model. For example, in this sample spectra the $\left[OIII\right]\lambda 5007\text{\AA}$ and $H\alpha$ lines can display uncertainties below $0.5\%$ of their integrated flux. Running synthetic test cases with similar likelihoods we learned that the algorithm found it hard to find a solution given the wide temperature and abundance priors provided. This uncertainty does not really represent the physical processes intrinsic randomness. Instead, it quantifies the emission line signal-to-noise ratio. Consequently, applying this uncertainty in the emission line likelihood implies very small sampling steps for the model parameters. As most of the proposed parameter values are rejected the simulation runs very slowly. To deal with this issue a minimum $2\%$ uncertainty was set on the emission lines likelihood. This is the value used in the test cases, which provided stable solutions and measurements with the expected uncertainty.

Fig.6 corresponds to the fitting of SHOC022. In most object fittings, the correlation between the ionic abundances and the electron temperatures seem to be more pronounced than the one displayed in Fig. 5 for synthetic inputs. Still, it can be established that the data is well represented by the theoretical model and the correlation between parameters is the one expected. An exception, however, is the optical depth $\tau$ . The surface plots in Fig.6 for this parameter imply a sampling limited by a mathematical boundary rather than a true fitting. The reason behind these results can be found in the available helium lines in the sample spectra. As it was explained in F2018, due to a mismatch between our data and the available Single Stellar Populations (SSP) synthesis models wavelength range: only three $HeI$ lines include a correction for the underlying stellar population ( $HeI\lambda 4471$ , $5876$ and $6678\text{\AA}$ ). These emission lines, however, are weakly affected by fluorescence. As shown in Benjamin et al. (2002), the correction $k_{\tau}$ in eq.1 for the $HeI6678\text{\AA}$ line accounts for less than $0.4\%$ even for the highest $\tau$ values. Therefore, it is very hard to fit $\tau$ in a parameter space, where the other dimensions have a much larger impact in the final line flux. Actually, the $\tau$ surface distributions in Fig.6 correspond to the input prior in Table 1: a log-normal distribution with $\mu=0$ and $\sigma=0.4$ . Therefore, the optical depth cannot be measured from the provided observation inputs. Nevertheless, the fluorescence excitation is still being taken into account for this parameter and the simulation behaves as a standard Monte Carlo algorithm for a user-declared parameter distribution. This is an improvement from the helium abundance analysis in F2018. The scatter plots for the complete HIIGs sample can be found in the online support material.

Once the simulation convergence quality has been confirmed, the next step consists in evaluating the fitting result. In a frequentist analysis this can be done via the $\chi^{2}$ test. This is the approach followed by Aver et al. (2015) and references therein, where galaxies with $\chi^{2}<5.5$ are excluded from the $Y_{P}$ regression. In the Bayesian paradigm, there is not an universal estimator which can quantify the fitting quality. However, in this particular analysis an efficient evaluation can be accomplished by comparing the output flux distribution with the input observed emission fluxes. This comparison can be found in the tables included in the supporting online material. This information can also be represented graphically and is shown in Fig.7: The fitted emission line flux distributions are plotted from bluest to reddest wavelength. The distributions are colour coded according to the ion responsible for each transition. Each plot cell includes the observed line flux as a vertical line and its uncertainty as a shaded area. These results correspond to the initial fitting of the object. In general, it can be appreciated that the flux distributions centre is very close to the observed value and their width is representative of the measurement uncertainty. There are, however, some flux discrepancies which require additional explanation:

•

$H_{\gamma}$ : The algorithm fails to fit this emission line. The reason can be found in the instrumental setup: For most objects this line lies right at the low wavelength edge. Consequently, the pixels belonging to this line have a greater noise, which is not being propagated by our current reduction pipeline. This was clearly observed during the $c\left(H\beta\right)$ calculation in F2018, where this line did not match the extinction calculated from the comparison between the $H\alpha$ to $H\beta$ flux ratio with the case B recombination theoretical value. Consequently, this line have been excluded in most objects fitting except for those with the higher redshift.

•

$HeI$ lines: not all helium lines are fitted with the same precision. For example in SHOC022 the fitted flux for $HeI\lambda 6678\text{\AA}$ is $6\%$ above the observed value while the $HeI\lambda 4471\text{\AA}$ and $HeI\lambda 5876\text{\AA}$ fluxes disagree only by a $0.5\%$ and a $2\%$ respectively. In the particular case of the helium lines this behaviour can be explained by the absorption from the underlying stellar population. As it was shown in Fig. 1, each line is affected differently by this effect. Additionally, in the SSP synthesis applied in F2018 the fitted continuum does not include the uncertainty. Consequently, our current methodology cannot quantify which helium line absorptions are better fitted (as indeed was also the case in F2018). The methodology presented in this paper, though, makes it much easier to check which are the lines contributing more to the helium abundance uncertainty.

•

$\left[OII\right]\lambda\lambda 7319,7330\text{\AA}$ lines: In order to include the recombination correction from Liu et al. (2001) our algorithm takes as an input the integrated flux from both lines even though our spectral resolution can separate the $\left[OII\right]$ doublet.

•

$\left[NII\right]$ lines: Some of our objects show a broad $H\alpha$ component, that made difficult the deblending of the narrow $H\alpha$ and $[NII]$ lines as explained in F2018. The observed ratio of the $[NII]$ doublet lines differ from the theoretical value for some of the objects and large uncertainty was measured. Since the $\left[NII\right]\lambda 6548\text{\AA}$ line is almost three times weaker than $\left[NII\right]\lambda 6583\text{\AA}$ and therefore harder to deblend, it was excluded in those fittings were a large mismatch was found in both lines fittings.

•

$\left[SIII\right]$ lines: Even though the $\left[SIII\right]\lambda 6312\text{\AA}$ line is fitted with good precision, there is a small mismatch for the infrared lines due to telluric contamination. As discussed in F2018, even with the extra observations to calibrate the sky features the ratio between the $\left[SIII\right]\lambda 9069\text{\AA}$ and $\left[SIII\right]\lambda 9531\text{\AA}$ lines is within $5\%$ of the expected value. Our theoretical model, however, does not include the sky contribution in the $\left[SIII\right]$ computation, hence the slight disagreement between the observed fluxes and the simulated ones. Our analysis follows the same approach as in F2018: both lines are included in the chemical analysis except in those cases where the telluric calibration observations were not available. For those objects, only the line less affected by the sky was used.

Once the discrepancies between the flux distributions generated by the simulation and the observed line fluxes have been dealt with, we may proceed with the next phase: evaluating the posterior distributions for the model parameters. This provides an opportunity to compare the results with those obtained classically. The results for SHOC022 are given in Table 9 in the appendix. For the rest of the sample, the results are given online. In these tables, column 1 corresponds to the parameter label, column 2 is the measurement published in F2018, columns 3 and 4 are the trace mean and standard deviation, columns 5, 6 and 7 are the $50^{th}$ (median), $16^{th}$ and $84^{th}$ percentiles and column 8 is the percentage difference between the F2018 measurement and the present one. The graphical simulations presented in section 4 are also available for each object. The one corresponding to the final fitting of SHOC022 can be found in Fig. 8 where the vertical line in the histograms corresponds to the parameter measurement in F2018 and the shaded area represents its uncertainty. It can be easily noticed that the $\tau$ posterior is actually the prior log-normal as previously discussed. In general, however, most parameters agree well, within the $5\%$ , with those from F2018. The following paragraphs discuss patterns encountered in most object fittings:

•

$T_{low}$ and $T_{high}$ : We have shown in F2018 that, when one of the temperature diagnostic lines $\left[OIII\right]\lambda 4363\text{\AA}$ or $\left[SIII\right]\lambda 6312\text{\AA}$ had larger uncertainty than the other, it was excluded from the analysis and the corresponding temperature was calculated using eq. 3. In order to duplicate the same procedure here, two actions were necessary. First to replace the corresponding prior by eq. 3. However, as discussed in section 4 second test case, once eq. 3 is added to the simulation both auroral lines, and their uncertainty, dominate the temperature fitting. Therefore, as done in F2018, the second step consisted in excluding that line from the fitting.

•

$n_{e}\left[SII\right]$ : In many objects there is a perfect match between the new $n_{e}\left[SII\right]$ value and the one published in F2018. This is expected since the density prior is well constrained from the density measured via the $\left[SII\right]\lambda\lambda 6717$ , $6731\text{\AA}$ ratio. In some cases, including SHOC022 in Fig. 8, we obtain lower densities (e.g. $n_{e}\left[SII\right]=37\,cm^{-3}$ vs $n_{e}\left[SII\right]=50\,cm^{-3}$ in F2018 37% smaller in this case). As it was described in F2018 the reported density values in cases were $n_{e}\left[SII\right]<75\,cm^{-3}$ the density applied in the chemical analysis was $n_{e}\left[SII\right]=50\pm 25\,cm^{-3}$ . Therefore these results are consistent with those previously published.

•

$c\left(H\beta\right)$ : We found that for some objects the extinction measured via the current algorithm was higher and with a greater uncertainty than those in F2018. In particular, MRK475 and PHL293,with very low extinction values in F2018 $\left(cH\beta\approx 0.01\right)$ the new extinction coefficient measurements are one order of magnitude larger. There are several reasons why this should be the case: to start with, in the current analysis there are thirteen additional parameters in the fitting. Secondly, while in the traditional methodology $H\beta$ and $H\alpha$ dominate the $c\left(H\beta\right)$ calculation, in the current algorithm all emission lines are contributing to its measurement. In this study, the number of collisionally excited lines is greater than the recombination ones. Consequently, they weight heavily in the extinction determination. Nevertheless, as long as the simulated line fluxes match the observed ones it can be concluded that the new value quantifies the extinction adequately.

•

$Ar^{+3}$ : This ion abundance measurement depends exclusively on the $\left[ArIV\right]\lambda 4741\text{\AA}$ observation. This is a weak line even when compared with the helium ones. In the extreme case of $IZw18$ , the uncertainty for this line is above $20\%$ of the observed flux. Due to these issues the sampler had problems to fit the $Ar^{3+}$ abundance in this galaxy. Consequently, the prior for this parameter was changed from the one in Table 1 to a normal distribution with $\mu=4.5$ and $\sigma=0.5$ . This prior provided a more informative sampling improving the fitting for this parameter.

5.1 Comparison between methods

The fitting results for the complete sample are tabulated as follows: the electron density $\left(n_{e}\left[SII\right]\right)$ , the low and high ionisation temperatures and the logarithmic extinction coefficient at $H\beta$ are shown in columns 2, 3, 4 and 5 in Table 3. The optical depth, defined as the mean value from the prior distribution, was found to be $\tau\approx 1.05$ for most objects. Table 4 includes the $y^{+}$ and $y^{2+}$ abundances in columns 2 and 3, the oxygen abundances, columns 4 and 5 and nitrogen $N^{+}$ abundance in column 6. Table 5 displays the sulfur ionic abundances in columns 2 and 3 and the argon ionic abundances in columns 5 and 6. The $ICF\left(S^{3+}\right)$ is shown in Column 4 and was obtained following F2018, calculated via the $S^{2+}$ , $Ar^{2+}$ , $Ar^{3+}$ abundances. Finally, Table 6 displays the element abundances: column 2 corresponds to the helium abundance while columns 3 and 4 correspond to the helium mass fractions computed using the oxygen and sulfur abundances respectively. columns 5, 6 and 7 are the oxygen, nitrogen and sulfur abundances as $12+log\left(X\right)$ . We calculate the total abundances as in F2018, by adding the ionic abundances for each element using the traces as in a standard Monte Carlo.

In general, these measurements agree very well with the ones published in F2018. The differences between both methods is below $3\%$ for most parameters except $n_{e}$ and $c\left(H\beta\right)$ as discussed above. A special case is the galaxy PHL293B. The comparison between the simulation fluxes and the observed ones is as good as for the rest of the sample. The parameter fittings, however, disagree considerably from those in F2018. The measured electron temperatures are $10\%$ and $5\%$ above the F2018 values for $T_{e}\left[SIII\right]$ and $T_{e}\left[OIII\right]$ respectively and the $c\left(H\beta\right)$ magnitude increased from $0.01$ to $0.23$ once all the lines are taken in consideration. Terlevich et al. (2014) presented WHT-ISIS and X-shooter spectra to study the broad $\left(FWHM=1000\,\nicefrac{{km}}{{s}}\right)$ and very broad $\left(FWHM=4000\,\nicefrac{{km}}{{s}}\right)$ components, as well as the blue shifted and absorption components in $H\beta$ for this object. Reviewing observations for this galaxy from 2005-2013 they concluded that this uncommon emission must be caused by the young ionising cluster wind. In a forthcoming paper we will try to adapt the emission model to better fit this object’s complex structure. FTDTR-10 also shows some disagreement with the helium abundance previously measured. In order to understand this discrepancy, it will be necessary to repeat the chemical analysis using published spectra from other sources such as the SDSS. For this paper, however, we exclude these galaxies from the $Y_{P}-Z$ regressions. The simulation results for these HIIGs are also available online.

The methodology applied in this paper presents the advantage of being able to include (or exclude) lines very easily and to quantify their impact on the target abundances. This includes the chlorine and neon abundances not measured in the current study. As discussed in section 4, this algorithm computes ionic transition emissivity via the traditional approach: parameterised equations as a function of the electron density and temperature. In contrast, the treatment in F2018 consisted in a bilinear interpolation in a very fine grid generated with PyNeb for each transition. For the hydrogen and metals considered both approaches provide results which agree within 1% for the complete temperature and density domains considered. This is not the case for helium, that depends on temperature and density in a remarkably more complex manner. Some of the parametric equations considered for the Porter et al. (2012) recombination coefficients include those from Porter et al. (2007), Olive & Skillman (2004) and Pérez-Montero (2017). None of them guarantee a precision better than 1% for the region of interest. In order to improve the fitting, three temperature domains were considered: the $8000\,K-12000\,K$ , $12000\,K-15000K$ and $15000K-22000K$ ranges. Even with this treatment, however, the discrepancy for some HeI emissivities goes up to $5\%$ at uneven $T_{e}$ and $n_{e}$ intervals. The ideal solution consists of introducing a bilinear interpolation for the emissivity calculation. Such an approach has been implemented by Foreman-Mackey et al. (2019) in the exoplanet library. This will be the approach to follow in the next algorithm iteration.

The simultaneous fitting of all the emission line fluxes provides a more realistic analysis on how each one behaves. In general, the helium lines were more likely to display a missmatch with their expected flux values than those from the oxygen or sulfur ions. As explained in section 2 the flux absorption from the underlying stellar population varies with wavelength. This is a plausible cause for the disagreement in the $HeI$ lines. A large discrepancy between the $HeI$ lines was observed in the two starforming regions in $MRK36$ . This galaxy is the only one in our sample where two bursts were observed within the same slit. It is possible that some cross talk in the spectra extraction is contributing to the mismatch between the helium lines measurement. Additionally, the scheme applied in F2018 to fit the HIIGs continua does propagate the uncertainty introduced in the $HeI$ lines intensity. Therefore, the current algorithm does not quantify which recombination lines are more contaminated by the underlying stellar continua. There are additionally processes contributing to the individual temperatures. For example, the $HeI\lambda 5876\text{\AA}$ intensity may be contaminated by the $NaI$ $D_{1}$ and $D_{2}$ lines. Finally, a fluorescence contribution could not be properly fitted by the available $HeI$ lines. These issues, however, are also affecting the F2018 and therefore the present $y$ and $y^{2+}$ measurements agree very well with those published in F2018.

To the best of our knowledge, this is the first algorithm capable of fitting both the recombination and collisionally excited line spectra in a parameter space including the electron density and temperature. There exist, however, many multi-dimensional chemical samplers based on strong lines. An example of these libraries are HII-CHI-MISTRY by Pérez-Montero (2014), IZI by Blanc et al. (2014), BOND by Vale Asari et al. (2016) or GAME by Ucci et al. (2018). Unlike the algorithm presented here, these tools fit the input fluxes from photoionisation model grids rather than a direct computation. These techniques provide their own advantages. For example, they make possible the chemical analysis of high redshift objects, they allow the modelling of more complex temperature and density structures or can help to solve the double value nature of some strong line diagnostics. There is a lot to be learned from the techniques applied by these tools to explore photoionisation grids. However, due to the high precision needed for determining $Y_{P}$ , strong line methods are not suitable for this work. The strong line diagnostic strategy will be explored in future developments of this algorithm. For example, instead of using the fluorescence correction from Benjamin et al. (2002), the algorithm could derive a better correction from photoionisation grids tailored for each object. Currently, the closest analogues to our algorithm are found in studies to determine $Y_{P}$ .

Izotov & Thuan (2004) and Izotov et al. (2007, 2013) and references therein developed a self-consistent helium abundance determination method. This scheme consists in a Monte Carlo process chain to account for the systematic effects on the $HeI$ lines. Izotov and colleagues vary stochastically the $N_{e}\left(He^{+}\right)$ , $T_{e}\left(He^{+}\right)$ and $\tau\left(3889\right)$ over a range of expected values to compute $y^{+}$ . Afterwards, the best solution for the helium abundance is fitted via a $\chi^{2}$ test. Each emission line is assigned a different underlying stellar absorption and the helium lines are included in the $\chi^{2}$ minimisation. The methodology introduced by Olive & Skillman (2001) and Olive & Skillman (2004) and later enhanced in Aver (2012); Aver et al. (2013) and references therein, also follows a frequentist paradigm applied to Izotov and collaborators data. Their algorithm fits the nine dimensions space simultaneously via a MCMC sampler. Both authors, however, emphasize how the $\chi^{2}$ fitting quality depends on the available $HeI$ lines sensitivity to model parameters. For example, the data in Izotov et al. (2014) reaches the $HeI\lambda 10830\text{\AA}$ line which, unlike most recombination features, is very sensitive to the electron density. Aver et al. (2015) concluded after using these observations in their algorithm, that the uncertainty in $Y_{P}$ decreased by more than $50\%$ .

As a final quality check, we may compare our results with the ones by Aver et al. (2015). This can be easily done by adding two additional dimensions to the model: $T_{e}\left(He\right)$ and $n_{e}\left(He\right)$ . These are the electron temperature and density contributing exclusively to the $HeI$ and $HeII$ emissivity. In this $5^{th}$ test case, the $HeI\lambda 3889$ and $10830\text{\AA}$ lines have been included in the analysis in order to make the simulations as similar as possible to each other. Table 7 compares the results from the synthetic test case in Aver et al. (2015) with the one solved using our HMC algorithm. This table does not include the metals electron density and temperatures nor the ionic abundances contributing to the collisionally excited lines since Aver et al. (2015) does not include them. Similarly, the absorption on hydrogen $\left(a_{H}\right)$ and helium $\left(a_{He}\right)$ are not tabulated since our algorithm does not cover them. Once this is done, the number of parameters for the recombination lines is very similar: the algorithm from Aver and collaborators includes a correction for collisional excitation on the hydrogen lines, which our current model does not. This correction consists in a parametrisation by Anderson et al. (2000) and Anderson et al. (2002), which depends on $T_{e}$ and the ratio of neutral to ionized hydrogen atoms, $\xi$ . On the other hand, our algorithm includes the $y^{2+}$ abundance which is measured from the $HeII4686\text{\AA}$ line. We conclude from the results that both methods have similar accuracy and precision. The precision on the helium temperature and density is not as good as for the collisionally excited lines. This is expected due to the low sensitivity of the recombination lines to these parameters. Indeed this is the reason why the helium abundance is measured to such a good accuracy despite the relatively large uncertainty on the electron temperature. This is, though, not good enough for the primordial helium abundance determination, where ideally one would like to reach accuracies better than $1\%.$ A more suitable prior for the helium temperature could be designed by modelling the $T\left(He\right)$ priors from the available $T_{e}\left[SIII\right]$ or $T_{e}\left[OIII\right]$ as it was done by Peimbert & Peimbert (2002). Additionally, the width in this temperature prior could also be added as a model dimension. This would provide a quantification of the temperature fluctuations, and thus, a direct comparison with the $Y_{P}$ regression by Peimbert et al. (2017).

5.2 $Y_{P}$ regression

The helium mass fractions displayed in Table 6 were calculated using the same procedure as in F2018. In summary, for the oxygen and nitrogen $Y_{P}$ regressions, the $Y_{P,\,O}=f\left(\nicefrac{{HeI}}{{H}},\,\nicefrac{{O}}{{H}}\right)$ formulation from Peimbert & Torres-Peimbert (1974) and Lequeux et al. (1979) was applied. For the regression with sulfur the classical relation was modified to:

[TABLE]

where the $\nicefrac{{O}}{{S}}$ ratio is computed from the sulfur to oxygen relation from the low metallicity BCDs (Blue Compact Dwarfs) sample published by Dors Jr. et al. (2016): $log\left(\nicefrac{{S}}{{O}}\right)=-1.53\pm 0.03$ . As it was done in F2018 the $Y_{P,\,O}$ value is used in the $Y-N$ regression since the oxygen abundance is necessary to compute the nitrogen one. A Monte Carlo algorithm described in F2018 was used to perform the $Y-O$ , $Y-N$ and $Y-S$ regressions, as well as the multivariable regression using the three metals.

Fig.9 shows the regressions using oxygen, nitrogen and sulfur. The results from the linear fitting are listed in Table 8. The first four rows correspond to the primordial helium abundance computations using the chemical abundances from the Bayesian algorithm. The following four rows correspond to the values published in F2018 using the traditional methodology. In the final rows in Table 8 we reproduce recent $Y_{P}$ determinations from the literature. It can be concluded from this comparison that the $Y_{P}$ determinations from both methodologies agree within their uncertainties. The largest disagreement may be found in the $Y_{P,\,N}$ measurement. This could be due to having excluded the $\left[NII\right]\lambda 6548\text{\AA}$ line in some of the objects due to discrepancies with the observed values. The current $Y_{P,\,N}$ result is closer to the oxygen and sulfur regressions than the one given in F2018. Slightly lower values for $Y_{P,\,O}$ and $Y_{P,\,S}$ can be appreciated also. This can be explained by the helium emissivity computation as discussed previously. Nevertheless, the results match very well those from the traditional determinations. Noticeably, the uncertainty has decreased for all the $Y_{P}$ determinations. This is a remarkable achievement, specially once it is considered how quickly and cleanly the new methodology can complete a full direct method chemical analysis. Since the new scheme fits all the ionic species simultaneously we take the $Y_{P,\,O-N-S}$ determination as the preferred one:

[TABLE]

This result is consistent with the determinations from the Planck experiment in a framework provided by a standard Big Bang cosmology.

6 Conclusions

The nebular abundances computed using the traditional direct method are compared with those from a new Bayesian algorithm which simultaneously fits a 14 parameters chemical model. The emission line fluxes belong to the 21 HIIGs sample presented in F2018 where the recombination lines have been corrected by the underlying stellar population. The model parameters are one electron density, two electron temperatures (for the low and high ionization species), the logarithmic extinction coefficient at $H\beta$ , the optical depth for the $HeI$ transitions and nine ionic abundances: $Ar^{2+}$ , $Ar^{3+}$ , $y^{+}$ , $y^{2+}$ , $O^{+}$ , $O^{2+}$ , $N^{+}$ , $S^{+}$ and $S^{2+}$ . The sampling of this relatively big parameter space was successful thanks to machine learning tools. The new HMC sampling shortens the simulation length from several hours to a couple of minutes while providing a stable convergence in contrast to the better known MCMC samplers. The main conclusions from this analysis are:

•

The direct method was adapted to a Bayesian paradigm by a careful design for the model parameter priors. Uninformative (wide) probability distributions were successfully applied for the electron temperatures, the ionic abundances and the extinction coefficients. In contrast, informative priors were necessary for the electron density and the optical depth. For $n_{e}$ this was necessary due to the fact that for these low density HIIGs sample (below the low density limit) none of the available collisionally excited ratios are good diagnostic tools. On the other hand, The $HeI$ lines available in our data are weakly affected by the fluorescence excitation. Consequently, the optical depth which parametrises it could not be fitted in these observations. We still took it into consideration while calculating the HeI fluxes as if it were a standard Monte Carlo simulation.

•

A set of synthetic test cases with increasing number of model parameters are presented to confirm the accuracy and precision of the sampler. Once the algorithm is applied on real observations the behaviour and correlation observed on the parameters remain constant. An issue, however, was encountered during the fitting of real observations due to the high precision in the strong lines such as $\left[OIII\right]\lambda 5007\text{\AA}$ and $H\alpha$ . A minimum accuracy of $2\%$ was imposed in the observed fluxes. The precision in the results is as good, though, as the one expected from these high quality observations.

•

In order to test the success of the simulation, we propose a graphical comparison between the input fluxes and the output flux distributions from the fitting. We found in general a very good match between the two. Some objects, however, displayed a sizeable mismatch in the $H\gamma$ and $\left[NII\right]\lambda 6548\text{\AA}$ fluxes which can be explained by technical constraints in the observations. In those cases, the lines were excluded from the fitting.

•

In general, we obtain a very good match with the abundances published in F2018. Nevertheless, the new values for the logarithmic extinction coefficient were higher for some objects and with a consistent larger uncertainty. This can be explained by the difference in procedure: while in the traditional methodology only $H\alpha$ and $H\beta$ end up contributing to the $c\left(H\beta\right)$ calculation in this new analysis all the lines are. In spite of this, the abundances obtained are similar, as the extinction for these objects is relatively small.

•

An additional synthetic test case was presented in order to reproduce the test case in Aver et al. (2015). The $Te\left(He\right)$ and $n_{e}\left(He\right)$ are added as additional parameters in the $HeI$ and $HeII$ fluxes calculation. Despite the difference in mathematical and computational schemes our algorithm can replicate the accuracy and precision in the measurement of the parameters in common. It is argued, however, that given the low dependency of the recombination lines on electron temperature and density, the input prior design may have a larger weight in the fitting results. In the near future we plan to explore new observations and theoretical models to better constrain the $T_{e}\left(He\right)$ and $n_{e}\left(He\right)$ measurement.

•

The primordial helium abundance measurement from the new chemical abundances is in good agreement with those in F2018. The biggest disagreement was found for the $Y-N$ regression which resulted in a lower $Y_{P,\,N}$ value obtained in the present work. This is a consequence of the excluded $\left[NII\right]\lambda 6548\text{\AA}$ lines due to a mismatch with the model fluxes. The present result, however, is in good agreement with $Y_{P,\,O}$ and $Y_{P,\,S}$ . Since this technique fits all the ionic abundances simultaneously, we take the multivariable linear regression as our chosen measurement giving $Y_{P,\,O-N-S}=0.243\pm 0.005$ . This result is consistent with Standard Big Bang Nucleosynthesis.

We obtain in general lower uncertainties than we did with the classical direct method. This new chemical analysis allows to increase the complexity of the theoretical model. In future work, we will focus in including third party data grids during the sampling process. This will improve the emissivity computation for the helium lines which was one of the challenges found in this paper. Finally, we hope to include the continua computation in the fitting. This task has two purposes: a better quantification of the stellar continuum absorption on the helium lines and a better estimation of $T\left(He\right)$ via the nebular continuum jumps. The algorithms developed in this work are not yet available in the standard python distribution channels. However, all the scripts can be found at the following github account until they are properly published (https://github.com/Vital-Fernandez).

Acknowledgements

Vital Fernández is indebted to Enrique Pérez and Daniel Miralles for their generous discussions on the nebular and stellar continua determination, respectively. Similarly, he is grateful to Erick Aver for his prompt replies on the physical and computational characteristics of his helium abundance determination algorithm.

We are indebted to an anonymous referee whose detailed comments contributed to the improvement of the paper.

This work could not have been accomplished without the work and generosity from the PyMC3 development team. V. F. also expresses his gratitude to Dan Foreman-Mackey for his insight on how to use data science algorithms in astronomical research.

We thank the Spanish allocation committee (CAT) for awarding observing time and the cheerful technical support from the observatory personnel. Vital Fernández is grateful to the Mexican research Council (CONACYT) for supporting this research through studentship 554031/300844 and Elena and Roberto Terlevich acknowledge CONACYT for research grant CB-2008-103365. This work has been supported by DGICYT grants AYA2013-47742-C4-3-P and AYA2016-79724-C4-1-P. Partial financial support came also from project SELGIFS: PIRSES-GA- 2013-612701-SELGIFS. Vital Fernández and Elena and Roberto Terlevich are grateful to the hospitality of the Departamento de Física Teórica at the Universidad Autónoma de Madrid, Spain during a visit to advance with this study.

Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the US Department of Energy, the Japanese Monbukagakusho and the Max Planck Society. The SDSS Web site is http://www.sdss.org.

This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Appendix A Extra material

This section includes the tabulated results for the galaxy SHOC022 fitting presented in the discussion, as well as the atomic data references and their parametrisations.

Bibliography69

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Aggarwal & Keenan (2000) Aggarwal K. M., Keenan F. P., 2000, Vizie R Online Data Catalog, 212
2Anderson et al. (2000) Anderson H., Ballance C. P., Badnell N. R., Summers H. P., 2000, Journal of Physics B: Atomic, Molecular and Optical Physics , 33, 1255 · doi ↗
3Anderson et al. (2002) Anderson H., Ballance C. P., Badnell N. R., Summers H. P., 2002, Journal of Physics B: Atomic, Molecular and Optical Physics , 35, 1613 · doi ↗
4Aver (2012) Aver E., 2012, Ph D thesis, University of Minnesota
5Aver et al. (2011) Aver E., Olive K. A., Skillman E. D., 2011, Journal of Cosmology and Astroparticle Physics , 2011, 043 · doi ↗
6Aver et al. (2013) Aver E., Olive K. A., Porter R. L., Skillman E. D., 2013, Journal of Cosmology and Astro-Particle Physics , 2013, 017 · doi ↗
7Aver et al. (2015) Aver E., Olive K. A., Skillman E. D., 2015, ar Xiv:1503.08146 [astro-ph]
8Benjamin et al. (2002) Benjamin R. A., Skillman E. D., Smits D. P., 2002, The Astrophysical Journal , 569, 288 · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

A Bayesian direct method implementation to fit emission line spectra: Application to the primordial He abundance determination

Abstract

keywords:

1 Introduction

2 The data

3 Bayesian direct method.

4 Model test cases

5 Results and discussion

5.1 Comparison between methods

5.2 YPY_{P}YP​ regression

6 Conclusions

Acknowledgements

Appendix A Extra material

5.2 $Y_{P}$ regression