Apparent cosmic acceleration from type Ia supernovae

TL;DR

This paper compares the timescape cosmology model to the standard Lambda-CDM model using supernova data, finding similar fits and suggesting the importance of inhomogeneous models in understanding cosmic acceleration.

Contribution

It demonstrates that the timescape model, an inhomogeneous cosmology without dark energy, fits supernova data as well as the standard model, highlighting the role of spatial curvature evolution.

Findings

Timescape model fits supernova data as well as Lambda-CDM.

Cosmic acceleration in the timescape model is marginal and non-zero.

Systematic effects influence supernova light curve parameters independently of cosmology.

Abstract

Parameters that quantify the acceleration of cosmic expansion are conventionally determined within the standard Friedmann-Lemaitre-Robertson-Walker (FLRW) model, which fixes spatial curvature to be homogeneous. Generic averages of Einstein's equations in inhomogeneous cosmology lead to models with non-rigidly evolving average spatial curvature, and different parametrizations of apparent cosmic acceleration. The timescape cosmology is a viable example of such a model without dark energy. Using the largest available supernova data set, the JLA catalogue, we find that the timescape model fits the luminosity distance-redshift data with a likelihood that is statistically indistinguishable from the standard spatially flat $\Lambda$ cold dark matter cosmology by Bayesian comparison. In the timescape case cosmic acceleration is non-zero but has a marginal amplitude, with best-fitting apparent deceleration parameter, $q_0=-0.043^{+0.004}_{-0.000}$. Systematic issues regarding standardization of supernova light curves are analysed. Cuts of data at the statistical homogeneity scale affect light curve parameter fits independent of cosmology. A cosmological model dependence of empirical changes to the mean colour parameter is also found. Irrespective of which model ultimately fits better, we argue that as a competitive model with a non-FLRW expansion history, the timescape model may prove a useful diagnostic tool for disentangling selection effects and astrophysical systematics from the underlying expansion history.

Tables

Table 1. Table 1: All nuisance parameters in each model have identical priors. In the timescape model Ω M ​ 0 subscript Ω 𝑀 0 \Omega_{\lower 1.0pt\hbox{$\scriptstyle M0$}} is defined in terms of f v0 subscript 𝑓 v0 f_{\rm v0} hence we take the latter to be the more ‘fundamental’ parameter and assign the prior to it.

Cosmological parameter	Prior distribution	Range
Timescape: $f_{\mathrm{v0}}$	Uniform	$[0.500,0.799]$
Flat $\Lambda$CDM: ${\Omega }_{M0}=1-{\Omega }_{\Lambda 0}$	Uniform	$[0.143,0.487]$
Nuisance parameters	Prior distribution	Range
$\alpha$	Uniform	$[0,1]$
$x_{1,0}$	Uniform	$[-20,20]$
${\sigma }_{x_{1,0}}$	Uniform on ${\log }_{10}{\sigma }_{x_{1,0}}$	$[-5,2]$
$\beta$	Uniform	$[0,4]$
$c_0$	Uniform	$[-20,20]$
${\sigma }_{c_0}$	Uniform on ${\log }_{10}{\sigma }_{c_0}$	$[-5,2]$
$M_{B,0}$	Uniform	$[-20.3,-18.3]$
${\sigma }_{M_{B,0}}$	Uniform on ${\log }_{10}{\sigma }_{M_{B,0}}$	$[-5,2]$
Additional stretch and colour	Uniform	$[-20,20]$
parameters for models II-VIII

Table 2. Table 2: Best-fitting MLE parameters corresponding to the likelihood ℒ ​ ( Data | M ) ℒ conditional Data M \mathcal{L}({\rm Data}|\mathrm{M}) with the model M M \mathrm{M} representing the cosmological model, the SALT2 procedure and the intrinsic distributions of SneIa parameters. SneIa at redshifts z < 0.033 𝑧 0.033 z<0.033 (statistical homogeneity scale) are excluded. Confidence limits are given for the one free cosmological parameter. In the timescape case this corresponds to f v0 = 0.778 − 0.056 ​ ( 1 ​ σ ) ​  0.104 ​ ( 2 ​ σ ) + 0.063 ​ ( 1 ​ σ ) ​  0.155 ​ ( 2 ​ σ ) subscript 𝑓 v0 subscript superscript 0.778 0.063 1 𝜎 0.155 2 𝜎 0.056 1 𝜎 0.104 2 𝜎 f_{\rm v0}=0.778^{+0.063\;(1\sigma)\;0.155\;(2\sigma)}_{-0.056\;(1\sigma)\;0.104\;(2\sigma)} . The value of M B , 0 subscript 𝑀 B 0 M_{\lower 1.0pt\hbox{$\scriptstyle\rm B,0$}} is obtained for h = 0.668 ℎ 0.668 h=0.668 for the timescape, and h = 0.7 ℎ 0.7 h=0.7 for the two FLRW models. The difference of parameters from NGS16 is principally due to the SHS cut at z min = 0.033 subscript 𝑧 min 0.033 z_{\rm min}=0.033 , the effect of which is seen in Fig.  2 .

Model	${\Omega }_{M0}$	$\alpha$	$x_{1,0}$	${\sigma }_{x_{1,0}}$	$\beta$	$c_0$	${\sigma }_{c_0}$	$M_{B,0}$	${\sigma }_{M_{B,0}}$
Timescape	${0.309}_{-0.088\left(1\sigma \right)\mathrm{\hspace{0.33em}0.210}\left(2\sigma \right)}^{+0.070\left(1\sigma \right)\mathrm{\hspace{0.33em}0.127}\left(2\sigma \right)}$	0.134	0.1050	0.899	3.13	-0.0211	0.0689	-19.1	0.104
Spatially flat $\Lambda$CDM	${0.365}_{-0.031\left(1\sigma \right)\mathrm{\hspace{0.33em}0.060}\left(2\sigma \right)}^{+0.033\left(1\sigma \right)\mathrm{\hspace{0.33em}0.066}\left(2\sigma \right)}$	0.134	0.1061	0.899	3.14	-0.0215	0.0688	-19.0	0.104
Empty universe	$-$	0.133	0.1013	0.900	3.13	-0.0204	0.0690	-19.0	0.106

Table 3. Table 3: Selected parameters fit for z min = 0.033 subscript 𝑧 min 0.033 z_{\rm min}=0.033 , with the following empirical models for light-curve parameters: (I) constant x 1 , 0 subscript 𝑥 1 0 x_{\lower 1.0pt\hbox{$\scriptstyle 1,0$}} , constant c 0 subscript 𝑐 0 c_{\lower 1.0pt\hbox{$\scriptstyle 0$}} ; (II) global linear x 1 , 0 subscript 𝑥 1 0 x_{\lower 1.0pt\hbox{$\scriptstyle 1,0$}} , constant c 0 subscript 𝑐 0 c_{\lower 1.0pt\hbox{$\scriptstyle 0$}} ; (III) split linear x I , 1 , 0 subscript 𝑥 𝐼 1 0 x_{\lower 1.0pt\hbox{$\scriptstyle I,1,0$}} , constant c 0 subscript 𝑐 0 c_{\lower 1.0pt\hbox{$\scriptstyle 0$}} ; (IV) constant x 1 , 0 subscript 𝑥 1 0 x_{\lower 1.0pt\hbox{$\scriptstyle 1,0$}} , global linear c 0 subscript 𝑐 0 c_{\lower 1.0pt\hbox{$\scriptstyle 0$}} ; (V) constant x 1 , 0 subscript 𝑥 1 0 x_{\lower 1.0pt\hbox{$\scriptstyle 1,0$}} , split linear c 0 , I subscript 𝑐 0 𝐼 c_{\lower 1.0pt\hbox{$\scriptstyle 0,I$}} ; (VI) global linear x 1 , 0 subscript 𝑥 1 0 x_{\lower 1.0pt\hbox{$\scriptstyle 1,0$}} , global linear c 0 subscript 𝑐 0 c_{\lower 1.0pt\hbox{$\scriptstyle 0$}} ; (VII) split linear x 1 , 0 subscript 𝑥 1 0 x_{\lower 1.0pt\hbox{$\scriptstyle 1,0$}} , split linear c 0 subscript 𝑐 0 c_{\lower 1.0pt\hbox{$\scriptstyle 0$}} ; (VIII) global linear x 1 , 0 subscript 𝑥 1 0 x_{\lower 1.0pt\hbox{$\scriptstyle 1,0$}} , split linear c 0 subscript 𝑐 0 c_{\lower 1.0pt\hbox{$\scriptstyle 0$}} . Notes: k ≡ 𝑘 absent k\equiv number of free parameters; quantities ⟨ Φ ⟩ ≡ ( ∑ N I ​ Φ I ) / ( ∑ I = 1 4 N I ) delimited-⟨⟩ Φ subscript 𝑁 𝐼 subscript Φ 𝐼 superscript subscript 𝐼 1 4 subscript 𝑁 𝐼 \langle\Phi\rangle\equiv(\sum N_{\lower 1.0pt\hbox{$\scriptstyle I$}}\,\Phi_{\lower 1.0pt\hbox{$\scriptstyle I$}})/(\sum_{I=1}^{4}N_{\lower 1.0pt\hbox{$\scriptstyle I$}}) denote an average over subsamples with I = 1 ​ … ​ 4 𝐼 1 … 4 I=1\dots 4 for x 1 , 0 , I subscript 𝑥 1 0 𝐼 x_{\lower 1.0pt\hbox{$\scriptstyle 1,0,I$}} , c 0 , I subscript 𝑐 0 𝐼 c_{\lower 1.0pt\hbox{$\scriptstyle 0,I$}} and I = 1 ​ … ​ 3 𝐼 1 … 3 I=1\dots 3 for x z , I subscript 𝑥 𝑧 𝐼 x_{\lower 1.0pt\hbox{$\scriptstyle z,I$}} , c z , I subscript 𝑐 𝑧 𝐼 c_{\lower 1.0pt\hbox{$\scriptstyle z,I$}} for split models or ⟨ Φ ⟩ ≡ Φ delimited-⟨⟩ Φ Φ \langle\Phi\rangle\equiv\Phi otherwise; BIC = = Bayesian Information Criterion; E = 𝐸 absent E= Bayesian evidence; C b = subscript 𝐶 𝑏 absent C_{b}= Bayesian complexity; Δ BIC = BIC model − BIC I subscript Δ BIC subscript BIC model subscript BIC I \Delta_{\lower 1.0pt\hbox{$\scriptstyle\rm BIC$}}=\hbox{BIC}_{\rm model}-\hbox{BIC}_{\lower 1.0pt\hbox{$\scriptstyle\rm I$}} and ln ⁡ B 1 = ln ⁡ ( E I / E model ) subscript 𝐵 1 subscript 𝐸 I subscript 𝐸 model \ln B_{\lower 1.0pt\hbox{$\scriptstyle 1$}}=\ln(E_{\lower 1.0pt\hbox{$\scriptstyle\rm I$}}/E_{\rm model}) are evaluated with cosmological model fixed; ln ⁡ B 2 = ln ⁡ ( E TS / E model ) subscript 𝐵 2 subscript 𝐸 TS subscript 𝐸 model \ln B_{\lower 1.0pt\hbox{$\scriptstyle 2$}}=\ln(E_{\lower 1.0pt\hbox{$\scriptstyle\rm TS$}}/E_{\rm model}) is evaluated with light-curve model fixed.

Timescape	$k$	${\Omega }_{M0}$	$⟨x_{1,0}⟩$	$⟨c_0⟩$	$⟨x_z⟩$	$⟨c_z⟩$	BIC	${\Delta }_{\mathrm{BIC}}$	$\ln E$	$C_b$	$\ln B_1$	$f_{\mathrm{v0}}$
I (NGS16)	9	0.309	0.105	-0.021			-185.0		80.38	8.53		0.778
II	10	0.278	-0.073	-0.021	0.511		-199.6	-14.6	81.67	9.28	-1.28	0.802
III	15	0.278	-0.183	-0.021	0.806		-169.9	15.1	66.46	14.11	13.93	0.801
IV	10	0.000	0.104	0.002		-0.065	-249.9	-64.9	78.52	9.19	1.86	1.000
V	15	0.010	0.092	0.054		-0.351	-157.7	27.3	81.29	14.15	-0.91	0.993
VI	11	0.000	-0.071	0.001	0.499	-0.062	-189.7	-4.7	79.35	10.21	1.94	1.000
VII (RH16)	21	0.000	-0.123	0.054	0.490	-0.348	-200.0	-15.0	65.85	19.47	14.53	1.000
VIII	16	0.000	-0.085	0.061	0.501	-0.348	-229.3	-44.3	82.03	15.30	-1.65	1.000
Flat $\Lambda$CDM	$k$	${\Omega }_{M0}$	$⟨x_{1,0}⟩$	$⟨c_0⟩$	$⟨x_z⟩$	$⟨c_z⟩$	BIC	${\Delta }_{\mathrm{BIC}}$	$\ln E$	$C_b$	$\ln B_1$	$\ln B_2$
I (NGS16)	9	0.365	0.106	-0.021			-192.5		80.30	8.93		0.08
II	10	0.353	-0.069	-0.021	0.503		-241.2	-48.7	81.64	10.01	-1.34	0.03
III	15	0.353	-0.186	-0.021	0.847		-159.8	32.7	66.62	14.60	13.68	-0.16
IV	10	0.303	0.106	-0.002		-0.057	-192.9	-0.4	79.60	9.98	0.70	-1.08
V	15	0.296	0.093	0.052		-0.354	-228.5	-36.0	83.77	14.87	-3.47	-2.47
VI	11	0.292	-0.069	-0.002	0.501	-0.057	-179.1	13.4	80.87	10.97	2.89	-1.52
VII (RH16)	21	0.286	-0.127	0.051	0.534	-0.352	-155.3	37.2	68.97	20.58	11.33	-3.12
VIII	16	0.286	-0.080	0.059	0.499	-0.354	-232.8	-40.3	84.95	15.89	-4.65	-2.92
Empty	$k$	${\Omega }_{M0}$	$⟨x_{1,0}⟩$	$⟨c_0⟩$	$⟨x_z⟩$	$⟨c_z⟩$	BIC	${\Delta }_{\mathrm{BIC}}$	$\ln E$	$C_b$	$\ln B_1$	$\ln B_2$
I (NGS16)	8	-	0.101	-0.020			-181.5		78.18	8.11		2.20
II	9	-	-0.078	-0.019	0.517		-190.1	-8.6	79.92	9.02	-1.74	1.75
III	14	-	-0.095	-0.020	0.749		-218.9	-37.4	64.43	13.75	13.76	2.03
IV	9	-	0.098	0.002		-0.054	-185.7	-4.2	78.56	9.05	-0.37	-0.03
V	14	-	0.087	0.054		-0.336	-180.4	1.1	79.85	14.19	-1.66	1.45
VI	10	-	-0.072	0.002	0.489	-0.051	-186.3	-4.8	79.62	10.17	0.23	-0.27
VII (RH16)	20	-	-0.122	0.054	0.460	-0.332	-198.7	-17.2	64.31	18.68	13.88	1.55
VIII	15	-	-0.081	0.061	0.482	-0.334	-221.4	-39.9	80.74	14.89	-2.55	1.30

Equations

$$\mu \equiv 25+5{\log }_{10}\left(\frac{d_{\text{lower}1.0ptL}}{Mpc}\right),$$

$${\mu }_{\text{lower}1.0pt\mathrm{SN}}=m_{\text{lower}1.0pt\mathrm{B}}^{\ast }-M_{\text{lower}1.0pt\mathrm{B}}+\alpha x_{\text{lower}1.0pt1}-\beta c,$$

$$\Delta {\mu }_{\text{lower}1.0pt\mathrm{B}}\equiv (m-M)-(m_{\text{lower}1.0pt\mathrm{B}}^{\ast }-M_{\text{lower}1.0pt\mathrm{B}})=\alpha x_{\text{lower}1.0pt1}-\beta c,$$

$$x_{\text{lower}1.0pt1,0}\to x_{\text{lower}1.0pt1,0,J}+x_{\text{lower}1.0ptz,J}z,and{c}_{\text{lower}1.0pt0}\to c_{\text{lower}1.0pt0,J}+c_{\text{lower}1.0ptz,J}z,$$

$$\mathcal{L}$$

$$=\int \mathrm{d}{M}_{\text{lower}1.0pt\mathrm{B}}^N\mathrm{d}{x}_{\text{lower}1.0pt1}^N\mathrm{d}{c}^N\mathcal{P}\left[(\hat{z},{\hat{m}}_{\text{lower}1.0pt\mathrm{B}}^{\ast },{\hat{x}}_{\text{lower}1.0pt1},\hat{c}{)}_i|(M_{\text{lower}1.0pt\mathrm{B}},x_{\text{lower}1.0pt1},c{)}_i,\mathrm{M}\right]$$

$$to85.35826pt\text{hfil}\times \mathcal{P}\left[(M_{\text{lower}1.0pt\mathrm{B}},x_{\text{lower}1.0pt1},c{)}_i|\mathrm{M}\right],$$

$$\mathcal{P}[(M_{\text{lower}1.0pt\mathrm{B}},x_{\text{lower}1.0pt1},c{)}_i|\mathrm{M}]=\prod _i^N\mathcal{P}[(M_{\text{lower}1.0pt\mathrm{B},i},x_{\text{lower}1.0pt1,i},c_{\text{lower}1.0pti})|\mathrm{M}]$$

$$to8.53581pt\text{hfil}=\prod _i^N\mathcal{N}(M_{\text{lower}1.0pt\mathrm{B},i};M_{\text{lower}1.0pt\mathrm{B},0},{\sigma }_{\text{lower}1.0pt{M}_{\mathrm{B},0}})$$

$$to42.67912pt\text{hfil}\times \mathcal{N}(x_{\text{lower}1.0pt1,i};x_{\text{lower}1.0pt1,0},{\sigma }_{x_{\text{lower}1.0pt1,0}})\mathcal{N}(c_{\text{lower}1.0pti};c_{\text{lower}1.0pt0},{\sigma }_{c_{\text{lower}1.0pt0}}),$$

$$\mathcal{L}$$

$$={\left[(2\pi {)}^{3N}\det \left({\Sigma }_d+A^{⊺}{\Sigma }_{l}A\right)\right]}^{-1/2}$$

$$to14.22636pt\text{hfil}\times \exp \left[-\frac{1}{2}\left(\hat{Z}-Y_{\text{lower}1.0pt0}A\right){\left({\Sigma }_d+A^{⊺}{\Sigma }_{l}A\right)}^{-1}{\left(\hat{Z}-Y_{\text{lower}1.0pt0}A\right)}^{⊺}\right]$$

$$Y_{\text{lower}1.0pt0}A=(M_{\text{lower}1.0pt\mathrm{B},0}-\alpha x_{\text{lower}1.0pt1,0}+\beta c_{\text{lower}1.0pt0},x_{\text{lower}1.0pt1,0},c_{\text{lower}1.0pt0},\dots ).$$

$$\mathcal{P}[(M_{\text{lower}1.0pt\mathrm{B}},x_{\text{lower}1.0pt1},c{)}_i|\mathrm{M}]$$

$$to8.53581pt\text{hfil}=\prod _{J=1}^4\prod _{i=1}^{N_{\text{lower}1.0ptJ}}\mathcal{N}(M_{\text{lower}1.0pt\mathrm{B},i},M_{\text{lower}1.0pt\mathrm{B},0},{\sigma }_{\text{lower}1.0pt{M}_{\mathrm{B},0}})$$

$$to12.80373pt\text{hfil}\times \mathcal{N}(x_{\text{lower}1.0pt1,i},x_{\text{lower}1.0pt1,0,J}+x_{\text{lower}1.0ptz,J}\hat{z},{\sigma }_{x_{\text{lower}1.0pt1,0}})\mathcal{N}(c_{\text{lower}1.0pti},c_{\text{lower}1.0pt0,J}+c_{\text{lower}1.0ptz,J}\hat{z},{\sigma }_{c_{\text{lower}1.0pt0}}).$$

$${\mu }_{\text{lower}1.0pt\Lambda \mathrm{CDM}}(0.286)-{\mu }_{\text{lower}1.0pt\Lambda \mathrm{CDM}}(0.365)=0.1287z-0.0085z^2$$

$$to22.76219pt\text{hfil}-0.0481z^3+0.0249z^4+0.0161z^5-0.0232z^6+\dots$$

$${\mu }_{\text{lower}1.0pt\mathrm{TS}}(1.0)-{\mu }_{\text{lower}1.0pt\mathrm{TS}}(0.778)=0.0674z+0.0444z^2$$

$$to22.76219pt\text{hfil}-0.0242z^3+0.0190z^4-0.0193z^5+0.0173z^6+\dots$$

$${\mu }_{\text{lower}1.0pt\mathrm{FLRW}}=25+5{\log }_{\text{lower}1.0pt10}\left(\frac{\mathtt{c}z}{H_{\text{lower}1.0pt0}Mpc}\right)$$

$$\displaystyle+{5\over\ln 10}\Bigl{\{}{\textstyle{1\over 2}}(1-q_{\lower 1.0pt\hbox{$\scriptstyle 0$}})z+{\textstyle{1\over 24}}\Bigl{[}9{q_{\lower 1.0pt\hbox{$\scriptstyle 0$}}}^{2}-10q_{\lower 1.0pt\hbox{$\scriptstyle 0$}}-7-4j_{\lower 1.0pt\hbox{$\scriptstyle 0$}}-\Omega_{\lower 1.0pt\hbox{$\scriptstyle k0$}}\Bigr{]}z^{2}$$

$$\displaystyle\hbox to36.98857pt{\hfil}+{\textstyle{1\over 24}}\Bigl{[}s_{\lower 1.0pt\hbox{$\scriptstyle 0$}}+5-10{q_{\lower 1.0pt\hbox{$\scriptstyle 0$}}}^{3}-16{q_{\lower 1.0pt\hbox{$\scriptstyle 0$}}}^{2}-9q_{\lower 1.0pt\hbox{$\scriptstyle 0$}}+j_{\lower 1.0pt\hbox{$\scriptstyle 0$}}\left(8q_{\lower 1.0pt\hbox{$\scriptstyle 0$}}+7\right)$$

$$\displaystyle\hbox to65.44133pt{\hfil}-4\Omega_{\lower 1.0pt\hbox{$\scriptstyle k0$}}\left(q_{\lower 1.0pt\hbox{$\scriptstyle 0$}}+1\right)\Bigr{]}z^{3}+\dots\Bigr{\}},$$

$$d_{\text{lower}1.0ptL}$$

$$\mathcal{H}(y)\equiv \sqrt{{\Omega }_{\text{lower}1.0ptR0}+{\Omega }_{\text{lower}1.0ptM0}y+{\Omega }_{\text{lower}1.0ptk0}{y}^2+{\Omega }_{\text{lower}1.0pt\Lambda 0}{y}^4},$$

$$sinn(x)\equiv \{\begin{array}{ll}\sinh (x),& {\Omega }_{\text{lower}1.0ptk0}>0\\ x,& {\Omega }_{\text{lower}1.0ptk0}=0\\ \sin (x),& {\Omega }_{\text{lower}1.0ptk0}<0\end{array},$$

$$d_{\text{lower}1.0ptL}=$$

$$d_{\text{lower}1.0ptA}=$$

$$=$$

$$\mathcal{F}(t)$$

$$to56.9055pt\text{hfil}+\frac{b^{1/3}}{\sqrt{3}}{\tan }^{-1}\left(\frac{2t^{1/3}-b^{1/3}}{\sqrt{3}{b}^{1/3}}\right),$$

$$z+1=\frac{(2+f_{\mathrm{v}}){f_{\mathrm{v}}}^{1/3}}{3f_{v0}^{1/3}{\bar{H}}_{\text{lower}1.0pt0}t}=\frac{2^{4/3}{t}^{1/3}(t+b)}{f_{v0}^{1/3}{\bar{H}}_{\text{lower}1.0pt0}t(2t+3b{)}^{4/3}},$$

$$f_{\mathrm{v}}(t)=\frac{3f_{v0}{\bar{H}}_{\text{lower}1.0pt0}t}{3f_{v0}{\bar{H}}_{\text{lower}1.0pt0}t+(1-f_{v0})(2+f_{v0})},$$

$$\tau =\frac{2}{3}t+\frac{2(1-f_{v0})(2+f_{v0})}{27f_{v0}{\bar{H}}_{\text{lower}1.0pt0}}\ln \left(1+\frac{9f_{v0}{\bar{H}}_{\text{lower}1.0pt0}t}{2(1-f_{v0})(2+f_{v0})}\right).$$

$$a\equiv {\bar{\gamma }}^{-1}\bar{a},$$

$${\vphantom{\tilde{a}}\bar{a}}={{\vphantom{\tilde{a}}\bar{a}}_{\lower 1.0pt\hbox{$\scriptstyle 0$}}\bigl{(}3{\vphantom{\tilde{H}}\bar{H}}_{\lower 1.0pt\hbox{$\scriptstyle\!0$}}t\bigr{)}^{2/3}\over 2+f_{\rm v0}}\left[3f_{\rm v0}{\vphantom{\tilde{H}}\bar{H}}_{\lower 1.0pt\hbox{$\scriptstyle\!0$}}t+(1-f_{\rm v0})(2+f_{\rm v0})\right]^{1/3}$$

$$\bar{\gamma }=\frac{1}{2}(2+f_{\mathrm{v}})=\frac{3(t+b)}{(2t+3b)}$$

$$q=\frac{-\left(1-f_{\mathrm{v}}\right)(8{f_{\mathrm{v}}}^3+39{f_{\mathrm{v}}}^2-12f_{\mathrm{v}}-8)}{{\left(4+f_{\mathrm{v}}+4{f_{\mathrm{v}}}^2\right)}^2}.$$

Full text

1 Inhomogeneous alternatives to dark energy

One of the most important discoveries in cosmology (Riess et al., 1998; Perlmutter et al., 1999) is the observation that the luminosity distances and redshifts of type Ia supernovae (SneIa) are well matched to the expansion history of a spatial homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) model only if the Universe began an epoch of accelerated expansion late in cosmic history. Since gravity with matter obeying the Strong Energy Condition is universally attractive, this demands a cosmological constant, $\Lambda$, or some other unknown source of spatially homogeneous dark energy with an equation of state, $P=w\rho c^{2}$, that subsequent tests find to be close to the cosmological constant case, $w=-1$. Independently of SneIa, since the late 1990s other data sets including the cosmic microwave background (CMB) and galaxy clustering statistics, have been found to independently require late epoch cosmic acceleration in the FLRW model.

Despite this success, the nature of dark energy remains a mystery for fundamental physics. Furthermore, a number of puzzles remain (Bull et al., 2016; Buchert et al., 2016) which range from significant anomalies (Cyburt et al., 2008) to lesser tensions (Delubac et al., 2015). While some puzzles may inevitably boil down to statistical sampling biases, systematic problems in data reduction or potentially unknown astrophysics, they may also in part be due to incorrect cosmological model assumptions. It is therefore important to consider alternatives to the standard cosmology, particularly any alternatives which do not change known gravitational physics but which reconsider assumptions that are not demanded theoretically in light of observations. The assumption of spatial homogeneity and isotropy falls in this category. While the remarkable isotropy of the CMB points to an initial state with a very high degree of smoothness, the late epoch Universe encompasses a complex cosmic web of structures (van de Weygaert et al., 2016). It is dominated in volume by voids that are threaded and surrounded by clusters of galaxies distributed in sheets, knots and filaments.

Spatial homogeneity is first encountered in some statistical sense in averages on scales \mathop{\hbox{{\lower 3.8pt\hbox{$>$}}\atop{\raise 0.2pt\hbox{$\sim$}}}}70–$120\,h^{-1}\hbox{Mpc}$ based on the 2-point galaxy correlation function (Hogg et al., 2005; Scrimgeour et al., 2012). This indicates that we can expect cosmology to be modelled by some averaged description on larger scales, but does not guarantee that such a description is given at all times by the FLRW geometry, however. Demanding FLRW geometry on arbitrarily large scales in space and time presupposes one particular answer to the fitting problem (Ellis, 1984; Ellis & Stoeger, 1987; Wiltshire, 2011). The fitting problem is an unanswered fundamental question in general relativity: on what scale(s) are matter and geometry dynamically coupled via Einstein’s equations? Since general relativity is a causal theory, the speed of light yields a natural upper limit. However, most of the energy density of the present epoch Universe in the energy-momentum tensor on the right hand side of Einstein’s equations is in non-relativistic matter subject to speeds much smaller than that of light. Consequently, the effective causal scale for density fluctuations is the matter horizon (Ellis & Stoeger, 2009), the comoving region which has significantly contributed matter to the local physical environment of a given observer.

The matter horizon gives a natural effective causal upper limit for direct applicability of the Einstein equations. Regardless of what the scale is, if Einstein’s equations apply only on small scales then generic large-scale averages do not lead to the Einstein equations for a prescribed averaged matter source. Instead one obtains averaged cosmological models with backreaction (Buchert, 2000, 2001; Buchert & Räsänen, 2012), whose evolution generically can differ significantly from that of the FLRW geometry. Since the cosmic expansion history is phenomenologically very well described by a FLRW model with a small number of parameters, then any physically relevant model of backreaction must be subject to some simplifying principles that restrict many possible mathematical choices for defining the average expansion history and its relation to local observables (Buchert & Carfora, 2002, 2003).

The timescape scenario (Wiltshire, 2007a, b, 2009; Duley, Nazer & Wiltshire, 2013) provides such principles. Within Buchert’s scalar averaging formalism (Buchert, 2000, 2001), the relationship between statistical volume averages and local observables is restricted by applying the Cosmological Equivalence Principle (Wiltshire, 2008), a generalization of the Strong Equivalence Principle. The clocks and rulers which best describe average cosmic evolution (the ‘bare’ or ‘volume–average’ quantities) in Buchert’s statistical formalism are reinterpreted, and differ systematically from local clocks and rulers of observers in environments where the matter density differs significantly from the volume average. Given the smooth initial state of the Universe, this only becomes significant late in cosmic history, when spatial curvature gradients grow between bound structures (which by definition are greater than critical density) and the voids which dominate by volume.

The timescape model has passed a number of independent observational tests (Leith, Ng & Wiltshire, 2008; Smale & Wiltshire, 2011; Smale, 2011; Duley, Nazer & Wiltshire, 2013; Nazer & Wiltshire, 2015). Its distance–redshift relation is very close to that of particular $\Lambda$ cold dark matter ($\Lambda$CDM) models over small redshift ranges, but effectively interpolates (Wiltshire, 2009, 2014) between spatially flat $\Lambda$CDM cosmologies with different values of \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} and \Omega_{\lower 1.0pt\hbox{\scriptstyle\Lambda 0}} over larger redshift ranges. In particular, when the timescape model is fit to the angular diameter distance of the sound horizon in the CMB then the spatially flat $\Lambda$CDM model with the closest comoving distance at redshifts z\mathop{\hbox{{\lower 3.8pt\hbox{$>$}}\atop{\raise 0.2pt\hbox{$\sim$}}}}100 has a value of \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} 15–27% lower than that of the spatially flat $\Lambda$CDM model with the closest comoving distance at redshifts $z<1.2$ (Wiltshire, 2009, 2014).

Geometric tests of the timescape expansion history are most developed (Wiltshire, 2009), and give rise to measures (Clarkson, Bassett & Lu, 2008; Räsänen, Bolejko & Finoguenov, 2015) which will definitely distinguish both the timescape model and other alternatives (Larena et al., 2009; Lavinto, Räsänen & Szybka, 2013) from the standard FLRW model using Euclid satellite data (Sapone, Majerotto & Nesseris, 2014). On the other hand, tests of the CMB anisotropy spectrum in the timescape model are at present limited by systematic uncertainties of 8–13% in parameters which relate to the matter content (Nazer & Wiltshire, 2015). This is a consequence of backreaction schemes having not yet been applied to the primordial plasma.

2 Supernova redshift-distance analysis

In the case of the redshift range probed by SneIa the difference between the timescape and $\Lambda$CDM cosmologies is comparable to the systematic uncertainties that arise between different methods for fitting the light curves of SneIa to obtain “standard candles”. In particular, in the last full analysis of the timescape model (using data available in 2010) Smale & Wiltshire (2011) found significant differences between data reduced by the MLCS2k2 (Multicolor Light-Curve Shape) fitter (Jha et al., 2007) and the SALT/SALT2 (Spectral Adaptive Light-curve Template) fitters (Guy et al., 2005, 2007). While the relative Bayesian evidence was sometimes ‘positive’ (but not very strong), the conclusion as to which cosmological model fitted better depended on the light-curve fitting method. Consequently the empirical nature of light-curve fitting may mask effects due to the underlying expansion history if this deviates from the FLRW geometry.

The significantly larger Joint Light-curve Analysis (JLA) SneIa catalogue (Betoule et al., 2014) now makes possible a renewed comparison of the timescape and $\Lambda$CDM models, as well as further investigation of the systematics of light-curve fitting. Recently, Nielsen, Guffanti & Sarkar (2016) (NGS16) have used the JLA catalogue to reinvestigate systematic issues associated with SneIa light-curve fitting within FLRW cosmologies using the SALT2 method. They adopted maximum likelihood estimators (MLE) that take into account the underlying Gaussianity of particular light-curve parameters (March et al., 2011). NGS16 concluded that the significance for cosmic acceleration, as compared to an empty Milne model (or any cosmology with constant expansion), is “marginal” (at \mathop{\hbox{{\lower 3.8pt\hbox{$<$}}\atop{\raise 0.2pt\hbox{$\sim$}}}}3\sigma significance). This conclusion was challenged by Rubin & Hayden (2016) (RH16), who introduced 12 additional light-curve parameters to allow for possible unaccounted systematics, concluding that the $2.8\sigma$ significance found by NGS16 increased to $3.7\sigma$ for a general FLRW model, or to $4.2\sigma$ for the spatially flat case. However, RH16 did not consider whether the increased model complexity was justified from a Bayesian standpoint.

In the SALT2 method each observed supernova redshift is used to determine a theoretical distance modulus,

[TABLE]

where d_{\lower 1.0pt\hbox{\scriptstyle L}} is the luminosity distance for each cosmological model. This is then compared to the observed distance modulus, which is related to the supernova light-curve by

[TABLE]

where m_{\lower 1.0pt\hbox{\scriptstyle\rm B}}^{*} is the apparent magnitude at maximum in the rest-frame $B$ band, M_{\lower 1.0pt\hbox{\scriptstyle\rm B}} is the corresponding absolute magnitude of the source, $x_{1}$ and $c$ are empirical parameters that describe the light-curve stretch and colour corrections for each supernova, while $\alpha$ and $\beta$ are parameters that are assumed to be constant for all SneIa.

The theoretical distance modulus (1) is determined for a bolometric flux, which is not directly measured. The SALT2 (Guy et al., 2005, 2007) relation (2) can thus be viewed as a model for a band correction, \Delta\mu_{\lower 1.0pt\hbox{\scriptstyle\rm B}}, that is linear in the variables x_{\lower 1.0pt\hbox{\scriptstyle 1}} and $c$,

[TABLE]

where $m$ and $M$ are the bolometric apparent and absolute magnitudes in the observer and emitter rest frames respectively.

In the SALT2 method, the light-curve parameters are simultaneously fit together with the free cosmological parameters on the entire data set.

NGS16 assumed that all SneIa in the JLA catalogue (Betoule et al., 2014) are characterized by parameters, M_{\lower 1.0pt\hbox{\scriptstyle\rm B}}, x_{\lower 1.0pt\hbox{\scriptstyle 1}} and $c$, drawn from the same independent global Gaussian distributions, with means M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}}, x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} and c_{\lower 1.0pt\hbox{\scriptstyle 0}}, and standard deviations \sigma_{\lower 1.0pt\hbox{\scriptstyle M_{\rm B,0}}}, \sigma_{\lower 1.0pt\hbox{\scriptstyle x_{\rm 1,0}}} and \sigma_{\lower 1.0pt\hbox{\scriptstyle c_{\rm 0}}} respectively. These 6 free parameters were then fitted along with the light-curve parameters $\alpha$, $\beta$ and the cosmological parameters.

RH16 claimed that the mean light-curve stretch and colour parameters, x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} and c_{\lower 1.0pt\hbox{\scriptstyle 0}}, of the Gaussian distributions analysed by NGS16 show some redshift dependence. This may be partially due to astrophysical effects in the host population, or – particularly for the colour parameter – may arise from the colour–luminosity relation combined with redshift–dependent detection limits. In other words, Malmquist type biases may not be completely corrected for in the JLA catalogue (Betoule et al., 2014). In the absence of a known astrophysical model for such corrections, RH16 introduced 12 additional empirical parameters by replacing the global Gaussian means according to

[TABLE]

where the index $J$ runs over the four independent subsamples in the JLA catalogue: (1) SNLS (SuperNova Legacy Survey); (2) SDSS (Sloan Digital Sky Survey); (3) nearby supernovae; (4) HST (Hubble Space Telescope), with x_{\lower 1.0pt\hbox{\scriptstyle z,4}}=0, c_{\lower 1.0pt\hbox{\scriptstyle z,4}}=0 on account of limited HST data. The widths \sigma_{\lower 1.0pt\hbox{\scriptstyle x_{\rm 1,0}}}, \sigma_{\lower 1.0pt\hbox{\scriptstyle c_{\rm 0}}} were still treated as global parameters.

We perform a Bayesian comparison of fits of the JLA catalogue (Betoule et al., 2014) to the luminosity distance-redshift relation for the spatially flat $\Lambda$CDM model, and for the timescape model (Wiltshire, 2007a, b, 2009; Duley et al., 2013). We first use the MLE methodology of NGS16 directly, and then investigate the effect of changes to light-curve fitting suggested by RH16. An empty universe with constant expansion rate is also analysed, as a convenient demarcation of accelerating from non-accelerating expansion in the FLRW case.

Details of the theoretical luminosity distances used in (1) are given in Appendix A. The model differences that we are testing are best appreciated by comparing the distance moduli of the timescape and spatially flat $\Lambda$CDM models relative to an empty universe, as shown in Fig. 1. The timescape distance modulus, \mu_{\lower 1.0pt\hbox{\scriptstyle\rm TS}}(z), is closer to $\Lambda$CDM than the empty case. Nonetheless, \mu_{\lower 1.0pt\hbox{\scriptstyle\rm TS}}(z) is always closer to \mu_{\lower 1.0pt\hbox{\scriptstyle\rm empty}}(z) than \mu_{\lower 1.0pt\hbox{\scriptstyle\Lambda{\rm CDM}}}(z) is, a consequence of cosmic acceleration being an apparent effect in the timescape model.

Further technicalities about systematic issues in implementing the SALT2 method are discussed in Appendix B.

3 Statistical methods

3.1 The likelihood construction

We adopt the likelihood construction (March et al., 2011) used by NGS16. The likelihood, $\mathcal{L}$, is the probability density of the observed data – here (\hat{z},{\hat{m}}_{\lower 1.0pt\hbox{\scriptstyle\rm B}}^{*},\hat{x}_{\lower 1.0pt\hbox{\scriptstyle 1}},\hat{c})_{i}, $i=1,2,\dots,N$ on $N$ supernovae – given a model, $\mathrm{M}$. The likelihood can be written as (March et al., 2011)

[TABLE]

where hatted quantities denote measured data values including all experimental noise, and unhatted quantities are intrinsic parameters that characterize the statistical distributions from which the supernovae are drawn. Only the intrinsic parameters satisfy the SALT2 relation (2). The empirical light-curve model (2) and the theoretical distance modulus (1) together constitute the model, $\mathrm{M}$.

The expansion in (3.1) allows one to explicitly model the intrinsic scatter of the supernovae. For the NGS16 model (I) we assume that the intrinsic parameters of each supernova are drawn from identical independent Gaussian distributions

[TABLE]

where \mathcal{N}(y\,;\,y_{\lower 1.0pt\hbox{\scriptstyle 0}},\sigma)\equiv(2\pi\sigma^{2})^{-1/2}\exp\left[-(y-y_{\lower 1.0pt\hbox{\scriptstyle 0}})^{2}/(2\sigma^{2})\right] for each triple \{y,y_{\lower 1.0pt\hbox{\scriptstyle 0}},\sigma\}, with $3N\times 3N$ diagonal covariance matrix \Sigma_{l}=\mathrm{diag}\left(\sigma^{2}_{\lower 1.0pt\hbox{\scriptstyle M_{\lower 1.0pt\hbox{$\scriptscriptstyle{\rm B},0$}}}},\sigma^{2}_{x_{1,0}},\sigma^{2}_{\lower 1.0pt\hbox{\scriptstyle c_{\lower 1.0pt\hbox{$\scriptscriptstyle 0$}}}},\sigma^{2}_{\lower 1.0pt\hbox{\scriptstyle M_{\lower 1.0pt\hbox{$\scriptscriptstyle{\rm B},0$}}}},\dots\right). The experimental part of the likelihood (3.1), \mathcal{P}\left[(\hat{z},{\hat{m}}_{\lower 1.0pt\hbox{\scriptstyle\rm B}}^{*},\hat{x}_{\lower 1.0pt\hbox{\scriptstyle 1}},\hat{c})_{i}\,|\,(M_{\lower 1.0pt\hbox{\scriptstyle\rm B}},x_{\lower 1.0pt\hbox{\scriptstyle 1}},c)_{i},\mathrm{M}\right] is assumed to be a Gaussian in the intrinsic supernova parameters, and the covariance matrix of experimental statistical and systematic uncertainties is denoted $\Sigma_{d}$. Performing the integral in (3.1) one obtains the final expression of the likelihood (Nielsen et al., 2016)

[TABLE]

where \hat{Z}=({\hat{m}}_{\lower 1.0pt\hbox{\scriptstyle{\rm B},1}}^{*}-\mu(\hat{z}_{\lower 1.0pt\hbox{\scriptstyle 1}}),\hat{x}_{\lower 1.0pt\hbox{\scriptstyle 1,1}},\hat{c}_{\lower 1.0pt\hbox{\scriptstyle 1}},{\hat{m}}_{\lower 1.0pt\hbox{\scriptstyle{\rm B},2}}^{*}-\mu(\hat{z}_{\lower 1.0pt\hbox{\scriptstyle 2}}),\dots) is a $3N$-dimensional row vector containing the distance modulus residual and light-curve data, Y_{\lower 1.0pt\hbox{\scriptstyle 0}}=(M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}},\,x_{\lower 1.0pt\hbox{\scriptstyle 1,0}},\,c_{\lower 1.0pt\hbox{\scriptstyle 0}},\,M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}},\dots) is a $3N$-dimensional row vector of the intrinsic Gaussian means, and $A$ is the block diagonal matrix that propagates Y_{\lower 1.0pt\hbox{\scriptstyle 0}} to

[TABLE]

Note that the cosmological model enters only explicitly through the conversion $\mu(\hat{z})$ of the observed redshift to a distance modulus. There can, however, be implicit model dependence in the experimental covariance matrix111The propagation of the error $\sigma_{z}$ to $\sigma_{\mu}$ depends on the model. However, by (23)–(25), to leading order for small $z$, $\sigma_{\mu}\simeq 5\sigma_{z}/(z\,\ln 10)$ for all cases. $\Sigma_{d}$ or in corrections made to data prior to the analysis.

To implement the RH16 parametrization (4) we replace (6) by

[TABLE]

We recover (3.1) with the one difference: in place of three repeated entries, the vector Y_{\lower 1.0pt\hbox{\scriptstyle 0}} is now partitioned into different pieces for each subsample,Y_{\lower 1.0pt\hbox{\scriptstyle 0}}=(M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}},\,x_{\lower 1.0pt\hbox{\scriptstyle 1,0,1}}+x_{\lower 1.0pt\hbox{\scriptstyle z,1}}{\hat{z}}_{\lower 1.0pt\hbox{\scriptstyle 1}},\,c_{\lower 1.0pt\hbox{\scriptstyle 0,1}}+c_{\lower 1.0pt\hbox{\scriptstyle z,1}}{\hat{z}}_{\lower 1.0pt\hbox{\scriptstyle 1}},\,\dots,\,M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}},\,x_{\lower 1.0pt\hbox{\scriptstyle 1,0,2}}+x_{\lower 1.0pt\hbox{\scriptstyle z,2}}{\hat{z}}_{\lower 1.0pt\hbox{\scriptstyle i}},\,c_{\lower 1.0pt\hbox{\scriptstyle 0,2}}+c_{\lower 1.0pt\hbox{\scriptstyle z,2}}{\hat{z}}_{\lower 1.0pt\hbox{\scriptstyle i}},\dots,\dots,M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}},x_{\lower 1.0pt\hbox{\scriptstyle 1,0,4}},c_{\lower 1.0pt\hbox{\scriptstyle 0,4}}).

From the likelihood (3.1) we can define frequentist confidence regions and goodness of fit measures or alternative Bayesian versions of these, following conventional statistical procedures summarized in Appendix C.

In practice, estimating the Bayesian evidence is a computationally intensive task, much more so than what is required to obtain parameter estimates. We use standard Markov Chain Monte Carlo (MCMC) methods to sample parameter space. We estimate the evidence using the publicly available MultiNest (Feroz et al., 2009) code,222This package is based on the Nested sampling algorithm (Skilling, 2004). with Python interface PyMultinest (Buchner et al., 2014), for the efficient evaluation of the evidence integral (32) with likelihood (3.1). The accuracy of the Bayesian evidence estimate is controlled by the number of ‘live’ points, n_{\lower 1.0pt\hbox{\scriptstyle\rm live}}, with an error \sigma\sim O(n^{-1/2}_{\lower 1.0pt\hbox{\scriptstyle\rm live}}). In our analysis we choose 1000 points for the 8 or 9 parameter base model and add 100 more points for each additional parameter.

3.2 Choice of priors

Given the sensitivity of the Bayes factor to priors it is important these are chosen as objectively as possible. The choice of priors are summarized in Table 1.

3.2.1 Nuisance parameters

All nuisance parameters are common to both timescape and $\Lambda$CDM models and we therefore assign the same priors to both models. Where possible,333Given the complications introduced by empirical changes (4) to x_{\lower 1.0pt\hbox{\scriptstyle 1,0}}, c_{\lower 1.0pt\hbox{\scriptstyle 0}}, we adopt uniform priors for for these parameters. we adopt priors that have been used in previous Bayesian studies of the SALT2 method (March et al., 2011; Shariff et al., 2016). The standard deviations \{\sigma_{\lower 1.0pt\hbox{\scriptstyle x_{1,0}}},\sigma_{\lower 1.0pt\hbox{\scriptstyle c_{0}}},\sigma_{\lower 1.0pt\hbox{\scriptstyle M_{\rm B,0}}}\} are ‘scale’ parameters (of the residuals) and so it is more appropriate to assign a log-uniform prior to these parameters. The priors for the nuisance parameters are wide to ensure the most likely regions of parameter space are supported, and provided they are wide enough, this will have no overall effect on the Bayes factor (as the evidence of each model will be similarly scaled).

3.2.2 Cosmological parameters

Only one free cosmological parameter can be constrained by supernovae: \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} for spatially flat $\Lambda$CDM or $f_{\rm v0}$ for the timescape model. Conventionally, the combination of \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} and H_{\lower 1.0pt\hbox{\scriptstyle 0}} for the standard cosmology is strongly constrained by the CMB acoustic peaks (Ade et al., 2016). Measurements of the Baryon Acoustic Oscillation (BAO) scale in galaxy clustering statistics (Aubourg et al., 2015; Alam et al., 2017) at low redshifts and the Lyman $\alpha$ forest (Delubac et al., 2015; Aubourg et al., 2015) provide independent constraints. In the case of the timescape model, however, our ability to model the CMB is still limited by systematic uncertainties of 8–13% (Nazer & Wiltshire, 2015).

We therefore determine priors for $f_{\rm v0}$ in the timescape model using best present knowledge. For the CMB we use results of a model-independent analysis of the acoustic peaks (Aghamousa et al., 2015) with Planck satellite data, and choose a prior from a 95% confidence fit of the angular scale of the sound horizon. To date BAO studies all implicitly assume the FLRW model, and do not yet provide an equivalent model independent constraint. We therefore adopt a prior using FLRW-model estimates of the angular diameter of the BAO scale, including the full range of values which are currently in tension (Delubac et al., 2015; Aubourg et al., 2015; Alam et al., 2017). We take generous 95% confidence limits determined by assuming that either the low redshift galaxy clustering results (Aubourg et al., 2015; Alam et al., 2017) or the $z=2.34$ Lyman-$\alpha$ results (Delubac et al., 2015) are correct. Priors for the spatially flat $\Lambda$CDM model are determined by an identical methodology. Further details are given in Appendix D.

4 Results

4.1 Analysis with supernova parameters drawn from global Gaussian

distributions

Since there is a degeneracy between the Hubble constant, H_{\lower 1.0pt\hbox{\scriptstyle 0}}, and the magnitude, M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}}, we fix H_{\lower 1.0pt\hbox{\scriptstyle 0}} for each model. The value of M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}} then depends on this choice. We are then left with one free cosmological parameter, the matter density parameter \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} in the spatially flat $\Lambda$CDM model, and the present epoch void fraction $f_{\rm v0}$ in the timescape model. We can alternatively define an effective “dressed matter density parameter” \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}={\textstyle{1\over 2}}(1-f_{\rm v0})(2+f_{\rm v0}) (Wiltshire, 2007a, 2009), which takes similar numerical values to the concordance $\Lambda$CDM model, allowing likelihood functions to be plotted on the same scale. (This parameter does not obey the Friedmann equation sum rule, however.) The 9 parameters \{\Omega_{\lower 1.0pt\hbox{\scriptstyle M0}},\alpha,x_{\lower 1.0pt\hbox{\scriptstyle 1,0}},\sigma_{\lower 1.0pt\hbox{\scriptstyle x_{\rm 1,0}}},\beta,c_{\lower 1.0pt\hbox{\scriptstyle 0}},\sigma_{\lower 1.0pt\hbox{\scriptstyle c_{\rm 0}}},M_{\lower 1.0pt\hbox{\scriptstyle\rm B,0}},\sigma_{\lower 1.0pt\hbox{\scriptstyle M_{\rm B,0}}}\} are then fit for each model by determining the likelihood function with all parameters other than \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} (or $f_{\rm v0}$) treated as nuisance parameters. The empty universe has 8 parameters since \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}=0.

4.1.1 Statistical homogeneity scale cuts

An important systematic issue in the timescape cosmology is the fact that an average expansion law only holds on scales greater than the statistical homogeneity scale (SHS) \mathop{\hbox{{\lower 3.8pt\hbox{$>$}}\atop{\raise 0.2pt\hbox{$\sim$}}}}70–$120\,h^{-1}\hbox{Mpc}$ (Hogg et al., 2005; Scrimgeour et al., 2012). This corresponds to a CMB rest frame redshift of order $z\mathop{\sim}\limits 0.023$–$0.04$. In fact, SneIa analyses using the MLCS method have typically excluded SneIa below a cutoff at $z=0.024$ (Riess et al., 2007). However, the JLA catalogue (Betoule et al., 2014) includes $53$ SneIa, with $z<0.024$.

Following Smale & Wiltshire (2011) we determine cosmological model distances in the CMB frame, but make a redshift cut at the SHS, taken at $\mathop{\sim}\limits 100\,h^{-1}\hbox{Mpc}$. Furthermore, to examine the effect of the SHS cut on the fit of light-curve parameters, we perform the entire analysis while progressively varying the minimum redshift in the range $0.01\leq z_{\rm min}<0.1$; i.e., up to a redshift 3 times larger than the SHS. Systematic effects associated with the SHS can then be revealed. Our key results will be quoted for a cut at z_{\lower 1.0pt\hbox{\scriptstyle\rm SHS}}=0.033 in the CMB rest frame. The best-fitting MLE parameters with $z_{\rm min}=0.033$ are presented in Table 2.

For the priors given in Table 1 the Bayesian evidence in favour of the timescape model relative to the spatially flat $\Lambda$CDM model is $\ln B=0.085\pm 0.01$ with $z_{\rm min}=0.033$, or $\ln B=0.600\pm 0.007$ with $z_{\rm min}=0.024$. Since $|\ln B|<1$ the two models444Both models have positive relative Bayesian evidence compared to the empty model. Although the evidence is not particularly strong, |\ln B_{\lower 1.0pt\hbox{\scriptstyle 2}}|\mathop{\sim}\limits 2.2 incorporates priors which demand standard recombination for both $\Lambda$CDM and timescape. By that criterion the empty model is simply ruled out. are statistically indistinguishable. This conclusion is insensitive to O(1$\sigma$) changes to the width of the uniform priors on $f_{\rm v0}$ and \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}, or to variations of the minimum redshift as shown in Fig. 2(e).

While the Bayes factors do not show significant variation with $z_{\rm min}$, the values of particular best-fitting light-curve parameters show a marked change at the SHS. As shown in Fig. 2, there is a marked 30% jump in c_{\lower 1.0pt\hbox{\scriptstyle 0}} as $z_{\rm min}$ is varied from $0.01$ up to $z\simeq 0.033$, when compared to the subsequent fluctuations if $z_{\rm min}$ is increased up to $0.1$. For x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} there is a similar jump, although a linear trend remains in the range $0.033<z_{\rm min}<0.1$. The parameter $\beta$ parameter shows a small (3%) jump up to the SHS followed by $\pm 1$% fluctuations, while $\alpha$ remains relatively constant, fluctuating by $\pm 2$% over the whole range.

Since the light-curve parameters are remarkably close for all three cosmologies while showing a jump as the SHS emerges (distinct from the residual c_{\lower 1.0pt\hbox{\scriptstyle 0}} trend for the empty model with z_{\rm min}\mathop{\hbox{{\lower 3.8pt\hbox{$>$}}\atop{\raise 0.2pt\hbox{$\sim$}}}}0.05) there is clear evidence for some systematic effect at precisely the scale where we expect it.

4.2 Analysis with linear redshift variation for mean stretch and

colour parameters

Although RH16 considered four distinct subsamples, the mean stretch parameter actually shows a global increasing trend in the $\Lambda$CDM case evident in (Rubin & Hayden, 2016, Fig. 1, left panels). Our procedure of varying the minimum redshift cut on the whole sample also isolates any global trend. Such a trend is indeed evident in Fig. 2(c) beyond the SHS, with $x_{1,0}$ increasing 40% as the minimum sample redshift increases from $z_{\rm min}=0.033$ to $z_{\rm min}=0.1$. Beyond $z_{\rm min}=0.034$ an equivalent global trend in the mean colour parameter, c_{\lower 1.0pt\hbox{\scriptstyle 0}}, is not evident in Fig. 2(d), however, except in the case of the empty universe, which shows a 13% decrease in $c_{0}$ between $z_{\rm min}=0.034$ and $z_{\rm min}=0.1$. A global shift in $x_{1,0}$ would seem more consistent with an astrophysical systematic in the source population, rather than sampling biases with different thresholds for different samples.

To fully understand the differences that arise on making the RH16 changes (4), we have also investigated the effect of adding a smaller number of free parameters, by considering linear $z$ relations in just one of the parameters x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} or c_{\lower 1.0pt\hbox{\scriptstyle 0}}, and the difference between global linear relations and a split by subsamples. The advantage of our fully Bayesian approach is that not only can we compare the relative Bayesian evidence for different cosmological models with the same light-curve parameters, but we can also compare the merits of different empirical light-curve models. The values of the Bayesian evidence are shown in Table 3, along with a selection of parameters. The changes to the parameters $\alpha$ and $\beta$ are negligible between models, and are not tabulated. We do not tabulate all additional (up to 12) parameters for the case of the split subsamples, but an average.

4.2.1 Stretch parameter x_{\lower 1.0pt\hbox{\scriptstyle 1,0}}

Consistent with remarks above, relative to the baseline NGS16 model I, light curve model II provides positive (but not strong) Bayesian evidence for a global linear trend in x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} independent of cosmological model, with \ln B_{\lower 1.0pt\hbox{\scriptstyle 1}}=1.28, $1.34$, $1.74$ for the timescape, $\Lambda$CDM and empty models respectively. The BIC evidence for the same conclusion is very strong (timescape, $\Lambda$CDM models) or strong (empty model). By contrast model III gives strong evidence |\ln B_{\lower 1.0pt\hbox{\scriptstyle 1}}|>13 against a split linear law in x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} independent of cosmological model. The Bayesian penalty for introducing new empirical parameters depends on the choice of the priors, but our conclusion is robust to reasonable changes. Furthermore, the frequentist BIC statistic \Delta_{\lower 1.0pt\hbox{\scriptstyle\rm BIC}} also strongly disfavours model III relative to models I, II in the $\Lambda$CDM and timescape cases.

4.2.2 Colour parameter c_{\lower 1.0pt\hbox{\scriptstyle 0}}

In contrast to the stretch parameter, results involving the colour parameter vary greatly with cosmological model. Despite model IV having the global minimum BIC statistic for timescape, \ln B_{\lower 1.0pt\hbox{\scriptstyle 1}} shows no significant evidence555The empty universe has marginal evidence, consistent with Fig. 2(d) for $z_{\rm min}>0.05$. for any global linear redshift law. Relative Bayesian evidence for a split linear law in c_{\lower 1.0pt\hbox{\scriptstyle 0}} is marginal for timescape, positive for the empty universe, and strong for $\Lambda$CDM, with \ln B_{\lower 1.0pt\hbox{\scriptstyle 1}}= -0.91, -1.66 and -3.47 respectively.

The original RH16 model VII suffers similar problems to model III in terms of Bayesian evidence, evidently on account of the split linear law in x_{\lower 1.0pt\hbox{\scriptstyle 1,0}}. However, model VIII has the strongest Bayesian evidence of all models. It adds a global linear redshift law in x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} to model V. The improvement in Bayesian evidence for model VIII relative to model V is marginal for timescape and the empty universe, and positive for $\Lambda$CDM, with \Delta\ln B_{\lower 1.0pt\hbox{\scriptstyle 1}}=-0.74, $-0.89$ and $-1.18$ respectively.

Although \ln B_{\lower 1.0pt\hbox{\scriptstyle 2}} for model VIII gives positive (but not strong) relative evidence for $\Lambda$CDM over timescape, any conclusion drawn from this depends on additional empirical light-curve parameters which now depend on the cosmological model666Some of the largest differences occur in the SNLS subsample: c_{\lower 1.0pt\hbox{\scriptstyle\rm 0,\hbox{$\scriptstyle\Lambda$}CDM}}=0.0483 and c_{\lower 1.0pt\hbox{\scriptstyle\rm 0,TS}}=0.0565, a 17% difference. For the NGS16 model, by contrast, differences are 2%.. Furthermore, since the maximum likelihoods are comparable, the difference in Bayesian evidence is primarily due to the timescape maximum likelihood being driven to the unphysical limit $f_{\rm v0}\to 1$ for any light-curve model with linear variations in c_{\lower 1.0pt\hbox{\scriptstyle 0}}, as is shown in Fig. 3, which compares likelihoods in \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} (or $f_{\rm v0}$) for the NGS16 and RH16 models for two choices of $z_{\rm min}$.

Very similar results were found by Smale & Wiltshire (2011) in applying SALT2 without the methodology of NGS16, leading to a large discrepancy in the predictions of the SALT2 and MLCS2k2 fitters for timescape. Since direct application of the NGS16 methodology to the JLA catalogue agrees with some previous MLCS2k2 fits to smaller data sets (Leith et al., 2008; Smale & Wiltshire, 2011), we conclude that systematics similar to linear redshift variations in c_{\lower 1.0pt\hbox{\scriptstyle 0}} may be the key to earlier discrepancies.

4.2.3 Cosmological model dependency of linear redshift changes to

SALT2 methodology

To understand the origin of such differences consider the Taylor series expansions (23)–(25) for the timescape, spatially flat $\Lambda$CDM and empty universe models, as given in Appendix A. Leading coefficients for (23) and (24) are shown graphically in Fig. 4 as a function of the free cosmological parameter.

All cosmologies show improvement to a global increase in x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} with redshift and piecewise decreases in c_{\lower 1.0pt\hbox{\scriptstyle 0}}, including the empty model which has no free parameter to adjust. However, if we incorporate linear corrections $x_{z}z$ to x_{\lower 1.0pt\hbox{\scriptstyle 1,0}}, or $c_{z}z$ to c_{\lower 1.0pt\hbox{\scriptstyle 0}}, in the SALT2 relation, then the difference of (2) and (23)–(25) gives a potential degeneracy between empirical parameters $x_{z}$ or $c_{z}$ and changes in the free cosmological parameter if the linear term in (23)–(25) can be changed without greatly altering the next most important O($z^{2}$) term. Such a possibility is admitted by $\Lambda$CDM but not timescape.

For $\Lambda$CDM, the O($z$) term in (24) is linear in \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}, and the O($z^{2}$) term is quadratic in \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} with a minimum at \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}={\textstyle{8\over 27}}=0.296. For model V with split linear redshift laws in c_{\lower 1.0pt\hbox{\scriptstyle 0}} only, the best-fitting \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} coincides precisely with this minimum. The decrease in \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} by adding a global linear $z$ dependence to x_{\lower 1.0pt\hbox{\scriptstyle 1,0}} is approximately the same, \Delta\Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}=-0.01, in going from model V to VIII, or from model I to II. The difference in (24) between models I and VII/VIII,

[TABLE]

is dominated by the linear redshift changes, with negligible changes in the O($z^{2}$) term.

By contrast the terms in the Taylor series (23) for timescape are very slowly varying monotonic functions of $f_{\rm v0}$ on the range $0.6<f_{\rm v0}\leq 1.0$ (as shown in Fig. 4), so changes in \mu_{\lower 1.0pt\hbox{\scriptstyle\rm TS}} are much more constrained. The difference in (23) between models I and VII/VIII, is

[TABLE]

A large change in $f_{\rm v0}$ is required make changes in \mu_{\lower 1.0pt\hbox{\scriptstyle\rm TS}} comparable to the $\Lambda$CDM case, and the effect of increasing $f_{\rm v0}$ increases both the O($z$) and O($z^{2}$) terms. As seen in Fig. 5, the likelihood is consequently peaked along narrow ridges in the (x_{\lower 1.0pt\hbox{\scriptstyle z}},$f_{\rm v0}$) and (c_{\lower 1.0pt\hbox{\scriptstyle z,J}},$f_{\rm v0}$) planes, with almost constant values of x_{\lower 1.0pt\hbox{\scriptstyle z}} and c_{\lower 1.0pt\hbox{\scriptstyle z,J}} and no upper bound on $f_{\rm v0}$.

5 Discussion

Our study has a number of important consequences. Firstly, the timescape and spatially flat $\Lambda$CDM model luminosity distance–redshift fits to the JLA catalogue are statistically indistinguishable using either the approach of NGS16, or with modifications to only the mean stretch parameter. As shown in Table 3 the Bayesian complexity, $C_{b}$, is lower (better) for timescape than for $\Lambda$CDM, for every choice of light-curve model.

This completely reframes a debate (Nielsen et al., 2016; Rubin & Hayden, 2016; Shariff et al., 2016; Haridasu et al., 2017; Tutusaus et al., 2017) about whether cosmic acceleration is marginal or not, within the confines of a FLRW expansion history. Current supernova data does not distinguish between the standard $\Lambda$CDM model and the non-FLRW expansion history of the timescape model, which has non-zero apparent cosmic acceleration but with a marginal amplitude. The apparent deceleration parameter (22) for the best-fitting value of Table 2 is q_{\lower 1.0pt\hbox{\scriptstyle 0}}\equiv q(f_{\rm v0})=-0.043^{+0.004}_{-0.000}.

Within the class of FLRW models the significance of cosmic acceleration is often assessed by comparison to the empty universe model. That model is unphysical, since standard nucleosynthesis and recombination can never occur in a universe with $a(t)\propto t$ regardless of its matter content.777In particular, the $R_{h}=ct$ model is unphysical for this reason (Lewis, Barnes & Kaushik, 2016). The timescape model has positive \ln B_{\lower 1.0pt\hbox{\scriptstyle 2}} compared to the empty universe.888This is true for the NGS16 model I and all light-curve models for which \ln B_{\lower 1.0pt\hbox{\scriptstyle 1}} shows an improvement independent of cosmology, viz. models II, V, VIII. Nonetheless, |\mu_{\lower 1.0pt\hbox{\scriptstyle\rm TS}}(z)-\mu_{\lower 1.0pt\hbox{\scriptstyle\rm empty}}(z)|<|\mu_{\lower 1.0pt\hbox{\scriptstyle\Lambda{\rm CDM}}}(z)-\mu_{\lower 1.0pt\hbox{\scriptstyle\rm empty}}(z)| (c.f. Fig. 1) at late epochs, for a simple physical reason. The timescape model is void dominated at $z<1$, and the expansion of individual voids is close to an empty universe. While the timescape model has apparent acceleration at late epochs, its expansion law is closer to that of an empty universe than that of the $\Lambda$CDM model.

The second important consequence of our study is that allowing linear changes with redshift to the mean colour parameter, c_{\lower 1.0pt\hbox{\scriptstyle 0}}, produces cosmological model dependency. Since the redshift–distance relation of the timescape model effectively interpolates (Wiltshire, 2009, 2014) between those of spatially flat $\Lambda$CDM models with different values999Note that the Planck best-fitting value \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}=0.3175 (Ade et al., 2016) is lower than the best-fitting value for the spatially flat $\Lambda$CDM model value \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}=0.365 from Table 2, consistent with the timescape expectation. of \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}, particular care must be taken with piecewise linear relations in redshift.

The improved 16 parameter model VIII (this being a better fit than the original 21 parameter RH16 model) has positive Bayesian evidence for $\Lambda$CDM relative to the timescape model. However, this is contingent on degeneracies in the likelihood function between the free cosmological parameter and additional empirical parameters. The RH16 parametrization allows the $\Lambda$CDM deceleration parameter q_{\lower 1.0pt\hbox{\scriptstyle 0}}=-1+{\textstyle{3\over 2}}\Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} contained in the O($z$) term of the Taylor series (24) to be adjusted101010For the NGS16 model I and models VII/VIII one has best fits q_{\lower 1.0pt\hbox{\scriptstyle 0}}=-0.453 and q_{\lower 1.0pt\hbox{\scriptstyle 0}}=-0.571 respectively. The respective spatially flat $\Lambda$CDM values quoted by Rubin & Hayden (2016), namely q_{\lower 1.0pt\hbox{\scriptstyle 0}}=-0.412, $-0.552$, (or \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}=0.392, $0.299$), differ mostly on account of our SHS cut at $z_{\rm min}=0.033$. near the global minimum \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}=0.296 of the O($z^{2}$) term in (24). However, the same procedure drives the timescape free parameter to an unphysical limit, $f_{\rm v0}\to 1$. No fundamental model underlies the empirical parametrization (4). Variations in c_{\lower 1.0pt\hbox{\scriptstyle 0}} are most plausibly related to selection effects, given we cannot fit them by a global law. However, selection effects would be more correctly modelled by removing the tail of a Gaussian distribution rather than shifting its mean linearly in redshift.

Our results show that NGS16 did not account for every possible selection bias that remains in the JLA catalogue, consistent with some comments of RH16. Nonetheless, NGS16 are correct to point out the possible pitfalls in fitting SneIa data when empirical light-curve parameters are mixed with the cosmological parameters of a single class of cosmological models. If SneIa are to be used to distinguish cosmological models, then systematic uncertainties and selection biases should be corrected in as model independent manner as possible before the data is reduced.

A related issue which remains to be explored is the extent to which the corrections for selection biases that have already been made in the JLA catalogue depend on the FLRW model. Betoule et al. (2014) follow a procedure of Mosher et al. (2014, Sec. 6.2), who used the SNANA package to estimate selection biases. Simulated data (using the FLRW model) is used in such estimates. While efforts have been made to consider different dark energy equations of state (Mosher et al., 2014), models which do not satisfy the Friedmann equation fall outside the scope of such analyses.

Whether or not the timescape model is ultimately a better fit than the standard FLRW model, it may provide a useful diagnostic tool in comparing methods for SneIa light-curve reduction purely at the empirical level. In particular, it has an analytic non-FLRW redshift–distance relation which is very close to that of the $\Lambda$CDM model, but which is considerably more constrained in the free parameter $f_{\rm v0}$ than the $\Lambda$CDM model is in \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}.

Finally, Figs. 2(b),(c),(d) show evidence for a $\simeq 100\,h^{-1}\hbox{Mpc}$ statistical homogeneity scale which has an effect on global fits of light-curve parameters – most notably a 30% shift of c_{\lower 1.0pt\hbox{\scriptstyle 0}} – independent of the cosmological model. These systematics, which occur at a scale relevant from independent observations (Hogg et al., 2005; Scrimgeour et al., 2012), must be explained irrespective of the cosmological model.

Acknowledgements

We thank Chris Blake, Tamara Davis, Ahsan Nazer, David Rubin, Subir Sarkar and Bonnie Zhang for helpful discussions and correspondence, and Thomas Buchert for hospitality at the ENS, Lyon, France. This work was supported by Catalyst grant CSG-UOC1603 administered by the Royal Society of New Zealand.

Code availability: The code and data used in this analysis are available at https://doi.org/10.5281/zenodo.831360 .

Appendix A Luminosity distances in the FLRW and timescape cosmologies

We compare SneIa observations to distance moduli (1) for theoretical luminosity distances determined from the FLRW and timescape models. Regardless of the matter content of the universe, the distance modulus for any general FLRW model can be expanded as a Taylor series of derivatives of the scale factor $a(t)$ for small redshifts, $z$. This leads to a distance modulus (Visser, 2004)

[TABLE]

where ${\tt c}$ is the speed of light, and H_{\lower 1.0pt\hbox{\scriptstyle 0}}, q_{\lower 1.0pt\hbox{\scriptstyle 0}}, j_{\lower 1.0pt\hbox{\scriptstyle 0}}, s_{\lower 1.0pt\hbox{\scriptstyle 0}} and \Omega_{\lower 1.0pt\hbox{\scriptstyle k0}} are the present values of the Hubble, deceleration, jerk, snap and spatial curvature parameters: $H(t)\equiv a^{-1}\partial_{t}a$; $q(t)\equiv-a^{-1}H^{-2}\partial_{t}^{2}a$; $j(t)\equiv a^{-1}H^{-3}\partial_{t}^{3}a$; $s(t)\equiv a^{-1}H^{-4}\partial_{t}^{4}a$; \Omega_{\lower 1.0pt\hbox{\scriptstyle k}}(t)\equiv-k{\tt c}^{2}(Ha)^{-2}.

The luminosity distance-redshift relation in the $\Lambda$CDM model is given exactly by

[TABLE]

where \Omega_{\lower 1.0pt\hbox{\scriptstyle R0}}, \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}} and \Omega_{\lower 1.0pt\hbox{\scriptstyle\Lambda 0}} are respectively the present epoch values of the radiation, non-relativistic matter and cosmological constant density parameters, which at all epochs obey the Friedmann equation sum rule \Omega_{\lower 1.0pt\hbox{\scriptstyle R}}+\Omega_{\lower 1.0pt\hbox{\scriptstyle M}}+\Omega_{\lower 1.0pt\hbox{\scriptstyle k}}+\Omega_{\lower 1.0pt\hbox{\scriptstyle\Lambda}}=1. Since \Omega_{\lower 1.0pt\hbox{\scriptstyle R0}}=4.15\times 10^{-5}h^{-2}, the radiation term can be neglected on the distance scales of supernovae. Furthermore, for FLRW models \Omega_{\lower 1.0pt\hbox{\scriptstyle k0}} is constrained to be close to zero by the angular diameter distance of the sound horizon in the CMB. Thus we will restrict attention to the spatially flat $\Lambda$CDM model, with two effective free parameters, H_{\lower 1.0pt\hbox{\scriptstyle 0}} and \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}\simeq 1-\Omega_{\lower 1.0pt\hbox{\scriptstyle\Lambda 0}}. We use eq. (12) with \Omega_{\lower 1.0pt\hbox{\scriptstyle k0}}=0, \Omega_{\lower 1.0pt\hbox{\scriptstyle R0}}=0, for computations but note that in the Taylor series (11), q_{\lower 1.0pt\hbox{\scriptstyle 0}}=-1+{\textstyle{3\over 2}}\Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}, j_{\lower 1.0pt\hbox{\scriptstyle 0}}=1, s_{\lower 1.0pt\hbox{\scriptstyle 0}}=1-{\textstyle{9\over 2}}\Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}.

We also consider the FLRW model with linear expansion $a(t)\propto t$. This solution is obtained by setting \Omega_{\lower 1.0pt\hbox{\scriptstyle k0}}=1, \Omega_{\lower 1.0pt\hbox{\scriptstyle R0}}=\Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}=\Omega_{\lower 1.0pt\hbox{\scriptstyle\Lambda 0}}=0 in (12) or \Omega_{\lower 1.0pt\hbox{\scriptstyle k0}}=1, q_{\lower 1.0pt\hbox{\scriptstyle 0}}=j_{\lower 1.0pt\hbox{\scriptstyle 0}}=s_{\lower 1.0pt\hbox{\scriptstyle 0}}=\dots=0 in (11). Following NGS16 we denote this the empty universe, but note any matter content is admissible as long as the luminosity distance is exactly d_{\lower 1.0pt\hbox{\scriptstyle L}}={\tt c}\,z(1+{\textstyle{1\over 2}}z)/H_{\lower 1.0pt\hbox{\scriptstyle 0}}.

The timescape model (Wiltshire, 2007a, b, 2009; Duley et al., 2013), does not evolve by the Friedmann equation, and its distance modulus does not yield a Taylor series that coincides with (11) beyond the leading Hubble term. Instead observables are determined by conformally matching radial null geodesics of the regional “finite infinity” geometry of observers in gravitationally bound systems to a statistical geometry determined by fitting a spherically symmetric metric to a solution (Wiltshire, 2007b, 2009; Duley et al., 2013) of the Buchert equations (Buchert, 2000, 2001).

For the purpose of supernova distance analysis, the radiation density parameter (though somewhat differently calibrated to the CMB (Duley et al., 2013)) is negligible at the present epoch. To an accuracy of 0.3% the expansion history at late epochs is then given analytically (Wiltshire, 2007b, 2009). The “dressed” luminosity distance, d_{\lower 1.0pt\hbox{\scriptstyle L}}, and angular diameter distance, d_{\lower 1.0pt\hbox{\scriptstyle A}}, are given by

[TABLE]

where the volume-average time parameter, $t$, is defined implicitly in terms of the redshift by

[TABLE]

b\equiv 2(1-f_{\rm v0})(2+f_{\rm v0})/(9f_{\rm v0}{\vphantom{\tilde{H}}\bar{H}}_{\lower 1.0pt\hbox{\scriptstyle!0}}), $f_{\rm v0}$ is the present epoch value of the void volume fraction,

[TABLE]

and {\vphantom{\tilde{H}}\bar{H}}_{\lower 1.0pt\hbox{\scriptstyle!0}} is the “bare Hubble constant” related to the observed Hubble constant by {\vphantom{\tilde{H}}\bar{H}}_{\lower 1.0pt\hbox{\scriptstyle!0}}=2(2+f_{\rm v0})H_{\lower 1.0pt\hbox{\scriptstyle 0}}/(4{f_{\rm v0}}^{2}+f_{\rm v0}+4). The parameter $t$ is related to the time parameter, $\tau$, measured by typical observers in bound structures by

[TABLE]

The effective dressed scale factor $a(\tau)$ is given by

[TABLE]

where ${\vphantom{\tilde{a}}\bar{a}}$ is the bare or volume-average scale factor and ${\vphantom{\tilde{\gamma}}\bar{\gamma}}$ is the phenomenological lapse function. These have simple analytic forms in terms of the volume-average time parameter, $t$, namely

[TABLE]

and

[TABLE]

respectively (Wiltshire, 2007b, 2009). The bare Hubble parameter, ${\vphantom{\tilde{H}}\bar{H}}\equiv\partial_{t}{\vphantom{\tilde{a}}\bar{a}}/{\vphantom{\tilde{a}}\bar{a}}$, and dressed Hubble parameter, $H\equiv\partial_{\tau}a/a$, are given respectively by ${\vphantom{\tilde{H}}\bar{H}}=(2+{f_{\rm v}})/(3t)$ and $H=\left(4{f_{\rm v}}^{2}+{f_{\rm v}}+4\right){\vphantom{\tilde{H}}\bar{H}}/[2(2+{f_{\rm v}})]$. The bare deceleration parameter, ${\vphantom{\tilde{q}}\bar{q}}\equiv-{\vphantom{\tilde{a}}\bar{a}}^{-1}{\vphantom{\tilde{H}}\bar{H}}{}^{-2}\partial_{t}^{2}{\vphantom{\tilde{a}}\bar{a}}$, is always positive. However, on account of the different time parameters (18) the dressed deceleration parameter inferred by observers in bound systems, $q\equiv-a^{-1}H^{-2}\partial_{\tau}^{2}a$, may change sign from positive to negative, indicating apparent acceleration. The dressed deceleration parameter is given by

[TABLE]

The onset of apparent acceleration is determined by a root of the cubic in ${f_{\rm v}}$ in the numerator of (22), and begins when ${f_{\rm v}}\simeq 0.587$.

One may substitute (13)–(15) in (1) and then invert (16) as a series in $z$ using an algebraic computing package to arrive at a low redshift Taylor series for the distance modulus, \mu_{\lower 1.0pt\hbox{\scriptstyle\rm TS}}, equivalent to (11) for the FLRW model. The first terms are given below, along with equivalent expressions for the spatially flat $\Lambda$CDM and empty universe models as determined from (11):

[TABLE]

Here the term \mu_{\lower 1.0pt\hbox{\scriptstyle 0}}(z)\equiv 25+5\log_{\lower 1.0pt\hbox{\scriptstyle\rm 10}}[{\tt c}z/(H_{\lower 1.0pt\hbox{\scriptstyle 0}}\,\hbox{Mpc})]=25+5\log_{\lower 1.0pt\hbox{\scriptstyle\rm 10}}(2997.9\,h^{-1})+5\log_{\lower 1.0pt\hbox{\scriptstyle\rm 10}}z is common to all models, the Hubble constant being H_{\lower 1.0pt\hbox{\scriptstyle 0}}=100\,h\,\hbox{km}\;\,\hbox{sec}^{-1}\,\hbox{Mpc}^{-1}.

Appendix B Implementation of the SALT2 method

The SALT relation (2) refers to the actual emitter (em) and observer (obs), but the luminosity distance relations (12) and (14) refer to ideal observers who determine an isotropic distance–redshift relation. Consequently, the theoretical relations (12) and (14) must be transformed to the frame involving the actually measured redshift $\hat{z}=(\lambda_{\rm obs}-\lambda_{\rm em})/\lambda_{\rm em}$ before the SALT relation is applied. The luminosity distance entering (1) is then

[TABLE]

where for each cosmological model, D(z)=d_{\lower 1.0pt\hbox{\scriptstyle L}}/(1+z)=(1+z)d_{\lower 1.0pt\hbox{\scriptstyle A}} is the (effective) comoving distance, and

[TABLE]

gives the measured redshift, $\hat{z}$, in terms of the cosmological redshift, $z$, the local Doppler redshifts of observer, $z^{\rm pec}_{\rm obs}$, and emitter, $z^{\rm pec}_{\rm em}$, and gravitational redshifts at the two locations, $z^{\phi}_{\rm obs}$ and $z^{\phi}_{\rm em}$.

For our observations, $\hat{z}$, is the heliocentric redshift as the Earth’s annual motion is averaged to the rest frame of the Sun. In the standard cosmology gravitational potential effects are assumed to be small, and the only relevant terms on the r.h.s. of (27) are assumed to be local boosts of order $v/c\mathop{\sim}\limits 10^{-3}$. This leads to 0.1% corrections to the luminosity distance which are often neglected. However, as noted by Calcino & Davis (2017) differences of 0.1% in d_{\lower 1.0pt\hbox{\scriptstyle L}} lead to order 1% corrections to cosmological parameters, which we have confirmed in our analysis.

In the timescape model, as in any inhomogeneous cosmology, expansion below the $\mathop{\sim}\limits 100\,h^{-1}\hbox{Mpc}$ SHS will generally differ from that of a global average geometry plus local boosts. Equivalently, very slowly varying time-dependent gravitational potentials also make a contribution to (27). Such terms encode non-kinematic differential expansion (Bolejko, Nazer & Wiltshire, 2016) from inhomogeneities below the SHS. Spatial variations in the term $z^{\phi}_{\rm em}$ may have significant consequences for interpreting the local “peculiar velocity field” for sources within \mathop{\hbox{{\lower 3.8pt\hbox{$<$}}\atop{\raise 0.2pt\hbox{$\sim$}}}}100\,h^{-1}\hbox{Mpc} of our location (Wiltshire et al., 2013; McKay & Wiltshire, 2016) but any net anisotropy on SneIa redshifts on larger scales should only make a small correction to the standard boost between the heliocentric and CMB frames. Indeed, it could be a source for a small systematic redshift uncertainty of the type considered by Calcino & Davis (2017). However, we do not investigate that possibility in the present paper as the RH16 empirical light-curve models we study are already very complex. The peculiar velocity and gravitational potential terms in (27) that we are unable to determine will be assumed to contribute to statistical uncertainties in measured redshifts only.

We therefore compute cosmological luminosity distances is the CMB rest frame, exclude data below the SHS, and study the effect of different choices for this cutoff.111111Since we do not constrain H_{\lower 1.0pt\hbox{\scriptstyle 0}}, we do not specifically investigate the relationship between light-curve parameters and determinations of the local Hubble constant, which have been discussed in the past (Jha et al., 2007; Conley et al., 2007; Hicken et al., 2009; Smale & Wiltshire, 2011; Wiltshire et al., 2013; McKay & Wiltshire, 2016). In the timescape model higher average values of H_{\lower 1.0pt\hbox{\scriptstyle 0}} are expected below the SHS. Furthermore, we apply the SALT2 relation in the heliocentric frame using the values tabulated in the JLA catalogue (Betoule et al., 2014), and calculate the corresponding CMB rest frame redshifts ourselves.121212We use the NASA/IPAC Extragalactic Database standard, $371\,\hbox{km}\;\,\hbox{sec}^{-1}$ in the direction $(\ell,b)=(264.14^{\circ},48.26^{\circ})$ (Fixsen et al., 1996).

We do not use the JLA tabulated CMB frame redshifts (Betoule et al., 2014) since in addition to our $z^{\rm pec}_{\rm obs}$ correction, these values also already include a correction, $z^{\rm pec}_{\rm em}$, for the peculiar velocity field (Hudson et al., 2004; Conley et al., 2011) of galaxies up to $z=0.071$, implicitly assuming the FLRW model.

Appendix C Model comparison

C.1 Frequentist approach

We are interested in the dependence of the likelihood (3.1) on the model parameters, $\Theta$. We write $\mathcal{L}(\Theta)\equiv\mathcal{L}(\mathrm{Data}|\Theta,\mathrm{M})$. We are interested in a subset of parameter-space \Theta_{\lower 1.0pt\hbox{\scriptstyle 1}}\subset\Theta, for which we construct a profile likelihood \mathcal{L}_{p}(\Theta_{\lower 1.0pt\hbox{\scriptstyle 1}})\equiv\max_{\Theta_{2}}[\mathcal{L}(\Theta)], where maximization is over the nuisance parameters \Theta_{2}=\Theta\setminus\Theta_{\lower 1.0pt\hbox{\scriptstyle 1}}. In our case \Theta_{\lower 1.0pt\hbox{\scriptstyle 1}} is usually the free parameter(s) of the cosmological model, and $\Theta_{2}$ the intrinsic supernovae parameters and the empirical parameters, $\alpha$, $\beta$, of the light-curve fitter.

Confidence regions for the parameters of interest are estimated from the coverage probability $p_{\rm cov}$ of a region in the $k$-dimensional slice of parameter space, k\equiv\mathop{\rm dim}\Theta_{\lower 1.0pt\hbox{\scriptstyle 1}}, given asymptotically by the integral

[TABLE]

where \mathcal{L}_{\rm max}\equiv\max_{\Theta_{\lower 1.0pt\hbox{\scriptscriptstyle 1}}}[\mathcal{L}_{p}(\Theta_{\lower 1.0pt\hbox{\scriptstyle 1}})], and $f_{\chi^{2}}(x,k)$ is the probability density function of a $\chi^{2}$ distributed variable with $k$ degrees of freedom. Having constructed confidence intervals from (28), one can compare nested models.

Since we wish to compare independent non-nested models,131313We note that only models II, IV, VI are extensions of the 9 parameter base model, i.e., model I is nested in II, IV and VI, while II and IV are nested models of VI. Model V is nested in model VIII. we need to go beyond the procedure of (28). The Akaike Information Criterion (AIC) (Akaike, 1974) and Bayesian Information Criterion (BIC) (Schwarz, 1978) are two widely used measures of the relative information loss for non-nested models, given respectively by

[TABLE]

where $k$ is the number of independent model parameters, and $N$ the number of data points fit. The AIC estimate of relative probability of minimal information loss for two models is P_{AIC}\equiv\exp\left[-{\textstyle{1\over 2}}\left({\rm AIC}_{\lower 1.0pt\hbox{\scriptstyle 1}}-{\rm AIC}_{\lower 1.0pt\hbox{\scriptstyle 2}}\right)\right], and similarly for BIC. The BIC gives a greater penalty for introducing additional parameters than AIC if $N\geq 8$. Differences of at least 2, 6 and 10 are considered to be respectively ‘positive’, ‘strong’ and ‘very strong’ evidence (Kass & Raftery, 1995) for the model with the lower IC value. Both tests reduce to a maximum likelihood comparison when $k_{1}=k_{2}$, as is the case for the timescape and spatially flat $\Lambda$CDM models.

C.2 Bayesian approach

The frequentist methods place emphasis on the maximum likelihood, which is of limited use. We therefore perform a fully Bayesian analysis of the JLA data set to determine the relative statistical support for each cosmological model, as well as for the introduction of additional redshift dependent light curve parameters.

Given data, $D$, and a model, $\mathrm{M}$, determined by a set of $n$ parameters $\Theta=(\theta_{1},\theta_{2},\ldots,\theta_{n})$, by Bayes theorem the posterior probability distribution, $p(\Theta|D,M)$, is given by

[TABLE]

where $\mathcal{L}(\Theta)\equiv p(D|\Theta,\mathrm{M})$ is the likelihood, $\pi(\Theta|\mathrm{M})$ is the prior distribution and $p(D|\mathrm{M})$ is the Bayesian evidence. The prior represents a subjective initial state of belief in the parameters based on external information or previous experiments, while the evidence is a normalization constant,

[TABLE]

to ensure the posterior is a probability distribution. It is independent of parameters and as such does not play a role in parameter estimation, but becomes important for model comparison.

Given two models, $\mathrm{M}_{1}$ and $\mathrm{M}_{2}$, for the same data $D$, the Bayes factor (Kass & Raftery, 1995)

[TABLE]

gives a measure for which model is more probable in view of the data. The Bayes factor implicitly applies the principle of Occam’s razor141414The AIC and BIC statistics also include a penalty using simple approximations to Bayesian methods which derive from different assumptions about the priors. The factor of two difference in the IC evidence scale (Kass & Raftery, 1995) reflects the factor of 2 multiplying $\ln\left(\mathcal{L}_{\rm max}\right)$ in the definitions (29), (30). with a penalty for adding extra parameters. This makes model selection natural in the Bayesian framework. Values of $B>1$ indicate preference for model 1, $B<1$ for model 2. On a standard scale, evidence with $|\ln B|<1$ is ‘not worth more than a bare mention’ (Kass & Raftery, 1995) or ‘inconclusive’ (Trotta, 2007), while $1\leq|\ln B|<3$, $3\leq|\ln B|<5$ and $|\ln B|\geq 5$ indicate ‘positive’, ‘strong’ and ‘very strong’ evidences respectively (Kass & Raftery, 1995).

In the Bayesian approach the nuisance parameters are marginalized over, i.e., integrated out from the posterior $p(\Theta|D,\mathrm{M})$. E.g., the marginal posterior of $\theta_{1}$ is obtained from the $n$-dimensional posterior by

[TABLE]

and from this 1-dimensional distribution parameter inferences can be made. The posterior mean value is given by

[TABLE]

and more generally

[TABLE]

for some parameter dependent quantity $f$. Credible intervals, or uncertainties in parameters, can also be obtained from the marginal posterior. E.g., a $68\%$ equal-tailed credible interval is defined in such a way that $(1-0.68)/2=0.16$ of the probability lies on either side of the interval.

In cases where the Bayes factor is close to unity giving no clear preference for either model, the Bayesian complexity (Spiegelhalter et al., 2002) can provide a secondary measure to the model selection process. It is defined as

[TABLE]

where

[TABLE]

is the Kullback-Leibler divergence measuring the information gain of the inference, and \widehat{D_{\lower 1.0pt\hbox{\scriptstyle\rm KL}}} is a point estimator evaluated at the posterior mean $\bar{\Theta}$ measuring the expected information gain:

[TABLE]

where we have used Bayes theorem in the second equality. As the data may not be able to constrain all parameters, the Bayesian complexity determines the effective number of parameters supported by the data. Thus for models with $|\ln B|<1$, the model with the lower $C_{b}$ indicates the simpler model and is therefore preferred. By defining the effective chi-squared $\chi^{2}(\Theta)\equiv-2\ln\mathcal{L}$ and invoking Bayes theorem (31), we can rewrite (37) as

[TABLE]

with $\overline{\chi^{2}}$ being the posterior mean of $\chi^{2}$.

Appendix D Cosmological model priors

We construct priors for the timescape model (Wiltshire, 2007a, b, 2009; Duley et al., 2013) based on CMB and BAO observations, to the best of our knowledge. We will also construct equivalent priors for the $\Lambda$CDM model based on the same assumptions. The resulting priors are wider than in conventionally assumed, but do not unfairly weight a Bayesian comparison by integrating the $\Lambda$CDM model likelihood function over a narrow parameter range as compared to the timescape case.151515If we were to use conventional narrower priors for $\Lambda$CDM then the timescape model is either unfairly advantaged or disadvantaged, depending on whether the maximum likelihood lies within the range of the narrower prior or not. For the NGS16 model, for example, this is not the case for the spatially flat $\Lambda$CDM model, and the timescape model is unfairly advantaged. For model VIII the situation is reversed.

D.1 CMB acoustic scale constraint

In the case of the CMB, a cosmology independent analysis of the angular scale and heights of the acoustic peaks has been undertaken by Aghamousa et al. (2015) from the Planck data. We use the information resulting from the angular scale of the acoustic peaks alone. The angular scale depends on the angular diameter distance of the sound horizon alone, which is constrained in the timescape model. By contrast, the relative peak heights are related to the baryon–to–photon ratio, \eta_{\lower 1.0pt\hbox{\scriptstyle B\gamma}}, and the spectral index, $n_{s}$, which are parameters with the largest systematic uncertainties in the timescape case.

A non-parametric fit of the acoustic scale alone gives 286\leq\ell_{\lower 1.0pt\hbox{\scriptstyle\rm A}}\leq 305 at 95% confidence (Aghamousa et al., 2015). Our CMB prior is then determined by demanding that the angular diameter distance of the sound horizon at decoupling matches the corresponding angular scale \theta_{\lower 1.0pt\hbox{\scriptstyle\rm A}}=\pi/\ell_{\lower 1.0pt\hbox{\scriptstyle\rm A}}; i.e., 0.01030\leq\theta_{\lower 1.0pt\hbox{\scriptstyle\rm A}}\leq 0.01098. In earlier work (Leith et al., 2008; Smale & Wiltshire, 2011; Duley et al., 2013), given that non-parametric fits had not been performed, we had demanded a match to the FLRW parametric estimate of the acoustic scale \theta_{\lower 1.0pt\hbox{\scriptstyle\rm A}}=0.01041 to within 2%, 4% or 6%. The non-parametric fit represents a considerable improvement, particularly since the FLRW model value is not in the mid-range of the non-parametric 95% confidence interval.

To constrain the angular diameter distance of the sound horizon d_{\rm A\,dec}={\bar{D}}_{s}(z_{\rm dec})/\theta_{\lower 1.0pt\hbox{\scriptstyle\rm A}} in the timescape model, we determine the redshift of decoupling, $z_{\rm dec}$, and the comoving distance of the sound horizon $\bar{D}_{s}$ at that epoch (Duley et al., 2013; Nazer & Wiltshire, 2015), which require the baryon–to–photon ratio to be specified. In the FLRW model this ratio is very tightly constrained by the ratio of CMB peak heights, as first measured by WMAP (Bennett et al., 2003). However, to achieve a similarly precise constraint in the timescape model we would need to include backreaction in the primordial plasma (Nazer & Wiltshire, 2015), which is beyond the scope of current investigations. In previous work (Duley et al., 2013; Leith et al., 2008; Smale & Wiltshire, 2011) we used a range of pre-WMAP baryon–to–photon ratios (Tytler et al., 2000), 4.6<10^{10}\eta_{\lower 1.0pt\hbox{\scriptstyle B\gamma}}<5.6, for which all light element abundance measurements are within $2\sigma$, i.e., with no primordial lithium abundance anomaly. In the present analysis, we wish to use the same priors on \eta_{\lower 1.0pt\hbox{\scriptstyle B\gamma}} for both the timescape and $\Lambda$CDM models, and thus need to include the standard model value \eta_{\lower 1.0pt\hbox{\scriptstyle B\gamma}}=2.736\times 10^{-8}\Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}h^{2}=(6.08\pm 0.07)\times 10^{-10} for which the primordial lithium abundance is problematic. We therefore adopt the more conservative pre-WMAP range given by Olive et al. (2000), namely 4.2<10^{10}\eta_{\lower 1.0pt\hbox{\scriptstyle B\gamma}}<6.3.

D.2 Baryon Acoustic Oscillation constraints

Determinations of the BAO scale from galaxy clustering at low redshifts and Lyman alpha forest statistics at $z=2.34$ provide complementary constraints on the expansion history. In previous work (Duley et al., 2013; Leith et al., 2008; Smale & Wiltshire, 2011) we simply demanded that the timescape effective comoving BAO scale match a single estimate determined from the FLRW cosmology to within $\pm 2$%, $\pm 4$% or $\pm 6$%, which was a crude method but the best available given the earlier precision of measurements. The number and precision of measurements has now improved.

For the present investigation, we have considered estimates of the BAO scale at different redshifts (Delubac et al., 2015; Aubourg et al., 2015; Alam et al., 2017) using the best available data from the BOSS survey. Unfortunately the standard FLRW cosmology plays an implicit role in the data reduction, and limits the extent to which bounds can be placed on non-FLRW models. The systematic issues can be most directly understood by noting that the BAO scale is determined separately in the angular and radial directions, by converting angular separations and redshift separations in the galaxy–galaxy correlation function into the displacements

[TABLE]

where $r_{d}$ is the present comoving scale of the sound horizon at the baryon drag epoch, $d_{H}(z)\equiv c/H(z)$, and the subscript “fid” refers to quantities computed in a fiducial FLRW model that is used to convert the raw angular and redshift displacements into 3-dimensional comoving space. (Here we neglect the effect of redshift–space distortions which are also often modelled with $N$-body Newtonian simulations based on the $\Lambda$CDM model.)

The conversion to 3-dimensional comoving space can be problematic for a non-FLRW model. While use of purely angular results should pose no problems for the timescape cosmology, the conversion of redshift increments to a radial comoving distance involves different assumptions about spatial curvature in the FLRW and timescape models. One could in principle use the values determined by a fiducial $\Lambda$CDM model (Delubac et al., 2015; Aubourg et al., 2015; Alam et al., 2017) to recompute the radial comoving distance except for an additional problem: in particular redshift ranges the relative Alcock–Paczyński factor \left[\alpha_{\perp}/\alpha_{\parallel}\right]_{\lower 1.0pt\hbox{\scriptstyle\rm\Lambda CDM}}/\left[\alpha_{\perp}/\alpha_{\parallel}\right]_{\lower 1.0pt\hbox{\scriptstyle\rm TS}}=\left[H(z)d_{\lower 1.0pt\hbox{\scriptstyle A}}(z)\right]_{\lower 1.0pt\hbox{\scriptstyle\rm\Lambda CDM}}/\left[H(z)d_{\lower 1.0pt\hbox{\scriptstyle A}}(z)\right]_{\lower 1.0pt\hbox{\scriptstyle\rm TS}} between a fiducial $\Lambda$CDM model and the timescape model varies over the redshift slices $\Delta z\mathop{\sim}\limits 0.2$ used in the BOSS survey (Alam et al., 2017) by an amount similar in magnitude to the uncertainty. Consequently, to have any confidence in radial measurements, one really needs to recompute the radial BAO scale from the raw data assuming a fiducial timescape model. That is beyond the scope of the present paper.

For the present analysis we will consequently restrict constraints on the BAO scale to 2$\sigma$ bounds obtained from the angular estimates of BOSS data (Alam et al., 2017) at low redshifts 0.38\mathop{\hbox{{\lower 3.8pt\hbox{$<$}}\atop{\raise 0.2pt\hbox{$\sim$}}}}z\mathop{\hbox{{\lower 3.8pt\hbox{$<$}}\atop{\raise 0.2pt\hbox{$\sim$}}}}0.61 and in the Lyman-$\alpha$ forest (Delubac et al., 2015) at $z=2.34$. In the former case, the radial and angular measurements are actually also somewhat correlated. Consequently, and also in view of the fact that the measurements at low and high redshifts are in tension in the $\Lambda$CDM model, we will take bounds that result from the union of the constraints at low and high redshifts, rather than their intersection. In practice, the bounds are mostly set by the Lyman-$\alpha$ measurement since it has a much larger uncertainty.

D.3 Joint constraints

The joint 2$\sigma$ confidence regions from applying the CMB constraint and either the low-$z$ or $z=2.34$ BAO constraint is shown in Fig. 6. Since a range of possible baryon–to–photon ratios are admitted, with no information from the relative heights of the acoustic peaks used in either model, the width of the allowed regions is larger than in conventional analyses for $\Lambda$CDM.

For timescape the confidence regions are $f_{\rm v0}\in(0.588,0.765)$ at 1$\sigma$, $f_{\rm v0}\in(0.500,0.799)$ at 2$\sigma$, $f_{\rm v0}\in(0.378,0.826)$ at 3$\sigma$. The corresponding effective dressed \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}={\textstyle{1\over 2}}(1-f_{\rm v0})(2+f_{\rm v0}) is \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}\in(0.325,0.534) at 1$\sigma$, \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}\in(0.281,0.625) at 2$\sigma$, and \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}\in(0.245,0.740) at 3$\sigma$. For spatially flat $\Lambda$CDM the corresponding confidence regions are \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}\in(0.162,0.392) at 1$\sigma$, \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}\in(0.143,0.487) at 2$\sigma$, and \Omega_{\lower 1.0pt\hbox{\scriptstyle M0}}\in(0.124,0.665) at 3$\sigma$. We adopt the 2$\sigma$ bounds as priors.