A blinded determination of $H_0$ from low-redshift Type Ia supernovae, calibrated by Cepheid variables

TL;DR

This study performs a blinded reanalysis of low-redshift Type Ia supernovae data calibrated by Cepheid variables to measure the Hubble constant, aiming to address discrepancies between local and CMB-based measurements.

Contribution

It introduces a blinded, end-to-end analysis incorporating covariance matrices and simultaneous fitting, providing an independent validation of previous $H_0$ estimates with a detailed error budget.

Findings

Measured $H_0$ as 72.5 km s$^{-1}$ Mpc$^{-1}$ with 4.4% total uncertainty.

Identified statistical errors as the dominant source of uncertainty.

Systematic errors mainly from Cepheid fit variations and host galaxy effects.

Abstract

Presently a ${>}3\sigma$ tension exists between values of the Hubble constant $H_0$ derived from analysis of fluctuations in the Cosmic Microwave Background by Planck, and local measurements of the expansion using calibrators of type Ia supernovae (SNe Ia). We perform a blinded reanalysis of Riess et al. 2011 to measure $H_0$ from low-redshift SNe Ia, calibrated by Cepheid variables and geometric distances including to NGC 4258. This paper is a demonstration of techniques to be applied to the Riess et at. 2016 data. Our end-to-end analysis starts from available CfA3 and LOSS photometry, providing an independent validation of Riess et al. 2011. We obscure the value of $H_0$ throughout our analysis and the first stage of the referee process, because calibration of SNe Ia requires a series of often subtle choices, and the potential for results to be affected by human bias is significant. Our analysis departs from that of Riess et al. 2011 by incorporating the covariance matrix method adopted in SNLS and JLA to quantify SN Ia systematics, and by including a simultaneous fit of all SN Ia and Cepheid data. We find $H_0 = 72.5 \pm 3.1$ (stat) $\pm 0.77$ (sys) km s$^{-1}$ Mpc$^{-1}$ with a three-galaxy (NGC 4258+LMC+MW) anchor. The relative uncertainties are 4.3% statistical, 1.1% systematic, and 4.4% total, larger than in Riess et al. 2011 (3.3% total) and the Efstathiou 2014 reanalysis (3.4% total). Our error budget for $H_0$ is dominated by statistical errors due to the small size of the supernova sample, whilst the systematic contribution is dominated by variation in the Cepheid fits, and for the SNe Ia, uncertainties in the host galaxy mass dependence and Malmquist bias.

Tables

Table 1. Table 1: Recent nearby SNe Ia and their host galaxies used in R11, along with observations of Cepheids in these galaxies.

Galaxy	SN Ia	$N_{\mathrm{Cepheids}}$
NGC 4536	SN 1981B	69
NGC 4639	SN 1990N	32
NGC 3370	SN 1994ae	79
NGC 3982	SN 1998aq	29
NGC 3021	SN 1995al	26
NGC 1309	SN 2002fk	36
NGC 5584	SN 2007af	95
NGC 4038	SN 2007sr	39
NGC 4258	-	165
Total		570

Table 2. Table 2: Summary of selected Cepheid fits to carry forward to global fit (i.e. rejection, anchors, and period cut used). The positions of the best fit values for b W subscript 𝑏 𝑊 b_{W} and Z W subscript 𝑍 𝑊 Z_{W} in the b W , Z W subscript 𝑏 𝑊 subscript 𝑍 𝑊 b_{W},Z_{W} -plane (represented in Figure  1 ) are also given, as well as the symbols for these fits in Figure  1 . The top half of the table (solid symbols) lists fits with the P < 60 𝑃 60 P<60  day cut, whilst the bottom half (empty symbols) contains fits without.

	Rej ($T$)	Anchora	Symbol
top left	2.25	n4258+LMC+MW	solid green diamond
top left	2.5	n4258+LMC+MW	solid turquoise diamond
top left	R11	n4258+LMC+MW	solid orange diamond
top left	R11	LMC+MW	solid orange $\Delta$
middle	2.5	n4258+LMC	solid turquoise hexagon
middle	R11	MW	solid orange star
lower left	2.25	n4258+MW	solid green square
lower	2.25	LMC	solid green $\nabla$
lower	2.25	n4258	solid green circle
top	2.25	n4258+LMC+MW	empty green diamond
top	2.5	n4258+LMC+MW	empty turquoise diamond
top	R11	n4258+LMC+MW	empty orange diamond
top	R11	LMC+MW	empty orange $\Delta$
top	2.25	LMC+MW	empty green $\Delta$
middle	2.25	n4258+MW	empty green square
right	R11	n4258	empty orange circle
right	2.5	n4258	empty turquoise circle
lower right	2.25	n4258	empty green circle

Table 3. Table 3: Observations of nearby SNe Ia in Table  1 , including sources of photometry, SALT2 instruments, magnitude systems (including filters) where available. Lightcurves of the two earliest supernovae were given as standard photometry only.

SN Ia	Photometry source	Magnitude system and filters
SN 1981B	Buta & Turner (1983)	Standard UBVR
SN 1990N	Lira et al. (1998)	Standard UBVRI
SN 1994ae	Riess et al. (2005)	AndyCama BVRI
SN 1998aq	Riess et al. (2005)	4Shooter/AndyCam UBVRI
SN 1995al	Riess et al. (2009a)	AndyCam UBVRI
SN 2002fk	CfA3b	4Shooter2 UBVRI
	LOSSc	KAIT3/NICKEL BVRI
SN 2007af	CfA3	Keplercam BVri
	LOSS	KAIT3/KAIT4 BVRI
SN 2007sr	CfA3	Keplercam BVri
	LOSS	KAIT3/4 BVRI

Table 4. Table 4: Results of preliminary SN-only fits for various cuts.

SN cut	$N_{\mathrm{SN}}$	$\alpha$	$\beta$	$ℳ$
default	171	0.164 (0.013)	3.07 (0.14)	-3.240 (0.036)
higher ${\chi }^2$	175	0.167 (0.013)	3.12 (0.13)	-3.244 (0.036)
lower ${\chi }^2$	163	0.158 (0.013)	3.04 (0.16)	-3.256 (0.038)
$z>0.0233$	96	0.163 (0.016)	2.73 (0.17)	-3.252 (0.038)
stricter $C$	160	0.158 (0.015)	2.93 (0.18)	-3.238 (0.037)
str. ${\sigma }_{X_1},{\sigma }_C$	164	0.171 (0.014)	3.10 (0.14)	-3.232 (0.038)
str. ${\sigma }_{X_1}$	165	0.171 (0.013)	3.10 (0.15)	-3.245 (0.037)
str. $E\left(B-V\right)$	166	0.167 (0.013)	3.06 (0.15)	-3.241 (0.036)
d	187	0.165 (0.013)	3.11 (0.14)	-3.234 (0.035)

Table 5. Table 5: Global fit results for supernova parameters α , β 𝛼 𝛽 \alpha,\beta for each SN cut, averaged over Cepheid fits. The default rejection (bolded) detailed in Section  5.3 is chosen as our reference fit.

SN cut	$\alpha$	$\beta$
default	0.165	3.09
higher ${\chi }^2$	0.167	3.134
$z>0.0233$	0.162	2.759
lower ${\chi }^2$	0.158	3.057
stricter $C$	0.156	2.974
stricter ${\sigma }_{X_1}$	0.171	3.106

Table 6. Table 6: Global fit results for Cepheid parameters { b W , Z W , M W } subscript 𝑏 𝑊 subscript 𝑍 𝑊 subscript 𝑀 𝑊 \{b_{W},Z_{W},M_{W}\} for each Cepheid fit, averaged over SN cuts. The bolded fit ( T = 2.25 𝑇 2.25 T=2.25 rejection, all three anchors, and no upper cut on the period) is chosen as our reference fit.

Rejection	Distance anchor	d	$b_W$	$Z_W$	$M_W$
2.25	alla	Y	-3.28	-0.19	-5.95
2.5	all	Y	-3.28	-0.21	-5.96
R11	all	Y	-3.26	-0.14	-5.95
2.5	n4258+LMC	Y	-3.23	-0.45	-6.11
2.25	n4258+MW	Y	-3.31	-0.52	-5.89
R11	LMC+MW	Y	-3.25	-0.11	-5.91
2.25	n4258	Y	-3.23	-0.55	-6.03
2.25	LMC	Y	-3.24	-0.55	-6.16
R11	MW	Y	-3.21	-0.35	-5.83
2.25	all	N	-3.17	-0.11	-5.95
2.5	all	N	-3.20	-0.10	-5.96
R11	all	N	-3.21	-0.06	-5.94
2.25	n4258+MW	N	-3.16	-0.41	-5.89
2.25	LMC+MW	N	-3.16	-0.07	-5.90
R11	LMC+MW	N	-3.20	-0.03	-5.89
2.25	n4258	N	-3.04	-0.44	-6.10
2.5	n4258	N	-3.06	-0.34	-6.11
R11	n4258	N	-3.09	-0.23	-6.08

Table 7. Table 7: Global fit results for degenerate parameters M B subscript 𝑀 𝐵 M_{B} and ℋ ℋ \mathcal{H} , averaged over SN cuts. The bolded fit ( T = 2.25 𝑇 2.25 T=2.25 rejection, all three anchors, and no upper cut on the period) is chosen as our reference fit.

Rejection	Distance anchor	d	$M_B$	$\mathscr{H}$
2.25	alla	Y	-18.953	-15.709
2.5	all	Y	-18.967	-15.722
R11	all	Y	-18.932	-15.687
2.5	n4258+LMC	Y	-19.061	-15.816
2.25	n4258+MW	Y	-18.892	-15.647
R11	LMC+MW	Y	-18.889	-15.644
2.25	n4258	Y	-18.988	-15.743
2.25	LMC	Y	-19.122	-15.877
R11	MW	Y	-18.759	-15.514
2.25	all	N	-18.929	-15.685
2.5	all	N	-18.953	-15.708
R11	all	N	-18.924	-15.679
2.25	n4258+MW	N	-18.859	-15.614
2.25	LMC+MW	N	-18.875	-15.631
R11	LMC+MW	N	-18.868	-15.623
2.25	n4258	N	-18.996	-15.751
2.5	n4258	N	-19.015	-15.770
R11	n4258	N	-18.970	-15.725

Table 8. Table 8: Summary of uncertainties in ℋ ℋ \mathcal{H} from Section  6.2.3 (Equation  26 and 27 ), converted to relative errors in H 0 subscript 𝐻 0 H_{0} using Equation  28 , and added in quadrature in line with R11. The statistical error below are those reported by MultiNest (Table  13 ) and include terms (i) and (ii) described at the start of Section  6.3 . The systematic error is from the variation between fit results with different choices of Cepheid fits, and secondarily SN cuts, i.e. term (iii).

Anchor	all	NGC 4258 only
$\mathscr{H}$	-15.698	-15.748
${\sigma }_{\mathscr{H}}$
Statistical	0.093	0.107
Systematic	0.023	0.023
Relative $H_0$ error (%)
Statistical	4.3	4.9
Systematic	1.1	1.1
Total	4.4	5.0

Table 9. Table 9: Relative contributions to the uncertainty in H 0 subscript 𝐻 0 H_{0} (i.e. the variance) from individual statistical and systematic sources uncertainties, calculated as described in Betoule et al. ( 2014 , Section 6.2) .

Source of	Relative contribution	Described in
uncertainty	to ${\sigma }^2\left(H_0\right)$ (%)	Appendix
Statistical:
Lightcurves	62.1	C.2
SALT2 training	1.2	C.2
Total statistical	63.3
Systematic:
Malmquist bias	13.7	B.2,C
Host galaxy	13.0	C.5
Lightcurve model	6.8	C
Calibration	3.1	C.4
Peculiar velocities	0.04	C.6
Milky Way extinction	0.03	C.3
Total systematic	36.7

Table 10. Table 10: Systematics in κ 𝜅 \kappa

Instrument	Filters	ZP index	${\lambda }_{\mathrm{eff}}$ index
MegaCam	griz	0-3	50-53
Standard	UBVRI	4-8	54-58
KeplerCam	UsBVria	9-13	59-63
4Shooter2	UsBVRI	14-18	64-68
Swope	ugriBV	19-24	69-74
SDSS	ugriz	25-29	75-79
KAIT1-4	BVRI	30-45	80-95
NICKEL	BVRI	46-49	96-99

Table 11. Table 11: Results of the Cepheid-only fits described in Section  4.3 from each combination of distance anchor, rejection and period cut for each Cepheid fit. The best fit Cepheid parameters { b W , Z W , M W } subscript 𝑏 𝑊 subscript 𝑍 𝑊 subscript 𝑀 𝑊 \{b_{W},Z_{W},M_{W}\} are given, as well as the number of Cepheids remaining after rejection and intrinsic scatter.

Rejection	Distance anchor	d	$N_{\mathrm{cepheids}}$	${\sigma }_{\mathrm{int},C}$	$b_W$	$Z_W$	$M_W$
$T=2.25$	n4258	Y	439	0.17	-3.23 (0.07)	-0.54 (0.13)	-6.03 (0.07)
$T=2.5$	n4258	Y	463	0.27	-3.22 (0.08)	-0.49 (0.14)	-6.06 (0.07)
R11	n4258	Y	379	0.21	-3.18 (0.07)	-0.32 (0.14)	-6.05 (0.08)
$T=2.25$	LMC	Y	439	0.17	-3.24 (0.05)	-0.54 (0.13)	-6.16 (0.07)
$T=2.5$	LMC	Y	464	0.27	-3.24 (0.05)	-0.50 (0.14)	-6.14 (0.08)
R11	LMC	Y	379	0.21	-3.22 (0.05)	-0.32 (0.14)	-6.07 (0.07)
$T=2.25$	MW	Y	439	0.17	-3.24 (0.07)	-0.54 (0.13)	-5.83 (0.05)
$T=2.5$	MW	Y	463	0.27	-3.24 (0.07)	-0.49 (0.15)	-5.83 (0.05)
R11	MW	Y	379	0.21	-3.20 (0.07)	-0.32 (0.15)	-5.82 (0.05)
$T=2.25$	n4258+LMC	Y	439	0.18	-3.23 (0.05)	-0.46 (0.11)	-6.10 (0.05)
$T=2.5$	n4258+LMC	Y	466	0.28	-3.23 (0.05)	-0.42 (0.12)	-6.10 (0.05)
R11	n4258+LMC	Y	379	0.21	-3.22 (0.05)	-0.29 (0.12)	-6.05 (0.05)
$T=2.25$	n4258+MW	Y	437	0.17	-3.31 (0.06)	-0.50 (0.12)	-5.89 (0.04)
$T=2.5$	n4258+MW	Y	464	0.27	-3.30 (0.07)	-0.46 (0.14)	-5.90 (0.04)
R11	n4258+MW	Y	379	0.21	-3.26 (0.06)	-0.30 (0.14)	-5.89 (0.04)
$T=2.25$	LMC+MW	Y	435	0.16	-3.27 (0.05)	-0.12 (0.10)	-5.91 (0.05)
$T=2.5$	LMC+MW	Y	464	0.28	-3.27 (0.05)	-0.15 (0.11)	-5.92 (0.06)
R11	LMC+MW	Y	379	0.21	-3.25 (0.05)	-0.08 (0.11)	-5.90 (0.06)
$T=2.25$	n4258+LMC+MW	Y	434	0.16	-3.28 (0.05)	-0.17 (0.10)	-5.95 (0.04)
$T=2.5$	n4258+LMC+MW	Y	463	0.27	-3.28 (0.05)	-0.19 (0.10)	-5.96 (0.04)
R11	n4258+LMC+MW	Y	379	0.21	-3.26 (0.05)	-0.12 (0.11)	-5.94 (0.04)
$T=2.25$	n4258	N	521	0.2	-3.04 (0.05)	-0.42 (0.12)	-6.10 (0.07)
$T=2.5$	n4258	N	540	0.26	-3.06 (0.06)	-0.32 (0.13)	-6.11 (0.07)
R11	n4258	N	444	0.21	-3.09 (0.06)	-0.21 (0.13)	-6.08 (0.07)
$T=2.25$	LMC	N	523	0.21	-3.11 (0.04)	-0.39 (0.12)	-6.12 (0.07)
$T=2.5$	LMC	N	544	0.28	-3.12 (0.04)	-0.26 (0.13)	-6.06 (0.07)
R11	LMC	N	444	0.21	-3.13 (0.04)	-0.20 (0.13)	-6.04 (0.07)
$T=2.25$	MW	N	521	0.2	-3.07 (0.05)	-0.42 (0.12)	-5.80 (0.05)
$T=2.5$	MW	N	539	0.26	-3.09 (0.05)	-0.30 (0.13)	-5.81 (0.05)
R11	MW	N	444	0.21	-3.12 (0.06)	-0.20 (0.13)	-5.81 (0.05)
$T=2.25$	n4258+LMC	N	523	0.21	-3.11 (0.04)	-0.37 (0.11)	-6.09 (0.05)
$T=2.5$	n4258+LMC	N	539	0.26	-3.12 (0.04)	-0.30 (0.11)	-6.08 (0.05)
R11	n4258+LMC	N	444	0.21	-3.13 (0.04)	-0.21 (0.12)	-6.05 (0.05)
$T=2.25$	n4258+MW	N	520	0.2	-3.16 (0.05)	-0.40 (0.12)	-5.89 (0.04)
$T=2.5$	n4258+MW	N	538	0.26	-3.16 (0.05)	-0.29 (0.13)	-5.90 (0.04)
R11	n4258+MW	N	444	0.21	-3.17 (0.06)	-0.17 (0.13)	-5.89 (0.04)
$T=2.25$	LMC+MW	N	519	0.2	-3.14 (0.04)	-0.05 (0.10)	-5.89 (0.05)
$T=2.5$	LMC+MW	N	546	0.29	-3.14 (0.04)	-0.01 (0.10)	-5.89 (0.06)
R11	LMC+MW	N	444	0.21	-3.16 (0.04)	-0.01 (0.10)	-5.88 (0.05)
$T=2.25$	n4258+LMC+MW	N	517	0.19	-3.15 (0.04)	-0.09 (0.09)	-5.95 (0.04)
$T=2.5$	n4258+LMC+MW	N	544	0.28	-3.16 (0.04)	-0.08 (0.10)	-5.96 (0.04)
R11	n4258+LMC+MW	N	444	0.21	-3.17 (0.04)	-0.04 (0.10)	-5.94 (0.04)

Table 12. Table 12: Global fit results for distance modulus offsets { Δ ​ μ i } Δ subscript 𝜇 𝑖 \{\Delta\mu_{i}\} for the eight galaxies.

Cepheid rejection	Distance anchor	d	$\Delta {\mu }_{4536}$	$\Delta {\mu }_{4639}$	$\Delta {\mu }_{3370}$	$\Delta {\mu }_{3982}$	$\Delta {\mu }_{3021}$	$\Delta {\mu }_{1309}$	$\Delta {\mu }_{5584}$	$\Delta {\mu }_{4038}$
2.25	n4258+LMC+MW	Y	1.56 (0.05)	2.34 (0.07)	2.77 (0.06)	2.43 (0.07)	2.86 (0.09)	3.30 (0.07)	2.37 (0.05)	2.31 (0.09)
2.5	n4258+LMC+MW	Y	1.58 (0.06)	2.41 (0.08)	2.79 (0.06)	2.47 (0.09)	2.89 (0.09)	3.30 (0.08)	2.36 (0.06)	2.35 (0.10)
R11	n4258+LMC+MW	Y	1.54 (0.06)	2.35 (0.08)	2.78 (0.06)	2.43 (0.08)	2.86 (0.09)	3.18 (0.08)	2.37 (0.06)	2.32 (0.10)
2.5	n4258+LMC	Y	1.59 (0.06)	2.44 (0.08)	2.79 (0.07)	2.49 (0.09)	2.89 (0.09)	3.31 (0.08)	2.38 (0.06)	2.39 (0.10)
2.25	n4258+MW	Y	1.52 (0.05)	2.37 (0.07)	2.73 (0.06)	2.34 (0.08)	2.87 (0.09)	3.26 (0.07)	2.33 (0.06)	2.32 (0.09)
R11	LMC+MW	Y	1.56 (0.06)	2.37 (0.08)	2.80 (0.06)	2.45 (0.08)	2.88 (0.09)	3.20 (0.08)	2.38 (0.06)	2.33 (0.10)
2.25	n4258	Y	1.56 (0.06)	2.39 (0.07)	2.74 (0.06)	2.35 (0.08)	2.88 (0.09)	3.26 (0.07)	2.35 (0.06)	2.34 (0.09)
2.25	LMC	Y	1.56 (0.06)	2.39 (0.07)	2.75 (0.06)	2.35 (0.08)	2.89 (0.09)	3.27 (0.07)	2.35 (0.06)	2.34 (0.09)
R11	MW	Y	1.53 (0.06)	2.36 (0.08)	2.76 (0.07)	2.39 (0.09)	2.85 (0.10)	3.17 (0.08)	2.36 (0.06)	2.33 (0.10)
2.25	n4258+LMC+MW	N	1.64 (0.05)	2.31 (0.07)	2.76 (0.05)	2.39 (0.07)	2.81 (0.08)	3.27 (0.06)	2.40 (0.05)	2.30 (0.07)
2.5	n4258+LMC+MW	N	1.62 (0.06)	2.41 (0.08)	2.79 (0.06)	2.44 (0.08)	2.89 (0.09)	3.30 (0.07)	2.41 (0.05)	2.35 (0.08)
R11	n4258+LMC+MW	N	1.57 (0.06)	2.34 (0.08)	2.78 (0.06)	2.41 (0.08)	2.87 (0.09)	3.21 (0.07)	2.39 (0.05)	2.32 (0.08)
2.25	n4258+MW	N	1.58 (0.05)	2.35 (0.07)	2.72 (0.05)	2.32 (0.08)	2.87 (0.08)	3.24 (0.07)	2.38 (0.05)	2.32 (0.07)
2.25	LMC+MW	N	1.67 (0.05)	2.34 (0.07)	2.79 (0.05)	2.42 (0.07)	2.84 (0.08)	3.30 (0.07)	2.43 (0.05)	2.32 (0.07)
R11	LMC+MW	N	1.59 (0.06)	2.36 (0.08)	2.80 (0.06)	2.43 (0.08)	2.89 (0.09)	3.22 (0.07)	2.41 (0.05)	2.33 (0.08)
2.25	n4258	N	1.64 (0.05)	2.37 (0.07)	2.73 (0.05)	2.34 (0.08)	2.89 (0.08)	3.25 (0.07)	2.38 (0.05)	2.32 (0.07)
2.5	n4258	N	1.61 (0.06)	2.41 (0.08)	2.77 (0.06)	2.40 (0.08)	2.88 (0.09)	3.27 (0.07)	2.40 (0.06)	2.34 (0.08)
R11	n4258	N	1.55 (0.06)	2.34 (0.08)	2.75 (0.06)	2.37 (0.08)	2.85 (0.09)	3.18 (0.07)	2.37 (0.06)	2.31 (0.08)

Table 13. Table 13: Results of all global fits described in Section  6.1 , from each combination of Cepheid fit and SN cut. The SN parameters { α , β } 𝛼 𝛽 \{\alpha,\beta\} , the Cepheid parameters { b W , Z W , M W } subscript 𝑏 𝑊 subscript 𝑍 𝑊 subscript 𝑀 𝑊 \{b_{W},Z_{W},M_{W}\} , and the zero points { ℋ , M B \{\mathcal{H},M_{B} , and μ 4258 } \mu_{4258}\} are displayed. The indicative fits from which statistical uncertainties for parameters are retrieved are bolded.

Cepheid	Distance	d	SN cut	$\alpha$	$\beta$	$\mathscr{H}$	$M_B$	$b_W$	$Z_W$	$M_W$	${\mu }_{4258}$
rejection	anchor
2.25	all	Y	default	0.165 (0.010)	3.09 (0.11)	-15.714 (0.094)	-18.955 (0.089)	-3.28 (0.05)	-0.19 (0.10)	-5.95 (0.04)	29.36 (0.04)
2.25	all	Y	higher ${\chi }^2$	0.167 (0.010)	3.14 (0.11)	-15.706 (0.092)	-18.952 (0.087)	-3.28 (0.05)	-0.19 (0.10)	-5.95 (0.04)	29.36 (0.04)
2.25	all	Y	$z>0.0233$	0.162 (0.012)	2.77 (0.13)	-15.705 (0.093)	-18.954 (0.087)	-3.28 (0.05)	-0.20 (0.10)	-5.95 (0.04)	29.36 (0.04)
2.25	all	Y	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.702 (0.092)	-18.953 (0.087)	-3.28 (0.05)	-0.19 (0.09)	-5.95 (0.04)	29.36 (0.04)
2.25	all	Y	stricter $C$	0.156 (0.011)	2.99 (0.14)	-15.717 (0.092)	-18.956 (0.087)	-3.28 (0.05)	-0.20 (0.10)	-5.95 (0.04)	29.36 (0.05)
2.25	all	Y	stricter ${\sigma }_{X_1}$	0.171 (0.010)	3.11 (0.11)	-15.706 (0.093)	-18.950 (0.088)	-3.28 (0.05)	-0.19 (0.09)	-5.95 (0.04)	29.36 (0.04)
2.5	all	Y	default	0.165 (0.010)	3.09 (0.10)	-15.722 (0.095)	-18.964 (0.090)	-3.28 (0.05)	-0.22 (0.10)	-5.96 (0.04)	29.35 (0.05)
2.5	all	Y	higher ${\chi }^2$	0.166 (0.010)	3.14 (0.11)	-15.722 (0.094)	-18.967 (0.089)	-3.28 (0.05)	-0.21 (0.10)	-5.96 (0.04)	29.35 (0.05)
2.5	all	Y	$z>0.0233$	0.162 (0.012)	2.77 (0.13)	-15.722 (0.096)	-18.971 (0.091)	-3.28 (0.05)	-0.21 (0.10)	-5.97 (0.04)	29.35 (0.05)
2.5	all	Y	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.717 (0.094)	-18.968 (0.089)	-3.29 (0.05)	-0.21 (0.10)	-5.96 (0.04)	29.35 (0.04)
2.5	all	Y	stricter $C$	0.156 (0.011)	2.99 (0.14)	-15.729 (0.095)	-18.968 (0.090)	-3.28 (0.05)	-0.21 (0.10)	-5.96 (0.04)	29.35 (0.05)
2.5	all	Y	stricter ${\sigma }_{X_1}$	0.171 (0.010)	3.11 (0.11)	-15.720 (0.094)	-18.964 (0.090)	-3.28 (0.05)	-0.21 (0.10)	-5.96 (0.04)	29.35 (0.04)
R11	all	Y	default	0.165 (0.010)	3.08 (0.11)	-15.688 (0.096)	-18.930 (0.091)	-3.26 (0.05)	-0.14 (0.10)	-5.95 (0.05)	29.35 (0.05)
R11	all	Y	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.10)	-15.684 (0.095)	-18.929 (0.090)	-3.26 (0.05)	-0.14 (0.11)	-5.95 (0.04)	29.35 (0.05)
R11	all	Y	$z>0.0233$	0.162 (0.013)	2.76 (0.13)	-15.684 (0.094)	-18.933 (0.089)	-3.26 (0.05)	-0.14 (0.11)	-5.95 (0.05)	29.35 (0.04)
R11	all	Y	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.682 (0.093)	-18.934 (0.088)	-3.26 (0.05)	-0.14 (0.10)	-5.95 (0.05)	29.35 (0.05)
R11	all	Y	stricter $C$	0.156 (0.011)	2.97 (0.14)	-15.696 (0.094)	-18.935 (0.089)	-3.26 (0.05)	-0.14 (0.10)	-5.95 (0.04)	29.35 (0.04)
R11	all	Y	stricter ${\sigma }_{X_1}$	0.172 (0.011)	3.11 (0.11)	-15.685 (0.094)	-18.929 (0.090)	-3.26 (0.05)	-0.14 (0.11)	-5.95 (0.04)	29.35 (0.04)
2.25	all	N	default	0.164 (0.010)	3.09 (0.11)	-15.689 (0.093)	-18.929 (0.088)	-3.18 (0.04)	-0.11 (0.09)	-5.95 (0.04)	29.34 (0.04)
2.25	all	N	higher ${\chi }^2$	0.166 (0.010)	3.13 (0.10)	-15.683 (0.092)	-18.928 (0.088)	-3.17 (0.04)	-0.11 (0.09)	-5.95 (0.04)	29.34 (0.04)
2.25	all	N	$z>0.0233$	0.160 (0.012)	2.76 (0.13)	-15.684 (0.091)	-18.933 (0.086)	-3.17 (0.04)	-0.11 (0.09)	-5.95 (0.04)	29.34 (0.04)
2.25	all	N	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.678 (0.092)	-18.929 (0.088)	-3.17 (0.04)	-0.11 (0.09)	-5.95 (0.04)	29.34 (0.04)
2.25	all	N	stricter $C$	0.155 (0.011)	2.98 (0.15)	-15.691 (0.093)	-18.929 (0.088)	-3.17 (0.04)	-0.11 (0.09)	-5.95 (0.04)	29.34 (0.04)
2.25	all	N	stricter ${\sigma }_{X_1}$	0.170 (0.010)	3.11 (0.11)	-15.683 (0.093)	-18.926 (0.088)	-3.17 (0.04)	-0.11 (0.09)	-5.95 (0.04)	29.34 (0.04)
2.5	all	N	default	0.165 (0.010)	3.09 (0.11)	-15.710 (0.093)	-18.951 (0.089)	-3.20 (0.04)	-0.09 (0.10)	-5.96 (0.04)	29.33 (0.04)
2.5	all	N	higher ${\chi }^2$	0.166 (0.010)	3.14 (0.11)	-15.707 (0.093)	-18.952 (0.089)	-3.20 (0.04)	-0.09 (0.10)	-5.96 (0.04)	29.33 (0.04)
2.5	all	N	$z>0.0233$	0.161 (0.012)	2.76 (0.13)	-15.707 (0.094)	-18.955 (0.089)	-3.20 (0.04)	-0.10 (0.10)	-5.96 (0.04)	29.33 (0.04)
2.5	all	N	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.702 (0.093)	-18.953 (0.089)	-3.20 (0.04)	-0.09 (0.10)	-5.96 (0.04)	29.33 (0.04)
2.5	all	N	stricter $C$	0.156 (0.011)	2.98 (0.14)	-15.715 (0.095)	-18.953 (0.090)	-3.20 (0.04)	-0.10 (0.10)	-5.96 (0.04)	29.33 (0.04)
2.5	all	N	stricter ${\sigma }_{X_1}$	0.171 (0.010)	3.11 (0.11)	-15.706 (0.092)	-18.950 (0.088)	-3.20 (0.04)	-0.10 (0.10)	-5.96 (0.04)	29.33 (0.04)
R11	all	N	default	0.165 (0.010)	3.09 (0.11)	-15.682 (0.095)	-18.923 (0.090)	-3.21 (0.04)	-0.06 (0.10)	-5.94 (0.04)	29.34 (0.04)
R11	all	N	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.11)	-15.677 (0.095)	-18.923 (0.090)	-3.21 (0.04)	-0.05 (0.10)	-5.94 (0.04)	29.34 (0.04)
R11	all	N	$z>0.0233$	0.162 (0.013)	2.76 (0.13)	-15.678 (0.094)	-18.927 (0.089)	-3.21 (0.04)	-0.06 (0.10)	-5.94 (0.04)	29.34 (0.05)
R11	all	N	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.674 (0.092)	-18.926 (0.087)	-3.21 (0.04)	-0.05 (0.10)	-5.94 (0.04)	29.34 (0.04)
R11	all	N	stricter $C$	0.156 (0.011)	2.97 (0.15)	-15.687 (0.094)	-18.925 (0.089)	-3.21 (0.04)	-0.06 (0.10)	-5.94 (0.04)	29.34 (0.04)
R11	all	N	stricter ${\sigma }_{X_1}$	0.171 (0.011)	3.11 (0.11)	-15.679 (0.095)	-18.923 (0.090)	-3.21 (0.04)	-0.06 (0.10)	-5.94 (0.04)	29.34 (0.04)
2.5	n4258+LMC	Y	default	0.165 (0.010)	3.09 (0.11)	-15.818 (0.097)	-19.059 (0.093)	-3.23 (0.05)	-0.45 (0.12)	-6.11 (0.06)	29.43 (0.05)
2.5	n4258+LMC	Y	higher ${\chi }^2$	0.167 (0.010)	3.14 (0.11)	-15.814 (0.097)	-19.059 (0.092)	-3.23 (0.05)	-0.45 (0.12)	-6.11 (0.06)	29.43 (0.05)
2.5	n4258+LMC	Y	$z>0.0233$	0.162 (0.013)	2.77 (0.13)	-15.816 (0.096)	-19.065 (0.091)	-3.23 (0.05)	-0.45 (0.12)	-6.11 (0.05)	29.43 (0.05)
2.5	n4258+LMC	Y	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.810 (0.096)	-19.061 (0.091)	-3.23 (0.05)	-0.45 (0.12)	-6.11 (0.05)	29.43 (0.05)
2.5	n4258+LMC	Y	stricter $C$	0.156 (0.011)	2.99 (0.15)	-15.822 (0.098)	-19.060 (0.092)	-3.23 (0.05)	-0.45 (0.12)	-6.11 (0.05)	29.43 (0.05)
2.5	n4258+LMC	Y	stricter ${\sigma }_{X_1}$	0.171 (0.011)	3.11 (0.11)	-15.816 (0.098)	-19.060 (0.093)	-3.23 (0.05)	-0.45 (0.12)	-6.11 (0.05)	29.43 (0.05)
2.25	n4258+MW	Y	default	0.165 (0.010)	3.08 (0.11)	-15.650 (0.097)	-18.891 (0.092)	-3.31 (0.06)	-0.52 (0.12)	-5.89 (0.04)	29.32 (0.04)
2.25	n4258+MW	Y	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.10)	-15.646 (0.097)	-18.892 (0.092)	-3.31 (0.06)	-0.52 (0.12)	-5.89 (0.04)	29.32 (0.04)
2.25	n4258+MW	Y	$z>0.0233$	0.162 (0.013)	2.76 (0.13)	-15.645 (0.098)	-18.893 (0.093)	-3.31 (0.06)	-0.52 (0.12)	-5.89 (0.04)	29.32 (0.04)
2.25	n4258+MW	Y	lower ${\chi }^2$	0.159 (0.010)	3.06 (0.12)	-15.642 (0.095)	-18.893 (0.091)	-3.31 (0.06)	-0.52 (0.12)	-5.89 (0.04)	29.32 (0.04)
2.25	n4258+MW	Y	stricter $C$	0.157 (0.011)	2.97 (0.15)	-15.654 (0.098)	-18.893 (0.094)	-3.31 (0.06)	-0.51 (0.12)	-5.89 (0.04)	29.32 (0.04)
2.25	n4258+MW	Y	stricter ${\sigma }_{X_1}$	0.172 (0.010)	3.11 (0.11)	-15.646 (0.097)	-18.890 (0.092)	-3.31 (0.06)	-0.52 (0.12)	-5.89 (0.04)	29.32 (0.04)
2.25	n4258+MW	N	default	0.165 (0.010)	3.08 (0.11)	-15.617 (0.093)	-18.858 (0.089)	-3.16 (0.05)	-0.40 (0.11)	-5.89 (0.04)	29.29 (0.04)
2.25	n4258+MW	N	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.10)	-15.613 (0.093)	-18.859 (0.089)	-3.16 (0.05)	-0.41 (0.12)	-5.89 (0.04)	29.29 (0.04)
2.25	n4258+MW	N	$z>0.0233$	0.162 (0.012)	2.75 (0.13)	-15.612 (0.094)	-18.861 (0.089)	-3.16 (0.05)	-0.41 (0.11)	-5.89 (0.04)	29.29 (0.04)
2.25	n4258+MW	N	lower ${\chi }^2$	0.158 (0.010)	3.05 (0.12)	-15.609 (0.095)	-18.861 (0.090)	-3.16 (0.05)	-0.40 (0.11)	-5.89 (0.04)	29.29 (0.04)
2.25	n4258+MW	N	stricter $C$	0.156 (0.011)	2.96 (0.14)	-15.621 (0.095)	-18.859 (0.090)	-3.16 (0.05)	-0.41 (0.11)	-5.89 (0.04)	29.28 (0.04)
2.25	n4258+MW	N	stricter ${\sigma }_{X_1}$	0.172 (0.010)	3.10 (0.11)	-15.610 (0.096)	-18.854 (0.091)	-3.16 (0.05)	-0.41 (0.11)	-5.89 (0.04)	29.28 (0.04)
R11	LMC+MW	Y	default	0.165 (0.010)	3.09 (0.11)	-15.646 (0.102)	-18.887 (0.098)	-3.25 (0.05)	-0.11 (0.11)	-5.91 (0.06)	29.29 (0.07)
R11	LMC+MW	Y	higher ${\chi }^2$	0.167 (0.010)	3.14 (0.11)	-15.640 (0.101)	-18.886 (0.097)	-3.25 (0.05)	-0.11 (0.11)	-5.91 (0.06)	29.29 (0.07)
R11	LMC+MW	Y	$z>0.0233$	0.162 (0.012)	2.76 (0.13)	-15.641 (0.102)	-18.890 (0.097)	-3.25 (0.05)	-0.11 (0.11)	-5.91 (0.06)	29.29 (0.07)
R11	LMC+MW	Y	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.640 (0.102)	-18.891 (0.098)	-3.25 (0.05)	-0.11 (0.11)	-5.91 (0.06)	29.29 (0.07)
R11	LMC+MW	Y	stricter $C$	0.156 (0.011)	2.97 (0.14)	-15.654 (0.101)	-18.893 (0.097)	-3.25 (0.05)	-0.11 (0.11)	-5.91 (0.06)	29.29 (0.07)
R11	LMC+MW	Y	stricter ${\sigma }_{X_1}$	0.171 (0.011)	3.11 (0.11)	-15.644 (0.102)	-18.888 (0.098)	-3.25 (0.05)	-0.11 (0.11)	-5.91 (0.06)	29.29 (0.07)
2.25	LMC+MW	N	default	0.164 (0.010)	3.09 (0.10)	-15.634 (0.098)	-18.875 (0.094)	-3.16 (0.04)	-0.07 (0.09)	-5.90 (0.05)	29.26 (0.06)
2.25	LMC+MW	N	higher ${\chi }^2$	0.166 (0.010)	3.14 (0.10)	-15.630 (0.099)	-18.874 (0.095)	-3.16 (0.04)	-0.07 (0.10)	-5.90 (0.05)	29.26 (0.06)
2.25	LMC+MW	N	$z>0.0233$	0.161 (0.012)	2.76 (0.13)	-15.631 (0.098)	-18.879 (0.093)	-3.16 (0.04)	-0.07 (0.09)	-5.90 (0.05)	29.26 (0.06)
2.25	LMC+MW	N	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.624 (0.101)	-18.875 (0.097)	-3.16 (0.04)	-0.07 (0.10)	-5.90 (0.06)	29.26 (0.07)
2.25	LMC+MW	N	stricter $C$	0.155 (0.011)	2.98 (0.14)	-15.637 (0.099)	-18.875 (0.095)	-3.16 (0.04)	-0.07 (0.10)	-5.90 (0.05)	29.26 (0.06)
2.25	LMC+MW	N	stricter ${\sigma }_{X_1}$	0.170 (0.010)	3.11 (0.11)	-15.629 (0.100)	-18.873 (0.096)	-3.16 (0.04)	-0.07 (0.10)	-5.90 (0.06)	29.26 (0.07)
R11	LMC+MW	N	default	0.165 (0.010)	3.08 (0.11)	-15.626 (0.102)	-18.867 (0.097)	-3.20 (0.04)	-0.02 (0.10)	-5.89 (0.06)	29.26 (0.07)
R11	LMC+MW	N	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.10)	-15.620 (0.101)	-18.866 (0.097)	-3.20 (0.04)	-0.03 (0.10)	-5.89 (0.06)	29.26 (0.07)
R11	LMC+MW	N	$z>0.0233$	0.162 (0.012)	2.76 (0.13)	-15.623 (0.102)	-18.872 (0.098)	-3.20 (0.04)	-0.03 (0.10)	-5.89 (0.06)	29.26 (0.07)
R11	LMC+MW	N	lower ${\chi }^2$	0.158 (0.010)	3.05 (0.12)	-15.618 (0.100)	-18.869 (0.095)	-3.20 (0.05)	-0.02 (0.10)	-5.89 (0.06)	29.26 (0.07)
R11	LMC+MW	N	stricter $C$	0.156 (0.011)	2.97 (0.15)	-15.630 (0.103)	-18.869 (0.098)	-3.20 (0.05)	-0.03 (0.10)	-5.89 (0.06)	29.26 (0.07)
R11	LMC+MW	N	stricter ${\sigma }_{X_1}$	0.171 (0.010)	3.10 (0.11)	-15.621 (0.103)	-18.865 (0.098)	-3.20 (0.04)	-0.02 (0.10)	-5.89 (0.06)	29.26 (0.07)
2.25	n4258	Y	default	0.165 (0.010)	3.08 (0.11)	-15.747 (0.109)	-18.989 (0.104)	-3.23 (0.07)	-0.55 (0.12)	-6.03 (0.07)	29.40 (0.06)
2.25	n4258	Y	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.11)	-15.740 (0.108)	-18.985 (0.104)	-3.23 (0.07)	-0.55 (0.12)	-6.03 (0.07)	29.40 (0.06)
2.25	n4258	Y	$z>0.0233$	0.162 (0.013)	2.76 (0.13)	-15.744 (0.107)	-18.993 (0.102)	-3.23 (0.07)	-0.55 (0.12)	-6.04 (0.07)	29.41 (0.06)
2.25	n4258	Y	lower ${\chi }^2$	0.158 (0.010)	3.05 (0.12)	-15.737 (0.110)	-18.989 (0.105)	-3.23 (0.07)	-0.55 (0.12)	-6.03 (0.07)	29.40 (0.06)
2.25	n4258	Y	stricter $C$	0.157 (0.011)	2.97 (0.15)	-15.752 (0.109)	-18.990 (0.105)	-3.23 (0.07)	-0.55 (0.12)	-6.04 (0.07)	29.41 (0.06)
2.25	n4258	Y	stricter ${\sigma }_{X_1}$	0.172 (0.011)	3.10 (0.11)	-15.741 (0.108)	-18.985 (0.104)	-3.23 (0.07)	-0.55 (0.12)	-6.03 (0.07)	29.40 (0.06)
2.25	n4258	N	default	0.165 (0.010)	3.08 (0.11)	-15.754 (0.107)	-18.995 (0.104)	-3.04 (0.05)	-0.44 (0.11)	-6.10 (0.07)	29.40 (0.06)
2.25	n4258	N	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.11)	-15.747 (0.107)	-18.992 (0.103)	-3.04 (0.05)	-0.43 (0.11)	-6.10 (0.07)	29.40 (0.06)
2.25	n4258	N	$z>0.0233$	0.162 (0.013)	2.75 (0.13)	-15.750 (0.106)	-18.999 (0.102)	-3.04 (0.05)	-0.44 (0.11)	-6.10 (0.07)	29.40 (0.06)
2.25	n4258	N	lower ${\chi }^2$	0.158 (0.010)	3.05 (0.12)	-15.747 (0.105)	-18.998 (0.102)	-3.04 (0.05)	-0.43 (0.11)	-6.10 (0.07)	29.41 (0.06)
2.25	n4258	N	stricter $C$	0.156 (0.011)	2.96 (0.15)	-15.759 (0.107)	-18.998 (0.103)	-3.04 (0.05)	-0.44 (0.11)	-6.10 (0.07)	29.40 (0.06)
2.25	n4258	N	stricter ${\sigma }_{X_1}$	0.171 (0.010)	3.10 (0.11)	-15.748 (0.106)	-18.992 (0.102)	-3.04 (0.05)	-0.44 (0.11)	-6.10 (0.07)	29.40 (0.06)
2.5	n4258	N	default	0.165 (0.010)	3.08 (0.11)	-15.775 (0.107)	-19.017 (0.103)	-3.06 (0.06)	-0.33 (0.13)	-6.12 (0.07)	29.41 (0.06)
2.5	n4258	N	higher ${\chi }^2$	0.166 (0.010)	3.13 (0.10)	-15.769 (0.108)	-19.015 (0.104)	-3.06 (0.06)	-0.34 (0.13)	-6.11 (0.07)	29.40 (0.06)
2.5	n4258	N	$z>0.0233$	0.162 (0.012)	2.76 (0.13)	-15.768 (0.108)	-19.017 (0.104)	-3.06 (0.06)	-0.34 (0.13)	-6.11 (0.07)	29.40 (0.06)
2.5	n4258	N	lower ${\chi }^2$	0.158 (0.010)	3.05 (0.12)	-15.762 (0.108)	-19.013 (0.104)	-3.06 (0.06)	-0.33 (0.13)	-6.11 (0.07)	29.40 (0.06)
2.5	n4258	N	stricter $C$	0.156 (0.011)	2.97 (0.15)	-15.779 (0.110)	-19.017 (0.106)	-3.06 (0.06)	-0.34 (0.13)	-6.11 (0.07)	29.40 (0.06)
2.5	n4258	N	stricter ${\sigma }_{X_1}$	0.172 (0.010)	3.10 (0.11)	-15.770 (0.107)	-19.014 (0.103)	-3.06 (0.06)	-0.34 (0.13)	-6.11 (0.07)	29.40 (0.06)
R11	n4258	N	default	0.165 (0.010)	3.08 (0.11)	-15.728 (0.109)	-18.969 (0.105)	-3.09 (0.06)	-0.23 (0.13)	-6.08 (0.07)	29.40 (0.06)
R11	n4258	N	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.10)	-15.725 (0.107)	-18.971 (0.103)	-3.09 (0.06)	-0.23 (0.13)	-6.08 (0.07)	29.40 (0.06)
R11	n4258	N	$z>0.0233$	0.162 (0.013)	2.75 (0.13)	-15.723 (0.109)	-18.971 (0.104)	-3.09 (0.06)	-0.24 (0.13)	-6.08 (0.07)	29.40 (0.06)
R11	n4258	N	lower ${\chi }^2$	0.159 (0.010)	3.05 (0.12)	-15.718 (0.110)	-18.970 (0.106)	-3.09 (0.06)	-0.23 (0.13)	-6.07 (0.07)	29.40 (0.06)
R11	n4258	N	stricter $C$	0.156 (0.011)	2.97 (0.14)	-15.733 (0.110)	-18.972 (0.106)	-3.09 (0.06)	-0.23 (0.13)	-6.08 (0.07)	29.40 (0.06)
R11	n4258	N	stricter ${\sigma }_{X_1}$	0.171 (0.011)	3.10 (0.11)	-15.723 (0.109)	-18.967 (0.105)	-3.09 (0.06)	-0.23 (0.13)	-6.08 (0.07)	29.41 (0.06)
2.25	LMC	Y	default	0.165 (0.010)	3.08 (0.11)	-15.880 (0.104)	-19.122 (0.100)	-3.24 (0.05)	-0.55 (0.12)	-6.16 (0.07)	29.54 (0.07)
2.25	LMC	Y	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.11)	-15.873 (0.102)	-19.119 (0.098)	-3.24 (0.05)	-0.55 (0.12)	-6.16 (0.07)	29.54 (0.07)
2.25	LMC	Y	$z>0.0233$	0.162 (0.013)	2.76 (0.13)	-15.878 (0.104)	-19.127 (0.098)	-3.24 (0.05)	-0.56 (0.12)	-6.16 (0.07)	29.54 (0.07)
2.25	LMC	Y	lower ${\chi }^2$	0.159 (0.010)	3.05 (0.12)	-15.869 (0.105)	-19.121 (0.101)	-3.24 (0.05)	-0.55 (0.12)	-6.16 (0.07)	29.54 (0.07)
2.25	LMC	Y	stricter $C$	0.157 (0.011)	2.97 (0.15)	-15.884 (0.103)	-19.123 (0.099)	-3.24 (0.05)	-0.55 (0.12)	-6.16 (0.07)	29.53 (0.07)
2.25	LMC	Y	stricter ${\sigma }_{X_1}$	0.172 (0.011)	3.10 (0.11)	-15.876 (0.103)	-19.119 (0.099)	-3.24 (0.05)	-0.55 (0.12)	-6.16 (0.07)	29.53 (0.07)
R11	MW	Y	default	0.165 (0.010)	3.08 (0.11)	-15.515 (0.111)	-18.757 (0.106)	-3.21 (0.07)	-0.35 (0.14)	-5.83 (0.05)	29.19 (0.07)
R11	MW	Y	higher ${\chi }^2$	0.167 (0.010)	3.13 (0.11)	-15.513 (0.112)	-18.758 (0.108)	-3.21 (0.07)	-0.35 (0.14)	-5.83 (0.05)	29.19 (0.07)
R11	MW	Y	$z>0.0233$	0.162 (0.013)	2.75 (0.14)	-15.511 (0.111)	-18.760 (0.106)	-3.20 (0.07)	-0.36 (0.14)	-5.83 (0.05)	29.19 (0.07)
R11	MW	Y	lower ${\chi }^2$	0.158 (0.010)	3.06 (0.12)	-15.508 (0.111)	-18.760 (0.107)	-3.21 (0.07)	-0.35 (0.14)	-5.83 (0.05)	29.19 (0.07)
R11	MW	Y	stricter $C$	0.156 (0.011)	2.97 (0.15)	-15.521 (0.112)	-18.759 (0.107)	-3.21 (0.07)	-0.35 (0.14)	-5.83 (0.05)	29.19 (0.07)
R11	MW	Y	stricter ${\sigma }_{X_1}$	0.171 (0.011)	3.10 (0.11)	-15.513 (0.110)	-18.757 (0.106)	-3.21 (0.07)	-0.35 (0.14)	-5.83 (0.05)	29.19 (0.07)

Equations

$$v(z)=H_{0}D(z)$$

$$D(z)=\frac{1}{1+z}{D}_L(z),$$

$$D_L(z)\equiv (1+z)R_0{S}_k(\chi ),$$

$$S_k(x)$$

$$v(z)=c{\int }_0^z\frac{d{z}^{\prime }}{E(z^{\prime })},$$

$$E(z)=\sqrt{{\Omega }_M(1+z{)}^3+{\Omega }_k(1+z{)}^2+{\Omega }_{\Lambda }}$$

$$v(z)=\frac{cz}{1+z}\left[1+\frac{1}{2}(1-q_0)z-\frac{1}{6}(1-q_0-3q_0^2+j_0)z^2\right].$$

$$H_0=\frac{v(z)(1+z)}{D_L(z)}\frac{S_k(\chi )}{\chi }.$$

$$H_0$$

$$=\frac{cz}{D_L(z)}\left[1+\frac{1}{2}(1-q_0)z-\frac{1}{6}(1-q_0-3q_0^2+j_0)z^2\right].$$

$$\mathcal{L}(\Theta )$$

$${\chi }^2(\Theta )$$

$$m_W=b_W({\log }_{10}P-1)+Z_W\Delta {\log }_{10}[O/H{]}_{ij}+M_W+\mu .$$

$$m_B=M_B-\alpha X_1+\beta C+\mu$$

$$m_B=M_B^{\ast }-\alpha X_1+\beta C+5{\log }_{10}\left(\frac{czf(z)}{H_0}\right)+25.$$

$$m_{Wij}$$

$$m_{Bi}$$

$$m_{Bk}$$

$$\mathcal{H}:=5{\log }_{10}{H}_0-25.$$

$${\chi }_c^2=\sum _{ij}\frac{({\hat{m}}_{Wij}-m_{Wij,\mathrm{mod}}{)}^2}{{\hat{m}}_{Wij,\mathrm{err}}^2+{\sigma }_{\mathrm{int},\mathrm{C}}^2}.$$

$$m_B^{†}$$

$$\text{with }{{\hat{m}}^{†}}_{B\mathrm{mod}}$$

$${\chi }_{\mathrm{SN}}^2$$

$${\chi }_{\mathrm{global}}^2$$

$${\chi }_{\mathrm{nearby}}^2$$

$$\alpha$$

$$\beta$$

$$b_W$$

$$Z_W$$

$$M_W$$

$$M_B$$

$$\mathcal{H}$$

$$M_B$$

$$\mathcal{H}$$

$$\frac{{\sigma }_{H_0}}{H_0}=\frac{\log (10)}{5}{\sigma }_{\mathcal{H}}.$$

Full text

††pubyear: 2017††pagerange: A blinded determination of $H_{0}$ from low-redshift Type Ia supernovae, calibrated by Cepheid variables–13

1 Introduction

The Hubble constant $H_{0}$ has proven difficult to measure since the discovery of the Universe’s expansion almost a century ago (Hubble, 1929), following the prediction of the latter in Friedmann’s equations (Friedmann, 1922). As given in the Hubble law $v=H_{0}D$ (first derived by Lemaître, 1927), $H_{0}$ sets the cosmic distance scale via the present expansion rate of the local Universe. The quest to make precise measurements of $H_{0}$ has been a continual challenge in observational cosmology, due to the difficulty of making accurate distance measurements.

Recently, discrepant values obtained from local and global measurements have propelled the Hubble constant back into the spotlight. Observations of cosmic microwave background (CMB) anisotropies with the Planck satellite found $H_{0}=67.3\pm 1.2$ km s*-1* Mpc*-1* (Planck Collaboration et al., 2014), assuming a standard $\Lambda$CDM cosmology. This value is ${\sim}2.7\sigma$ lower than in Riess et al. (2011, hereafter R11), who measure $H_{0}=73.8\pm 2.4$ km s*-1* Mpc*-1* from observations of type Ia supernovae (SNe Ia) in the more local Universe. While the Planck measurement is dependent on an underlying cosmological model, the SN Ia-based measurement is model-independent. The precision of these values highlights the importance of the tension between the two modes of measurements, which has increased to over $3\sigma$ significance in the updated analyses in Riess et al. (2016, hereafter R16) (finding $H_{0}=73.0\pm 1.8$ km s*-1* Mpc*-1*), and Planck Collaboration et al. (2016) (finding $H_{0}=67.8\pm 0.9$ km s*-1* Mpc*-1*).

Numerous reanalyses of the SN Ia-based measurement have followed, many of which have focussed on the methods for the rejection of Cepheid outliers. Efstathiou (2014, hereafter E14) questions and revises the outlier rejection algorithm in R11, concluding $H_{0}=72.5\pm 2.5$ km s*-1* Mpc*-1* assuming a null metallicity dependence of the Leavitt law. Recently, Cardona et al. (2017) uses Bayesian hyper-parameters to down-weight portions of the Cepheid data for both R11 and R16 data sets, finding $H_{0}=73.75\pm 2.11$ km s*-1* Mpc*-1* for the R16 data. Moreover, the dependence of the intrinsic magnitude of SNe Ia on host galaxy properties has been explored in recent years (e.g. Sullivan et al., 2010). Rigault et al. (2013, 2015) find a relationship between peak brightness and star formation rate, and infer an overestimate of ${\sim}3$ km s*-1* Mpc*-1* in the R11 value of $H_{0}$ arising from the fact that the calibration set of SNe Ia exist in galaxies which necessarily contain Cepheids, hence are likely to be late-type galaxies. However, Jones et al. (2015) repeat the same analysis, with an increased sample size and the R11 selection criteria applied, and find no significant difference in the brightness of SNe Ia in star-forming and passive environments.

The CMB data in Planck has been reanalysed in Spergel et al. (2015), who find a similar value to Planck Collaboration et al. (2014), of $H_{0}=68.0\pm 1.1$ km s*-1* Mpc*-1*. Bennett et al. (2014) provides a CMB-based measurement which is independent of Planck, by combining data from WMAP9, the South Pole Telescope (SPT) and Atacama Cosmology Telescope (ACT), and baryon acoustic oscillation (BAO) measurements from BOSS, finding a value of $H_{0}=69.3\pm 0.7$ km s*-1* Mpc*-1* (with a slight increase to $H_{0}=69.7\pm 0.7$ km s*-1* Mpc*-1* if SN Ia data from R11 are included), which is slightly less discrepant with SN Ia-based values. Strong lensing provides an independent but model-dependent local measurement of $H_{0}$: the Suyu et al. (2017, (H0LiCOW)) program studies time delays between multiple images of quasars in strong gravitational lens systems, and find $H_{0}=71.9^{+2.4}_{-3.0}$ km s*-1* Mpc*-1* (Bonvin et al., 2017) in flat $\Lambda$CDM. It is noteworthy that the H0LiCOW analysis was performed blind to derived cosmological parameters Bonvin et al. (e.g.  2017, section 4.4); we discuss the importance of blinding in our analysis in Section 2.4.

One of the biggest open questions in cosmology today is whether the tension in $H_{0}$ signifies new physics – where inconsistencies between results from supernovae and the CMB arise from the model-dependence of the measurement, and disappear when the correct model is used – or is the result of some systematic error in one or both measurements that has yet to be accounted for. Possible theoretical modifications to standard $\Lambda$CDM to reconcile the tension in $H_{0}$ include an increased neutrino effective number (the existence of dark radiation), and/or a more negative dark energy equation-of-state parameter $w$ at late times. Di Valentino et al. (2016) explore these scenarios in a higher-dimensional parameter space, with their findings supporting phantom dark energy with $w{\sim}-1.3$, while Wyman et al. (2014); Dvorkin et al. (2014); Leistedt et al. (2014) focus on the implications of an additional massive sterile neutrino species. Meanwhile, Bernal et al. (2016) examine the model dependence of the Universe’s distance scale (anchored by $H_{0}$ and by the scale $r_{S}$ of the sound horizon at radiation drag, at late and early times respectively) by reconstructing its expansion history with minimal cosmological assumptions.111This is possible as the combination of SNe Ia and BAO as probes constrains the product $r_{S}H_{0}$ in a model-independent way. They conclude that the tension in $H_{0}$ translates to a mismatch in the normalisations provided by $H_{0}$ and $r_{s}$ at two opposite ends of the distance ladder.

A genuine inconsistency in the value of the Hubble constant at low and high redshifts would have profound consequences. Therefore it is imperative to fully understand uncertainties in the measured values of $H_{0}$, and to preclude possible human biases on the result. The most effective way of achieving the latter is to blind the value of $H_{0}$ throughout the analysis.

The use of data from R11 is for a proof of concept, necessary for our blind analysis technique, and to be followed shortly with the same analysis applied to R16 data. Numerous improvements over R11 have been made in R16, in the analysis as well as the size and quality of data. Changes to the outlier rejection and the Cepheid metallicity calculations have addressed some of the concerns raised in E14. However, our analysis involves both a simultaneous fit to all data sets, and the accepted methodology of recent supernova cosmology analyses (Conley et al., 2011; Betoule et al., 2014) for considering SN Ia systematics. Both of these points carry significant differences from the R11 and R16 analysis chains, and have yet to be included in a reanalysis. Nor has the supernova data been revisited in its entirety, starting from the available photometry. Thus, we are motivated by the desire to provide such a validation of the supernova data, and by the current relevance and importance of the Hubble constant, to produce in this work an independent, blinded, end-to-end reanalysis of the R11 data to determine $H_{0}$ and its uncertainty.

In summary, we combine the framework for calibrating a SN Ia Hubble diagram with Cepheid variables, with the best estimates of supernova systematics via covariance matrices. We determine $H_{0}$ using the magnitude-redshift relation (i.e. a Hubble diagram) of low-redshift SNe Ia, with their zero point set by Cepheid variables in host galaxies of eight nearby SNe Ia, which are in turn calibrated by very long baseline interferometry (VLBI) observations of megamasers in NGC 4258, and other geometric distances to the Large Magellanic Cloud (LMC) and Cepheids in our Galaxy.

This paper is structured as follows. First we present an overview of our methods in §2, followed by the distinct sets of data with the equations relating them in §3. In §4 we focus on the Cepheid variables and perform a fit to the Cepheid data only, comparable to the E14 reanalysis of R11. Next in §5 we discuss type Ia supernovae, including details of fitting SNe Ia on a Hubble diagram and results of a preliminary SN-only fit. This fit relies on computations of individual supernova systematic terms in the form of covariance matrices, which are examined in depth in Appendix C. Finally, §6 ties together the Cepheid and SN Ia information into a combined and simultaneous fit of all data; we conclude with a discussion of these results in §7.

2 Methods

In this section we paint a broad picture of our approach to measuring $H_{0}$, postponing specific details of and equations relating to the data to Section 3. We begin with the theory and mathematics of finding $H_{0}$ in the cosmology analysis, followed by the astronomy that enables this: distance measurements with Type Ia supernovae and Cepheid variables as standard candles. Equally important are the Bayesian statistics that underpin the analysis, and the method for blinding the result.

2.1 Theory of extracting $H_{0}$

In its traditional formulation Hubble’s law states that the recession velocity of objects is proportional to their distance:

[TABLE]

where the constant of proportionality $H_{0}$ represents the present expansion rate of the Universe, scaled by its size (i.e. $H_{0}=\frac{\dot{a}}{a}$ where $a$ is the scale factor and overdot indicates differentiation with respect to time, $t$). Methods of determining $H_{0}$ typically involve taking the ratio of the two sides of Equation 1. We expand on the subtleties of this below.

The distance in Hubble’s Law is related to the luminosity distance by

[TABLE]

assuming a flat Universe.222To include curvature note that the present distance to an object at redshift $z$ is given by $D(z)=R_{0}\chi$, with $\chi$ being the comoving coordinate and $R_{0}$ the scale factor at the present day with dimensions of distance (in the equation for $H_{0}$ above $a(t)\equiv R(t)/R_{0}$). Then luminosity distance is defined as

$\displaystyle D_{L}(z)\equiv(1+z)R_{0}S_{k}(\chi),$

(3)

with $R_{0}\equiv c/(H_{0}\sqrt{|\Omega_{k}|})$ and

$\displaystyle S_{k}(x)$ $\displaystyle=\begin{cases}\sin x&k=1,\\ x&k=0,\\ \sinh x&k=-1,\end{cases}$

(4)

so $D(z)=\frac{1}{1+z}\frac{\chi}{S_{k}(\chi)}D_{L}(z)$.

The luminosity distance $D_{L}(z)$ can be determined observationally (i.e. with no knowledge of cosmological parameters) using standard candles. These have known absolute magnitudes $M$, so taking the difference between $M$ and the apparent magnitude $m$ gives the distance modulus $\mu\equiv m-M$ and hence the luminosity distance $D_{L}\equiv 10^{\frac{\mu-25}{5}}$ Mpc. In practice the process of measuring distances is far from straightforward, and is outlined in Section 2.2.

On the left hand side, $v(z)$ is the predicted velocity due to expansion for a galaxy at redshift $z$.333For simplicity, we do not distinguish here between redshifts in different reference frames, and only use one redshift $z$. We distinguish between different redshifts, particularly in Equation 2, in Appendix B.4. The exact expression for $v(z)$ is given by integrating the Universe’s expansion up to redshift $z$:

[TABLE]

where $E(z)\equiv H(z)/H_{0}$ is a function of cosmological parameters, as defined in Peebles (1993),444In Friedmann-Lemaître-Robertson-Walker (FLRW) cosmologies, $E(z)$ is given by (e.g. Carroll et al., 1992)

$\displaystyle E(z)=\sqrt{\Omega_{M}(1+z)^{3}+\Omega_{k}(1+z)^{2}+\Omega_{\Lambda}}$

(6)

where $\Omega_{M}$ and $\Omega_{\Lambda}$ are respectively the densities of normal matter and the cosmological constant (relative to the critical density), $k$ is the curvature, and $\Omega_{k}\equiv 1-\Omega_{M}-\Omega_{\Lambda}$ (zero in a flat Universe). If dark energy is something other than a cosmological constant, with a generic equation of state $w$, replace $\Omega_{\Lambda}$ with $\Omega_{\rm de}(1+z)^{3(1+w)}$. and $v(z)$ is independent of $H_{0}$.555It is interesting to note that $v(z)$ is independent of $H_{0}$; it depends only on redshift and cosmological parameters such as $\Omega_{M}$ and $\Omega_{\Lambda}$. That may seem unintuitive, but it is velocity as a function of distance $v(D)$ that is function of $H_{0}$ (things that are moving faster have gone further). Velocity as a function of redshift $v(z)$ works differently since redshift is not proportional to distance. A galaxy’s redshift is determined by how much the Universe has expanded since the light was emitted. That depends on the travel time, which does depend on the densities that cause the Universe to accelerate or decelerate (and thus for the light to take longer or shorter times to reach us), but not on $H_{0}$.

At low redshifts the cosmological dependence of $v(z)$ is very weak and it is a good approximation to use a second order Taylor expansion in terms of the deceleration and jerk parameters $q_{0}$ and $j_{0}$666We assume a standard $\Lambda$CDM cosmology with $\Omega_{M}{\sim}0.3,\Omega_{\Lambda}{\sim}0.7$, fixing $q_{0}=-0.55$ and $j_{0}=1$.. Thus we follow R11 and use,

[TABLE]

At low redshift Equations 5 and 7 both reduce to the familiar $v(z)\approx cz$. At moderate redshifts ($z<0.1$), Equation 7 closely approximates most observationally reasonable cosmological models. We explored the uncertainty associated with assuming Equation 7 and the cosmology stated in Footnote 6, finding the impact to be small: varying either $\Omega_{M}$ or $w$ by 0.1 changes $H_{0}$ by 0.015 km s*-1* Mpc*-1* or 0.1 km s*-1* Mpc*-1*, respectively, in the sense that an increase in $\Omega_{M}$ or $w$ causes an increase in $H_{0}$. The maximal difference in $\mathcal{M}$ induced by varying $q_{0},j_{0}$ within values allowed by 1$\sigma$ contours in Betoule et al. (2014) is an order of magnitude smaller than its statistical uncertainty.

Rearranging Equations 1, 2, and 7 gives us the equation for $H_{0}$ as a function of observables, $z$ and $D_{L}$,777For non-zero curvature, Equation 9 becomes

$\displaystyle H_{0}=\frac{v(z)(1+z)}{D_{L}(z)}\frac{S_{k}(\chi)}{\chi}.$

(8)

[TABLE]

Thus determining $H_{0}$ amounts to comparing the velocity in Equation 7 – derived from the measured redshift – to the observed luminosity distance, measured with standard candles. The equations encapsulating this process are detailed in Section 3.

2.2 Measuring distance

Astronomical distances can be measured using standard candles: standardisable objects with known absolute magnitude which, combined with an apparent magnitude, give the distance modulus. These distances are often relative rather than absolute. Since each mode of measurement is useful only over a limited range of distances, multiple standard candles are tied together to form a so-called distance ladder. At the bottom of the ladder are absolute distances determined from geometric methods (i.e. trigonometric parallax), only accurate at relatively small distances. Then nearby standard candles (i.e. Cepheid variables) give distances relative to this geometric scale; similarly, each rung of the ladder is calibrated on the previous.

Standard candles (a distance scale) provide one approach to measuring cosmological parameters including $H_{0}$. Alternatively, standard rods (a length scale) in the form of baryon acoustic oscillations (BAO; e.g. Eisenstein et al., 2005) provide complementary (and for some parameters, orthogonal) constraints, most recently in Planck Collaboration et al. (2016). Weinberg et al. (2013) provides a review of cosmological probes; we refer the interested reader to its section 4 for a review of BAO.

In our determination of $H_{0}$ we rely on two standard candles: type Ia supernovae and Cepheid variables. These together prescribe a relative distance scale for the low-$z$ SNe Ia. The absolute calibration is given by the geometric maser distance of NGC 4258 from Humphreys et al. (2013). The Cepheid variables lie in this galaxy and eight other galaxies containing nearby SNe Ia, calibrating the supernovae. The absolute distances and measured redshifts of the low-$z$ supernovae are combined to determine $H_{0}$ as described in Section 2.1, through equations detailed in Section 3. Next we briefly describe each standard candle.

Cepheid variables are pulsating supergiants with periods of days to hundreds of days, well-characterised by their luminosity via the empirical Leavitt law (Leavitt, 1908; Leavitt & Pickering, 1912) – also commonly known as the Period-Luminosity relation. The brightness and regular pulsation of Cepheid variables as well as their ease of discovery and classification make Cepheids reliable distance indicators in the nearby Universe, and the basis of the cosmic distance ladder (Freedman & Madore, 2010). Some difficulties and systematics include crowding and confusion (which necessitate outlier rejection), metallicity, and extinction; these are discussed further in Section 4.2.

SNe Ia are thought to be thermonuclear explosions of accreting white dwarfs, with two qualities which recommend them as excellent distance indicators up to high redshift: they are intrinsically very bright, and highly standardisable in terms of their apparent peak brightness, lightcurve shape, and colour (Phillips, 1993; Hamuy et al., 1996; Riess et al., 1996) – see Section 3.2.2 for more details. Indeed SNe Ia have played a pivotal role in recent observational cosmology, particularly in the discovery of the accelerating Universe (Riess et al., 1998; Schmidt et al., 1998; Perlmutter et al., 1999). In the past decade SN Ia samples have greatly expanded, reducing statistical uncertainty. However, observations of SN Ia are still subject to numerous systematics, which can be significant and correlated between SNe. These systematics include: calibration uncertainties, dust, and corrections for peculiar velocities, and host galaxy mass, and will be discussed in Appendix C.

2.3 Bayesian statistical methods

We estimate the fit parameters $\Theta$ (given in Section 3) in a Bayesian framework, relying on the principle of sampling the likelihood $\mathcal{L}(\Theta)$ over the parameter space to determine the posterior distribution function (PDF). The generalised likelihood is determined from the $\chi^{2}$ statistic, a function of $\Theta$:

[TABLE]

Here $\boldsymbol{\hat{m}}$ and $\boldsymbol{m_{\rm{mod}}}$ are the observed and theoretical magnitude vectors (over all data) respectively,888We retain this convention where it is necessary to explicitly distinguish the data from the model. and $\mathbf{C}$ is the covariance matrix of uncertainties in $\boldsymbol{\hat{m}}$. The model $\boldsymbol{m_{\rm{mod}}}$ is implicitly a function of $\Theta$. In each fit outlined in Section 3.4 an expression for $\chi^{2}$ will be given explicitly, i.e. Equations 20, 22, and 23. When uncertainties are uncorrelated (i.e. $C$ is diagonal) Equation 11 reduces to the more familiar

[TABLE]

2.3.1 PDF estimation

In higher dimensional parameter spaces the computational expense of calculating and integrating the likelihood necessitates Monte Carlo techniques to statistically sample the parameter space, the most common being Markov Chain Monte Carlo (MCMC). These techniques are useful for parameter estimation or model selection (see e.g. Davis & Parkinson, 2016). Nested sampling (Skilling, 2004) is another such technique, in which the likelihood is evaluated at sample ‘live’ points drawn from an iteratively replaced distribution until convergence, where the posterior is recovered. The MultiNest algorithm (Feroz & Hobson, 2008; Feroz et al., 2009; Feroz et al., 2013) is a robust nested sampling tool for retrieving posterior samples from distributions which may have multiple peaks or ‘modes’. We use the implementation PyMultiNest described in Buchner et al. (2014) for most fits (details in Section 3.4). For some lower-dimensional fits (Section 5.4) we use emcee (Foreman-Mackey et al., 2013), a Python implementation of MCMC.

Each algorithm takes as inputs the data, a prior distribution within the parameter space (which live points are drawn from), and the likelihood as a function of the data and parameters. In addition we select the sampling efficiency for parameter estimation and the number of live points (walkers). MultiNest outputs include the best fit (maximal likelihood) parameters and the marginalised posterior distribution for each parameter. In our fits the marginalised PDFs appear symmetrical and Gaussian (e.g. Figs 3 and 4), so we take our best estimates of values and uncertainties of each parameter from the mean and standard deviation of the marginalised PDF.

2.4 Blind analysis

To perform a blind analysis is to obscure the principal aspects of the result until the analysis is complete. The overarching motivation for blinding is to eliminate the impact of human biases on the result, including confirmation bias. Preconceptions about the ‘correct’ value for a result are irrelevant to the validity of the analysis and can only reduce the value of the findings. Conversely a blind analysis has all the more bearing for having reached its conclusion blind. Croft & Dailey (2011) find evidence of confirmation bias in recent measurements of cosmological parameters and recommend blinding; similarly Maccoun & Perlmutter (2015) argue for its necessity. In recent years the practice of blind analysis has become standard in particle physics, and is increasingly adopted in cosmology.

Our priority is to hide the value of $H_{0}$ so as to not influence its result, so we blind the parameter $\mathcal{H}$ which contains equivalent information.999We fit for the parameter $\mathcal{H}:=5\log_{10}H_{0}-25$, which is linear in magnitude (Equation 19), instead of $H_{0}$. We also blind the SN Ia magnitude zero point $M_{B}$ which has the most interaction with $\mathcal{H}$, and is the best-constrained in the literature, relative to other parameters in $\Theta$. We implement these blinds in the analysis and data respectively. For any likelihood function containing $\mathcal{H}$ (i.e. involving the low-$z$ SNe) we make the shift $\mathcal{H}\mapsto\mathcal{H}+\mathfrak{o}_{H}$ for an offset $\mathfrak{o}_{H}$. Meanwhile we effectively shift $M_{B}$ by adding another offset $\mathfrak{o}_{M}$ to all SN magnitudes $m_{B}$. Both offsets $\mathfrak{o}_{H}$ and $\mathfrak{o}_{M}$ are unknown real numbers, randomly drawn from normal distributions and never printed. These are seeded by distinct known numbers to ensure that the offsets are constant and can be retrieved. Our method allows the recovery of the true unblinded values by simply subtracting the offsets once the blind is lifted.

We choose to not blind the other parameters which appear in the preliminary Cepheid- or SN-only fits, primarily because these parameters do not have strong enough priors from the literature to introduce human bias. Moreover, the variation we observe in the preliminary values of the nuisance parameters $\{b_{W},Z_{W},\alpha,\beta\}$ is useful for informing which preliminary fits to carry forward to the global fits. Knowing the preliminary nuisance parameters will not bias our results because ‘best’ versions of the preliminary fits are not chosen; instead we select a representative sample of these fits and use the scatter to quantify the systematic uncertainties.

3 Data and Analysis Techniques

This section describes our Cepheid and SN Ia data, using equations for the apparent magnitude of each data set to demonstrate the relationships between them. These are followed by an outline of the steps of the fit.

3.1 Data samples

Our analysis uses three sets of data:

Cepheid variables: 570 spread between nine nearby galaxies (see Table 1), namely:

• •

165 in the distance anchor NGC 4258, and

• •

405 in eight galaxies that host recent nearby SNe Ia. 2. 2.

Anchor (‘nearby’) supernovae: 8 recent SNe Ia in the nearby galaxies (also in Table 1). 3. 3.

Low-$z$ SNe Ia: 280 low-redshift ($z<0.06$) SNe Ia from the CfA3 (Hicken et al., 2009a) and LOSS (Ganeshalingam et al., 2010) samples.

Together these three data sets allow us to calibrate our distance ladder. The galaxy NGC 4258 hosts the water masers that give us a precise absolute local distant measurement (Humphreys et al., 2013), and allows us to calibrate the Cepheids. As in R11, we also use the LMC and Milky Way (MW) as distance anchors in combination with NGC 4258, relying on independent distances measured from detached eclipsing binaries (Pietrzyński et al., 2013) to Cepheids in the LMC, and Hipparcos and HST parallax measurements of Cepheids in our Galaxy (van Leeuwen et al., 2007). The Cepheids in turn enable us to calibrate the absolute magnitudes of the eight supernovae that occurred in nearby galaxies with quality Cepheid measurements. These then allow us to calibrate the whole supernova sample, which ultimately gives most of the constraining power for our $H_{0}$ measurement. In practice we perform a global fit to all of these samples together. In the next section we outline the equations needed to relate all of these standard candles and extract a measurement of $H_{0}$ following the theory in Section 2.1.

Since the purpose of this paper is to provide an independent analysis of the data in R11, we adopt an identical sample in order to make a faithful comparison. Our aim is to use the same framework to analyse newer data sets including SNe Ia in the CfA4 survey (Hicken et al., 2012) and Cepheids in R16 at a later stage.

3.2 Equations for apparent magnitude

3.2.1 Cepheids

Our first data set, the Cepheid variables, allow us to infer distances to the nearby galaxies via the Leavitt law (also commonly known as the period-luminosity relation):

[TABLE]

Equation 13 relates the apparent ‘extinction-free’ (Wesenheit) magnitude $m_{W}$,101010We use the quantity $M_{W}\equiv V-R_{V}(V-I)$ constructed in Madore (1982) from the Wesenheit function (van den Bergh, 1975), from V- and I-band absolute magnitudes. Assuming constant ratio $R_{V}$ of total to selective absorption, $M_{W}$ is independent of extinction. We fix $R_{V}=A_{V}/E(V-I)=3.1$ as in R11. period ($P$; in days), and metallicity of a Cepheid at distance modulus $\mu$. The slopes $b_{W}$ and $Z_{W}$ represent the dependence of the magnitude on period and metallicity; the zero point $M_{W}$ physically represents the Wesenheit magnitude of a Cepheid in our Galaxy (at a distance of 10 pc), with a period of 10 days. We use relative values of the metallicity ($\Delta\log_{10}[O/H]_{ij}:=\log_{10}[O/H]-8.9$) and period to pivot the fit near the data.

3.2.2 Type Ia supernovae

Type Ia supernovae comprise our remaining data. A spectroscopically normal SN Ia has a lightcurve parametrized by its brightness (hence distance), observed colour and decline rate. These measures are represented by different quantities in various SN Ia frameworks; in SALT2 (Guy et al., 2007) these are the apparent magnitude $m_{B}$ at time of B-band maximum, ‘stretch’ $X_{1}$ and colour $C$ (roughly corresponding to $B-V$ at maximum), related by:

[TABLE]

where $M_{B}$ is the canonical SN Ia absolute magnitude, and $\alpha,\beta$ are SALT2 nuisance parameters for the stretch and colour dependences.

SNe Ia in more massive galaxies are brighter after these standard corrections for colour and stretch, as discussed in Appendix B.3. To account for this we replace $M_{B}$ in Equation 14 with the corrected absolute magnitude $M_{B}^{*}$, which can take two discrete values depending on the host galaxy mass: $M_{B}$ or $M_{B}+\Delta M_{B}$. We will fix $\Delta M_{B}$ (see Appendix B.3) and fit for the three global parameters $\{\alpha,\beta,M_{B}\}$.

Our second data set contains the eight ‘nearby’ SNe Ia in Table 1, with apparent magnitudes given by Equation 14 (with $M_{B}^{*}$ instead of $M_{B}$). A SN Ia and Cepheid in the same galaxy have common distance modulus $\mu$ in Equations 14 and Equation 13; thus, the Cepheids calibrate the nearby SNe, which in turn determine the SN Ia magnitude zero point $M_{B}$.

The much larger sample of 280 SNe Ia makes up our third data set. These ‘low-$z$’ supernovae originate from CfA3 and LOSS, with details to follow in Section 5.1. Once we have calibrated their absolute magnitudes using the eight ‘nearby’ supernovae, we can use the theory derived in Section 2.1 to relate their measured magnitudes to the value of $H_{0}$. Assuming Equation 9 and writing $f(z)\equiv 1+\frac{(1-q_{0})z}{2}-\frac{(1-q_{0}-3q_{0}^{2}+j_{0})z^{2}}{6}$, we have in place of Equation 14

[TABLE]

3.3 Global fit

We will fit Equations 13, 14, 15 simultaneously for a combined fit to all Cepheid and SN Ia data. We rewrite these equations, making explicit the indexing: $i$ varies over the eight nearby galaxies (and the SNe Ia they contain), $j$ varies over Cepheids in these galaxies and NGC 4258, $k$ varies over the low-$z$ SNe.

[TABLE]

In Equation 18 we separate the intercept of Equation 15 into parameters $M_{B}^{*}$ (also appearing in Equation 17) and a constant term $\mathcal{H}$, which contains the same information as $H_{0}$:

[TABLE]

We fit for all 16 parameters appearing in Equations 16–18; explicitly these are $\Theta=\{\alpha,\beta,\mathcal{H},M_{B},b_{W},Z_{W},M_{W},\mu_{4258},\Delta\mu_{i}\}$ where $i$ varies over the eight nearby galaxies. Note that we fit for $M_{B}$ instead of $M_{B}^{*}$ as the latter is not a constant. The distance moduli in Equations 16 and 17 are expressed as offsets $\Delta\mu_{i}\equiv\mu_{i}-\mu_{4258}$, relative to NGC 4258.

Equations 16–18 assume a distance anchor of NGC 4258. The use of the LMC and MW as alternate or additional anchors is explored, and discussed in Appendix A.2. We impose a strong Gaussian prior $\mu_{4258}=29.404\pm 0.066$ on the distance, measured from VLBI observations of megamasers in Humphreys et al. (2013)111111This distance is slightly higher than the older value $\mu_{4258}=29.31$ assumed in Riess et al. (2012); this increase acts to decrease $H_{0}$ relatively. whenever NGC 4258 is used as an anchor, and similarly $\mu_{\rm{LMC}}=18.494\pm 0.049$ if the LMC is included.

3.4 Steps in fitting process

We break the process of fitting all data to Equations 16–18 into three steps to streamline the process: the data and parameters are separated into spheres of influence so that results from the Cepheid- and SN-only fits – in particular their dependences on factors such as rejection, cuts, and distance anchors – can be isolated, inspected, and selectively carried forward to the global fit.

The three steps are as follows. First we fit all Cepheid data simultaneously for parameters $\left\{b_{W},Z_{W},M_{W},\{\mu_{i}\}\right\}$ to Equation 16. Separately, we fit only the low-$z$ SNe Ia to Equation 18. The parameters $M_{B}$ and $\mathcal{H}$ are degenerate when constrained by only the low-$z$ data, so we fit for their difference $\mathcal{M}:=M_{B}-\mathcal{H}$, as well as SN Ia parameters $\alpha,\beta$. Finally, a global fit is performed (independent of the first two steps) of both data sets and the nearby SNe Ia to Equations 16, 17, 18 simultaneously. This step is similar to the Cepheid-only fit but also includes $\mathcal{H}$ and the SN parameters $\{\alpha,\beta,M_{B}\}$. Final values for all parameters including $H_{0}$ are extracted from this global fit. Each preliminary fit is described in detail in Sections 4.3 and 5.4, and the global fit in Section 6.1.

The Bayesian methods for parameter estimation (MultiNest for the high-dimensional Cepheid-only and global fits, and MCMC for the SN-only fit – described in Section 2.3) require priors, which may be uniform, on each parameter in $\Theta$. While some other parameters are predominantly influenced by a subset of the data (namely the nuisance parameters $b_{W}$ and $Z_{W}$ which only appear in the Cepheid-only fit, and $\alpha$ and $\beta$ which are predominantly determined by low-$z$ SNe), it would be statistically invalid to place Gaussian priors on these parameters in the global fit based on results of either preliminary fit. However non-uniform priors based on external data are allowed; our priors on $\mu_{4258}$ (and $\mu_{\rm{LMC}}$) are Gaussian if these galaxies are included as calibrators, and we constrain $b_{W}$ with a Gaussian prior informed by the LMC Cepheids in fits which are not anchored on the LMC (discussed in Appendix A.4). For the remaining parameters in $\Theta$ we set uniform priors over generous intervals.

Our approach differs from the R11 and E14 analyses, which both perform two independent steps: (i) using only the low-$z$ SN data, determine and fix $a_{V}$ (the intercept of the SN Ia $m-z$ relation equivalent to $0.2\mathcal{M}$ in our analysis) and (ii) from the Cepheids only, determine the Leavitt law parameters $b_{W},Z_{W}$, and $zp_{4258}$ (a zero point comparable to our $M_{W}$). The Cepheid parameters are combined with the nearby SNe Ia lightcurves to find the quantity $m^{0}_{v,4258}$ which signifies the fiducial peak apparent magnitude of a SN Ia in NGC 4258; this quantity is then combined with $a_{V}$ and $\mu_{4258}$ to give $H_{0}$ (R11, equation 4). We emphasize that, in contrast, our final global fit is truly simultaneous in that it allows each parameter in $\Theta$ to be influenced by all Cepheid and supernova data in the nearby galaxies and low-$z$ sample. Consequently we allow the data sets to interact freely with each other, enabling us to capture covariances between parameters.

4 Cepheid Leavitt law fit

Here we describe an initial simultaneous fit of the Cepheids in all nine galaxies to the Leavitt law (Equation 16). This has two purposes: to estimate the parameters $\left\{b_{W},Z_{W},M_{W},\{\mu_{i}\}\right\}$ for each fit (which uniquely define a Cepheid data set), and to examine the dependence of these parameters (particularly the period and metallicity coefficients $b_{W}$ and $Z_{W}$) on factors explored in R11 and E14 – namely the rejection algorithm and threshold, distance anchor, and inclusion of longer-period Cepheids, discussed in Appendices A.1-A.3. Some of these fits, with associated Cepheid data sets, are selected to be carried forward to the global fit.

We emphasize that the process of choosing these fits is motivated by the desire to capture and quantify variation that arises in results when different (but also valid) choices are made in the fitting process, rather than by the aim of choosing a ‘best’ fit; this will become clear in Fig. 1 and its discussion. Thus we do not blind this part of the analysis (the Cepheid-only fit), because the results are not final, and also because they do not directly reveal or affect the value of $H_{0}$.

4.1 Observations

The Cepheids in the nine galaxies in Table 1 were discovered or reobserved in the Supernovae and $H_{0}$ for the Equation of State (SH0ES) project (Riess et al., 2009b) on the HST, from Cycle 15. Infrared (F160W) observations of the SN Ia host galaxies were made using the Wide Field Camera 3 (WFC3). We refer the reader to R11, section 2 for descriptions of observations and data reduction. Our initial data set consists of 570 Cepheids from R11, table 2, excluding those marked ‘low P’; this number is reduced to 488 if we adopt the $P<60$ day cut on Cepheids, following E14.

We supplement the sample of 157 Cepheids in NGC 4258 with LMC and MW Cepheids, used as alternative anchors (discussed in Appendix A.2). Persson et al. (2004) presents near-infrared photometry of 92 Cepheids, of which 53 have optical measurements in Sebo et al. (2002), which we use for determining Wesenheit magnitudes. Two of these 53 Cepheids have period greater than 60 days, which we exclude if we impose the period cut on the Cepheids in the supernova host galaxies. We also make use of 13 Cepheids in the Milky Way from van Leeuwen et al. (2007, table 2) (excluding Polaris, an overtone pulsator), which have combined parallaxes from Hipparcos and HST data.

4.2 Cepheid systematics

Cepheid variables are powerful distance indicators to nearby galaxies, however they are subject to systematics. We briefly mention those that affect our method, and refer the reader to Freedman & Madore (2010, section 3.1) and references therein for further discussion of Cepheid systematics. In Appendices A.1 to  A.3 we test and report the dependence of the Leavitt law parameters on aspects of the Cepheid fit, namely outlier rejection, distance anchor, and cut on Cepheid period.

Careful treatment of Cepheids starts with their discovery and identification, where crowding and confusion can lead to misidentification. Light from a Cepheid can be blended with nearby or background sources, and aliasing or sampling problems can cause the wrong period to be inferred. Thus, outliers from the Leavitt law fit must be identified and rejected. Moreover the intrinsic scatter in the Leavitt law must be taken into account in assessing the goodness-of-fit; outliers that are rejected should lie well outside the so-called instability strip.

The secondary dependence of Cepheid luminosities on atmospheric metallicity is an ongoing area of research, and remains contentious. This effect arises from changes in the atmospheres and structure of Cepheids with their chemical composition, which impacts colours and magnitudes. There is evidence of a mild metallicity dependence at optical wavelengths (Kennicutt et al., 1998; Sakai et al., 2004; Macri et al., 2006; Scowcroft et al., 2009), which is weaker in the infrared. In the LMC, using spectroscopic [Fe/H] measurements, Freedman & Madore (2011) find that $Z_{H}$ (the metallicity dependence in the H-band) is close to zero. Efstathiou (2014, section 3.2) argues that these LMC observations, along with theoretical considerations, give cause to applying an external prior on the metallicity dependence centred at $Z_{W}{\sim}0$. We discuss this prior, which we find is inconsistent with the R11 data, in Appendix A.4.

Historically the zero point of the Leavitt law has proven difficult to measure, due to uncertainties in parallax measurements. To circumvent this, more accurate absolute distances have been pursued, including VLBI measurements of water megamasers in NGC 4258 (Humphreys et al., 2013). Multiple distance anchors are also tested and combined to reduce the impact of any single distance anchor. The effects of varying and combining anchors is explored in this analysis in Appendix A.2, following R11 and E14.

4.3 Cepheid-only fit

Our fit to all Cepheid data is based on E14 with the difference that we do not assume the SN Ia zero point (the quantity $a_{V}$ in R11 and E14) or indeed any SNe data. This is because we intend to fit the Cepheids separately from the SN data, whereas E14 calculates values of $H_{0}$ from the Cepheid fits, assuming SN Ia data from R11. All Cepheid data are fit to the Leavitt law (Equation 16) with MultiNest (Section 2.3.1). The 12 parameters of fit include the three nuisance parameters $\{b_{W},Z_{W},M_{W}\}$, the strongly constrained distance $\mu_{4258}$, and the eight distance modulus offsets $\{\Delta\mu_{i}\}$. We set an external Gaussian prior on $\mu_{4258}$, and by default place uniform priors for all other parameters over generous intervals. The $\chi^{2}$ function for the Cepheid fit is a function of $\{b_{W},Z_{W},M_{W},\mu_{4258}\}$, and $\{\Delta\mu_{i}\}$, and takes the form

[TABLE]

Here $m_{Wij,\mathrm{mod}}(b_{W},Z_{W},M_{W},\mu_{4258},\{\Delta\mu_{i}\})$ is the model magnitude of the $j$-th Cepheid in galaxy $i$ (given by Equation 16) and $\sigma_{\mathrm{int,C}}$ is the intrinsic scatter in Cepheid magnitude, from the width of the instability strip. For clarity, measured quantities are denoted with hats to distinguish them from model quantities. The logarithm of the likelihood $\mathcal{L}=e^{-\chi^{2}_{c}/2}$ and the priors on the fit parameters are inputs for MultiNest. We use 1000 live points in MultiNest and confirm that the precision is sufficient.121212For selected fits we repeat the outlier rejection and fitting steps, and find that the scatter in final parameters within ten runs is $<1\%$ of the statistical uncertainty.

4.4 Results of Cepheid-only fit

The results of all Leavitt law fits, for all combinations of distance anchor, outlier rejection, and upper period limit, are presented in Table 11. The details of these choices are given in Appendices A.1–A.3, along with the effect they have on fit results. The variation in the fits is visualised in Fig. 1 in $b_{W},Z_{W}$-space. The choice of these two parameters is obvious as they characterise the Leavitt law and are solely influenced by the Cepheid sample – all other parameters in $\Theta$ are influenced by the SN data, even the zero point $M_{W}$. Fig. 1 allows us to identify which of the Cepheid fits lie at the edges of the parameter space. The resultant scatter observed in Fig. 1 far exceeds the statistical uncertainties reported in Table 11. Therefore it is paramount that the systematic associated with varying the choices made in Appendices A.1 through A.3 is propagated carefully through the entire analysis process.

The choice of whether or not to apply the upper period limit of $P<60$ days has the most effect on the parameters, especially $b_{W}$. Fig. 1 reveals clearly the impact of including the longer-period Cepheids on parameters $b_{W}$ and $Z_{W}$, most notably splitting Fig. 1 down the middle vertically, i.e. by Leavitt law slope. Both parameters are smaller in magnitude by ${\sim}0.1$ when the longer-period Cepheids are included, indicating a weaker dependence of Cepheid magnitude on both period and metallicity. For the slope $b_{W}$, this difference dominates the statistical uncertainty and any other variation in $b_{W}$, whereas for $Z_{W}$ the resultant change from changing the period cut is comparable in size to the dependence on rejection algorithm, and the statistical uncertainty.

When the longer-period Cepheids are included, each of $b_{W}$ and $Z_{W}$ is better-constrained by the distance anchor, and the rejection algorithm, respectively: this is reflected in the vertical lines of empty markers with the same shape, and near-horizontal lines of markers with the same colour. That is, when the $P<60$ day cut is applied, the fit results are more sensitive to the choice of rejection algorithm and distance anchor. However, even without the cut, there remain strong dependences of $Z_{W}$ on rejection, and of $b_{W}$ on distance anchor.

Within each choice of period cut, the slope $b_{W}$ varies systematically with distance anchor: the NGC 4258 + MW and NGC 4258 anchors result in the lowest and highest $b_{W}$ respectively, with results form the other anchor combinations lying in between. The fits with both the LMC and MW in the anchors (upward triangles and diamonds) have the least spread in both parameters. With the exception of these fits, the data suggest a reasonably strong metallicity dependence with $-0.5<Z_{W}<-0.2$. As noted above, the results are sensitive to rejection algorithm, with the R11 rejection resulting in less negative values for $Z_{W}$ (and for $b_{W}$ with the $P<60$ day cut), followed by the E14 rejection with $T=2.5$ to a smaller extent.

We observe (Table 11) that there is little difference in values for $M_{W}$ between fits with and without the $P<60$ day cut, with the difference decreasing to zero for fits anchored on both n4258 and the Milky Way. However we defer further comment on $M_{W}$ (as well as $\{\Delta\mu_{i}\}$) to the discussion of global fit results. As $M_{W}$ is a magnitude zero point, and the $\Delta\mu_{i}$ are affected by the nearby SNe, the values of these parameters have potential to be influenced by the SN Ia data, and are expected to change with their inclusion.

4.4.1 Comparison to R11 and E14

We compare our fits in Table 11 to equivalent results in R11 and E14: our fits with R11 rejection and no period cut are compared to bolded fits in R11, table 2, and we compare our fits with the $P<60$ day cut to the results in E14, tables 2–4 without priors on $b_{W}$ and $Z_{W}$. Relative to R11, our $b_{W}$ values with LMC-only or MW-only anchors are slightly lower in magnitude (${\sim}-3.12$ instead of $-3.19$ in R11, a ${\sim}1\sigma$ difference). Moreover our fits with LMC + MW anchors result in a lesser metallicity dependence ($-0.2<Z_{W}<-0.1$ instead of $Z_{W}{\sim}-0.3$ in E14) – however uncertainties in $Z_{W}$ in these E14 fits exceed 0.1, and our $Z_{W}$ values (without the period cut) are supported by R11. Aside from these differences our results are in good agreement with R11 and E14, lying well within ranges allowed by statistical uncertainties. We retain 444 Cepheids when adopting the rejection flagged in R11, table 2 (close to the minimum of 448 reported in table 4 of R11) and only 379 with the $P<60$ day restriction. Applying the E14 rejection algorithm, our fits consistently result in lower numbers of remaining Cepheids by $10-20$, and consequently slightly lower $\sigma_{\mathrm{int,C}}$. It is worth noting that our methodology differs from E14 (and R11) in that we do not involve any SNe in the fit (omitting the third term in equation 14 of E14), whilst E14 includes the SN fit results by assuming a value of $a_{V}$ taken from R11. Presumably the complex ways of probing the multi-dimensional parameter space are leading to differences, albeit slight, between this work, R11, and E14, that cannot be easily reconciled. We believe a solution for the future is for authors to provide code and data sets used for calculations as part of publication, that can be used to better understand differences.

4.4.2 Selection for global fits

The choice of Cepheid fits to carry forward to the global fit is informed by their results, i.e. Leavitt law slope $b_{W}$ and metallicity dependence $Z_{W}$, as these parameters are only influenced by the Cepheid sample and are very minimally affected by the SN data. We are interested in the effect the choice of Cepheid sample (through varying aspects of the fit such as distance anchor, rejection, and upper period limit) has on these parameters in the global fit. In particular it is essential to quantify the systematic uncertainty in $\mathcal{H}$ with varying these choices.

We select 18 fits in total, summarised in Table 2. To span the full range of uncertainty induced by various Cepheid fits, we select fits at extremes of the parameter space (Fig. 1), with a selection of anchors and rejection algorithms. The combination of all three distance anchors has the most constraining power, so we include all of these fits to quantify the uncertainty within them.

Each fit has an associated set of best-fitting parameters with uncertainties, as well as (unless using the R11 rejection) the values of the intrinsic scatter and rejection threshold, which together uniquely define a set of Cepheids remaining after outlier rejection. These then make up the Cepheid data and priors for some parameters in $\Theta$, going into the global fit (Section 6).

5 Supernova fit

We now focus on the Type Ia supernovae. First we outline the data set and discuss systematics, then we describe various cuts on the SNe Ia and present the preliminary SN-only fit. This section is supplemented by details provided in Appendices B.1, B.2, B.3, B.4 on the lightcurve fitting method, and the corrections applied for Malmquist bias, host galaxy mass and peculiar velocities respectively. Also central to the subject are the computations of SN systematics in covariance matrices, which are also relegated to Appendix C for detailed discussion.

5.1 Observations

Our supernova data are identical to R11, consisting of eight ‘nearby’ SNe Ia in the galaxies hosting Cepheids (Table 1), and 280 unique ‘low-$z$’ SNe Ia from the 185 CfA3 (Hicken et al., 2009a) and 165 LOSS (Ganeshalingam et al., 2010) samples.131313There are 69 SNe in common between the samples; however SN 1998es was discarded because the lightcurve quality was so poor that the SALT2 lightcurve fit failed. Details of sources of photometry for the nearby supernovae are presented in Table 3. Natural photometry was not available for the oldest two, SN 1981B and SN 1990N. The most recent SNe are already in both CfA3 and LOSS, so we used combined photometry from both sources as described in Appendix B.1.1. The remaining three (SN 1994e, SN 1995al, SN 1998aq) were observed on the FLWO 1.2 m telescope with a variety of CCDs; we construct SALT2 instruments (including transmissions and zero points) using data from Jha et al. (2006).

CfA3 ran from 2001 to 2008 on the 1.2m telescope at FLWO almost entirely with the CfA3 4Shooter2 and Keplercam imagers (in UBVri and UBVRI filters respectively), while LOSS took place on the NICKEL and KAIT telescopes from 1998 to 2008 (in BVRI). Unlike more recent magnitude-limited surveys, CfA3 and LOSS targeted known galaxies and include SNe discovered by other sources, resulting in a more complex selection function and generally resulting in higher host galaxy masses (Appendix B.3). We refer the reader to the above works for further details of observations. Newer low-$z$ SNe Ia samples have since been published, notably CfA4 (Hicken et al., 2012), Carnegie Supernova Project (CSP-II; Contreras et al., 2010), Pan-STARRS (Rest et al., 2014), Palomar Transient Factory (PTF; Law et al., 2009), and La Silla-QUEST Supernova Survey (LSQ; Walker et al., 2015). However we retain the older CfA3-LOSS sample for this analysis to more faithfully compare our results to R11 and E14.

Photometry for the low-$z$ sample is sourced from Hicken et al. (2009b) and Ganeshalingam et al. (2013) in the natural systems of each filter set, with the exception of the CfA3 4Shooter2 and Keplercam U filters for which reliable measurements do not exist – we use photometry in the standard Johnson-Cousins UBVRI system as presented in Bessell (1990) for these passbands only, as well as the nearby SN 1981B and SN 1990N. We use SALT2 (Guy et al., 2007) to fit these SN Ia lightcurves for the quantities $m_{B},X_{1},C$, which are used to derive distances via Equation 14. Details of the lightcurve fitting are given in Appendix B.1.

One reason for our choice of SALT2 as a lightcurve fitter is that our framework for assessing SN Ia systematic uncertainties with covariance matrices (Section 5.2 and Appendix C) follows that in the SNLS-SDSS Joint Lightcurve Analysis (hereafter JLA; Betoule et al., 2014), which relies on the SALT2 model. In addition, SALT2 is the most modern fitter and used ubiquitously in cosmology analyses; thus our use allows for easier comparison and greater consistency. While R11 test the effects of fitting lightcurves with both SALT2 and MLCS2k2 (Jha et al., 2007)141414SALT2 differs from MLCS2k2 substantially in its treatment of extinction: instead of prescribing a reddening parameters $R_{V}$, all of the colour information (including the SN’s intrinsic colour and host extinction) is included in the single colour parameter $C$. lightcurve fitters, we use SALT2 only. This is justified, as the latest version SALT2.4 (described in Betoule et al., 2014) was released in parallel with simulations in Mosher et al. (2014) which assess and quantify the uncertainty associated with the choice of lightcurve fitter (and the lightcurve model itself) in covariance matrices (Appendix C.1). Hence it is unnecessary to use of multiple fitters to assess the the aforementioned systematic uncertainty.

5.2 Supernova systematics

As a statistical sample, type Ia supernovae are high fidelity standard candles. However as astronomical objects, SNe Ia are diverse and subject to systematics, with their measurable quantities (absolute brightness, observed colour and decline rate) dependent on factors which correlate with their progenitors and environments. Countless investigations into these correlations and their origins are partly motivated by the need to reduce residual scatter from these intrinsic SN Ia variations. Observations of supernovae are also influenced by factors such as galactic extinction, misclassifications, and differing telescope magnitude systems. Most of these effects are not sufficiently well-understood or accurately modelled to correct for them entirely. It is therefore essential to quantify the size of systematics; even when efforts have been made to apply corrections we still wish to estimate the uncertainty in the correction.

Our approach to accounting for SN Ia uncertainties follows methods in JLA, which are largely based on those in the Supernova Legacy Survey (hereafter SNLS; Conley et al., 2011). These use individual covariance matrices for each systematic, tracking correlated uncertainties between different SN quantities (i.e. $m_{B},X_{1},C$), between different supernovae. Advantages of the covariance matrix method over the more traditional method of adding systematics in quadrature are discussed in Conley et al. (2011, section 4); these include the ability to fully capture correlations in uncertainties, and the ease of including or reproducing uncertainties in further analyses. Details are our computations are provided in Appendix C.

5.3 Cuts on supernova sample

We make quality cuts on our SN Ia sample to eliminate potential biases from poorly constrained lightcurves and peculiar events, and to remain within the bounds of the SALT2 model. With the intent of replicating the sample in R11 as closely as possible, we broadly follow the cuts described in CfA3 (Hicken et al., 2009b) and LOSS (Ganeshalingam et al., 2013), also using cuts in SNLS and JLA – described in Guy et al. (2010, section 4.5), Conley et al. (2011, section 2.1), Betoule et al. (2014, section 4.5) – as guidance or as alternate cuts. In summary our criteria are as follows:

• •

low Milky Way extinction $E(B-V)<0.2$

• •

exclude local SNe Ia not in the Hubble flow $z>0.01$

• •

goodness-of-fit from SALT2 $\chi^{2}/\rm{DoF}<8$

• •

first detection by +5 days, relative to B-band maximum

• •

exclude stretch outliers $|X_{1}|<3$

• •

exclude colour outliers $|C|<0.5$

• •

well-constrained stretch $\sigma_{X_{1}}<0.8$

• •

well-constrained colour $\sigma_{C}<0.1$.

The above encompass cuts in CfA3 and LOSS, with stricter cuts on the date of first detection and lightcurve goodness-of-fit (originally at +10 days and $\chi^{2}/\rm{DoF}=15$ in CfA3 respectively), and with additional cuts on the uncertainties in $X_{1}$ and $C$ to further exclude supernovae which have large uncertainties in their stretch or colour. Our cuts are also informed by visual inspection of individual lightcurves and their SALT2 fits, particularly in placing boundaries for the lightcurve goodness-of-fit, uncertainties in stretch and colour, and date of first detection. In summary, we exclude supernovae at very low redshift (i.e. not yet in the Hubble flow), significantly extinguished by Milky Way dust, detected too late, with poorly constrained stretch and colour. We also exclude SNe Ia with poor SALT2 fits, and SNe that are too blue or red or have very fast or slow decline to exclude peculiar objects and ensure our sample fit within the SALT2 model.

Furthermore we test some alternate cuts, including some suggested in JLA and original CfA3/LOSS cuts which we have changed above. We repeat the SN-only fit with these cuts to test the effect on the SN fit parameters, carrying some through to the global fit. In particular, we follow R11 in raising the low-redshift cut to $z=0.0233$,151515This is to reduce possible bias from local coherent flows, or a possible local underdensity (a so-called ‘Hubble bubble’). The latter is discussed in R11 and Conley et al. (2011); however there is no conclusive evidence for its existence. and test strengthening or relaxing the lightcurve goodness-of-fit threshold to $\chi^{2}/\rm{DoF}<5$ or $\chi^{2}/\rm{DoF}<15$, and relaxing the date of first detection to +10 days. Following JLA we examine the effects of imposing a stricter bound on the colour ($|C|<0.3$), the uncertainty on the stretch ($\sigma_{X_{1}}<0.5$), and Milky Way extinction. These tests are important as the influence of these alternate cuts on the fit results is not straightforward or obvious; moreover no particular cut is necessarily more valid than the others. We discuss these results and their significance in Section 5.5. Histograms showing $X_{1}$ and $C$ distributions for several cuts are included in Appendix B.5.

5.4 SN-only fit

Analogous to the Cepheid-only fit in Section 4.3, we perform a preliminary fit of only the low-$z$ SNe Ia to Equation 18 using the MCMC routine emcee (Section 2.3.1), to identify the dependence of the SN parameters on the different cuts in Section 5.3. To clearly separate the data and model in Equation 18 we define the quantity $m_{B}^{\dagger}$ for the apparent SN magnitude corrected for stretch and colour:

[TABLE]

Explicitly the $\chi^{2}$ function for the low-$z$ SN fit is

[TABLE]

where the covariance matrix $\mathbf{C_{m_{B}^{\dagger}}}$ is derived from covariances in all SN parameters $\{m_{B},X_{1},C\}$, as given in Equation 38 in Appendix C along with detailed explanations of statistical and systematic contributions.

It is evident from Equations 18 and 21 that the SN-only fit is degenerate: we cannot constrain both $M_{B}$ and $\mathcal{H}$ simultaneously; the nearby SNe are necessary to constrain $M_{B}$. Instead we fit for the difference $\mathcal{M}:=M_{B}-\mathcal{H}$, adopting the blinds for each $M_{B}$ and $\mathcal{H}$ noted in Section 2.4 i.e. with the transformations $m_{B}\mapsto m_{B}+\mathfrak{o}_{M_{B}}$ in Equations 17 and 18 and $\mathcal{H}\mapsto\mathcal{H}+\mathfrak{o}_{\mathcal{H}}$ in the likelihood incorporating Equation 22 (a function of both $M_{B}$ and $\mathcal{H}$ through Equation 21). The marginalised posterior distributions (mean and 1$\sigma$ width) for $\alpha$ and $\beta$ are presented in Table 4 and plotted in Fig. 2; these results are dependent on the choice of quality cuts on the SN sample described in Section 5.3.

5.5 Results of SN-only fit

The results of the SN-only $\chi^{2}$-minimizing fit are presented in Table 4, while Fig. 2 shows the differences in fits with various SN cuts lie in the $\alpha,\beta$-plane (this is analogous to Fig. 1, which displays the numerous Cepheid fits in $b_{W},Z_{W}$-space). We discuss the dependence of the fit results on the various cuts, and select cuts with results spanning the parameter space to carry forward to the global fit to assess the associated systematic uncertainty.

The notable outlier is the higher low-redshift cut ($z>0.0233$), effecting a much lower value of $\beta$ than the other cuts. This cut, along with the stricter colour and stricter goodness-of-fit cuts, results in lower $\alpha$ also. The lowest and highest values of $\alpha$ correspond to lower $\chi^{2}$, and stricter $\sigma_{X_{1}}$ respectively. Fig. 13 in Appendix B.5 shows normalised $X1$ and $C$ distributions for the $z>0.0233$ cut: there are marginally slower-declining SNe compared to the default, but overall the distributions appear similar. It does not appear that the discrepant fit results from this cut are the result of a change in the colour or stretch distribution of the sample; indeed our tests with jackknifed samples (described below) indicate this is likely the result of removing a large portion of the sample (over 40% relative to the default). Disregarding the $z>0.0233$ cut, the variation in $\alpha$ and $\beta$ with the different cuts we test appears only slightly larger than the typical statistical uncertainties in these parameters (Fig. 2).

We use jackknife resampling to assess the statistical significance of the dependence of results on the cuts in Table 4. For several cuts (the lower lightcurve $\chi^{2}$/DoF, higher redshift cut and stricter cuts on colour or uncertainties in stretch and colour) we draw subsamples of size $N_{\rm{SN}}$ (Table 4) of the 171 SNe selected by the default cut. For each cut we compare the systematic change in fit results (parameters $\alpha,\beta,\mathcal{M}$) from the new cut to results from repeated jackknifed subsamples of size $N_{\rm{SN}}$ and their scatter. These reveal a systematic variation of $1-3\sigma$ from the default for almost all combinations of parameters and cuts (where $\sigma$ is the scatter within the numerous jackknifed subsamples). Thus the differences between rows of Table 4 cannot be solely attributed to shot noise, and the variation due to different cuts must be propagated to the global fit (Section 6) and treated as a contribution to the total systematic uncertainty. However, we will find that the variation from the choice of SN cut is dwarfed by the analogous source of uncertainty from the choice of Cepheid fits.

6 Global Fit Results

This section contains our final simultaneous fits to all Cepheid and SNe Ia data. We set out parameters and equations for this fit, and present fit results for all parameters, including the dependence of results on choices within the individual Cepheid and supernova data sets. We summarise our uncertainties, and discuss their increase compared to other analyses of the same data. Finally, we break down the statistical and systematic contributions to the uncertainty budget.

6.1 Global fit

We fit all Cepheid and supernova data simultaneously to Equations 16, 17, 18 as described in Section 3.3. We minimize a global $\chi^{2}$ function (a function of $\Theta=\left\{\alpha,\beta,M_{B},\mathcal{H},b_{W},Z_{W},M_{W},\mu_{4258},\{\Delta\mu_{i}\}\right\}$), which has contributions from the Cepheids and low-$z$ SNe remaining after cuts (given in Equations 20 and 22 respectively), and also an equivalent term to $\chi^{2}_{\rm{low-}z}$ for the eight nearby SNe:

[TABLE]

The bolded quantities in Equation 24 are vectors, over the eight nearby SNe Ia. The terms contributing to the nearby covariance matrix $\mathbf{C_{m_{B}^{\dagger},n}}$ are covariances between SALT2 quantities $m_{B},X_{1}$, and $C$, and the diagonal intrinsic scatter $\sigma_{\mathrm{int,SN}}$.

The global simultaneous fit is 16- or 17-dimensional (without and with the LMC included as a distance anchor respectively), and performed using MultiNest as described in Section 2.3.1. We are ultimately interested in $\mathcal{H}$, which contains the value of $H_{0}$. However to demonstrate degeneracies and correlations between parameters, we display in Figs 3 and 4 marginalised contour plots of the posterior distribution of an example fit (with $T=2.25$, NGC 4258+LMC+MW anchor, $P<60$ day cut Cepheid fit and default SN cuts).161616Figs 3, 4, and 9 were created with the ChainConsumer package (Hinton, 2016). The former posterior distribution is marginalised over the eight $\Delta\mu_{i}$, while the latter is also marginalised over $\mu_{4258}$ and the SN and Cepheid parameters which are strongly constrained by initial fits: $\{\alpha,\beta,b_{W},Z_{W}\}$. Fig. 3 shows a strong positive correlation between $M_{B}$ and $\mathcal{H}$ as expected from their degeneracy in the low-$z$ SN sample (Equation 18), and less apparent correlations between the ‘zero point-like’ parameters $\{\mathcal{H},M_{B},M_{W},\mu_{4258}\}$. In contrast, the five other parameters $\{\alpha,\beta,b_{W},Z_{W},\mu_{4258}\}$ each are largely independent of the other parameters (Fig. 3).

We repeat the global fit for each of 18 Cepheid fits in Table 2 and six supernova cuts determined in Section 5.5 from Fig. 2. Each Cepheid fit and SN cut corresponds to a subset of the total sample to use in the global fit, and associated values of best-fitting parameters, as well as $\sigma_{\mathrm{int,C}}$ for the Cepheids. In total there are 108 fits; the analysis of these results and the variation therein follows.

6.2 Results of global fit

The best-fitting values and uncertainties of parameters in $\Theta$ are given in Table 13 for each of 108 fits. Fig. 5 displays these fits in various subspaces of the 16- or 17-dimensional space spanned by $\Theta$, focussing on parameters $\{\alpha,\beta,b_{W},Z_{W},M_{W},M_{B},\mathcal{H}\}$. We discuss the dependences that this figure shows (which motivate the averaged tables and figures later), then present results for the nuisance parameters and the parameters of interest: the SN Ia peak absolute magnitude $M_{B}$ and (proxy for the) Hubble constant $\mathcal{H}$, which are degenerate with each other.

In the remainder of the section we depart from the distinction we make between statistical and systematic uncertainties in Appendix C.2: the uncertainties returned by MultiNest, reported in Table 13, simply the 1$\sigma$ widths of the PDFs, do not distinguish between the statistical and systematic components of covariance matrices input into the fit in the likelihood. Henceforth, we refer to this uncertainty from the MultiNest fit as statistical, and the variation observed in e.g. Fig. 5 between global fits with differing supernova cuts or Cepheid fits as systematic.

6.2.1 Dependence of parameters

Fig. 5 highlights the following dependence of parameters on properties of the global fit:

The Cepheid parameters $\{b_{W},Z_{W},M_{W}\}$ depend only on the choice of Cepheid fit (carried forward from Section 4.4.2), reflecting the variation observed in Fig. 1. Thus there is negligible scatter in values for these parameters between fits with the same Cepheid data, regardless of the SN cut (Fig. 5(a), (b)). 2. 2.

Similarly, the SN parameters $\{\alpha,\beta\}$ depend most strongly on cuts, and minimally on Cepheid fit, although there is more scatter than in $\{b_{W},Z_{W}\}$. On average fits without an upper period cut on the Cepheids result in slightly lower $\alpha$ by ${\sim}0.01$, for each SN cut (Fig. 5(c)). 3. 3.

The Cepheid and SN zero points $M_{W}$ and $M_{B}$ both depend predominantly on the Cepheid fit (Fig. 5(b)), reflecting the fact that the SNe Ia are calibrated on the Cepheids. While $M_{W}$ depends directly on the Cepheid data (Equation 16), the influence on $M_{B}$ is through its interaction with $M_{W}$ via the distance modulus offsets $\{\Delta\mu_{i}\}$ (Equations 16 and 17). We note that $M_{W}$ has negligible dependence on SN cut, whereas $M_{B}$ varies slightly with the choice of SN cut (with a spread of ${\sim}0.01$ within each choice of Cepheid fit). 4. 4.

As mentioned in Section 3.4, $\mathcal{H}$ is degenerate with $M_{B}$. Fig. 5(d) shows this degeneracy between the parameters, and that the difference $\mathcal{M}=M_{B}-\mathcal{H}$ lies on a straight line. Within each choice of Cepheid fit there is slight systematic dependence only on the choice of SN cut. There is no systematic difference between these parameters from fits with and without a cut on Cepheid period. 5. 5.

$\{\Delta\mu_{i}\}$: the values of the distance modulus offsets from the global fit depend significantly on the Cepheid fit, as shown in Figs 7 and 7.

In summary, it is expected that the SN cuts determine parameters $\{\alpha,\beta\}$, and the Cepheid fits determine parameters $\{b_{W},Z_{W},M_{W}\}$. However the interaction of the ‘zero point-like’ parameters is more subtle, and emerges from the simultaneous fit of the three data samples, most obvious in Fig. 5(b). Even though the parameter $M_{B}$ only appears in the SN apparent magnitudes (Equations 17, 18), it is most strongly influenced by the Cepheid data via $M_{W}$, as the two parameters are tied to their respective data sets through the distance modulus offsets $\{\Delta\mu_{i}\}$. Furthermore, $M_{B}$ and $\mathcal{H}$ are degenerate with their difference determined by the low-$z$ SNe. Thus the resultant value of $\mathcal{H}$, hence $H_{0}$, is sensitive both to the choice of SN cut (via the $M_{B}-\mathcal{H}$ degeneracy) and to the choice of Cepheid fit (via the influence of $M_{W}$ on $M_{B}$). Unsurprisingly, the most extreme values of $M_{B}$ and $\mathcal{H}$ (both driven by $M_{W}$, as seen in Fig. 5(b)) arise from Cepheid fits anchored on only the LMC or MW (most and least negative, represented by dark purple and pink symbols, respectively). It is clear from Fig. 5(d) that the variation with Cepheid fits (anchor and rejection) is at least an order of magnitude larger than the variation with SN cuts, even when the fits anchored on the LMC or MW only are excluded.

6.2.2 Nuisance parameter results

Tables 5 and 6 contain results for the supernova and Cepheid nuisance parameters, averaged over the Cepheid fits and SN cuts respectively. We choose to average over these aspects of the fit that have minimal effect on the parameters, as shown in Fig. 5: the SN parameters in (c) predominantly depend on shape (SN cut) and not on colour (Cepheid fit), while the Cepheid parameters in (a) depend entirely on colour and not on shape. We omit statistical uncertainties of parameters in these tables as they can be obtained from the full set of results in Table 13. For the nuisance parameters we select a single best fit (bolded in Table 13 and indicated in Fig. 5). This is preferable to averaging over results in Tables 5 and 6, which are asymmetric, based on different premises (e.g. different distance anchors), and include more questionable fits (e.g. those SN cuts that reject a larger fraction of the total). Thus we use the maximal variation in these values to inform our systematic error budget, but not to influence the best fit.

Final values for the nuisance parameters are taken from the bolded reference fits, which have the default SN cut and the Cepheid fit with all three anchors, $T=2.25$ rejection, and no cut on Cepheid period. We have chosen this fit because the results are representative and centred amongst the different choices. The Cepheid fit here also aligns with fits selected in E14 and R11. As in R11, we choose to not impose a cut on Cepheid period, and note the effects of including this cut on nuisance parameters described in Section 4.4: both $b_{W}$ and $Z_{W}$ are more negative with the $P<60$ day cut, while there is no difference in $M_{W}$ when all three distance anchors are used. The global SN results differ slightly from the initial results in Table 4, and we again note that the most deviant (lowest) values of $\alpha$ or $\beta$ are where a large number of SNe have been rejected; the remaining cuts are in agreement with values derived from the default cut. In summary, the fit parameters and uncertainties from Tables 5 and 6 are:

[TABLE]

These statistical uncertainties are found from Table 13. We generally take the maximal variation measured from the reference fits in Tables 5 and 6 as the systematic uncertainty, with the following exceptions. We disregard the higher low-redshift cut (associated with a large fraction of the SNe being discarded) in estimating the systematic uncertainty in $\beta$ – see discussion in Section 5.5. For the uncertainty in $Z_{W}$, we only consider the variation over fits which include both the LMC and MW in the distance anchor: the constraints on the metallicity dependence provided by different distance anchors are inconsistent with each other, so we only consider these fits for estimating the uncertainty for the nuisance parameter $Z_{W}$ alone (i.e. the other anchors are considered for estimating uncertainties on $M_{B}$ and $\mathcal{H}$, in Section 6.2.3.) From Fig. 5(c) it is clear that the statistical uncertainties in the SN parameters are around double the systematic uncertainty if we disregard the higher low-redshift SN cut. The opposite is true for the Cepheid parameters, where the statistical uncertainties are dwarfed by systematic variation with differing fits. If we restrict our analysis to only Cepheid fits anchored on all three galaxies, the statistical and systematic uncertainties are comparable in size.

The systematic errors are asymmetric for most parameters, especially for $\beta$ (due to the outlying $z>0.0233$ cut) and $Z_{W}$. This can be observed in Fig. 5(a)–(c), where it is evident our reference fits do not lie centrally within the parameter subspaces. Fig. 5(b) shows that the MW as a distance anchor drives $M_{W}$ up, while the LMC (and to a lesser extent NGC 4258) drives $M_{W}$ down, an effect which propagates to $M_{B}$ and $\mathcal{H}$ (Fig. 5(d)). Fits anchored on all three distance anchors lie centrally. Our Cepheid nuisance parameters remain consistent with R11 and E14 as initially found in Section 4.4.1.

We note that our best-fitting value for $\alpha$ is significantly higher than found in JLA (Betoule et al. (2014), table 10) and LOSS (Ganeshalingam et al., 2013), by ${\sim}0.02$ (around double the total uncertainty in $\alpha$). This difference occurs consistently over a range of SN cuts. While the JLA analysis always determines $\alpha$ from the low-$z$ sample in conjunction with a higher-redshift sample, Ganeshalingam et al. (2013) finds $\alpha=0.146\pm 0.007$ from the LOSS sample, which overlaps with ours considerably and is over a similar redshift range. Our results for $\beta$ are consistent with the literature with the exception of the $z>0.0233$ SN cut, which results in a value ${\sim}1\sigma$ below the other cuts (the triangles in Fig. 5(c)). The impact of this cut on $\mathcal{H}$ can be seen in Fig. 5(d): the triangles (higher low-$z$ cut) have higher $\mathcal{H}$ than the other shapes (cuts) for each colour/fill (Cepheid fit). This effect is much smaller than the differences from varying the Cepheid fit. Nevertheless, it is in agreement with the increase of $H_{0}$ with increasing low-$z$ observed in R16, fig. 13.

The remaining nuisance parameters are the distance modulus offsets $\{\Delta\mu_{i}\}$, which, like all other nuisance parameters, are eventually marginalised over. Their values depend primarily on the Cepheid fits. The full table of fit values is left to Table 12 in Appendix D. The $\Delta\mu_{i}$ are visualised in Figs 7 and 7 with different colour/fill representing Cepheid fit. Fig. 7 gives some insight into the interplay and correlations between distance moduli of different galaxies, while Fig. 7 shows the scatter and relative values of the $\Delta\mu_{i}$ from different fits. The statistical uncertainties in $\Delta\mu_{i}$ from individual fits range from 0.05 to 0.1, and is comparable to the scatter over different fits.

6.2.3 Results for $M_{B}$ and $\mathcal{H}$

We now consider the parameters $M_{B}$ and $\mathcal{H}$ which, together, directly reveal $H_{0}$. The degeneracy between them is apparent in Fig. 5(d), which also shows that their primary dependence is on the Cepheid fits. Thus in Table 7 we present the global fit results averaged over the SN cuts.171717We report $M_{B}$ and $\mathcal{H}$ to 3 decimal places, unlike most other parameters which have been truncated to 2 decimal places (but not rounded in the analysis). These two quantities are of particular interest, and it is desirable to retain precision in both their values and uncertainties throughout this section. Given that the fits in Table 7 anchored on all three galaxies are spread out, we average these fits rather than choose a best fit, and take the maximal variation in these fits as the systematic uncertainty. There is a slight systematic difference between fits in Table 7 with and without the upper period limit (on average, $\mathcal{H}$ is decreased by 0.015 mag where the $P<60$ day cut is applied). From a theoretical standpoint, we have no reason to preference one cut over the other. Thus our best estimates for $M_{B}$ and $\mathcal{H}$ are averaged over all fits anchored on all three galaxies (including fits both with and without the upper period limit), represented by solid and empty navy, green and dark purple markers in Figs. 5–8.

Our best estimates are

[TABLE]

Here the statistical uncertainties are found from relevant fits in Table 13, in which a representative fit is bolded (with default SN cut, $T=2.25$ rejection, and no upper period limit). The above systematic uncertainties are given by the maximal variation in values with the combined NGC 4258+LMC+MW distance anchor. We impose this constraint on the anchor so that we can fairly assess the systematic uncertainty when all available distance information is used, and to allow better comparison with R11 and E14 who primarily report errors with all three anchors. In Section 6.3 we investigate and discuss uncertainties in $\mathcal{H}$, including converting from an absolute error in the logarithmic quantity $\mathcal{H}$ to a relative error in $H_{0}$.

We next consider fits anchored on NGC 4258 only, to estimate the uncertainty in this case, and for the sake of comparison with R11 and E14. These fits are represented by the empty turquoise and pink markers, and by all red markers in Figs. 8 and 5. We average these results from Table 7 to find Equation 27, as with Equation 26. The systematic uncertainties are given by the maximal variation in values derived from these fits, and the statistical uncertainties are found from NGC 4258-anchored fits in Table 13, with one representative fit bolded.

[TABLE]

The resultant value of $\mathcal{H}$ in Equation 27 is 0.05 mag lower (corresponding to a 2.3% decrease in $H_{0}$) compared to where all three anchors are used (Equation 26). Moreover, $M_{B}$ (which is largely degenerate with $\mathcal{H}$) is also 0.05 mag lower (brighter). The systematic uncertainty (i.e. the spread in values between different fits) is the same, while the statistical uncertainties are larger, reflective of the fact that a distance scale is anchored on a smaller set of data.

Our best estimate of the peak SN Ia brightness $M_{B}$ (in Equation 26, from the three-galaxy anchor) appears mildly higher (dimmer) than values reported in JLA (assuming $H_{0}=70$ km s*-1* Mpc*-1*), which are $M_{B}=-19.05\pm 0.02$ from all SN Ia data, or $M_{B}=-19.02\pm 0.03$ from a lower-redshift subsample consisting of low-$z$ and SDSS supernovae (table 10 of Betoule et al. (2014)). However, the supernova-only fit in JLA cannot constrain both $M_{B}$ and $H_{0}$, which are degenerate. As they have assumed a value for $H_{0}$ (while we have fitted separately using a distance ladder), our numerical values for $M_{B}$ are not directly comparable, but merely reflect the influence of different values of $H_{0}$.

Returning to $H_{0}$, Equation 26 corresponds to a value of $H_{0}=72.5\pm 3.1(\rm{stat})\pm 0.77(\rm{sys})$ km s*-1* Mpc*-1* (total uncertainty of 4.4%) from the combined NGC 4258+LMC+MW anchors. If we assume the older distance $\mu_{4258}=29.31$ in R11 (Footnote 11), our best estimate increases to $H_{0}=73.8\pm 3.2(\rm{stat})\pm 0.78(\rm{sys})$. These central values agree with R11 ($H_{0}=73.8\pm 2.4$) and E14 ($H_{0}=72.5\pm 2.5$), which respectively assume $\mu_{4258}=29.31$ and 29.404. Using only NGC 4258 as a distance anchor (and the new Humphreys et al. (2013) value of $\mu_{4258}=29.404$) gives $H_{0}=70.9\pm 3.5(\rm{stat})\pm 0.75(\rm{sys})$ km s*-1* Mpc*-1*, which is $2.3\%$ lower than with the three anchors. The uncertainties in $\mathcal{H}$ are broken down in Section 6.3.1 and summarised in Table 8. We next discuss the uncertainties in $H_{0}$; their size informs the significance of the tension between values of the Hubble constant measured using different probes, so they are of equal interest to the values.

6.3 Uncertainties

We have emphasized the importance of quantifying and incorporating the scatter in parameters arising from varying aspects of the SN and Cepheid fits, and indeed we use this overall variation in results to gauge the systematic uncertainty in these parameters. However, we have also seen that the statistical uncertainty dominates for the supernova parameters $\alpha$ and $\beta$ (Fig. 2 and 5(c)), as well as for $M_{B}$ and $\mathcal{H}$ when only considering the systematic variation between fits with all three anchors (Fig. 5(d)). This dominance reflects the fact that the SN Ia samples, especially the nearby sample (i.e. in Cepheid hosts), are relatively small with large errors when compared to the Cepheids. Hence the SNe are statistically limited while the Cepheids are not.

For clarity we divide the contributions to the total uncertainty in the parameters into three classes:

The statistical (in the usual sense) portion of the uncertainty reported by MultiNest, which is dominated by noise in the nearby and low-$z$ supernovae. 2. 2.

The systematic elements of $\mathbf{C_{\eta}}$, which make up remainder of the uncertainty reported by MultiNest, listed in Table 9. 3. 3.

The systematic uncertainty associated with varying aspects of the fit between reasonable alternatives is dominated by the variation in the choice of Cepheid fit, as shown in Figs 5 and 8. For our final value and uncertainty of $H_{0}$ we focus on fits with all three anchors only (with some consideration of fits with only NGC 4258 as an anchor for the sake of comparison to R11 and E14). Then in effect we are only considering the variation with the rejection algorithm and the cut on Cepheid period.

6.3.1 Uncertainties in $H_{0}$

We now address the uncertainty in the Hubble constant $H_{0}$ explicitly, using results in Section 6.2.3 (Equations 26 and 27). As the quantity $\mathcal{H}$ is related to the logarithm of $H_{0}$, its absolute error informs the relative error in $H_{0}$, via

[TABLE]

Table 8 summarises our calculations of the final uncertainty in $H_{0}$ from Equations 26 and 27. We find using Equation 28 relative errors in $H_{0}$ of 4.3% statistical and 1.1% systematic (corresponding to terms (i) and (ii) combined, and (iii) respectively as described at the start of Section 6.3) from all three distance anchors. From using only NGC 4258 as an anchor, these errors are 4.9% statistical, 1.1% systematic, 5.0% total. The final uncertainty in $H_{0}$ (the bottom row of Table 8) is the quadrature sum of the above statistical and systematic terms. Table 8 is comparable to the lower portion of table 5 of R11 (and subsequently table 7 of R16), which lists all systematic and statistical uncertainties contributing in quadrature to the uncertainty in $H_{0}$.

6.3.2 Increase in error compared to R11 and E14

Our final uncertainty in $H_{0}$ is 4.4% total (4.3% statistical and 1.1% systematic, with the statistical term inclusive of contributions from SN Ia covariance matrices) for the NGC 4258+LMC+MW distance anchor, which is significantly larger than previously found for the same data set (by 1% absolutely, or a $\sim$20% increase): R11 and E14181818E14 adopts the error in R11 for the SN Ia-side of the analysis. report total uncertainties of 3.3%191919The errors reported in table 5 of R11 are: 2.9% statistical, 1.0% systematic, 3.1% total. However the final error given with all three distance anchors conservatively includes the larger statistical error associated with using two distance anchors instead of three, resulting in a total of 3.3%. and 3.4% respectively. If NGC 4258 alone is used as a distance anchor, the above errors increase to 5.0% total (4.9% statistical and 1.1% systematic) for our fit, 4.1% (4.0% statistical, 1.0% systematic) from R11, and 4.7% total from E14. The difference between our errors and those found in E14 is smaller with the NGC 4258 anchor compared to when all three anchors are used – however, both are significantly larger than found in R11. For the remainder of Section 6.3.2, our discussion of errors pertains to fits with all three distance anchors.

Although it appears that the increase in our error lies in the statistical term (with the systematic term remaining the same), it is important to note the significant differences in how these terms are derived and defined in this work (given in points (i)–(iii) at the start of Section 6.3), compared to R11. Explicitly, the covariance matrices which quantify our SN Ia systematic terms directly contribute to the statistical errors in our global fits (i.e. increase the widths of PDFs) via the likelihood, while our systematic term contains the variation in parameters resulting from changing features of the fits. In comparison, the errors in each part of the calibration chain from the distance calibrators to the SNe Ia are separately given in R11, table 5. The total uncertainty is a quadrature sum of these individual terms, and the systematic variation described in R11, section 4.

The two major differences in our analysis which potentially contribute to the increased error are the inclusion of the supernova covariance matrices, detailed in Appendix C.1, and the simultaneous fit to all parameters, described in Section 3.3. As outlined above, it is not possible to make a direct comparison between contributions to our error and errors given in R11, with the aim of isolating the source of the discrepancy. However, we attempt to separate the influences of the covariance matrices and simultaneous fit, replicating the quadrature sum in R11 as closely as we can below.

First, we isolate the effect of the supernova covariance matrices alone on the size of uncertainties, by considering the error in the intercept of the SN Ia $m-z$ relation: this is $\mathcal{M}$ determined from our SN-only fit (Table 4 in Section 5.4), and is equivalent to $5a_{V}=3.485\pm 0.010$ in R11. Our error in $\mathcal{M}$ is ${\sim}0.036$, over three times larger than in R11. This error is roughly halved to 0.019 if we only consider the strictly statistical covariance matrix, i.e. $\mathbf{C_{\rm{stat}}}$ in Equation 40.202020Neglecting the uncertainty from the finiteness of the SALT2 training sample reduces the error slightly to 0.017, which reflects the statistical error in the observed SNe Ia only (i.e. the tridiagonal matrix $\mathbf{C_{\rm{stat,diag}}}$). For the same supernova data, our statistical-only error exceeds the total error in R11. Including the SN Ia systematic covariance matrices doubles the error again. We infer that the increase in error in this analysis compared to R11 is attributable to both the covariance matrix method of accounting for correlated SN Ia uncertainties, and to the individual systematic covariance matrices this method comprises.

Next, we attempt to replicate the quadrature summation of terms in R11, table 5 (most of which unfortunately do not have equivalent terms in our analysis) using projected uncertainties from our global fit. It is important to note that this comparison is not directly equivalent, because we are marginalising simultaneously over all nuisance parameters. With this caveat, we break down the uncertainty in the overall value of $H_{0}$ into three components: the uncertainties in the distance anchor, in the calibration of the SN Ia absolute magnitude $M_{B}$ via Cepheids, and in the measurement of the local expansion rate via SNe Ia (given in the intercept $\mathcal{M}$).212121This decomposition essentially follows equation 4 of R11. These can be determined separately from three disjoint data sets, as follows. The anchor distance is constrained by external data: the megamaser distance to NGC 4258 has a 0.066 mag uncertainty, corresponding to 3.0% in $H_{0}$. Only the low-$z$ SNe Ia are used to constrain $\mathcal{M}$ (or $5a_{V}$), with a 0.036 mag or 1.7% uncertainty. The calibration transfer from the Cepheids to the SNe Ia occurs in the simultaneous fit of the nearby supernova and Cepheid data222222For the uncertainty in $M_{B}$ to be independent of the error in $\mathcal{M}$, only these data can be included. to Equations 16 and 17. The resultant uncertainty in $M_{B}$ is 0.103 mag (with only the NGC 4258 anchor) and incorporates both the uncertainty in the SN Ia-Cepheid calibration and the uncertainty in the distance anchor; thus the former is 0.079 mag or 3.6% in $H_{0}$.232323The same calculation with all three anchors results in the same number. In the setup of R11, equation 4, this SN Ia-Cepheid calibration uncertainty is the error in $m^{0}_{v,4258}$, which is equivalent to $M_{B}+\mu_{4258}$. In quadrature, these three terms sum to 5.0% in $H_{0}$ using the NGC 4258 anchor, and 4.3% with all three anchors. This decomposition, whilst approximate, indicates that a quadrature sum of uncertainties in independent components results in similar uncertainties to our simultaneous fit. Thus, the simultaneous fit does not by itself result in the increase in statistical error.

6.3.3 Relative size of SN Ia uncertainties

We now examine the breakdown of uncertainties contributing to the statistical error, which include the multiple statistical and systematic uncertainties in SN Ia parameters making up $\mathbf{C_{\eta}}$ as constructed in Appendix C.

To visually assess the impact on confidence contours we compare results from MultiNest with different covariance matrix inputs. For an example global fit (with Cepheid fit $T=2.5$, NGC 4258 anchor, no priors, default SN cuts) we test each systematic, and compare their results from MultiNest. The full expression for the covariance matrix $\mathbf{C_{\eta}}$ for observed SN Ia quantities is described in Appendix C.1. As entries of $\mathbf{C_{\eta}}$ in MultiNest we try the following: only the statistical contribution $\mathbf{C_{\rm{stat}}}$ (described in Appendix C.2), each single systematic term added to $\mathbf{C_{\rm{stat}}}$, and all systematics added i.e. $\mathbf{C_{\rm{stat}}}+\mathbf{C_{\rm{sys}}}$ (the default for all global fits). The confidence contours with statistical uncertainties only and with all systematics are easily distinguishable in Fig. 9, but the contours with individual systematics are not. Thus for clarity we only show in Fig. 9 the systematic term from the uncertainty in host mass correction corrections ($\mathbf{C_{\rm{host}}}$ in Equation 40, described in Appendix C.5), in addition to $\mathbf{C_{\rm{stat}}}$ and $\mathbf{C_{\rm{stat}}}+\mathbf{C_{\rm{sys}}}$. The difference between the contours is slight, indicating that the uncertainties in the parameters only increase slightly when covariance matrices for different systematics are added to the statistical term $\mathbf{C_{\rm{stat}}}$.

Following the method in JLA (Betoule et al., 2014, section 6.2), we quantify the relative contributions, replacing the parameters $\{\Omega_{M},w,\alpha,\beta,M_{B},\Delta M\}$ with our parameters $\{H_{0},M_{B},\alpha,\beta\}$ (the only parameters in $\Theta$ which can be influenced by the low-$z$ SN Ia covariance matrices). The breakdown of relative contributions to the variance in $H_{0}$ from each term (the purely statistical term $\mathbf{C_{\rm{stat}}}$, and each systematic) are reported in Table 9. We emphasize that each of these numbers represents a proportion of the uncertainty (terms (i) and (ii) in Section 6.3 combined) from the systematic or statistical term in question alone, rather than reflecting an uncertainty in $H_{0}$.

From Table 9 and Fig. 9 it is clear that $\mathbf{C_{\rm{stat}}}$ is the largest component of $\mathbf{C_{\eta}}$. Even though the contributions to $\mathbf{C_{\eta}}$ from SN systematics are included in the statistical uncertainty, all of these systematics together are smaller than the SN statistical uncertainties: the relative contributions to the variance in $H_{0}$ are dominated by $\mathbf{C_{\rm{stat,diag}}}$ (Table 9), and the contours with and without systematic covariance matrices added to $\mathbf{C_{\rm{stat}}}$ in Figure 9 are similar to each other. Of the systematic terms (Appendices C.3 –C.6), the most significant in their impact on $H_{0}$ are from uncertainties in the Malmquist bias correction (including the selection function) and the host mass correction (Appendices B.2 and C.5), followed by the uncertainty in lightcurve model. While JLA had found the photometric calibration (Appendix C.4), especially from low-$z$ SNe, to be the dominant uncertainty for $\Omega_{M}$ and $w$, its effect on $H_{0}$ is decidedly smaller. It is interesting to note that despite conservative estimates of both the uncertainties in Milky Way extinction and peculiar velocity correction (Appendices C.3 and C.6), their effects on the error in $H_{0}$ are negligible.

7 Conclusions and future work

Our independent analysis of the Riess et al. (2011) data complements the R11 and E14 analyses in understanding the local measurement of the Hubble constant from type Ia supernovae. This work occupies the unique position of combining the precise Cepheid calibration of nearby SNe Ia (R11) with the sophisticated, thorough treatment of supernova lightcurves and systematics within a SALT2 framework (Betoule et al., 2014). In the context of the present tension in $H_{0}$ we present the first blinded SN Ia-based determination of $H_{0}$, eliminating confirmation and other biases. This work is intentionally applied to the well understood historical work of R11 and E14, as a proof of concept. It is our intent to extend this analysis to the sample of Riess et al. (2016).

Our best estimate from R11 data is $H_{0}=72.5\pm 3.1(\rm{stat})\pm 0.77(\rm{sys})$ km s*-1* Mpc*-1* using a three-galaxy (NGC 4258+LMC+MW) anchor. The central value is in excellent agreement with both R11 (after correcting for the lower value of $\mu_{4258}$ adopted – see Footnote 11) and the E14 reanalysis. Our above value and uncertainty imply tension with the Planck value at ${\sim}1.5\sigma$ significance, which is smaller than found in previous analyses of the R11 data, due to our increased uncertainties. However, our blinded affirmation of the central value in R11 signifies that the discrepancy between SN Ia- and CMB-derived values of the Hubble constant continues to exist. While this discrepancy is less significant in our analysis than in the original analysis of R11 data, it has potential to be magnified by the improved data set in R16 (which has smaller statistical uncertainties compared to R11), and hence remains of interest. It is thus necessary to apply the techniques in this paper to the R16 data, in order to make a contemporary assessment of the significance of the tension in the Hubble constant.

Incidentally, we find a higher stretch coefficient $\alpha=0.165\pm 0.010$ for our low-$z$ supernovae compared to LOSS (which find $\alpha=0.147\pm 0.007$ over a similar redshift range) and JLA. This discrepancy at ${\sim}2\sigma$ is surprising, and prompts further investigation. While our SN Ia zero point $M_{B}=-18.94\pm 0.09$ appears higher than in JLA, this difference arises because the parameters are degenerate using supernovae only: the JLA analysis assumes $H_{0}=70$ km s*-1* Mpc*-1*, whereas we have used a distance ladder to constrain both $H_{0}$ and $M_{B}$.

We find a larger relative uncertainty in $H_{0}$ (4.4% total) compared to R11 and E14 analyses of the same data (3.3% and 3.4% total respectively), based on the NGC 4258+LMC+MW distance scale. The difference appears in the statistical error; our systematic term is similar to that in R11 (1.1% and 1.0% respectively), with the caveat that the separation of our total uncertainty into statistical and systematic components is not directly comparable to R11, as described in Section 6.3.2. Our larger error primarily arises from our use of covariance matrices to estimate SN Ia systematic uncertainties. Other significant differences in our analysis, which potentially contribute to the increased uncertainty, are our simultaneous fit of all three sets of data, allowing all parameters to interact, and our use of variation in results to quantify systematic error. These distinctions are in our view justified and desirable. Given the increase in uncertainty they produce compared to previous works, they are important to consider in future analyses.

As found in R11, our results are limited by statistics in the supernova samples. Steps to reduce this statistical uncertainty have been implemented in R16, namely increasing the number of nearby galaxies to 18 and improving the SN Ia photometry, to reduce the total uncertainty to 2.4%. We envisage a similar relative increase in precision when the techniques in this work are applied to the same data set. R16 also includes important changes to data analysis of the Cepheids. Other contributions to our error budget are the systematic uncertainty, which is dominated by the variation in the different Cepheid fits, and the SN Ia systematic terms in $\mathbf{C_{\eta}}$, the largest of which are $\mathbf{C_{\rm{bias}}}$ and $\mathbf{C_{\rm{host}}}$.

Foremost, we find that both the use of covariance matrices and the simultaneous fit of data from different rungs of a distance ladder will be important in future analyses in order to wholly account for uncertainties. Furthermore, our findings recommend more sophisticated techniques for quantifying host galaxy dependence of SN Ia magnitudes, and modelling of Malmquist bias – both of which have the potential to diminish the systematic error in $H_{0}$. These techniques are continually improved in supernova analyses, particularly in the pursuit of more precise measurements of dark energy, for example in the Dark Energy Survey (DES; Dark Energy Survey Collaboration et al., 2016). Meanwhile a uniform, non-targeted low-$z$ sample (e.g. the SkyMapper Transient Survey (Scalzo et al., 2017), or the Pan-STARRS1 Survey (Rest et al., 2014)) will simplify photometric calibration and the selection function, reducing associated uncertainties, and will avoid peculiar velocity biases from coherent flows. Adopting these changes will benefit future SN Ia-based $H_{0}$ measurements.

Acknowledgements

This research was supported by the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and ARC Laureate Fellowship Grant FL0992131. We also acknowledge support of a University of Southampton Diamond Jubilee International Visiting Fellowship. BZ thanks the Australian Astronomical Observatory and the John Shaw Foundation for their financial support.

We are grateful to the SNLS collaboration for providing us with the snpca code and the recalibrated SALT2 training surfaces, used in Sections C.2 and C.4 respectively. We thank Matthew Colless and Christina Magoulas for suggesting the 2M++ peculiar velocity model, Dan Scolnic for helpful discussions, and Anais Möller for useful comments. This manuscript has greatly benefitted from thorough and constructive feedback from Adam Riess, to whom we are very grateful. Finally, we express our sincere thanks to the anonymous referee who has made particularly helpful suggestions and comments, which have substantially improved this work. Furthermore, we acknowledge their support and participation in the unique blinding process which contributed considerably both to the paper and their workload.

Appendix A Dependence of Cepheid-only fit

A.1 Outlier rejection

We perform fits in two ways: either assuming the outlier rejection in Riess et al. (2011), or following the rejection method in Efstathiou (2014). The R11 algorithm rejects Cepheids from each galaxy (rather than from the global fit), based on their deviation from the best Leavitt law fit. This rejection does not take into account the size of the Cepheid uncertainties, so that points with small residuals but large uncertainties are selectively accepted (E14, section 3.1). Consequently a large fraction of the total number of Cepheids is rejected, including a set of subluminous low-metallicity Cepheids (later corrected in R16, as discussed in Appendix A.4). Moreover, the intrinsic scatter is overestimated, resulting in a low reduced $\chi^{2}$.

These limitations in the R11 rejection are addressed in the upgraded algorithm in E14, which rejects a Cepheid from the global fit if its magnitude deviates from the global fit by more than the quantity $T\sqrt{m_{W,\mathrm{err}}^{2}+\sigma_{\mathrm{int,C}}^{2}}$ for a threshold $T$ (set to 2.25 or 2.5), where $m_{W,\mathrm{err}}$ and $\sigma_{\mathrm{int,C}}$ are the uncertainty in the individual Cepheid’s measurement and the intrinsic scatter $\sigma_{\mathrm{int,C}}$ respectively. This process is iterative, with $\sigma_{\mathrm{int,C}}$ recalculated at each step (such that $\chi^{2}_{c}$ per degree of freedom ${\sim}1$) with increments of 0.1, where the sum in $\chi^{2}_{c}$ is always over only the Cepheids in NGC 4258 and SN hosts (i.e. not the LMC or MW). The rejection at each iteration is based on the best fit determined in the previous iteration, i.e. the mean and 1$\sigma$ uncertainty of the posterior distribution.

Initially $\sigma_{\mathrm{int,C}}$ is set to 0.3. Then in each iteration we perform the following steps:

perform a MultiNest fit to all remaining Cepheids, to find marginalised posterior distributions; 2. 2.

find and remove outliers based on these parameters; 3. 3.

compute the new $\sigma_{\mathrm{int,C}}$ for these parameters and the updated Cepheid sample.

These steps are repeated until convergence, i.e. until no Cepheids are rejected in the second step.

The variation in fit results from different outlier rejection is presented in Fig. 1 and Section 4.4. In general the R11 rejection results in less negative values of both $b_{W}$ and $Z_{W}$, attributable to the aforementioned subluminous and low-metallicity subsample that it rejects.

The fit is forced to be good for all three rejection algorithms: $\sigma_{\mathrm{int,C}}$ is engineered to result in $\chi^{2}/\rm{DoF}{\sim}1$. Thus the algorithms cannot be compared statistically; the outlier rejection method has the drawback of not allowing the uncertainty on $\sigma_{\mathrm{int,C}}$ to be estimated, and the related consequence that we (by construction) cannot assess goodness-of-fit. Alternative statistical methods used in recent SN Ia analysis can surpass these limitations, including Bayesian hierarchical models (March et al., 2011; Shariff et al., 2016), the alternate Bayesian framework in Rubin et al. (2015), and Approximate Bayesian Computation Jennings et al. (2016). Notably, these have been applied to determining $H_{0}$ from the R11 and R16 data sets in Cardona et al. (2017).

A.2 Distance anchors

Our equations in Section 3 assume NGC 4258 is the only distance calibrator. We can generalise these equations to allow for combinations of the three anchor galaxies in R11, adding Cepheids in the LMC and MW (data described in Section 4.1). As these additional Cepheids do not have metallicity measurements, we adopt the mean values from E14 of $12+\log_{10}[O/H]$ of 8.5 and 8.9 for LMC and MW Cepheids respectively. Here, we test the dependence of the Cepheid parameters on the distance anchor. For Cepheids in the LMC and MW the Leavitt law (Equation 13) takes the forms:

[TABLE]

We consider combinations of NGC 4258, LMC, and MW (seven in total) as distance calibrators. If NGC 4258 is not included, then no prior for $\mu_{4258}$ is imposed in MultiNest. However, the likelihood $\mathcal{L}$ in Equation 10 still depends on $\mu_{4258}$, which is indirectly constrained through the other anchors and $M_{W}$, and hence remains a fit parameter in $\Theta$. If the LMC is used as an anchor then it is necessary to include $\mu_{\rm{LMC}}$ as a parameter in $\Theta$; this always has a (Gaussian) prior set to reflect the Pietrzyński et al. (2013) measurement from eclipsing binaries of $\mu_{\rm{LMC}}=18.494\pm 0.049$. The likelihood $\mathcal{L}$ is affected as a term $\chi^{2}_{\rm{LMC}}$ (Equation 31) is added to Equation 20; similarly if the MW is used as an anchor then $\chi^{2}_{\rm{MW}}$ (Equation 32) is added. We assume $\sigma_{\mathrm{int,C}}=0.113$ and $0.1$ for the LMC and MW respectively following E14.

[TABLE]

A modification to the above is necessary if both the LMC and MW are used as distance anchors, to account for the calibration uncertainty between ground-based and HST photometry. We do this using the covariance matrix $\mathbf{C}_{\mathrm{LMC+MW}ij}=(\hat{m}_{Wi}^{2}+\sigma_{\mathrm{int,C}}^{2})\delta_{ij}+\sigma_{\mathrm{cal}}^{2}$ with $\sigma_{\mathrm{cal}}=0.04$ (R11). Instead of $\chi^{2}_{\rm{LMC}}+\chi^{2}_{\rm{MW}}$, we add the term

[TABLE]

to the $\chi^{2}_{c}$ term going into $\mathcal{L}$. Here, bolded quantities represent vectors over all LMC and MW Cepheids.

The results of varying the distance anchor are discussed in Section 4.4. Briefly, the inclusion of both the LMC and MW anchors constrains both $b_{W}$ and $Z_{W}$ more tightly.

A.3 Longer-period cepheids

The data include Cepheids with period ranging from ${\sim}10$ days, to over 200 days, and Cepheids of all periods are included in Leavitt law fits in R11 (except for those Cepheids marked ‘low P’ in R11, table 2). Bird et al. (2009) examine longer-period ($P>80$ day) Cepheids and find that these Cepheids obey a flatter Leavitt law, with a less negative period dependence $b_{W}$. Accordingly, recent studies of the Leavitt law (e.g. Freedman et al., 2011; Scowcroft et al., 2011) have excluded Cepheids with period greater than 60 days. Similarly, E14 in their reanalysis of the R11 data have imposed the same upper limit on Cepheid period because of the observed change in slope. It is pragmatic to follow these examples in only using Cepheids over a period range where the slope remains constant; however, it is also useful the full range of periods to accommodate the change in slope and for the sake of comparison with R11. Rather than making an argument to include the $P<60$ day cut or not, we perform fits with and without an upper limit on the period.

E14 reasons that while including longer-period Cepheids decreases the magnitude of the Leavitt law slope $b_{W}$, there is little impact on $H_{0}$ (Efstathiou, 2014, appendix A), so they only include $P<60$ day Cepheids in their fits. Our priors on $b_{W}$ differ slightly from E14 (discussed in Appendix A.4), and we are interested in the variation of Cepheid parameters with the choice of period cut (as with distance anchor and rejection algorithm in previous sections), so we test the dependence of fit results on the inclusion of an upper limit on period. Results of including longer-period Cepheids are lesser dependence on Cepheid slope and metallicity dependence (less negative $b_{W}$ and $Z_{W}$), as described in Section 4.4.

A.4 Slope and metallicity priors

We test and discuss the Gaussian priors on $b_{W}$ and $Z_{W}$ described in E14 (but not mentioned in R11), and explain our choices for our fits. E14 performs Cepheid fits with and without Gaussian priors centred at $b_{W}=-3.23$ and $Z_{W}=0$, motivated by expectations of what the slope and metallicity dependence should be. We test the same priors in our fits but ultimately decide to not include these different priors as one of the variables in our fit, for reasons which follow.

Out of all the Cepheid data, the LMC Cepheids constrain the slope $b_{W}$ most tightly. Given the relative paucity of data on the Leavitt law, we always include this information on the Leavitt law in all our fits, independent of whether the LMC is used as a distance anchor. For the fits where the LMC is not included as an anchor, we impose the same Gaussian prior on the slope as in E14: $\langle b_{W}\rangle=-3.23$, $\sigma_{b_{W}}=0.10$. If the LMC is used as an anchor, there is already a contribution to the likelihood from these Cepheids, so it is inappropriate to reuse this information as a prior. Then, the inclusion of the prior on $b_{W}$ is prescribed by the distance anchor.

The metallicity priors in E14 are motivated by the observed strong dependence of the Cepheids’ period on metallicity, in contrast with expectations that $Z_{W}{\sim}0$, based on theoretical considerations and measurements in the LMC (Freedman & Madore, 2011; Efstathiou, 2014, section 3.2, and references therein). However, the R11 sample of Cepheids demonstrates a strong metallicity dependence, with values of $Z_{W}$ around -0.3 or -0.5 for the R11 and E14 rejection algorithms respectively. The difference between values for $Z_{W}$ from the two approaches to outlier rejection (detailed in Appendix A.1) can be traced to a sample of low metallicity Cepheids that are rejected by cuts in R11 but not E14. This systematic difference (discussed in E14, section 3.2) arose from the erroneous extrapolation of metallicity gradients to large radii, and was later corrected in R16. Including both the LMC and MW as distance anchors reduces the magnitude of the metallicity dependence $Z_{W}$. As we have observed that the R11 Cepheid data do not support the $Z_{W}{\sim}0$ priors (weak or strong) in E14, it is most appropriate to exclude these Gaussian priors in our analysis.

Appendix B Supernova lightcurves and data

B.1 SALT2 lightcurve fits

For each supernova we use SALT2 to fit SN Ia lightcurves (i.e. determine parameters $m_{B},X_{1},C$ in Equation 14). The SALT2 model, based on its precursor SALT (Guy et al., 2005), is described in Guy et al. (2007) along with details of its training. Two newer versions, SALT2.2 and SALT2.4, with additional training samples, have been released with the SNLS and JLA analyses, respectively in Guy et al. (2010) and Betoule et al. (2014, hereafter B14). Notably these include supernovae from the SDSS-II survey (Sako et al., 2014) and high-$z$ SNe which have constrained the model better in the rest-frame ultraviolet region, eliminating the need to exclude U-band data as in R11.

Primary inputs for the SALT2 lightcurve-fitting routine snfit are photometry in each filter, heliocentric redshifts, and Milky Way extinction (obtained from Schlegel et al. (1998) dust maps) for each supernova. In addition, zero points and filter transmissions for each passband of each instrument are necessary. As CfA3 supernovae are included in the JLA analysis and the SALT2.4 release, we only need to create LOSS instruments KAIT1–4 and NICKEL. Adding these instruments involves including the filter transmissions provided by the Berkeley group242424http://heracles.astro.berkeley.edu/sndb/info, and determining zero points for BD+17∘4708252525This is the fundamental SDSS standard star that the Vega-based magnitude system @VEGA2 in SALT2.4 is calibrated on. using colour transformations in Ganeshalingam et al. (2010, table 4). Our photometry and the instruments used are briefly mentioned in Sections 5.1 and C.4; for more details see the survey papers Hicken et al. (2009a); Ganeshalingam et al. (2010).

B.1.1 Photometry consistency checks

We are able to compare our SALT2 fit results to published values for a subsample of our SNe Ia, namely some of the CfA3 SNe which appear in SNLS (Conley et al., 2011, hereafter C11) and JLA (B14). Hicken et al. (2009b) and Ganeshalingam et al. (2013) also report SALT2 fits of their samples, albeit with an older version of SALT2. We thoroughly check for consistency between these results and find agreement within quoted uncertainties, with no systematic differences.

Furthermore we examine both the photometry and lightcurve fits for the 69 SNe in the CfA3-LOSS overlap, taking into account the different magnitude systems. For the directly comparable passbands (the Bessell-like BVRI filters) the photometry is consistent, while differences from the SDSS ri filters in Keplercam are in line with expectations. We also compare results of lightcurve fits using only CfA3 photometry, only LOSS photometry, and both combined in a single lightcurve. We find that SALT2 parameters $\{\alpha,\beta\}$, determined from each survey separately, generally average to the parameters determined from the combined fit (which lie well within reported uncertainties of from either CfA3 or LOSS). Occasionally one survey dominates in its influence over the SALT2 parameters; this occurs equally often with each survey and only when the lightcurve quality is discernibly superior in terms of number of points, sampling frequency, and size of uncertainties. Furthermore we test for systematic offsets in the three SALT2 fit parameters $m_{B},X_{1},C$, and find none. The comparison of combined lightcurves and their SALT2 fits supports the consistency of CfA3 and LOSS and favours neither over the other; thus we concatenate photometry from CfA3 and LOSS instruments to obtain the highest quality lightcurves available for these SNe.

B.2 Malmquist bias correction

In magnitude limited surveys, intrinsically brighter objects are preferentially detected, leading to Malmquist bias: a skewed estimate of the absolute magnitude distribution. The Malmquist bias can be estimated by modelling the selection efficiency (i.e. the rate of successful detection as a function of magnitude) to match observed distributions (in redshift, stretch, and colour), then simulating the survey with SNANA (Kessler et al., 2009) to obtain the bias $\delta\mu$ in distance modulus. This procedure, described in e.g. Scolnic & Kessler (2016), is outside the scope of this work. Thus we adopt the estimate of the bias (for low-$z$ supernovae) in Betoule et al. (2014, section 5.3), which adopts a magnitude limited selection function, and uses the difference between the resultant bias and an unbiased regime as the uncertainty in the correction (the covariance matrix $\mathbf{C_{\rm{bias}}}$ in Appendix C.1). The targeted discovery of supernovae in CfA3 and LOSS means they should not be magnitude limited; however, as observed in the JLA low-$z$ sample, the colour distribution grows more blue with redshift, suggesting that some selection effect is at play. Using the JLA approximation is justified as our supernova sample is similarly distributed to the low-$z$ sample in JLA, and the bias correction is inherently approximate and has a miniscule impact on $H_{0}$.262626The difference between correcting for Malmquist bias and no correction is less than 0.01 in $\mathcal{H}$, or ${\sim}0.4\%$ in $H_{0}$.

B.3 Host galaxy dependence

The dependence of the intrinsic SN Ia brightness on properties of their host galaxies is well established, with numerous studies finding that more massive galaxies (correlated with higher metallicity and lower specific star formation rates) host SNe which are on average ${\sim}0.08$ mag brighter (e.g. Sullivan et al., 2010; Lampeitl et al., 2010; Kelly et al., 2010). To mitigate the systematic error that this effect introduces to the cosmological analysis, we follow Sullivan et al. (2011) and subsequent analyses (C11, B14) in adopting two discrete values for the SN Ia absolute magnitude, using the variable

[TABLE]

We fit for the parameter $M_{B}$ as indicated in Section 3, and fix the offset $\Delta M_{B}$ based on analyses in C11; B14; Sullivan et al. (2010); Childress et al. (2013), which determine $\Delta M_{B}=-0.08$ from SN Ia samples greater in size and redshift range than ours. We consider fitting for the magnitude offset using our data, and find a larger absolute difference $\Delta M_{B}=-0.15\pm 0.07$ (with some degeneracy with $\mathcal{M}$), which is still consistent with the established value in the literature. Given the large uncertainty on our value we adopt the more reliable reference value.

The host galaxy masses for our SNe are obtained from the literature where available, with 77 from JLA and 71 from a combination of Sako et al. (2014); Childress et al. (2013); Neill et al. (2009); Kelly et al. (2010). The stellar masses of nearby galaxies are all given in Neill et al. (2009). We were able to derive mass estimates for 72 of the remaining galaxies using SDSS photometry, following standard methods in Sullivan et al. (2006); Smith et al. (2012). We refer the reader to descriptions therein of the method, which relies on the ZPEG photometric redshift code (Le Borgne & Rocca-Volmerange, 2002, 2010) based on spectral energy densities from the PÉGASE.2 spectral synthesis code. Where possible, we check for consistency between multiple sources and our estimates. The distribution of the host masses of the CfA3 and LOSS SNe Ia (Fig. 10) clearly shows that they predominantly exist in more massive galaxies, with 241 out of 280 SNe lying in the high mass bin. This is a consequence of the targeted nature of these surveys, in contrast to the magnitude-limited surveys SNLS and SDSS in JLA. Thus we assign the remaining 60 SNe with unknown masses to the high mass bin with a large associated uncertainty, even though unknown hosts in JLA are assigned to the low mass bin (B14, section 5.2). The propagation of uncertainties in this correction through to SN parameters is later described in Appendix C.5.

B.4 Peculiar velocity corrections

Peculiar velocities arise from motion other than from cosmological expansion, such as dipole or bulk flows, local galaxy infall, and higher order coherent flows. These perturb the observed redshifts via the Doppler effect,272727A supernova’s peculiar motion changes not only its redshift but also its observed luminosity (Hui & Greene, 2006; Davis et al., 2011) as it experiences relativistic beaming. This in turn induces a deviation in the supernova’s peak magnitude; however this is approximately an order of magnitude smaller than the change in redshift. and can impact cosmological analyses. Hui & Greene (2006) show that neglecting correlations between peculiar velocity uncertainties at low-redshift results in a greatly underestimated zero point uncertainty, and degrades the precision of the dark energy equation of state parameter $w$. Moreover, correlated SN peculiar velocities can bias cosmological results: Davis et al. (2011) show that neglecting coherent flows results in a shift of $\Delta w=0.02$.

Thus, an effort to quantify the uncertainty induced by correlated motions is an essential part of any modern SN Ia cosmological analysis. Approaches to this include the addition of large ($300-400$ km s*-1*) uncertainties in redshifts to account for peculiar velocities (Ganeshalingam et al., 2013; Hicken et al., 2009b), and attempts to correct for peculiar velocities. The latter first appeared in SNLS (C11), which corrects redshifts on a supernova by supernova basis for the (line-of-sight) peculiar velocity at the location of the SN, as determined from a velocity field.

C11 uses the velocity field in Hudson et al. (2004), derived from the galaxy density field from the IRAS PSCz redshift survey (Branchini et al., 1999). While B14 adopts the same correction, we apply the same method using updated density and velocity fields from the 2M++ redshift compilation.282828Data available at http://cosmicflows.iap.fr/. In each case the velocity field is derived from the respective density fields under a linear biasing approximation.292929That is, the mass density and galaxy number density are proportional via the linear bias factor $b$, i.e. $\delta_{g}=b\delta$. In this regime, peculiar velocities are proportional to gravitational attraction:

$\displaystyle\boldsymbol{v}=\frac{\beta^{*}}{4\pi}\int_{0}^{R_{\rm{max}}}\delta_{g}(\boldsymbol{r}^{\prime})\frac{(\boldsymbol{r}^{\prime}-\boldsymbol{r})}{|\boldsymbol{r}^{\prime}-\boldsymbol{r}|^{3}}\ d^{3}\boldsymbol{r}^{\prime}+\boldsymbol{U},$

(35)

where $\boldsymbol{U}$ represents a residual dipole (in 2M++ this is the dipole of the Local Group), with $\beta^{*}=\frac{f(\Omega_{M})}{b}=0.43$ (Carrick et al., 2015) where $f(\Omega_{M})=\Omega_{M}^{0.55}$ for $\Lambda$CDM (Wang & Steinhardt, 1998).

In correcting supernova redshifts for peculiar velocities the aim is to isolate the redshift due purely to expansion. The several redshifts of interest are: the heliocentric redshift $z_{h}$, the CMB frame redshift $z_{\rm{cmb}}$, and the cosmological redshift $\bar{z}$. The latter two differ in that $\bar{z}$ is corrected for peculiar velocities from coherent flows; it is intended to reflect a velocity derived only from the expansion of space and therefore this is the redshift that should be used in $v(z)$ in Hubble’s law. The peculiar motions to consider in converting $z_{h}$ to $\bar{z}$ are the motions of the solar system ($\boldsymbol{v}_{\odot}^{\rm{pec}}$), and of the SN ($\boldsymbol{v}_{\rm{SN}}^{\rm{pec}}$), both relative to the CMB. Many SNe at low redshifts share some of the Local Group’s motion; by converting to a heliocentric frame (i.e. correcting for the Sun’s motion relative to the CMB) we are also overcorrecting for the motion of nearby SNe, necessitating the second correction. For a SN at position $\boldsymbol{n}$ from the Sun, these redshifts and velocity are related by303030The minus sign in front of $\boldsymbol{v}_{\odot}^{\rm{pec}}$ arises because we have defined it is the motion of the Sun relative to the CMB, rather than the other way around. (Davis et al., 2011)

[TABLE]

For our analysis of the low-$z$ SNe Ia, we use $\bar{z}$ derived in this way as the CMB-frame redshift. Unless otherwise specified, this is the redshift meant by $z$.

The exact form of the luminosity distance introduced in Equation 2 actually requires both the heliocentric and cosmological redshifts:

[TABLE]

This is because the factors affecting the $(1+z_{h})$ prefactor (redshifting and beaming) depend on the total relative velocities, whereas the cosmological distance only depends on $\bar{z}$, the redshift due to expansion (Calcino & Davis, 2017, Appendix A). The difference resulting from using $\bar{z}$ for both is negligible so we do not differentiate in our analysis. Calcino & Davis (2017, Section 4.2–4.3) quantify the effect of possible redshift systematic errors on the derivation of $H_{0}$ and find that a systematic redshift error as small as ${\sim}2.6\times 10^{-4}$ can result in a ${\sim}0.3$% bias in $H_{0}$.

The peculiar velocity corrections we make here are reliant on predicted velocity fields, which are intrinsically approximate. We discuss and quantify these uncertainties in Appendix C.6, and propagate them to the SN fit parameters. Moreover, we ensure that our corrections do not bias our results: if the results of our SN-only fit (Section 5.4) varied significantly with the introduction of the correction, then this effect would need to be explored and quantified. In this scenario, the impact of performing a velocity correction would greater than the uncertainty contribution in $\mathbf{C_{\rm{pecvel}}}$ (Appendix C.6). However, we find a negligible effect on the zero point $\mathcal{M}$ (less than 10% of the statistical uncertainty) in the SN-only fit when peculiar velocity corrections are omitted. Consequently the velocity correction cannot bias $H_{0}$ (as $\mathcal{H}=M_{B}-\mathcal{M}$).

B.5 Histograms for SN cuts

We include in Figs 13, 13, and 13 the stretch and colour distributions of our low-$z$ SNe Ia with various alternate cuts relative to the default, all described in Section 5.3. We observe that the $X_{1}$ distribution changes marginally with stricter cuts on $\sigma_{X_{1}}$ and a higher low-redshift cut. Naturally, the $C$ distribution is affected by a stricter cut on $C$. Otherwise, there is no significant impact on the stretch and colour distributions from alternate cuts, and in particular no evidence that the variation observed in Fig. 2 in Section 5.5 arises from such biases.

Appendix C Computing supernova systematics

This section is devoted to the construction of the SN Ia covariance matrices. We break down systematic and statistical terms that contribute to the error budget, and describe their propagation to the supernova magnitudes.

C.1 Overview of covariance matrices

Each covariance matrix tracks uncertainties in the vector $\eta=\{m_{Bi},X_{1i},C_{i}\}_{1\leq i\leq N}$, which contain the SALT2 quantities for all $N=280$ low-$z$ SNe. These matrices sum to the $3N\times 3N$ matrix $\mathbf{C_{\eta}}$ which encompasses all covariances in $\eta$, and are independent of $\alpha$ and $\beta$. For fitting the low-$z$ SN data we require covariances in supernova magnitudes in the form of $\mathbf{C_{m_{B}^{\dagger}}}$ (Equation 22), which is derived from $\mathbf{C_{\eta}}$ by conjugation with the $N\times 3N$ matrix $\mathbf{A}$ (with $\mathbf{A}_{ij}=\delta_{3i,j}+\alpha\delta_{3i+1,j}+\beta\delta_{3i+2,j}$):

[TABLE]

The remaining terms are diagonal uncertainties, which affect each SN individually, ascribed to uncorrelated uncertainties in redshift due to peculiar velocity uncertainties (distinct from the uncertainty in their corrections, described in Appendix C.6), and perturbances in SN Ia magnitudes caused by gravitational lensing and intrinsic scatter. We adopt the values for these used in C11 and B14, of $c\sigma_{z}=$ 150 km s*-1*, $\sigma_{\rm{lens}}=0.055z$, and $\sigma_{\mathrm{int,SN}}=0.12$.

To understand $\mathbf{C_{\eta}}$, we first separate it into statistical and systematic components, and later explain the distinction in Appendix C.2. The contributions to $\mathbf{C_{\rm{sys}}}$ we consider are from uncertainties in the following sources: peculiar velocity corrections, Milky Way extinction, host galaxy mass dependence, photometric calibration, Malmquist bias correction, and lightcurve model.313131This is equation 11 of B14 without the contamination term $\mathbf{C_{\rm{nonIa}}}$, which only concerns higher redshift SNe. For the nearby SNe (i.e. those in galaxies containing Cepheids), we only include $\mathbf{C_{\rm{diag}}}$ and $\mathbf{C_{\rm{host}}}$. The host mass correction has the potential to shift the magnitude scale by up to ${\sim}0.08$ mag, and is important in the context of the dependence of SN Ia magnitude on host galaxy properties (Appendix C.5). The other correction terms, for Malmquist bias and peculiar velocities, are redshift-dependent effects hence irrelevant for this sample. The remaining covariance matrices are not tied to the corrections in Appendices B.2–B.4, and are more precise than warranted, given the inhomogeneity and larger uncertainties in these data, so we neglect them.

[TABLE]

We follow standard techniques to compute each covariance matrix, which is to enfold partial derivatives of SN parameters with respect to each systematic, with the typical size of systematics:

[TABLE]

Here the sum is over all systematics $k$, each of size $\Delta k$. Equation 41 is applied directly to compute $\mathbf{C_{\rm{dust}}}$ and $\mathbf{C_{\rm{cal}}}$. These calculations are intrinsically approximate, yet even as estimates they are invaluable for gauging the contribution of each systematic term affecting SNe Ia, and affirming that we sufficiently account for each effect. Section 6.3.3 presents our assessment of these uncertainties.

We first digress to discuss the statistical term $\mathbf{C_{\rm{stat}}}$, then address the construction of each systematic in turn. We describe our calculations of the first four systematics from first principles. Computations of the bias and model uncertainties, as well as the non-diagonal part of $\mathbf{C_{\rm{stat}}}$ require estimates of the sample’s selection function (as discussed in Section B.2), thorough end-to-end simulations with several lightcurve models, and in-depth deconstruction and analysis of the SALT2 model. These have been performed in section 5.3 of B14, Mosher et al. (2014), and Guy et al. (Appendix A3 2010) respectively. We obtain our best estimate of these matrices for our low-$z$ SN Ia sample, and refer the reader to the aforementioned references for details.

C.2 Statistical uncertainties

The distinction between statistical and systematic errors blurs, as many uncertainties have sources for which both descriptors are appropriate. We adopt the separation used in C11, which defines statistical uncertainties as those that can be reduced by increasing the size of some data set. In this case the data sets are the measured low-$z$ supernovae, and the training set used to define the SALT2 parameters (Guy et al. (2010), updated in B14). We separate these two terms into matrices $\mathbf{C_{\rm{stat,diag}}}$ and $\mathbf{C_{\rm{stat,model}}}$ respectively.323232In B14 these two terms are combined as $\mathbf{C_{\rm{stat}}}$, while C11 sums $\mathbf{C_{\rm{stat,diag}}}$ and the three diagonal terms in our Equation 38 to their $\mathbf{D_{\rm{stat}}}$.

The matrix $\mathbf{C_{\rm{stat,diag}}}$ arises from uncertainties in the measurement of lightcurves, constructed directly from correlated uncertainties in $m_{B},X_{1},C$ (a $3\times 3$ covariance matrix for each SN) reported in SALT2 outputs. These terms are uncorrelated between different supernovae, so $\mathbf{C_{\rm{stat,diag}}}$ is tridiagonal (i.e. only a diagonal strip of width 3 is nonzero).

The latter matrix $\mathbf{C_{\rm{stat,model}}}$ contains the uncertainty in the SALT2 model from the finiteness of the training sample, which could theoretically be decreased by training SALT2 on more supernovae. Its calculation requires propagating statistical uncertainties in the lightcurve model through to supernova fit parameters $\eta$, as described in Guy et al. (Appendix A3 2010) and implemented in the snpca package.333333Private SNLS communication. We use the code salt2_stat from this package to directly compute $\mathbf{C_{\rm{stat,model}}}$.

C.3 Milky Way extinction

The calculation of our first systematic covariance matrix directly follows Equation 41. This contains the uncertainty in $\eta$ due to the uncertainty in the Milky Way extinction. The sole systematic $k$ is the uncertainty in the $E(B-V)$ value from dust maps (Schlegel et al., 1998); we follow the conservative estimate in B14 of a 20% relative uncertainty. Perturbing the value of the extinction (encoded in the dust keyword @MWEBV in SALT2 inputs) and refitting lightcurves give the partial derivatives in Equation 42:

[TABLE]

We verify that the partial derivatives of SN parameters $\eta$ with respect to Milky Way extinction are independent of the size of perturbation over a satisfactory range, and that our resultant $\mathbf{C_{\rm{dust}}}$ is identical to the same matrix reported in B14 for the 60 SNe Ia in common.

C.4 Calibration

B14 and C11 emphasize the significant contribution of uncertainties in the calibration of individual surveys to the total error budget. We follow the methods therein and in Betoule et al. (2013) to reproduce the calibration covariance matrix relevant to our SN sample and the instruments used to observe them. Computing the calibration matrix $\mathbf{C_{\rm{cal}}}$ relies on the same principle as in Appendix C.3, but over multiple systematics. Calibration uncertainties are grouped into two types of systematics: uncertainties in the magnitude zero point (shifting the overall flux scale) and in the effective wavelength (shifting the transmission function in wavelength space), for each filter. These are contained in the vector $\kappa$ and enumerated in Table 10.

The instruments and passbands to consider in $\kappa$ are those used for observing the low-$z$ SNe: 4Shooter2 and Keplercam for CfA3, and KAIT1–4 and Nickel for LOSS, and those involved in the training of SALT2 (i.e. used to observe the SNe in the training sample). The latter, and sizes of systematics in these passbands, are given in B14, table 5. It is necessary to include these training instruments and passbands as they influence measured magnitudes of training SNe hence the SALT2 model.

We directly consider the covariance matrix of calibration systematics $\mathbf{C_{\kappa ij}}=\langle\sigma_{\kappa_{i}}\sigma_{\kappa_{j}}\rangle$, which captures the correlations between systematics in different instruments and passbands. Then Equation 41 is equivalent to $\mathbf{C_{\rm{cal}}}=\mathbf{J\cdot C_{\kappa}\cdot J}^{T}$ where $\mathbf{J}_{ij}=\frac{\partial\eta_{i}}{\partial\kappa_{j}}$ is the Jacobian matrix, encoding partial derivatives of SN parameters. Then finding $\mathbf{C_{\rm{cal}}}$ amounts to constructing $\mathbf{C_{\kappa}}$, and calculating $\mathbf{J}$ from first principles. For the LOSS instruments we achieve the latter by either perturbing an element of $\kappa$ (i.e. shift the zero point of effective wavelength). For the other instruments, which were involved in training SALT2, we change to a different SALT2 model altogether i.e. one that was trained with the systematic shift in question applied. Each SN lightcurve is fitted again to find the difference, and the resultant Jacobian is smoothed in accordance with footnote 9 of B14.

To find $\mathbf{C_{\kappa}}$ we start with the same matrix from JLA and reindex it according to Table 10, appending the LOSS instruments and removing HST instruments NICMOS and ACSWF (which do not contribute to the SALT2 training sample). We approximate the elements of $\mathbf{C_{\kappa}}$ for LOSS instruments as diagonal: this is exact for the $\lambda_{\rm{eff}}$ elements, and a good approximation for the zero point. Using Ganeshalingam et al. (2010) as a guide, we take the zero point and $\lambda_{\rm{eff}}$ uncertainties to be 0.03 mag and 10 Å respectively. As LOSS observations (with KAIT1–4 and Nickel) of SNe Ia were not used for SALT2 training, only SNe in the sample with LOSS measurements have nonzero partial derivatives with respect to these instruments.

C.5 Host galaxy mass

The SN Ia magnitude zero point $M_{B}$ is corrected for the magnitude offset $\Delta_{M}$ between high and low host galaxy (stellar) mass bins, as described in Appendix B.3. The uncertainty in this correction is propagated to SN parameters in $\mathbf{C_{\rm{host}}}$. As in B14, we treat the systematic associated with this correction as having two components: from potentially having attributed an individual supernova to the wrong host mass bin, and from the arbitrariness of the $10^{10}M_{\odot}$ cut. Both effects are discrete, so the computation of $\mathbf{C_{\rm{host}}}$ differs from those of $\mathbf{C_{\rm{dust}}}$ and $\mathbf{C_{\rm{cal}}}$ which take partial derivatives with respect to continuous quantities.

Our calculation follows identically the method in B14. As discussed in Appendix B.3, our data do not justify fitting for $\Delta_{M}$, and instead we adopt a fixed value from literature. Only supernova magnitudes are affected: the only indices of nonzero components of $\mathbf{C_{\rm{host}}}$ correspond to the $m_{B}$ components of $\eta$, but for compatibility for the other matrices we make $\mathbf{C_{\rm{host}}}$ the same size ($3N\times 3N$). The vectors $\boldsymbol{H_{h}}$ and $\boldsymbol{H_{l}}$ are indicator functions of (the magnitudes of) the SNe with masses within an order of magnitude higher and lower respectively than the mass boundary $10^{10}M_{\odot}$, while $\boldsymbol{B}$ is and indicator function for the SNe whose host mass estimates and uncertainties in combination imply that they could be be assigned to either bin. Then the covariance matrix for the host mass correction is:

[TABLE]

For the overlap of our SNe with B14, results are very similar, with differences arising only from SNe with masses newly obtained or updated (Appendix B.3).

C.6 Peculiar velocities

As described in Appendix B.4 we have corrected individual SN redshifts for peculiar motion, using the 2M++ velocity field corrections. However, there is intrinsic uncertainty in these models, with variation between velocity fields generated from different galaxy density fields, and in some cases limited agreement between predicted and measured velocities (Springob et al., 2014; Scrimgeour et al., 2016).

Thus the significant contribution in the correction model itself must be taken into account. Below, we adopt the approach in C11 and B14, which is to use the uncertainty in the velocity field to inform $\mathbf{C_{\rm{pecvel}}}$, the contribution to $\mathbf{C_{\eta}}$ from peculiar velocities. We emphasize that this is distinct from the diagonal term $\sigma_{z}$ in Equation 38.

For a given density field $\delta_{g}$, the velocity field derived through Equation 35 can be parametrized by $\beta^{*}$, the ratio of the growth rate of density perturbations to the linear bias factor. In C11, $\beta^{*}$ is the systematic which encompasses the uncertainty in the peculiar velocity model; that is, $\mathbf{C_{\rm{pecvel}}}$ is derived through Equation 41 with $k=\beta^{*}$. As this treatment of uncertainty lies within only one density field and model (that is, it doesn’t account fully for velocities derived from different realisations/measurements of galaxy densities) we are conservative in using it; like C11 we perturb $\beta^{*}$ by five times its uncertainty.343434Pike & Hudson (2005) find $\beta^{*}=0.49\pm 0.04$ so C11 vary $\beta^{*}$ between 0.3 and 0.7. Likewise we adopt $\beta^{*}=0.43\pm 0.02$ (Carrick et al., 2015) in the correction. To compute $\mathbf{C_{\rm{pecvel}}}$, we measure the shift in $z_{\rm{cmb}}$ when $\beta^{*}$ is set to $0.33$ or $0.53$ instead. The resultant difference in $z_{\rm{cmb}}$ is propagated to $m_{B}$ using the derivative of Equation 15:

[TABLE]

This has no impact on the stretch and colour of SNe, so only the $m_{B}$ elements of $\mathbf{C_{\eta}}$ have non-zero entries from $\mathbf{C_{\rm{pecvel}}}$.

Appendix D Full tables of fit results

This appendix supplements Sections 4 and 6 with full tables of the Cepheid and global fits, and the averaged results for $\{\Delta\mu_{i}\}$ from the global fit.

D.1 Results of Cepheid-only fit

D.2 Full results for $\{\Delta\mu_{i}\}$

D.3 Full results of global fit