# Cosmological foundations revisited with Pantheon+

Zachary G. Lane1, Antonia Seifert1,2, Ryan Ridden-Harper1, Jenny Wagner3, and
David L. Wiltshire1
1School of Physical and Chemical Sciences — Te Kura Matū, University of
Canterbury, Private Bag 4800, Christchurch 8140, New Zealand
2Institut für Theoretische Physik, Universität Heidelberg, Philosophenweg 12,
D-69120 Heidelberg, Germany
3Bahamas Advanced Study Institute and Conferences, 4A Ocean Heights, Hill View
Circle, Stella Maris, Long Island, The Bahamas
Corresponding author
ZGL and AS contributed equally to this work. 0009-0003-8380-4003
0009-0005-9892-3667 0000-0003-1724-2885 0000-0002-4999-3838
0000-0003-1992-6682

(Accepted XXX. Received YYY; in original form 2 November 2023)

###### Abstract

We reanalyse the Pantheon+ supernova catalogue to compare a cosmology with
non-Friedmann-Lemaître-Robertson-Walker (FLRW) evolution, the “timescape
cosmology”, with the standard spatially flat $\Lambda$CDM cosmology. To this
end, we consider the Pantheon+ supernova catalogue, which is the largest
available Type Ia supernova dataset for a geometric comparison between the two
models. We construct a covariance matrix to be as independent of cosmology as
possible, including independence from the FLRW geometry and peculiar velocity
with respect to FLRW average evolution. Within this framework, which goes far
beyond most other definitions of “model independence”, we introduce new
statistics to refine Type Ia supernova (SNe Ia) light-curve analysis. In
addition to conventional $n$-point galaxy correlation functions used to define
the scale of statistical homogeneity, denoted SHSgn statistics, we introduce
the empirical statistic SHSα which enables a refined analysis of the
distribution biases of SNe Ia light-curve parameters $\beta c$ and $\alpha
x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$. We find that the stretch parameter,
$x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$, is subject to significant changes
below the SHSα convergence scale, found at $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=0.075^{+0.007}_{-0.009}$. Beyond this
scale, we conclude that a systematic error present in the JLA catalogue – and
the subject of previous debates – has been successfully controlled in the
Pantheon+ analysis.

For lower redshifts, the Bayesian analysis highlights important features
attributable to the increased number of low-redshift supernovae, the artefacts
of model-dependent light-curve fitting and the cosmic structure through which
we observe supernovae. This indicates the further need for surveys with more
uniform redshift distributions and cosmology-independent data reduction to
conduct a stronger investigation of the emergence of statistical homogeneity
and to compare alternative cosmologies in light of recent challenges to the
standard model.

“Dark energy” is generally invoked as a place–holder for “new physics”. Our
from-first-principles reanalysis of the Pantheon+ catalogue supports future
deeper studies of the interplay of matter and nonlinear spacetime geometry, in
a data-driven setting. For the first time in 25 years, we find evidence that
the Pantheon+ catalogue already contains such a wealth of data that with
further reanalysis, a genuine “paradigm shift” may soon emerge.

###### keywords:

cosmological parameters — dark energy — cosmology: observations — cosmology:
theory

††pubyear: 2023††pagerange: Cosmological foundations revisited with
Pantheon+–3

##  1 Introduction

Modern cosmology is entering a new era characterised by orders of magnitude
increase in the volume, variety and precision of available data, from both
ground–based and space–based telescopes, from one–off missions to new
dedicated observatories operating over decades. We are now learning to grapple
not only with the entire electromagnetic spectrum, but also with multi-
messenger astronomy.

Our new windows on the Universe – high–energy neutrinos and gravitational
waves from binary mergers of neutron stars and black holes – offer glimpses of
potential new physics which cannot be easily guessed at from consideration of
only electromagnetism, nuclear/particle physics or gravity in isolation. Vital
clues to unravelling new mysteries occur not only when all three
multimessenger channels are detected, but also when some channels are detected
in coincidence but others are not. The absence of directly detected neutrinos
from GW170817 (Abbott et al., 2017) – to date the most closely studied single
astronomical event ever – is a case in point.

In contemplating new physics we face the challenge that unlike
electromagnetism which is linear, non–Abelian gauge theories (nuclear physics)
and general relativity (gravity) are each fundamentally nonlinear. This paper
presents a data–driven detailed statistical reanalysis of the best available
Type Ia supernovae (SNe Ia) data (Scolnic et al., 2022). However, rather than
seeking new physics amongst the nonlinearities of either particle gauge
theories or of general relativity (GR) alone, we will compare standard
cosmology with a model that combines the nonlinearities of both key
theoretical approaches.

This latter model, the timescape cosmology (Wiltshire, 2007a, 2009b),
considers the fundamental questions, such as “what is the origin of scale?”,
“where is infinity?”, and “how do we compare standard to standardisable rulers
and clocks?”. Through this lens, timescape can then be viewed as
simultaneously addressing both the physics of the “dark sector” and the
physics of “modified gravity” based on theoretical assumptions which are far
more conservative than those commonly proposed to resolve cosmological
tensions (Di Valentino et al., 2021).

SNe Ia were the crucial dataset that first convinced (Riess et al., 1998;
Perlmutter et al., 1999) a majority of the community of the (apparently)
accelerating rate of present epoch cosmic expansion. After more than two
decades SNe Ia are still one of the most widely investigated and crucial
astronomical sources for determining geometric distances in the Universe. At a
time when significant improvements in precision are now turning “tensions”
into anomalies (Peebles, 2022), we present a re-examination of SNe Ia.

In the present analysis, we revisit the relationship of the foundations of
geometry to the calibration of measurements in observational cosmology. In
particular, a common usage of the phrase “model independence” is broadened to
include independence from geometrical assumptions. Here we present our
technical results alongside the underlying theories, to address numerous
research communities engaged in studying the evolution of the Universe. We
will examine how model-dependent assumptions potentially impact inherent
biases in the standardisation of astrophysical data and model selection. This
paper thereby introduces critical philosophical and scientific questions that
the scientific community must address.

###  1.1 $\Lambda$CDM, Timescape and Inhomogeneous Cosmology

The Friedmann-Lemaître-Robertson-Walker (FLRW) metric is built on two main
assumptions, spatial isotropy and homogeneity, and it has been at the
forefront of modern cosmology since the 1930s. Isotropy and homogeneity are
supported by measurements of the Cosmic Microwave Background (CMB). The CMB is
observed to have a near-perfect blackbody spectrum up to fluctuations of order
$|\delta T|/T\mathop{\sim}\limits 10^{-3}$ (Corey & Wilkinson, 1976; Smoot et
al., 1977; Cheng et al., 1979). The leading dipole fluctuation, with $|\delta
T|=(3362.08\pm 0.99)\,$ $\mu$K (Planck Collaboration et al., 2020), is usually
interpreted kinematically as being purely due to our local motion with respect
to a rest-frame.

With the leading dipole subtracted, the remaining anisotropies, with amplitude
$|\delta T|/T\mathop{\sim}\limits 10^{-5}$ (Planck Collaboration et al.,
2020), are due to intrinsic density fluctuations on the surface of last
scattering, or to inhomogeneities from structures on our line of sight. The
small amplitude of the anisotropies is taken as a support for the Cosmological
Principle. While this is well justified on the surface of last scattering
itself, the late epoch Universe contains a vast cosmic web of complex
structures, requiring the fundamental understanding of average homogeneity and
isotropy to be more nuanced (Aluri et al., 2023).

The standard cosmology is based on a FLRW background metric whose spatial
curvature is found to be negligible (Planck Collaboration et al., 2020;
Pierpaoli et al., 2000), i.e., the line element

$\mathrm{d}s^{2}_{\lower 2.0pt\hbox{$\scriptstyle\rm
FLRW$}}=-\mathrm{d}t^{2}+a^{2}(t)\left(\mathrm{d}r^{2}+r^{2}\
\mathrm{d}\Omega^{2}\right)$ (1)

is adapted to comoving observers at fixed $(r,\theta,\phi)$ who measure cosmic
time, $t$, as their proper time. Here $\mathrm{d}\Omega^{2}\equiv{\rm
d}\theta^{2}+\sin^{2}\theta\,{\rm d}\phi^{2}$ is the metric on the celestial
2-sphere. The word background contains an implicit assumption of a top–down
ontology (Wiltshire, 2011): the average evolution is assumed to contain the
particular split of space and time implicit in the FLRW geometry regardless of
the scale to which such an average is applied.

For averages over very large spatial scales, linear perturbation theory on a
global FLRW background is assumed to still apply today, at late cosmic epochs.
When gravitational instability results in small-scale spatial structures
forming, the inhomogeneities of standard cosmology enter the nonlinear regime,
which must be treated by numerical simulations. Although simulations using
full GR have been developed recently (Giblin et al., 2016; Bentivegna & Bruni,
2016; Macpherson et al., 2017), most large numerical simulations apply
Newtonian gravity to a simulation volume whose expansion is fit to a
background spatially flat FLRW model.

In the standard approach, which envisages Hubble expansion as a global Hubble
flow, all deviations from its average can always and everywhere be described
by local Lorentz boosts, or peculiar velocities. Such a situation is not
typical of generic models in general relativity, including inhomogeneous dust
cosmologies, which may exhibit non-kinematic differential expansion (Bolejko
et al., 2016).

With the standard kinematic interpretation of the CMB dipole, deviations from
isotropy on the largest distance scales – e.g., in galaxy and quasar
distributions – are best characterised in terms of the special relativistic
aberration and modulation effects due to our local boost relative to the CMB
rest-frame (Ellis & Baldwin, 1984).

The kinematic interpretation is now challenged, however, at more than five
standard deviations (Secrest et al., 2022, 2021). Building on an earlier
result,111Secrest et al. (2021) use a sample of 1.36 million very distant
quasars (with a mean redshift of 1.2) alone to find 4.9$\sigma$ tension. While
the size of the quasar catalogue has increased, statistically the most
significant aspect of the Secrest et al. (2022) result are the correlated
features in the smaller independent radio galaxy catalogue. Since the Secrest
et al. (2022) catalogue extends to fainter fluxes, the median redshift should
be slightly greater than the mean stated by Secrest et al. (2021), and is not
relevant for this analysis. The radio galaxies of Secrest et al. (2022) have a
median redshift of 0.8 (Condon et al., 1998). Secrest et al. (2022) combine a
larger quasar catalogue of 1.6 million source with a sample of 508,144
not–quite–so distant radio galaxies to find a 5.1$\sigma$ disagreement with
the kinematic interpretation. This anomaly now exceeds the $5\sigma$ threshold
used to demarcate signals worthy of “new physics” attention in high–energy
physics accelerator experiments – as has been noted by careful readers
(Peebles, 2022).

To address the challenges of the Secrest et al. (2022) result, in this paper
we will broaden the definition of model–independence to include global
non–FLRW models which agree with the standard cosmology to the extent that:

  1. (i)

a notion of an isotropic average Hubble–Lemaître expansion law exists on
scales $z\gtrsim z_{\rm min}$, where the redshift $z_{\rm min}$ is determined
by a revised notion of statistical homogeneity scale (SHS);

  2. (ii)

on scales $z\lesssim z_{\rm min}$ the directly observed small–scale
inhomogeneous matter distribution is produced via gravitational instability
from primordial density fluctuations. Viability requires consistency with all
observed phenomenology of the CMB radiation on one hand, and with nuclear
physics models on the other. The latter have been well–established in
terrestrial accelerator experiments at energies below that of the quark–gluon
transition (Rafelski, 2020).

Criterion (i) is addressed by general models that incorporate backreaction of
inhomogeneities (Buchert, 2000, 2001; Buchert et al., 2020),222For general
reviews see Buchert & Räsänen (2012); Wiltshire (2014). namely additional
terms in the averaged Einstein equations governing cosmic evolution on large
scales. Such terms arise from the bottom up when coarse–graining Einstein’s
equations for structure on small scales.

Criterion (ii) is addressed by phenomenological models which calibrate
radiation and matter transfer functions from the early Universe to non-FLRW
evolution in the nonlinear regime333In contrast to electromagnetism, which is
linear, non-Abelian gauge theories and general relativity are both nonlinear.
As a gauge theory, electromagnetism involves a single complex scalar
potential, $\psi$, and the internal symmetry group $U(1)$. Furthermore, in GR
with the absence of torsion, covariant derivatives in Maxwell’s equations
reduce to ordinary exterior derivatives: ${\rm d}{\mathbf{F}}=0$, ${\rm
d}{*}{\mathbf{F}}=\mu_{\lower 2.0pt\hbox{$\scriptstyle 0$}}{\mathbf{J}}$,
where $\mathbf{F}={\rm d}{\mathbf{A}}$, $\mathbf{A}$ is the 4-vector
potential, ${\mathbf{A}}={\rm d}\psi$, $\mathbf{J}$ the 4-current, and ${*}$
denotes the Hodge dual (Szekeres, 2004). at late epochs (Nazer & Wiltshire,
2015; Rácz et al., 2017).

The timescape model uses a particular implementation of Buchert’s averaging
scheme with regional average metrics444Inspired by Ellis (1984), in timescape
cosmology finite infinity regions are defined as (marginally expanding)
boundaries of spatially flat regions containing bound structures. An
alternative related quantity is the matter horizon (Ellis & Stoeger, 2009).
The quantities $\tau_{\rm w,v},\eta_{\rm w,v},a_{\rm w,v}$ are not defined
globally but can be adapted to the clocks and rulers of ideal regional
observers. Quantities in finite infinity region, $\tau_{\rm w},\eta_{\rm
w},a_{\rm w}$, differ from their counterparts in void domains, ${\cal
D}_{\lower 2.0pt\hbox{$\scriptstyle\rm v$}}$, since finite infinity regions
average over the many sources of gravitational energy within the largest bound
structures (binding, rotational, thermal, …) whereas voids are dominated by
gradients in the kinetic energy of expansion. The timescape model assumes that
“dark energy” results from a cumulative misidentification and miscalibration
of gradients in the quasilocal kinetic energy of cosmic expansion as
structures grow. See Wiltshire (2011, 2014) for further introduction to these
ideas, and Dam et al. (2017); Buchert et al. (2020) for later updates.

$\mathrm{d}s^{2}_{\mathop{\hbox{\it fi}}}=-\mathrm{d}\tau_{\rm w}^{2}+a_{\rm
w}^{2}(\tau_{\rm w})\left(\mathrm{d}\eta_{\rm w}^{2}+\eta_{\rm w}^{2}\
\mathrm{d}\Omega^{2}\right),$ (2)

for finite infinity regions within the walls and filaments of the cosmic web,
and

$\mathrm{d}s^{2}_{{\cal D}_{\lower 2.0pt\hbox{$\scriptstyle\rm\rm
v$}}}=-\mathrm{d}\tau^{2}_{\rm v}+a_{\rm v}^{2}(\tau_{\rm
v})\left(\mathrm{d}\eta_{\rm v}^{2}+\sinh^{2}(\eta_{\rm v})\
\mathrm{d}\Omega^{2}\right),$ (3)

at void centres (Wiltshire, 2007a, 2008, 2009b).

Both of these metrics have FLRW form: Eq. 2 is identical to Eq. 1 and the void
metric Eq. 3 is a negative spatial curvature FLRW metric. However, in the
timescape model these metrics are not exact solutions of the Einstein
equations in general. Rather, they are regional best-fit effective average
metrics obtained in one particular solution (Wiltshire, 2007a, b) to the
fitting problem (Ellis, 1984; Ellis & Stoeger, 1987). The effective metrics
apply to small regions within which the Einstein equations hold to a good
approximation, and which are in turn embedded in the global average geometry
of the complex cosmic web.

Applying Buchert’s scalar averaging formalism to the geometry, Eqs. 2 and 3,
of an ensemble of wall and voids yields the bare cosmological parameters
(Wiltshire, 2007a, 2009b)

$\displaystyle{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\rm M}$
$\displaystyle\equiv\frac{8\pi G\bar{\rho}_{\lower 2.0pt\hbox{$\scriptstyle\rm
M0$}}{\bar{a}_{0}}^{3}}{3{\vphantom{\tilde{H}}\bar{H}}^{2}\bar{a}^{3}}\,,$
$\displaystyle{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\rm R}$
$\displaystyle\equiv\frac{8\pi G\bar{\rho}_{\lower 2.0pt\hbox{$\scriptstyle\rm
R0$}}{\bar{a}_{0}}^{4}}{3{\vphantom{\tilde{H}}\bar{H}}^{2}\bar{a}^{4}}\,,$ (4)
$\displaystyle{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\rm K}$
$\displaystyle\equiv\frac{-k_{\lower 2.0pt\hbox{$\scriptstyle\rm v$}}f_{{\rm
v}i}^{2/3}f_{\lower 2.0pt\hbox{$\scriptstyle\rm
v$}}^{1/3}}{\bar{a}^{2}{\vphantom{\tilde{H}}\bar{H}}^{2}}\,,$
$\displaystyle{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\cal Q}$
$\displaystyle\equiv\frac{-\dot{f}_{\rm v}^{2}}{9f_{\lower
2.0pt\hbox{$\scriptstyle\rm v$}}(1-f_{\lower 2.0pt\hbox{$\scriptstyle\rm
v$}}){\vphantom{\tilde{H}}\bar{H}}^{2}}\,,$ (5)

in terms of a volume-average comoving time, $t$; the average expanding volume,
${\cal V}\propto{\bar{a}}^{3}$ with $\bar{a}^{3}=f_{\lower
2.0pt\hbox{$\scriptstyle{\rm v}i$}}a_{\rm v}^{3}+(1-f_{\lower
2.0pt\hbox{$\scriptstyle{\rm v}i$}})a_{\rm w}^{3}$ on effective spatial
hypersurfaces, $t=\mathop{\hbox{\,const}}.$; the bare Hubble parameter,
$\bar{H}=\dot{\bar{a}}/\bar{a}$; the void volume fraction $f_{\lower
2.0pt\hbox{$\scriptstyle\rm v$}}(t)=f_{\lower
2.0pt\hbox{$\scriptstyle\mathrm{v}i$}}a_{\rm v}^{3}/\bar{a}^{3}$ with initial
value $f_{\lower 2.0pt\hbox{$\scriptstyle{\rm v}i$}}$ at last scattering;
$k_{\rm v}=\mathop{\hbox{\,const}}$; and the subscript $0$ denotes present
epoch values of the average energy densities of nonrelativistic matter,
$\bar{\rho}_{\lower 2.0pt\hbox{$\scriptstyle\rm M$}}$; radiation,
$\bar{\rho}_{\lower 2.0pt\hbox{$\scriptstyle\rm R$}}$; etc.

An average global Hubble-Lemaître expansion law parameterised by
${\vphantom{\tilde{H}}\bar{H}}(t)$ only applies on scales $\gtrsim
70$–$120\,h^{-1}\mathop{\hbox{\,Mpc}}$ (Hogg et al., 2005; Scrimgeour et al.,
2012), and the transition to homogeneity is conventionally defined via density
contrasts, $\delta_{\rho}\equiv(\rho-\bar{\rho})/\bar{\rho}$, relative to a
global FLRW background (Gabrielli et al., 2005). To circumvent this for non-
FLRW models, the Statistical Homogeneity Scale (SHS) is introduced as a
neologism555Inventing modern neologisms concerning the nature of space and
time involves new and different challenges to even the ones Einstein faced.
Further philosophical discussions, aside from the Conclusion (Section 6), are
beyond the scope of this paper and will be published elsewhere (Wiltshire in
prep.). Since extraordinary claims require extraordinary evidence, when faced
with a genuine paradigm shift, the responsibility of any scientist to actively
seek to “kill” their own theory is even greater. The directive given by DLW to
the coauthors of this article was to apply precisely such an approach. for the
scale determined by $n$-point galaxy correlation functions.

In Buchert averages, volume-average parameters defined as in Eqs. 5 and 4 obey
a sum-rule666An inhomogeneous ${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\Lambda$}}$ parameter, induced by a global
cosmological constant, can also be added mathematically (Buchert et al.,
2020). While of formal interest, it is not part of timescape phenomenology and
is not considered here. In addition, Buchert et al. (2020) also introduce
terms for dynamical backreaction,
${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\cal P$}}$, and stress-energy backreaction,
${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\cal T$}}$, in the full generalisation of Eq. 6.
Such terms are likely to be important in the very early Universe, but not at
late cosmic epochs when non-relativistic matter dominates, and when
$|{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\cal
P,T$}}|\ll|{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\cal Q$}}|<0.044$. The latter inequality is a
strict upper bound in the timescape model but not in generic backreaction
schemes. We also use ${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm K$}}$ rather than the alternative
${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\cal R$}}$ to avoid confusion with the radiation
density parameter, ${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm R$}}$, which is significant in the early Universe.

${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\rm
M$}}+{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\rm
R$}}+{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\rm
K$}}+{\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\cal Q$}}=1,$ (6)

similarly to the Friedmann equation. By contrast,
${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower 2.0pt\hbox{$\scriptstyle\rm
K$}}$ does not scale as a simple power of $\bar{a}$. Instead, average spatial
curvature is directly related to the void fraction:
${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower 2.0pt\hbox{$\scriptstyle\rm
K$}}\propto f_{\rm v}^{1/3}/(\bar{a}^{2}\bar{H}^{2})$. Here $f_{\rm v}\to 0$
at early times but dominates at late epochs. The kinematical backreaction term
${\vphantom{\tilde{\Omega}}\bar{\Omega}}_{\lower
2.0pt\hbox{$\scriptstyle\rm\cal Q$}}$ gives rise to non-FLRW evolution if
$\dot{f}_{\rm v}\neq 0$.

Volume–average time, $t$, only coincides with proper time for ideal isotropic
observers whose regional spatial curvature coincides with the global volume
average.

Since the Universe is dominated in volume by voids, such locations are
necessarily within voids, though not at their centres. Phenomenological lapse
functions,

${{\vphantom{\tilde{\gamma}}\bar{\gamma}}_{\rm w}}\equiv\frac{{\rm d}t}{{\rm
d}\tau_{\rm w}},\qquad{{\vphantom{\tilde{\gamma}}\bar{\gamma}}_{\rm
v}}\equiv\frac{{\rm d}t}{{\rm d}\tau_{\rm v}}$ (7)

relate the proper times $\tau_{\rm w,v}$ to volume-average time, $t$. For
observers and sources in bound structures that have broken away from the
Hubble expansion only the wall lapse,
${{\vphantom{\tilde{\gamma}}\bar{\gamma}}_{\rm w}}$, is directly relevant.
Structure formation thus gives rise to a selection effect: our mass-biased
view of the Universe differs systematically from the volume-average.

Timescape phenomenology is based theoretically on an extension of the Strong
Equivalence Principle to regional averages (Wiltshire, 2008), via its uniform
quasilocal Hubble expansion condition. This is used to conformally match the
regional geometry (Eq. 2) to the non-FLRW geometry along the past light cone.
One then arrives at effective dressed cosmological parameters, including,
e.g., the dressed matter density with present epoch value

$\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm{\rm
M}0$}}\equiv\left.{{\vphantom{\tilde{\gamma}}\bar{\gamma}}_{\rm
w}}^{3}\bar{\Omega}_{\lower 2.0pt\hbox{$\scriptstyle\rm\rm
M$}}\right|_{0}=\frac{1}{2}(1-f_{\lower 2.0pt\hbox{$\scriptstyle\rm
v0$}})(2+f_{\lower 2.0pt\hbox{$\scriptstyle\rm v0$}}).$ (8)

Dressed parameters do not obey a sum rule akin to Eq. 6 for the bare
parameters, and cannot be interpreted in the conventional FLRW manner.
Nonetheless, since $\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm{\rm M}0$}}$
takes similar numerical values to the FLRW matter density parameter, its use
is convenient for plotting timescape parameter constraints alongside those for
spatially flat $\Lambda$CDM.

Among non-FLRW models with backreaction of inhomogeneities, the timescape
model is observationally the most well-developed. Most significantly, in 16
years since its expansion history was first constrained by SNe Ia data (Leith
et al., 2007), the timescape still stands as a viable alternative to
$\Lambda$CDM, offering a foundational explanation for “dark energy”. In
addition to being well-supported in early studies of gamma-ray-burst distance
indicators (Smale, 2011) and ongoing SNe Ia analyses (Dam et al., 2017),
constraints have also been derived from the CMB (Nazer & Wiltshire, 2015)
using Planck data (Planck Collaboration et al., 2014), and the Baryon Acoustic
Oscillation signal has been extracted directly (Heinesen et al., 2019) from
Baryon Oscillation Spectroscopic Survey (BOSS) data (Alam et al., 2017).

Since spatially flat $\Lambda$CDM has been a good phenomenological fit to
independent datasets for two decades, any viable contender to displace the
standard cosmology must necessarily produce many results which are very close.
An early study (Wiltshire, 2009b) thus concentrated on quantifying the
precision needed to distinguish timescape from $\Lambda$CDM: via the Clarkson-
Bassett-Lu (CBL) test (Clarkson et al., 2008), the $Om(z)$ test (Sahni et al.,
2008; Shafieloo et al., 2009), the Alcock–Paczyński test (Alcock & Paczynski,
1979) and the redshift-time drift test (Sandage, 1962). Among these
observational tests, the CBL test is perhaps the most significant as it yields
a test statistic $\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm k$}}(z)$ of the
Friedmann equation: $\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm k$}}(z)$ is
constant with redshift, and highly constrained by the CMB, for any FLRW
expansion history whether embedded in an exotic matter model or in a modified
gravity model. Any non-FLRW expansion history will differ from this. Timescape
yields differences of $1$–$3$% in the distance–redshift relation at any finite
$z$, leading to large cumulative correlated differences over a long lever arm
in the expansion history. Timescape was used for independent projections for
the Euclid mission in 2014 (Sapone et al., 2014). Consequently, now that
Euclid has finally been launched, in the next six years timescape itself will
be tested as will any FLRW expansion history, including theories routinely
studied by the majority of the community (Di Valentino et al., 2021).

The projections made a decade ago (Sapone et al., 2014) yield a potential
falsification of the Friedmann equation based on robust BAO extraction by
Euclid, and independent distances and redshifts for 1000 SNe Ia. The biggest
systematic uncertainties that the small inhomogeneous cosmology community has
always found the most daunting relate to the use of the FLRW expansion history
in data reduction pipelines777Adapting matter and radiation transfer functions
from the early Universe to the present epoch requires a proxy for the small
but nontrivial violation of the Friedmann equation at last scattering.
Timescape provides a quantitative framework for this, with initial values
$f_{{\rm v}i}\mathop{\sim}\limits 2\times 10^{-5}$,
$\bar{\Omega}_{Ki}\mathop{\sim}\limits 7\times 10^{-5}$ and
$\bar{\Omega}_{{\cal Q}i}\mathop{\sim}\limits-1\times 10^{-5}$ at decoupling
in Eq. 6. Nonetheless, by adapting class (Lesgourgues, 2011; Blas et al.,
2011) to match perturbed FLRW models to timescape Nazer & Wiltshire (2015)
found that the range of physically viable replacements for
$\bar{\Omega}_{{\cal Q}i}$ led to $8$–$13$% systematic uncertainties in
present epoch timescape model parameters (Duley et al., 2013, Table 1).. With
the most recent Pantheon$+$ data release (Scolnic et al., 2022), we are now in
a position to directly test the SNe Ia expansion history projection made a
decade ago by Sapone et al. (2014) by bootstrap resampling. In doing so, we
will also directly perform Bayesian tests that not only compare the spatially
flat $\Lambda$CDM model with the timescape model, but which have a very wide
impact in broadening criteria for model independence in SNe Ia analysis.

###  1.2 Supernova Redshift-Distance Analysis

Dam et al. (2017) investigated the redshift – luminosity distance data of Type
Ia supernovae (SNe Ia) using the Joint Light-curve Analysis (JLA) SNe Ia
catalogue (Betoule et al., 2014) to determine the likelihood of spatially flat
$\Lambda$CDM and timescape. This analysis found no statistical preference
favouring either model. The JLA sample consists of 740 unique supernovae and
was the largest SNe Ia catalogue available at the time of the Dam et al.
(2017) analysis. The JLA catalogue has since been superseded by Pantheon
(Scolnic et al., 2018) and most recently Pantheon$+$ (Scolnic et al., 2022).
The full Pantheon$+$ sample contains a dataset of 1708 supernova from 20
different surveys and 1550 unique supernovae. The Pantheon$+$ and JLA samples
have 584 supernovae in common.

The remaining 176 supernovae of the JLA catalogue have not been included in
Pantheon$+$ as Scolnic et al. (2022) imposed quality cuts on this survey,
partly to satisfy stricter requirements on pre-maximum luminosity data.
Specifically, Brout et al. (2022a) identify significant differences in the
distance modulus for SNLS and DES supernovae with redshifts exceeding 0.8.
These disparities are contingent upon the inclusion of the $U$ band data at
lower redshifts, which is incorporated in their SALT2 training set (Scolnic et
al., 2022). The SNe Ia relying on this data are removed from the dataset due
to the incompatibility for cross-calibration.

Both the JLA and Pantheon$+$ use the SALT2 algorithm (Guy et al., 2007) to
standardise SNe Ia light-curves for cosmological analysis. For each light-
curve SALT2 fits the time stretch parameter ($x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$), colour ($c$), and amplitude $x_{\lower
2.0pt\hbox{$\scriptstyle 0$}}$, which is then used to calculate the rest-frame
apparent magnitude of the $B$ band filter ($m_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}=-2.5\log_{10}(x_{\lower
2.0pt\hbox{$\scriptstyle 0$}})$). These parameters depend on the SALT2
surfaces used, which were updated between JLA and Pantheon$+$. The differences
arising from the update are discussed in Taylor et al. (2021) and result in
differences between the common subsample of SNe Ia.

The SALT2 parameters can be combined into the distance modulus $\mu$ via the
Tripp formula (Tripp, 1998),

$\mu=m-M=m_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}}^{*}-M_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}+\alpha x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}-\beta c,$ (9)

where $m_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}}^{*}$ denotes the apparent
magnitude at maximum in the rest-frame $B$ band and $M_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}$ is the corresponding absolute magnitude of
the source. $\alpha$ and $\beta$ are assumed to be constant, Scolnic et al.
(2022) adopt values of $0.148$ and $3.112$, respectively, for their nominal
fit. Since Scolnic et al. (2022) calculate their nominal fit with $w$CDM and
we wish to conduct a model-independent analysis similar to that of JLA, we do
not use these values.

The observational distance modulus serves as a model for band correction

$\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}}\equiv(m-M)-(m_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}^{*}-M_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}})=\alpha x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}-\beta c,$ (10)

to the theoretical distance modulus, shown by

$\mu\equiv 25+5\log_{10}\left(\frac{d_{\lower 2.0pt\hbox{$\scriptstyle
L$}}}{\rm Mpc}\right),$ (11)

which is determined using the bolometric flux. The luminosity distance $d_{L}$
can be calculated using the redshift and suitable model parameters
$\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm M0$}}$ and $\Omega_{\lower
2.0pt\hbox{$\scriptstyle\rm\Lambda 0$}}=1-\Omega_{\lower
2.0pt\hbox{$\scriptstyle\rm M0$}}$ for the spatially flat $\Lambda$CDM model
and $f_{\lower 2.0pt\hbox{$\scriptstyle\rm v0$}}$ for the timescape cosmology.
For details on this calculation, see Appendix A and Dam et al. (2017, Appendix
A). Evaluating the distance moduli relative to an empty Milne universe (Fig.
1) reveals the key differences between the models tested here. Since the
observed cosmic acceleration is an apparent phenomenon in the timescape model,
the distance modulus values, denoted as $\mu_{\lower
2.0pt\hbox{$\scriptstyle\rm TS$}}$, fall between those of the $\mu_{\lower
2.0pt\hbox{$\scriptstyle\rm\Lambda CDM$}}$ model and the empty case. However,
it is worth noting that they are positioned closer to the $\mu_{\lower
2.0pt\hbox{$\scriptstyle\rm\Lambda CDM$}}$, as illustrated in Fig. 1 and Dam
et al. (2017, Fig. 1).

Figure 1: The residual distance moduli for the P$+1690$ sample following
$\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm\Lambda CDM$}}(z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}})-\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm
Empty$}}(z_{\lower 2.0pt\hbox{$\scriptstyle\rm CMB$}})$ and $\mu_{\lower
2.0pt\hbox{$\scriptstyle\rm TS$}}(z_{\lower 2.0pt\hbox{$\scriptstyle\rm
CMB$}})-\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm Empty$}}(z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}})$, calculated with $H_{\lower
2.0pt\hbox{$\scriptstyle 0$}}=66.7\,\rm kms^{-1}\,Mpc^{-1}$ and assuming the
best-fitting parameters at $z_{\lower 2.0pt\hbox{$\scriptstyle\rm
CMB$}}=0.033$: $\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm M0$}}=0.433$ for
spatially flat $\Lambda$CDM and $f_{\lower 2.0pt\hbox{$\scriptstyle\rm
v0$}}=0.6751$ for timescape (TS). For this dataset, the maximum difference of
the distance moduli in the $\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm
TS$}}(z_{\lower 2.0pt\hbox{$\scriptstyle\rm CMB$}})<\mu_{\lower
2.0pt\hbox{$\scriptstyle\rm\Lambda CDM$}}(z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}})$ region, i.e. $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}<0.8850$, is $\mu_{\lower
2.0pt\hbox{$\scriptstyle\rm\Lambda CDM$}}(z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}})-\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm
TS$}}(z_{\lower 2.0pt\hbox{$\scriptstyle\rm CMB$}})=0.0222$ at $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=0.3375$. The median and the standard
deviation for random subsamples (referred to as RS) of size $N=1000$ (obtained
as described in Section 4.3) are shown for reference.

Another approach used extensively in the 1990s–2000s to ‘standardise’ SNe Ia
light-curves is the Multi-colour Light Curve Shape (MLCS) method, along with
other empirical methods. These methods were pivotal in the breakthrough
discovery by Riess et al. (1998) and Perlmutter et al. (1999). MLCS is an
elaborate empirical model designed to correct the absolute magnitude and
colour of SNe Ia light-curves, rendering them suitable as “standardisable
candles”.

One notable distinction between SALT and MLCS2k2 (an updated version of MLCS)
lies in their implementation of $K$-corrections (Jha et al., 2007). In the
case of SALT, these corrections are integrated into the model, while MLCS2k2
determines the $K$-corrections more empirically through coefficient modelling,
with uncertainties derived from measurement errors (Guy et al., 2005). The
approach in MLCS2k2 involves modelling multi-colour light-curves separately to
calculate luminosity distances, as opposed to estimating model parameters.

Several authors have observed that SALT2 and MLCS2k2 yield different results
in the computed distance modulus. The differing results are associated with
how each method corrects for and handles systematic uncertainties (Kessler et
al., 2009b; Hicken et al., 2009). SALT2 and MLCS2k2 introduce the dependence
on the cosmological model at different stages of the analysis and data
reduction. In MLCS2k2, the cosmological parameters are not entangled with the
distance modulus, whereas SALT2 utilises cosmological parameters to estimate a
distance modulus, thereby maintaining a degree of model dependence (Kessler et
al., 2009b; Hicken et al., 2009; Smale & Wiltshire, 2011).

As demonstrated in the study by Smale & Wiltshire (2011), differences between
various datasets reduced using MLCS2k2 and SALT2 are significant, especially
in the context of the timescape model. Even when SALT2 is applied to the
timescape model using simple_cosfitter, the light-curve parameters provided
have still been fitted with the standard cosmology. However, with the rigorous
approach that Pantheon$+$ takes with SALT2 we are able to extract the relevant
parameters before cosmological biases are introduced, which is crucial for
this analysis.

In this paper, we use the Pantheon$+$ supernova catalogue to conduct a fair
and more detailed comparison of both the spatially flat $\Lambda$CDM and
timescape models. To this end, we construct suitable input and covariance
files as described in Section 2. The methods used for the statistical analysis
are explained in Section 3. In Section 4, we analyse the light-curve
parameters found from the Pantheon$+$ sample by our statistical approach and
compare the different cosmological models. We continue with a discussion of
the most relevant findings in Section 5 and present our conclusions in Section
6.

##  2 Pantheon+ Data

The Pantheon$+$ data release (Scolnic et al., 2022) includes the heliocentric
and CMB corrected redshifts, the light-curve parameters $m_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}$, $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$,
and $c$ alongside the covariances of the light-curve parameters for all
observed supernovae. The CMB corrected redshifts given by Pantheon$+$ include
a correction for the peculiar velocity field, which implicitly assumes the
standard cosmological model based on the FLRW metric and Newtonian $N$-body
simulations. As we test a fundamentally different cosmology, we apply a boost
from the heliocentric rest-frame to calculate the redshift in the CMB rest-
frame without such corrections. In particular, we adopt the Fixsen et al.
(1996) values: $(\ell,b)=(264.14\degree,48.26\degree)$, $v_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=371.0\,$km s-1, $\rm RA_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=168.01\degree$, and $\rm DEC_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=-6.983\degree$. Consequently, we do not
consider any covariance terms associated with the peculiar velocity treatment
by Pantheon$+$ either.

Brout et al. (2022a) provide the results and the covariance matrix in terms of
the distance modulus, which is determined by the SALT2mu algorithm from the
SNANA package (Kessler et al., 2009a). This algorithm requires the input of a
reference and fiducial cosmological model to the light-curve SALT2 fitting
parameters888The results and the covariance are given in terms of $x_{0}$,
$x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ and $c$, and their covariance can be
transformed as ${\rm Cov}(m_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}},x_{\lower
2.0pt\hbox{$\scriptstyle 1$}},c)=\mathbb{J}\>{\rm Cov}(x_{\lower
2.0pt\hbox{$\scriptstyle 0$}},x_{\lower 2.0pt\hbox{$\scriptstyle
1$}},c)\>\mathbb{J}^{\intercal}$. to calculate the distance modulus. As the
reference cosmology used was $w$CDM, any analysis performed with these biases
would inherently favour models based on the FLRW framework. Therefore, we do
not account for several bias-corrections from the Hubble diagram analysis such
as Malmquist bias or peculiar velocity for a fair model comparison. The
reconstruction of the covariance matrix considering these changes with respect
to the one given in Brout et al. (2022a, Section 2.2) is presented here in
Sections 2.1 and 2.2. In doing so, we also consider the three light-curve
parameters taken directly from SALT2 as separate input parameters, instead of
combining them into a single distance modulus parameter, $\mu$. We extract
input and covariance files in terms of $m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}}$, $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ and $c$, to replicate the
methods performed by Dam et al. (2017) and Nielsen et al. (2016) with the JLA
catalogue. In this paper, our reanalysis of JLA differs from the approach
taken in Dam et al. (2017) as we calculate the $\sigma^{2}_{\lower
2.0pt\hbox{$\scriptstyle\rm lens$}}$ and $\sigma^{2}_{z}$ terms using the non-
peculiar velocity corrections.

###  2.1 Statistical Covariance

Considering the input from the fitted light-curve parameters yields a $3\times
3$ covariance block for each supernova. Combining all the $\mathcal{N}=1708$
supernovae gives a $3\mathcal{N}\times 3\mathcal{N}$ block-diagonal matrix,
which accounts for the SALT2 measurement uncertainty term $\sigma_{\lower
2.0pt\hbox{$\scriptstyle\rm meas$}}^{2}$ in the covariance matrix by Brout et
al. (2022a). The redshift covariance $\sigma^{2}_{z}$ and the lensing term
$\sigma^{2}_{\lower 2.0pt\hbox{$\scriptstyle\rm lens$}}=0.055\,z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}$ are calculated from the boosted CMB
redshift and added to the $(m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B,i$}},m_{\lower 2.0pt\hbox{$\scriptstyle\rm B,i$}})$ matrix elements. The
$\sigma^{2}_{\lower 2.0pt\hbox{$\scriptstyle\rm floor$}}$ from the SALT2mu
algorithm and $\sigma^{2}_{\lower 2.0pt\hbox{$\scriptstyle\rm vpec$}}$ terms
are not accounted for, as they introduce a bias towards the standard cosmology
based on the FLRW metric. The full statistical covariance is thus a
$3\mathcal{N}\times 3\mathcal{N}$ block-diagonal matrix build from
$\mathcal{N}$ blocks of the form

$\left(\begin{array}[]{ccc}(\Delta m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}})^{2}&\Sigma_{m_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}},x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}}&\Sigma_{m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}},c}\\\ \Sigma_{m_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}},x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}}&(\Delta x_{\lower 2.0pt\hbox{$\scriptstyle
1$}})^{2}&\Sigma_{x_{\lower 2.0pt\hbox{$\scriptstyle 1$}},c}\\\
\Sigma_{m_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}},c}&\Sigma_{x_{\lower
2.0pt\hbox{$\scriptstyle 1$}},c}&(\Delta
c)^{2}\end{array}\right)+\left(\begin{array}[]{ccc}\sigma^{2}_{z}+\sigma^{2}_{\lower
2.0pt\hbox{$\scriptstyle\rm lens$}}&&\\\ &0&\\\ &&0\end{array}\right),$ (12)

one for each supernova.

However, some of the underlying blocks from the light-curve fitting (left term
of Eq. 12) are not positive semi-definite. For the Pantheon$+$ analysis this
is not a problem, since their full statistical covariance matrix includes more
terms which create a positive semi-definite matrix. As our formulation of the
covariance matrix differs, the complete covariance matrix calculated is not
positive semi-definite. To rectify this, we drop 15 supernovae that have non-
positive semi-definite statistical $3\times 3$ covariance matrices (see Table
2) to get only non-negative real eigenvalues for the full statistical
covariance. The resulting nominal sample consists of $\mathcal{N}=1690$
supernovae, $N=1535$ of which are unique, and will be called P+1690. For use
in future catalogues (Section 4.7 below) the maximum number of unique SNe Ia
events will be denoted $N_{\rm u}$; here $N_{\rm u}=1535$.

###  2.2 Systematic Covariance

The systematic covariance by Brout et al. (2022a),

$C_{\lower 2.0pt\hbox{$\scriptstyle\rm syst$}}^{ij}=\sum_{\psi}\sigma_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}\mu^{i}\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}\mu^{j},$ (13)

accounts for the different calibration uncertainties of the light-curve
parameters and models by iterating over all different systematics ($\psi$)
mentioned in Brout et al. (2022a, Table 4) and assigning respective scale
factors $\sigma_{\lower 2.0pt\hbox{$\scriptstyle\psi$}}$. Here,
$\partial_{\lower 2.0pt\hbox{$\scriptstyle\psi$}}\mu=\mu_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}-\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm
nominal$}}$ is the deviation caused by each calibration. This gives a
$\mathcal{N}\times\mathcal{N}$ matrix for the systematic covariance associated
with the distance modulus $\mu$.

As the same uncertainties apply to our dataset, we calculate our systematic
covariance similarly; but do not account for any calibration associated with
peculiar velocities. Furthermore, we do not account for deviations in $\mu$
only, but consider the three parameters $m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}}$, $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ and $c$ separately999Brout et
al. (2022a) consider $\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}\Delta\mu=\Delta\mu_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}-\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm
nominal$}}$ with $\Delta\mu=\mu-\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm
ref$}}$ in Eq. 13, where $\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm ref$}}$ is
the result of their SALT2mu algorithm and thus assumes a standard FLRW-based
cosmology framework implicitly. However, as $\mu_{\lower
2.0pt\hbox{$\scriptstyle\rm ref$}}$ does not depend on the calibration
systematics $\psi$ except for peculiar velocity calibrations (which we
neglect), it drops out of the equation in our case. We do not consider
reference values for any of the three light-curve parameters for this reason..
Thus, instead of running over the supernova events only, $\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}\mu^{i}\in\\{\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}\mu^{1},...,\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}\mu^{\cal N}\\}$, we consider $\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}\tilde{\mu}^{i}\in\\{\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}}^{1},\partial_{\lower 2.0pt\hbox{$\scriptstyle\psi$}}x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}^{1},\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}c^{1},...,\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}}^{\cal N},\partial_{\lower 2.0pt\hbox{$\scriptstyle\psi$}}x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}^{\cal N},\partial_{\lower
2.0pt\hbox{$\scriptstyle\psi$}}c^{\cal N}\\}$ and obtain a $3\mathcal{N}\times
3\mathcal{N}$ matrix in this way.

###  2.3 Duplicate Supernovae

Observations of the same supernovae event from different surveys were included
in Pantheon$+$ as differing SNe Ia. To account for covariances between surveys
for a given supernova event, we determine a duplication covariance,

$C_{{\rm dupl},i}^{rs}=\delta\tilde{\mu}^{r}\delta\tilde{\mu}^{s},$ (14)

for all duplicated SNe Ia. The indices $r,s$ run over all surveys
$1,...,n_{\lower 2.0pt\hbox{$\scriptstyle i$}}$ that observed this particular
supernova event. This matrix is constructed similar to Eq. 13, but considers
the deviation from the mean over all $n_{\lower 2.0pt\hbox{$\scriptstyle i$}}$
surveys,
$\delta\tilde{\mu}^{r}=\tilde{\mu}^{r}-\langle\tilde{\mu}\rangle_{\lower
2.0pt\hbox{$\scriptstyle\rm{\rm surveys}$}}$. Again, we include
$\delta\tilde{\mu}^{r}\in\\{\delta m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}}^{1},\delta x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}^{1},\delta
c^{1},...,\delta m_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}}^{n_{i}},\delta
x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}^{n_{i}},\delta c^{n_{i}}\\}$, thus
$C_{{\rm dupl},i}^{rs}$ is a $3{n_{\lower 2.0pt\hbox{$\scriptstyle i$}}}\times
3{n_{\lower 2.0pt\hbox{$\scriptstyle i$}}}$ matrix.

We also add the standard deviation of the differences in the brightness of
duplicate supernovae (found by Scolnic et al. (2022) to be $0.102$ mag) to the
off-diagonal $(m_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}}^{r},m_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}^{s})$ matrix elements of these duplication
matrices $C_{{\rm dupl},i}$.

Combining the duplication covariances of all $N$ unique supernovae, we obtain
a block diagonal matrix

$C_{\lower 2.0pt\hbox{$\scriptstyle\rm
dupl$}}=\left(\begin{array}[]{ccc}C_{{\rm dupl},1}&&\\\ &\ddots&\\\ &&C_{{\rm
dupl},N}\end{array}\right)$ (15)

which is of dimension $3\mathcal{N}\times 3\mathcal{N}$ since
$\mathcal{N}=\sum_{i=1}^{N}n_{\lower 2.0pt\hbox{$\scriptstyle i$}}$.

The full covariance matrix is then given by the sum of the statistical,
systematic and duplication covariance. As all three of these matrices are
constructed to be positive semi-definite, so is the full covariance matrix.

##  3 Statistical Methods

Bayesian methods offer a complementary tool to frequentist methods. Given the
extent to which the advantages of each approach have been debated in all
scholarly enquiry, even much further back in academic publishing history than
the much–cited review of Kass & Raftery (1995), it is fair to say both
approaches have their uses, depending on circumstances. Bayesian methods
require prior knowledge that includes models or theories as well as data. The
crux of this paper, and Section 4.7 in particular, is just such a Bayesian
comparison.

Obtaining the Bayesian evidence is computationally intensive. Therefore, a
frequentist statistical approach was used to estimate the parameter values
(having a faster optimisation time), and a Bayesian sampler being used to
numerically compare the two models. Kass & Raftery (1995) provide
classifications for the strength of the evidence; with $|\ln{B}|<1$ having no
statistical preference, $1\leq|\ln{B}|<3$ indicating positive evidence,
$3\leq|\ln{B}|<5$ and $5\leq|\ln{B}|$ representing strong and very strong
evidence respectively.

We follow in the work of Dam et al. (2017), using the framework of Nielsen et
al. (2016) for the purpose of our analysis. There was a debate where Rubin &
Hayden (2016) criticised the methods of Nielsen et al. (2016) for not
considering Malmquist biases. However, Dam et al. (2017) determined that the
additional parameters introduced in Rubin & Hayden (2016) did not
significantly improve Bayesian test statistics. Therefore, for simplicity, we
restrict ourselves to the framework of Nielsen et al. (2016).

The same priors chosen by Dam et al. (2017) were adopted and are summarised in
Table 3. The parameters can be organised in terms of three different groups,
cosmological, nuisance101010Nuisance parameters are those additional to the
model parameters under consideration. Many empirical parameters are in this
class. However, for SNe Ia analysis, $H_{\lower 2.0pt\hbox{$\scriptstyle 0$}}$
is also treated as a nuisance parameter as it is degenerate with $M_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}$ in Eq. 10., and nuisance standard
deviations.

We adopt the Bayesian hierarchical likelihood construction implemented by
March et al. (2011), as used by Nielsen et al. (2016). This can be expressed
in the form

$\displaystyle\mathcal{L}$
$\displaystyle\equiv\mathrm{Pr}\left[\left.(\hat{m}_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}^{*},\hat{x}_{\lower 2.0pt\hbox{$\scriptstyle
1$}},\hat{c})_{\lower 2.0pt\hbox{$\scriptstyle i$}}\right|H\right]$
$\displaystyle=\int\mathrm{Pr}\left[\left.(\hat{m}_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}^{*},\hat{x}_{\lower 2.0pt\hbox{$\scriptstyle
1$}},\hat{c})_{\lower 2.0pt\hbox{$\scriptstyle i$}}\right|(M_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}},x_{\lower 2.0pt\hbox{$\scriptstyle
1$}},c)_{\lower 2.0pt\hbox{$\scriptstyle i$}},H\right]$ $\displaystyle\hbox
to22.76219pt{\hfil}\times\mathrm{Pr}\left[\left.(M_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}},x_{\lower 2.0pt\hbox{$\scriptstyle
1$}},c)_{\lower 2.0pt\hbox{$\scriptstyle
i$}}\right|H\right]\mathrm{d}M_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}}^{\mathcal{N}}\mathrm{d}x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}^{\mathcal{N}}\mathrm{d}c^{\mathcal{N}},$ (16)

where the quantities which are denoted with a hat are considered to be
observed values, and the true values (the quantities which are not denoted by
a hat) are drawn from Gaussian distributions and account for fluctuations from
various sources of noise (Nielsen et al., 2016; March et al., 2011; Dam et
al., 2017). The true data represents the intrinsic parameters utilised
explicitly in the SALT2 relation.

Nielsen et al. (2016) follow the analysis of March et al. (2011) and adopt
global, independent Gaussian distributions for $M$, $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ and $c$ to determine the probability density of
the true parameters (Nielsen et al., 2016; March et al., 2011). This is
determined by $\mathrm{Pr}[(M,x_{\lower 2.0pt\hbox{$\scriptstyle
1$}},c)|\Theta]=\mathrm{Pr}(M|\Theta)\times\mathrm{Pr}(x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}|\Theta)\times\mathrm{Pr}(c|\Theta)$; where for
$\tilde{\mu}\in\\{M,x_{\lower 2.0pt\hbox{$\scriptstyle 1$}},c\\}$, the
probability can be expressed as
$\mathrm{Pr}(\tilde{\mu}|\Theta)\equiv(2\pi\sigma_{\tilde{\mu}}^{2})^{-1/2}\times\exp\left(-\frac{1}{2}\left[(\tilde{\mu}-\tilde{\mu}_{0})/\sigma_{\tilde{\mu}_{0}}\right]^{2}\right)$.
By creating a vector of the true parameters, $Y=[M_{\lower
2.0pt\hbox{$\scriptstyle\rm B,1$}},x_{\lower 2.0pt\hbox{$\scriptstyle
1,1$}},c_{\lower 2.0pt\hbox{$\scriptstyle 1$}},...,M_{\lower
2.0pt\hbox{$\scriptstyle\rm B,\mathcal{N}$}},x_{\lower
2.0pt\hbox{$\scriptstyle 1,\mathcal{N}$}},c_{\lower
2.0pt\hbox{$\scriptstyle\cal N$}}]$ and a diagonal matrix of dispersion
uncertainties $\Sigma_{\lower 2.0pt\hbox{$\scriptstyle
l$}}=\mathrm{diag}(\sigma_{M_{\lower 2.0pt\hbox{$\scriptstyle
0$}}}^{2},\sigma_{x_{\lower 2.0pt\hbox{$\scriptstyle
1,0$}}}^{2},\sigma_{c_{\lower 2.0pt\hbox{$\scriptstyle 0$}}}^{2},...)$, this
can be constructed by

$\mathrm{Pr}(Y|\Theta)\equiv\det\left(2\pi\Sigma_{\lower
2.0pt\hbox{$\scriptstyle
l$}}\right)^{-1/2}\exp\left(-\frac{1}{2}\left(Y-Y_{\lower
2.0pt\hbox{$\scriptstyle 0$}}\right)\Sigma_{\lower 2.0pt\hbox{$\scriptstyle
l$}}^{-1}\left(Y-Y_{\lower 2.0pt\hbox{$\scriptstyle
0$}}\right)^{\intercal}\right),$ (17)

where $Y_{\lower 2.0pt\hbox{$\scriptstyle 0$}}=[M_{\lower
2.0pt\hbox{$\scriptstyle\rm B,0$}},x_{\lower 2.0pt\hbox{$\scriptstyle
1,0$}},c_{\lower 2.0pt\hbox{$\scriptstyle 0$}},M_{\lower
2.0pt\hbox{$\scriptstyle\rm B,0$}},...]$ is the zero-point vector of the
parameters, drawn from the intrinsic Gaussian averages.

To transform the absolute magnitude to an apparent magnitude minus the
distance moduli, the block diagonal matrix $A$ given by

$A=\left(\begin{array}[]{cccc}1&0&0&\\\ -\alpha&1&0&0\\\ \beta&0&1&\\\
&0&&\ddots\\\ \end{array}\right),$ (18)

is introduced, which transforms as $Y_{\lower 2.0pt\hbox{$\scriptstyle
0$}}A=[M_{\lower 2.0pt\hbox{$\scriptstyle\rm B,0$}}-\alpha x_{\lower
2.0pt\hbox{$\scriptstyle 1,0$}}+\beta c_{\lower 2.0pt\hbox{$\scriptstyle
0$}},x_{\lower 2.0pt\hbox{$\scriptstyle 1,0$}},c_{\lower
2.0pt\hbox{$\scriptstyle 0$}},M_{\lower 2.0pt\hbox{$\scriptstyle\rm
B,0$}}-\alpha x_{\lower 2.0pt\hbox{$\scriptstyle 1,0$}}+\beta c_{\lower
2.0pt\hbox{$\scriptstyle 0$}},...]$. Furthermore, the estimated experimental
covariance matrix of statistical and systemic errors, $\Sigma_{d}$, and the
new vector $\hat{Z}_{\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle 0$}}}$ are
considered. $\hat{Z}_{\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle 0$}}}$ is
constructed from the observational values for $m_{\lower
2.0pt\hbox{$\scriptstyle\rm B$}}^{*}$, $x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}$ and $c$ by $\hat{Z}_{\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle
0$}}}=[\hat{m}_{\lower 2.0pt\hbox{$\scriptstyle\rm B,1$}}^{*}-\mu_{\lower
2.0pt\hbox{$\scriptstyle 1$}}-\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle
0$}},\hat{x}_{\lower 2.0pt\hbox{$\scriptstyle 1,1$}},\hat{c}_{\lower
2.0pt\hbox{$\scriptstyle 1$}},...,\hat{m}_{\lower 2.0pt\hbox{$\scriptstyle\rm
B,\mathcal{N}$}}^{*}-\mu_{\lower 2.0pt\hbox{$\scriptstyle\rm\cal
N$}}-\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle 0$}},\hat{x}_{\lower
2.0pt\hbox{$\scriptstyle\rm 1,\mathcal{N}$}},\hat{c}_{\lower
2.0pt\hbox{$\scriptstyle\rm\cal N$}}]$ and transforms as $YA$. Together, this
yields the likelihood

$\displaystyle\mathcal{L}_{\Delta\mu_{0}}$
$\displaystyle=\int\mathrm{Pr}(\hat{Z}_{\Delta\mu_{0}}|Y,\Theta)\
\mathrm{Pr}(Y|\Theta)\ \mathrm{d}Y$
$\displaystyle=\det\left[2\pi\left(\Sigma_{\lower 2.0pt\hbox{$\scriptstyle
d$}}+A^{\intercal}\Sigma_{\lower 2.0pt\hbox{$\scriptstyle
l$}}A\right)\right]^{-1/2}$ $\displaystyle\hbox
to8.53581pt{\hfil}\times\exp\left[-\frac{1}{2}\left(\hat{Z}_{\Delta\mu_{\lower
2.0pt\hbox{$\scriptstyle 0$}}}-Y_{0}A\right)\left(\Sigma_{\lower
2.0pt\hbox{$\scriptstyle d$}}+A^{\intercal}\Sigma_{\lower
2.0pt\hbox{$\scriptstyle l$}}A\right)^{-1}\left(\hat{Z}_{\Delta\mu_{\lower
2.0pt\hbox{$\scriptstyle 0$}}}-Y_{0}A\right)^{\intercal}\right],$ (19)

for given zero-point offset $\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle 0$}}$.

This parameter must be included to account for the possibility that the Hubble
constant varies between subsamples as different redshift cuts are made. In
standard analyses where the whole sample is considered without implementing a
sequence of redshift cuts, this question is only considered once in a single
marginalisation step. The offset indicates a universal degeneracy between the
choice of the Hubble constant and the absolute magnitude of the supernovae.
Naturally, this “zero-point calibration” issue is intimately connected to the
Hubble tension. Our aim here is to compare two cosmological models
independently of the zero-point issue, using Bayesian methods. Consequently,
while the zero-point issue is of considerable interest, it is temporarily set
aside in making redshift cuts by marginalising over $\Delta\mu_{\lower
2.0pt\hbox{$\scriptstyle 0$}}$. A similar approach was adopted in past
analyses of the subsamples Gaur (2019); Lane (2022) of the Pantheon and
Pantheon+ surveys common to JLA, in which only distance moduli fit with FLRW
geometry assumptions were available. Considering $\hat{Z}_{\Delta\mu_{\lower
2.0pt\hbox{$\scriptstyle 0$}}}-Y_{0}A$ and thus
$\mathcal{L}_{\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle 0$}}}$ as a function
of $y+\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle 0$}}$, the reduced likelihood
obtained by the marginalisation is given by

$\displaystyle\mathcal{L}$ $\displaystyle=\frac{1}{6\sigma_{\lower
2.0pt\hbox{$\scriptstyle y$}}}\int_{-3\sigma_{\lower 2.0pt\hbox{$\scriptstyle
y$}}}^{3\sigma_{\lower 2.0pt\hbox{$\scriptstyle
y$}}}\mathcal{L}_{\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle
0$}}}(y+\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle
0$}})\mathrm{d}\Delta\mu_{\lower 2.0pt\hbox{$\scriptstyle 0$}},$ (20)

where $y=M_{\lower 2.0pt\hbox{$\scriptstyle\rm B,0$}}-\alpha x_{\lower
2.0pt\hbox{$\scriptstyle 1,0$}}+\beta c_{\lower 2.0pt\hbox{$\scriptstyle 0$}}$
is the zero-point value of $m_{\lower 2.0pt\hbox{$\scriptstyle\rm
B$}}^{*}-\mu$ and
$\sigma_{y}=\sqrt{\sigma_{M_{0}}^{2}+\alpha^{2}\sigma_{x_{\lower
2.0pt\hbox{$\scriptstyle 1,0$}}}^{2}+\beta^{2}\sigma_{c_{0}}^{2}}$ is its
standard deviation. As the marginalisation slows down the convergence
significantly, it is performed for the frequentist analysis only. We analyse
its effects for this approach but calculate the Bayesian results based on the
original likelihood from Eq. 19.

We utilised a nested Bayesian sampler PyMultiNest (Buchner et al., 2014),
which interacts with the MultiNest (Feroz & Hobson, 2008; Feroz et al., 2009,
2019) code. This was used to compute the evidence and determine the Bayes
factor in comparing the spatially flat $\Lambda$CDM and timescape models. The
accuracy of PyMultiNest is determined by the number of ‘live’ points used in
the sampling, with the error of order $\sigma\sim\mathcal{O}(n_{\lower
2.0pt\hbox{$\scriptstyle\rm live$}}^{-1/2})$; and the tolerance for the nested
sampling. We chose a tolerance of $10^{-3}$ and $n_{\lower
2.0pt\hbox{$\scriptstyle\rm live$}}=1000$ for nine parameters.

##  4 Results

Following Dam et al. (2017), we perform a likelihood analysis for subsets of
the P+1690 sample determined by redshift cuts $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}\geq z_{\lower 2.0pt\hbox{$\scriptstyle\rm
min$}}$ for varying $z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}$ (Fig. 4).
We investigate the light-curve parameters obtained in this way and compare the
results for analyses started from different subsamples instead of the full
P+1690 sample.

Since the light-curve parameters exhibit variations across different redshift
cuts, we will refine the concept of the Scale of Statistical Homogeneity
(SHS). We first denote SHS measures derived empirically from $n$–point galaxy
correlation functions as the SHSgn statistics. Thus conventional analyses
(Hogg et al., 2005; Sylos Labini et al., 2009; Scrimgeour et al., 2012) based
on 2-point galaxy correlation functions investigate the SHSg2 statistic. Our
analysis here will not be based on higher order galaxy correlation functions,
but on the empirical quantities $\beta c$ (colour) and $\alpha x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ (stretch) that appear in the Tripp relation Eq.
9. We therefore introduce the statistics SHSα (stretch) and SHSβ (colour) to
to explore this phenomenon in Section 4.4. The SHSα statistic, derived from
the light-curve parameters, serves as an indicator of the point at which the
results of SNe Ia free from biases and subject to numerical uncertainty,
approach a statistically homogeneous Universe.

The analysis of the light-curve parameters, the random subsamples and the SHSα
is based on frequentist methods, involving marginalised likelihoods if needed.
However, for the evaluation and comparison of the two models, we consider both
the frequentist and Bayesian approach as described in Section 3, and presented
in Section 4.7.

###  4.1 Pantheon+/JLA Common Subsample

Figure 2: The differences ($\Delta$) between the $\alpha x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ and $-\beta c$ Tripp equation terms calculated
for the Pantheon$+$ and JLA common subsample (Eq. 9). The differences between
the two datasets asymptotically approach $0.00360\pm 0.00013$ ($0.00378\pm
0.00012$) for $\alpha x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ and $0.0405\pm
0.0022$ ($0.0381\pm 0.0021$) for $-\beta c$ based on the timescape (spatially
flat $\Lambda$CDM) model. These differences arise from the change from SALT2
for the JLA analysis to SALT2-2021 for Pantheon$+$ (Taylor et al., 2021). The
vertical dotted line represents the statistical homogeneity scale (SHSα) we
define to be $z_{\lower 2.0pt\hbox{$\scriptstyle\rm CMB$}}=0.075$.

To compare the Pantheon$+$ results to the ones by Dam et al. (2017) based on
the JLA catalogue, we compute the common subsample of the two supernovae
catalogues. This common subset consists of 584 SNe Ia, for which both the JLA
and Pantheon$+$ data are accessible. Within this set of 584 supernovae, it was
necessary to exclude four of them in order to ensure that the covariance
matrix remains positive semi-definite (as detailed in Section 2.1 and Appendix
B). As a result, we will refer to this sample as P+580.

Due to the presence of duplicate supernova events, P+580 comprises 635
supernova observations. For the sake of direct comparison, we focus on the JLA
data encompassing the same 580 unique supernovae as those found in P+580,
denoting this dataset as JLA580. Our analysis encompasses both datasets,
revealing significant disparities in the light-curve parameters. Notably, the
asymptotic difference in the $-\beta c$ term, as described by Eq. 9, amounts
to $0.0381\pm 0.0021$ when assuming a spatially flat $\Lambda$CDM cosmology
(as shown in Fig. 2). This difference can be attributed to the alterations in
the SALT2 surface, as documented by Taylor et al. (2021).

###  4.2 Low-Redshift Supernovae

Figure 3: The distribution of the SNe Ia redshift in the CMB frame (upper) and
its influence on the $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ parameter for
the subsamples (lower, for spatially flat $\Lambda$CDM model). The spikes in
the low redshift regime for the P+1690 sample can be avoided by choosing a
reduced low-redshift sample or using the marginalising approach. The vertical
dotted line represents the SHSα found at $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=0.075$.

The light-curve parameters (e.g., the stretch $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$, shown in Fig. 3) have a large scatter for low
redshift cuts when fitted to the full Pantheon$+$ dataset. This behaviour is
not observed when analysing the full JLA sample (Dam et al., 2017), nor does
it appear when considering the P+580 subsample. To determine if this feature
is redshift dependent, we construct a new subsample P+1690, reduced z, which
mimics the P+580 distribution of supernovae for $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}<0.05$. As this dataset has minimal scatter
(Fig. 3), we confirm that this artefact is due to the high number of low
redshift supernovae with a wide range of $x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}$ and $c$ parameters reported by Pantheon$+$. As discussed in Section 3,
different subsamples of supernova observations within the full dataset might
rely on different values for the Hubble constant due to varying redshift cuts.
This is a known source for variations throughout the full sample and is
accounted for by marginalisation over the zero-point offset by choosing a
different likelihood (Eq. 20). As shown in Fig. 3, this approach also resolves
the artefacts in the low redshift regime.

###  4.3 Random Subsamples

Figure 4: The SALT2 parameters for JLA (left), Pantheon$+$ (middle) and the
marginalised results (right). The shaded regions give one standard deviation
around the median of the parameters at a given $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ calculated from our set of random
subsamples (referred to as RS). The right column shows results which have been
integrated with Eq. 20 to account for the choice of Hubble constant. For the
offset between Pantheon$+$ data (middle column) and JLA data (left column) see
Fig. 2. For plots of the products $\alpha x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}$ and $-\beta c$ that contribute to the distance modulus $\mu$, see Fig.
12.

The asymptotic results for the light-curve parameters are not numerically
stable, but instead, depend on the selected subsample of the Pantheon$+$
dataset. $50$ random subsamples were chosen of size $N=1000$, where $N$ refers
to the number of unique supernovae chosen (the actual random subsamples may
include more individual supernova observations). The median and standard
deviation of the light-curve parameters for varying $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ computed from these samples are shown in
Fig. 4 (middle column) alongside the full sample and the Pantheon$+$/JLA
common subsample. Wagner & Meyer (2019) use a similar method of taking random
subsamples to systemically investigate the sampling issues in the Pantheon
catalogue.

We conduct a similar analysis with the JLA data, choosing the size $N$ of the
random samples to match that of the common subsample ($N=580$). The results
(left column of Fig. 4) show that the full sample is within $2\sigma$ of the
randomised results.

For the Pantheon$+$ common dataset (P+580), we observe important differences
when comparing it to the full sample and the random subsamples of $N=1000$
(middle column of Fig. 4). The values of $\alpha$, $\beta$, $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ and $c$ are within 2–3$\sigma$ bounds of the
full sample, and within 1–2$\sigma$ of P+580. These variations required
further analysis to identify if the value changes are a product of taking
random subsamples of differing sizes.

Figure 5: The change in the light-curve parameters for increasing number of
unique SNe Ia ($N$) for the random samples for the timescape model (left) and
tests of the effects of distribution on the sample (right). The medians of the
samples were taken of all values after the SHSα ($z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=0.075$ up to $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}=0.1$) for both the left (50 random
subsamples) and right columns (containing 10 different random subsamples). The
star and cross markers indicate the values of the full (non-marginalised)
P+1690 and marginalised P+1690 evaluations respectively. Note that the P+1690
sample consists of $N=1535$ unique SNe Ia. The dashed lines indicate the
prospective advancement with sample size.

As the common subsample is much smaller than $N=1000$, we investigate the
possible effects of the sample size on the scatter of random subsamples. In
particular, we calculated the light-curve parameters of subsamples of P+1690
for111111The lowest value for $N$ was chosen for comparison with P+580.
$N=580,\,660,\,750,\,875,\,1000,\,1125,\,1250,\,1375$ at121212We choose
$z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}=0.033$ as an estimate for the
SHSα that Dam et al. (2017) used for their results. $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}=0.01,\,0.033,\,0.05,\,0.075,\,0.1$. Up to
statistical fluctuations, all of these random subsamples have the same
redshift distribution as the full P+1690 sample.

The model parameters, which vary with sample size ($N$), are depicted in Fig.
5. Our analysis reveals that for the non-marginalised calculations using the
likelihood from Eq. 19, the parameter $\beta$ demonstrates a clear dependence
on the sample size, whereas $\alpha$, $x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}$, $c$, and $\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm M0$}}$ fall
comfortably within the uncertainty bounds of the entire sample.

$\beta$ converges in parameter-space for each of the various samples across
differing redshift cuts; however, these parameters do not converge to a
uniform value for all sample sizes as expected. Interestingly, increasing the
size of each sample appears to increase the $\beta$ value for both the
marginalised and non-marginalised samples, therefore this phenomenon is not
attributable to a specific choice of the Hubble parameter.

To investigate what influence the redshift distribution has on the $\beta$
convergence, we compared various samples of size $N=580$ with differing
redshift distributions. For the samples highlighted in green in Fig. 5, we
selected all SNe Ia with an equal probability, thereby resulting in a redshift
distribution nearly identical to that of the full sample. For the ‘high’
sample, we considered the 290 SNe Ia from the P+1690 dataset with the highest
redshift (or with the lowest redshift for the ‘low’ sample) and selected the
remaining 290 SNe Ia from the remainder of the P+1690 dataset using the same
methodology as the uniform sampling.

As displayed in Fig. 5, the ‘low’ sample is consistent with the full sample
for all parameters, including $\beta$. This confirms that the behaviour
observed for the uniform subsamples is caused by the redshift distribution.
For further analysis of the parameter evolution, a full sample with a uniform
distribution would be needed.

Figure 6: The convergence of the $\alpha x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}$ light-curve parameter for the spatially flat $\Lambda$CDM model across
various redshift cuts. In the upper plot, a power–law model has been fit to
the data, and the green shaded region represents within 5% of the median value
within the range $0.1\leq z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}<0.14$
indicating when the model converges. The vertical dotted line represents the
SHSα found at $z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}=0.075$. The lower
plot provides an enlarged view of the data to illustrate finer details. Figure
7: Light-curve parameters $\beta$ and $c$ and their degenerate product for the
spatially flat $\Lambda$CDM model at higher redshifts to analyse how they
evolve with redshift cuts. The low-redshift scatter shows the power–law model
lower bound for our SHSα analysis. Here we see the effects of Malmquist bias
on our parameters, making us unable to use them for determining an SHSβ as an
alternative to the α. For the analysis of $\alpha$, $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ and $\alpha x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}$, see Fig. 13.

###  4.4 Scale of Statistical Homogeneity

Dam et al. (2017) report significant changes in the light-curve parameters in
the presence of an emerging SHS. For the JLA sample, a change in the
parameters $\beta$, $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ and $c$ across
the differing redshift cuts is clearly visible in Dam et al. (2017, Fig. 2)
and in the corresponding panels of Fig. 4. Specifically, the analysis of
randomly subsampled data validates the conclusions drawn from the entire JLA
dataset.

Our definition of the SHSα is exclusively rooted in empirical light-curve
fitting parameters. This approach is adopted to ensure the independence of
cosmological models and eliminate any assumptions associated with the
perturbed FLRW geometry. We achieved this by fitting a power–law model to the
$\alpha x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ parameter, represented by
$-az^{-n}+C$, where $a$, $n$, $C$ are parameters subject to fitting131313These
parameters were fit with the Markov Chain Monte-Carlo sampling tool emcee
(Foreman-Mackey et al., 2013) with values of $a=(23.2\pm 3.5)\times 10^{-6}$,
$n=1.76\pm 0.04$, and $C=(21.61\pm.10)\times 10^{-3}$., and $z$ represents the
redshift cut.

Figure 8: Normalised likelihood at $z_{\lower 2.0pt\hbox{$\scriptstyle\rm
CMB$}}=0.033$ of the model parameters, $\Omega_{\lower
2.0pt\hbox{$\scriptstyle\rm M0$}}$ for the spatially flat $\Lambda$CDM and
$f_{\lower 2.0pt\hbox{$\scriptstyle\rm v0$}}$ for the timescape model for
Pantheon$+$ and for the 50 $N=1000$ random subsamples (RS). The maximum
likelihood values are $f_{\lower 2.0pt\hbox{$\scriptstyle\rm v0,MLE$}}=0.675$,
and $\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm M0,MLE$}}=0.433$,
respectively, for the P+1690 sample. The integrated non-normalised likelihood
of the P+1690 sample is $2.73\times 10^{167}$ for spatially flat $\Lambda$CDM
and $2.06\times 10^{167}$ for timescape. Figure 9: Model parameter
$\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm M0$}}$ for spatially flat
$\Lambda$CDM and $f_{\lower 2.0pt\hbox{$\scriptstyle\rm v0$}}$ for the
timescape model, evaluated with Bayesian and frequentist statistical methods.
The dotted line represents the SHSα.

We consider the redshift cut to be within the SHSα when the power–law falls
within $5\%$ of the limit of $\alpha x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$
(with the uncertainties coming from 3 and 7% intervals). This limit is
estimated by the median value of $\alpha x_{\lower 2.0pt\hbox{$\scriptstyle
1$}}$ in the range $0.10\leq z_{\lower 2.0pt\hbox{$\scriptstyle\rm
min$}}<0.14$, which we chose as to when the parameter has obviously converged,
as seen in Fig. 6.

We employed the power–law model within the range of $0.0225\leq z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}<0.14$, where the bounds were chosen such
that they avoid low-redshift noise and reduced SNe sampling at higher redshift
(for additional details, refer to Section 4.2). From this fit we find a
parameter convergent SHSα of $z_{\lower 2.0pt\hbox{$\scriptstyle\rm
CMB$}}=0.075^{+0.007}_{-0.009}$, which is significantly greater than values
estimated in earlier studies (Dam et al., 2017; Hogg et al., 2005; Scrimgeour
et al., 2012). However, as opposed to the statistical homogeneity scale SHSg2
considered by (Hogg et al., 2005; Scrimgeour et al., 2012), our definition of
the SHSα does not use the two-point galaxy correlation function but relies on
the convergence of the light-curve parameters and the distance modulus. This
approach respects general features in the light-curve parameters arising from
unaccounted large scale structures that are only now apparent using the
Pantheon$+$ dataset (further discussed in Section 4.7). It thus leads to a
higher value of the SHSα.

To check the impact of the high redshift cut on the fit, we extend the upper
bound of $z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}$ to be 0.25. Parameter
values at $z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}=0.25$ have nearly 20%
less unique SNe Ia than $z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}=0.14$
and therefore are more prone to numerical fluctuations with increasing
redshift cut. For the extended range we find a SHSα of $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=0.086^{+0.009}_{-0.013}$, which is within
1$\sigma$ of our stated value. We therefore use our initial fit over the range
$0.0225\leq z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}<0.14$ for the SHSα.

To conduct our convergence analysis, we focus on the $\alpha x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ parameter to mitigate numerous potential issues
stemming from Malmquist bias, which falls beyond the scope of this current
analysis. When examining the $\beta$ or $c$ parameters there is no convergence
within the tested range, thus a SHSβ can not be defined based on the P+1690
sample.

The $\beta$ ‘constant’ parameter exhibits a pattern of fluctuation, peaking at
a value exceeding 4.30 just after $z_{\lower 2.0pt\hbox{$\scriptstyle\rm\rm
min$}}=0.07$ and subsequently decreasing to 3.77 at $z_{\lower
2.0pt\hbox{$\scriptstyle\rm\rm min$}}=0.25$ (as shown in Fig. 7). However,
this trend is counterbalanced by the supernovae appearing ‘bluer’ as we extend
our observations to higher redshifts. When we investigate the degeneracy
between $\beta$ and the colour parameter $c$, we are able to discern and
elucidate the influences of Malmquist bias on our data as we progress through
different redshift cuts (Fig. 7).

We also did not solely utilise the $\alpha$ parameter for our analysis. This
was to be in alignment with previous studies (e.g., Dam et al. (2017)), which
reported minimal variations in its value within the region of an emerging
SHSα. Furthermore, it exhibits a significant decline at $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}\sim 0.17$, rendering it an unreliable
indicator for recovering the SHSα (see Fig. 13). Nevertheless, it does play a
valuable role in constraining the scatter associated with the lower bound.

The decrease in the $\alpha$ parameter, however, is effectively balanced by an
increase in $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$, resulting in the
$\alpha x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ degenerate parameter
remaining nearly constant within the tested high-redshift range (Fig. 6 and
Fig. 13).

We explored various convergence methods, including statistical analyses of the
parameter spaces and gradients based on high-redshift cuts. However, the
‘power–law’ methodology demonstrated the advantage of requiring minimal human
intervention, thereby allowing for a more standardised approach to the
analysis.

Table 1: Median and standard deviation of the model and light-curve parameters beyond the scale of statistical homogeneity, $z=0.075$, for the different samples, models and statistical methods. |  |  | $\Omega_{\rm M0}$ | $\alpha$ | $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ | $\beta$ | $c$
---|---|---|---|---|---|---|---
JLA580 | TS | Freq. | 0.293 $\pm$ 0.011 | 0.1355 $\pm$ 0.0008 | 0.146 $\pm$ 0.004 | 3.192 $\pm$ 0.015 | -0.01849 $\pm$ 0.00029
|  | Bay. | 0.305 $\pm$ 0.018 | 0.1358 $\pm$ 0.0010 | 0.148 $\pm$ 0.009 | 3.188 $\pm$ 0.023 | -0.0176 $\pm$ 0.0008
|  $\Lambda$CDM | Freq. | 0.361 $\pm$ 0.006 | 0.1354 $\pm$ 0.0009 | 0.146 $\pm$ 0.004 | 3.197 $\pm$ 0.016 | -0.0184 $\pm$ 0.0003
|  | Bay. | 0.362 $\pm$ 0.007 | 0.1365 $\pm$ 0.0009 | 0.153 $\pm$ 0.009 | 3.208 $\pm$ 0.012 | -0.0182 $\pm$ 0.0008
JLA740 | TS | Freq. | 0.352 $\pm$ 0.004 | 0.1314 $\pm$ 0.0010 | 0.137 $\pm$ 0.005 | 3.151 $\pm$ 0.015 | -0.02198 $\pm$ 0.00020
|  | Bay. | 0.338 $\pm$ 0.009 | 0.1315 $\pm$ 0.0013 | 0.135 $\pm$ 0.004 | 3.147 $\pm$ 0.028 | -0.0219 $\pm$ 0.0007
|  $\Lambda$CDM | Freq. | 0.3769 $\pm$ 0.0027 | 0.1314 $\pm$ 0.0011 | 0.138 $\pm$ 0.005 | 3.154 $\pm$ 0.015 | -0.02191 $\pm$ 0.00023
|  | Bay. | 0.373 $\pm$ 0.007 | 0.1319 $\pm$ 0.0015 | 0.137 $\pm$ 0.003 | 3.162 $\pm$ 0.022 | -0.0219 $\pm$ 0.0008
P+580 | TS | Freq. | 0.289 $\pm$ 0.009 | 0.1490 $\pm$ 0.0009 | 0.148 $\pm$ 0.004 | 3.765 $\pm$ 0.028 | -0.02516 $\pm$ 0.00016
|  | Bay. | 0.291 $\pm$ 0.008 | 0.150 $\pm$ 0.0009 | 0.147 $\pm$ 0.007 | 3.76 $\pm$ 0.03 | -0.0254 $\pm$ 0.0009
|  $\Lambda$CDM | Freq. | 0.355 $\pm$ 0.005 | 0.1491 $\pm$ 0.0009 | 0.150 $\pm$ 0.004 | 3.778 $\pm$ 0.029 | -0.02439 $\pm$ 0.00018
|  | Bay. | 0.359 $\pm$ 0.007 | 0.1489 $\pm$ 0.0009 | 0.153 $\pm$ 0.005 | 3.79 $\pm$ 0.03 | -0.0244 $\pm$ 0.0014
P+1690 | TS | Freq. | 0.467 $\pm$ 0.003 | 0.1710 $\pm$ 0.0004 | 0.1172 $\pm$ 0.0011 | 4.185 $\pm$ 0.025 | -0.03304 $\pm$ 0.00016
|  | Bay. | 0.462 $\pm$ 0.009 | 0.1688 $\pm$ 0.0013 | 0.119 $\pm$ 0.005 | 3.992 $\pm$ 0.005 | -0.0334 $\pm$ 0.0004
|  $\Lambda$CDM | Freq. | 0.438 $\pm$ 0.003 | 0.1709 $\pm$ 0.0004 | 0.1182 $\pm$ 0.0011 | 4.217 $\pm$ 0.026 | -0.03276 $\pm$ 0.00014
|  | Bay. | 0.440 $\pm$ 0.007 | 0.1706 $\pm$ 0.0015 | 0.1168 $\pm$ 0.0026 | 3.9956 $\pm$ 0.0024 | -0.0323 $\pm$ 0.0008

###  4.5 Profile Likelihood

The (non-marginalised) profile likelihood at141414Again, this value was chosen
for comparison with Dam et al. (2017). $z_{\lower 2.0pt\hbox{$\scriptstyle\rm
CMB$}}=0.033$ gives the maximum likelihood value for the cosmological
parameter as $f_{\lower 2.0pt\hbox{$\scriptstyle\rm v0,MLE$}}=0.675$, and
$\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm M0,MLE$}}=0.433$ (Fig. 8) when
maximising the nuisance parameters, as seen by

$\mathcal{L}_{\lower 2.0pt\hbox{$\scriptstyle\rm
p$}}(\Theta)=\max_{\phi}\mathcal{L}(\Theta,\phi),$ (21)

where $\mathcal{L}$ is the likelihood as defined in Eq. 19 in terms of the
nuisance parameters $\phi$, and $\Theta$ is constrained to take a range of
values. The profile likelihoods have similarities to the JLA likelihoods from
Dam et al. (2017, Fig. 3b) in shape, however, the relative sizes of the
likelihoods differ between the JLA and Pantheon$+$ surveys. The JLA likelihood
seems to favour the timescape cosmology. while the Pantheon$+$ likelihood, by
contrast, favours the spatially flat $\Lambda$CDM cosmology. The integrated
likelihood of the spatially flat $\Lambda$CDM is 1.32 times higher than
timescape for the P+1690 sample with $z_{\lower 2.0pt\hbox{$\scriptstyle\rm
min$}}=0.033$ and 1.18 times as high for $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}=0.075$. In particular, the likelihood ratio
decreases when restricting to the SNe Ia beyond the SHSα.

The MLE value for the cosmological parameters is in the broad 1–2 $\sigma$
bounds chosen for the Bayesian analysis. Fig. 8 shows the likelihood for the
P+1690 sample and the random subsamples which are normalised for comparison.
When we consider the entire sample, the likelihood exhibits a narrower
Gaussian profile when contrasted with the subsamples of size $N=1000$.

###  4.6 Statistical Approaches

In our analysis, we considered the frequentist approach using the marginalised
(Eq. 20) and non-marginalised (Eq. 19) likelihoods as well as the Bayesian
evidence. Because of the computational resource requirements associated with
marginalising the likelihood, we exclusively performed the latter based on Eq.
19. A comparison of the Bayesian and frequentist approach based on the
cosmological model parameters is shown in Fig. 9, while the differences
between the marginalised and non-marginalised frequentist results are depicted
in Fig. 10. The values beyond the SHSα for all parameters, models and
statistical methods are given in Table 1.

In the frequentist approach, the marginalised likelihood Eq. 20 maintains the
same shape as the non-marginalised treatment, with different numerical values
(as shown in the right columns of Figs. 4, 10 and 5). In particular, the
results found for $\alpha$ and $\beta$ are significantly lower in the
marginalised case. Furthermore, the marginalised likelihood effectively
manages the dispersion within the low-redshift range (see Section 4.2), except
with low values of sample size ($N$).

For a quantitative analysis of the statistical approaches, we consider the
median and standard deviation of the model and light-curve parameters beyond
the SHSα (Table 1). For the light-curve parameter values $\alpha,x_{\lower
2.0pt\hbox{$\scriptstyle 1$}},\beta$ and $c$, there are no discernible
differences between the two statistical approaches nor for the two models.

Beyond the light-curve parameters, however, it is essential to note that the
timescape model depends on the parameter $f_{\lower
2.0pt\hbox{$\scriptstyle\rm v0$}}$ according to Eq. 8, while spatially flat
$\Lambda$CDM considers $\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm M0$}}$. The
marginalised results for these exhibit minimal deviation from the non-
marginalised approach, as visually depicted in the top row of Fig. 10. While
the cosmological parameter values are consistent between each statistical
approach and the samples, we find expected deviations between timescape and
spatially flat $\Lambda$CDM for each respective sample. Timescape utilises an
effective dressed parameter, calculated from $f_{\lower
2.0pt\hbox{$\scriptstyle\rm v0$}}$, and is not directly related to the
$\Omega_{\lower 2.0pt\hbox{$\scriptstyle\rm M0$}}$ of the spatially flat
$\Lambda$CDM model. In particular, because of various model limitations, we
emphasise that these values are not used to constrain the spatially flat
$\Lambda$CDM parameters.

Furthermore, these parameters consistently agree within 1–2$\sigma$ for the
frequentist and Bayesian values across different values of $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ (Fig. 9). This emphasises that these two
statistical methods are consistent. A more comprehensive statistical
comparison of the two cosmological models is carried out in Section 4.7 in
terms of a Bayesian analysis.

Figure 10: Median and standard deviation of 50 random subsamples of various
sizes from the Pantheon$+$ catalogue for the timescape model, based on the
non-marginalised (Eq. 19, left column) and the marginalised (Eq. 20, right
column) likelihoods to show how the different parameters change with sample
size. In the first row, we display the effective dressed parameter
$\Omega_{\rm M0}=\frac{1}{2}(1-f_{\lower 2.0pt\hbox{$\scriptstyle\rm
v0$}})(2+f_{\lower 2.0pt\hbox{$\scriptstyle\rm v0$}})$ from Eq. 8.

###  4.7 Model Comparison

Figure 11: Bayesian evidence $\ln{B}$ of the timescape model relative to the
spatially flat $\Lambda$CDM model when data $z<z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ is progressively removed and parameters
are refitted for the entire remaining sample $z\geq z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$. On the Jeffrey’s scale (Kass & Raftery,
1995), $1\leq|\ln{B}|<3$ represents “positive” relative Bayesian evidence,
while $|\ln{B}|<1$ is “not worth more than a bare mention”, $3\leq|\ln{B}|<5$
constitutes “strong” evidence and $|\ln{B}|\geq 5$ “very strong” evidence.
Here the evidence for either model is positive, rather than strong, with the
upper (blue) regions favouring timescape and the lower (red) regions
$\Lambda$CDM. The observed variation with $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ strongly suggests the presence of
additional systematic effects as expected in generic inhomogeneous
cosmologies, which will be presented in a forthcoming analysis. The vertical
dotted line represents the SHSα found at $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}=0.075$. The overall evidence is $0.80\pm
0.04$ for JLA580, $0.72\pm 0.04$ for JLA740 and $1.02\pm 0.04$ for P+580, all
favouring the timescape model. In contrast, the integrated evidence for the
P+1690 yields $-0.46\pm 0.04$, which favours the spatially flat $\Lambda$CDM
model.

We now present the synthesis of all previous results in terms of a Bayesian
model comparison, shown in Fig. 11. A Bayes factor for each pair of models is
determined by a joint refit of the parameters $\\{\alpha,x_{\lower
2.0pt\hbox{$\scriptstyle 1$}},\beta,c\\}$ on the r.h.s. of Eq. 10 for light-
curves and the energy–density fractions $\\{\Omega_{\lower
2.0pt\hbox{$\scriptstyle\rm M0$}}\\}$ ($\Lambda$CDM) or $\\{f_{\rm v0}\\}$
(timescape) that feed into $d_{\lower 2.0pt\hbox{$\scriptstyle L$}}$ on the
r.h.s. of Eq. 11. These correspond to empirical observations and cosmological
models151515The empty Milne universe has no free cosmological parameter, given
that the absolute magnitude $M_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}}$ and
Hubble constant $H_{\lower 2.0pt\hbox{$\scriptstyle 0$}}$ are marginalised
over as nuisance parameters. As discussed by Dam et al. (2017) this
complicates any Bayesian comparison. The Milne universe is of mathematical
interest as a limiting case but is not physically viable and is omitted in
this paper. respectively.

Within each panel, parameters (Eqs. 10 and 11) for all $z>z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ SNe Ia are refit as $z<z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ events are discarded, with $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ treated as a running variable.

The P+1690 panel has significant differences from the other panels. The
differences between the panels are nonetheless entirely consistent with the
expectations of the timescape cosmology, once all statistical distribution
biases are accounted for. In discussing Bayes factors, we must bear in mind
that the timescape model includes a deep principle – the Cosmological
Equivalence Principle (CEP)161616For an introductory essay on the CEP see
Wiltshire (2009a). (Wiltshire, 2008) – which provides an explanation for the
spatial isotropy of distant sources in large full sky surveys when averaged on
our celestial sphere. The CEP provides a strong contrast with the standard
cosmology in which average spatial isotropy (and homogeneity) are assumed.

Since both $\Lambda$CDM and timescape are compared here with precisely the
same SNe Ia data, the Bayes factor depends predominantly on how the
theoretical models depart from average isotropy and homogeneity. As discussed
in Section 1.1 in $\Lambda$CDM – or indeed in any FLRW model – one has
globally defined peculiar velocities, whereas for viable inhomogeneous
cosmologies, one must introduce the notion of non-kinematic differential
expansion (Bolejko et al., 2016). In our empirical Bayesian approach we could
compare any pair of FLRW models and/or non-FLRW general relativistic
cosmologies.

The specific feature of the timescape model which is compared here is the
uniform quasilocal Hubble expansion condition which automatically dispenses
with the Hubble tension (Wiltshire, 2007a, 2008, 2009b). It is the
mathematical model for cosmic structure and emergence of statistical
homogeneity which is the key to understanding the explanatory significance of
Fig. 11. With this insight in hand, a set of distribution effects are visible
across the four panels of Fig. 11, which we can group as follows.

####  4.7.1 Sensitivity

The P+1690 Pantheon+ full sample (lowest panel) contains significantly more
$z\lesssim z_{\rm min}$ than in the upper 3 panels. Since these additions to
the P+1690 distribution are at smallest cosmological distances where the
variation of the Hubble–Lemaître law is the greatest, we can expect to see a
greater sensitivity of the Bayes factors in this region. The results all
consistently favour the quasilocal uniform Hubble expansion hypothesis as
follows:

  1. (i)

Including the closest $z_{\rm min}\lesssim 0.01$ SNe Ia drives $\ln B$
evidence for the full sample from “not worth more than a bare mention” to
positive evidence in favour of timescape, and just below the $\ln
B\mathop{\sim}\limits 3$ threshold at which the Bayesian evidence is strong by
Jeffrey’s scale (Kass & Raftery, 1995). It is the value of the $z_{\rm min}=0$
intercept, $\ln B=2.93$, which drives the vertical scale of Fig. 11,
consistent with the expectation for a typical galaxy – like the Milky Way – in
a thin filament in a void.

  2. (ii)

The peak in the P+580 panel ($\ln B=2.55$) is significantly higher than in the
JLA740 panel ($\ln B=1.60$), both occurring within the range $0.02\lesssim
z_{\rm min}\lesssim 0.03$ or $60\,h^{-1}\mathop{\hbox{\,Mpc}}\lesssim
d_{\lower 2.0pt\hbox{$\scriptstyle L$}}\lesssim
95\,h^{-1}\mathop{\hbox{\,Mpc}}$. Since exactly the same SNe Ia events are
involved, this result is consistent with:

    1. a)

the more stringent colour selection criteria in Pantheon+;

    2. b)

the $z_{\rm min}=0.023$ cutoff used in the Riess et al. (2007) Gold07 MLCS2k2
analysis and results of the first Bayesian $\Lambda$CDM / timescape comparison
using that sample (Leith et al., 2007);

    3. c)

convergence (see Section 4.7.2 below);

    4. d)

the expected transition scale from the nearby nonlinear regime to the larger
scale linear regime;

    5. e)

the hypothesis that this feature is correlated with the Centaurus–Hydra
complex, i.e., the “Great Attractor” (Wiltshire et al., 2013; McKay &
Wiltshire, 2016).

  3. (iii)

The minimum $\ln B$ value in the P+1690 panel at $z_{\rm min}=0.04$ is
significantly lower than in the P+580 panel, with a Bayes factor $\ln B=-2.82$
that could naïvely be interpreted as favouring $\Lambda$CDM. However, the
P+580 subsample has a Bayes factor $\ln B=0.98$ at the same redshift cut in
favour of timescape, indicative of a distribution sampling bias. This
confirmed by the ordering $\ln B|_{\lower 2.0pt\hbox{$\scriptstyle\rm
P+580$}}>\ln B|_{\lower 2.0pt\hbox{$\scriptstyle\rm JLA580$}}>\ln B|_{\lower
2.0pt\hbox{$\scriptstyle\rm JLA740$}}>\ln B|_{\lower
2.0pt\hbox{$\scriptstyle\rm P+1690$}}$. These results are consistent with:

    1. a)

the more stringent colour selection criteria in Pantheon+;

    2. b)

the hypothesis that the feature is a statistical distribution bias resulting
from the degeneracy of $\beta c$ with $d_{\lower 2.0pt\hbox{$\scriptstyle
L$}}$ when $M_{\lower 2.0pt\hbox{$\scriptstyle\rm B$}}$ and $H_{\lower
2.0pt\hbox{$\scriptstyle 0$}}$ are marginalised over in Eq. 9 as nuisance
parameters;

    3. c)

convergence (see Section 4.7.2 below);

    4. d)

the hypothesis that this feature is correlated with structures leading to a
puzzling bulk flow in the $\Lambda$CDM/FLRW framework (Watkins et al., 2023;
Whitford et al., 2023). Structures along the
$\\{\ell,b\\}=\\{(297/117)\degree\pm 4\degree,(-6/6)\degree\pm 3\degree\\}$
axis, and in particular both those at “fine structure” redshifts
$z\mathop{\sim}\limits 0.033$–$0.06$ as well as the $z\lesssim 0.033$
foreground structures through which they are observed, are consistent with a
need for reinterpretation in the timescape framework.

Identifying structures along this axis is a challenge due to obscuration by
our own galaxy (Zone of Avoidance) compounding difficulties associated with
the identification of voids. The Chamaeleon overdensity in the South Pole Wall
centred at a redshift of $\sim 0.04$ (Pomarède et al., 2020) is one
overdensity matching what we denote a fine structure distance tension.

  4. (iv)

The $\ln B$ values for $z_{\rm min}>0.075$ to the right of the vertical dotted
line are close to $\ln B=0$ for the full P+1690 sample, but are consistently
at the marginally positive boundary $\ln B=1$ for all other cases. This is
consistent with Items (iii)a), (iii)b) and (iii)c) of Item (iii) above.

####  4.7.2 Convergence

In Figs. 5 and 10 bootstrap methods were applied to test the convergence of
the random subsamples of $N$ SNe Ia of the full P+1690 Pantheon+ as $N\to
1535$ SNe Ia (the P+1690 sample consists of 1535 unique SNe Ia). As was
already discussed in Section 4.4 this appears to be due to distribution biases
arising from structures in the range $0.04\lesssim z_{\rm min}\lesssim 0.06$.
In Fig. 11 these effects are associated with

  1. (i)

the height of the $\ln B$ peak value in the JLA580 and P+580 panels;

  2. (ii)

the magnitude of the minimum $\ln B$ in the P+1690 and P+580 panels;

  3. (iii)

$\ln B$ values for $z_{\rm min}>0.075$ (see Section 4.7.2, Item (iv) above).
We note that in the timescape cosmology it is expected that for very
large171717Here “very large” is determined by the radius of convergence of at
Taylor series expansion of the distance–redshift law relative to the precision
of datasets used. $z_{\rm min}\gg 0.033$, $\ln B\to 0$. The fact that the
upper 3 panels of Fig. 11 have $\ln B\to 1$ is consistent with:

    1. a)

$z_{\rm min}\mathop{\sim}\limits 0.1$ being a scale at which very small
corrections to average expansion law from inhomogeneities can still be
detected;

    2. b)

the fact that Malmquist distribution biases have thus far not been subtracted
in our analysis.

The key finding that our confirmation of convergence in bootstrap resampling
is associated with the value of $\beta c$ in Eq. 9 is deferred to Section 5.

####  4.7.3 Novel Selection Criteria

Sections 4.7.1 and 4.7.2 have features that suggest further refinement of our
empirical definition of the SHSα statistic is likely to improve the
subtraction of peculiar velocities/expansion in very large distance–redshift
surveys. The clearest signal for this are the systematic differences of $\ln
B$ values for $z_{\rm min}>0.075$ and the points made in Section 4.7.2, Item
(iii) above.

In particular, now that an old debate (Nielsen et al., 2016; Rubin & Hayden,
2016; Dam et al., 2017; Colin et al., 2020) is resolved, we see that the most
significant aspects of each viewpoint were correct within the context of
assumptions about unknown systematics. Each researcher in the debate brought
different perspectives – whether instrumentational, phenomenological or
theoretical. Systematics due to the interplay of model assumptions about
cosmic expansion and calibration of the entire available SNe Ia catalogue are
key to understanding the convergence of $\alpha x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}\to 0.0204\pm 0.0010$, which defines SHSα. Its
convergence can be seen in Fig. 6, in contrast with previous debates that
focussed on SHSg2, which fails to converge for the samples considered. In
particular, the $\beta$ parameter from the light-curve analysis does not
converge as $N\to N_{\rm u}=1535$ (for the full P+1690 catalogue) as seen in
Fig. 5 which leads to a lack of convergence for the SHSg2.

Inhomogeneity from cosmic structure underlies the key systematic interpreted
as peculiar velocities in the $\Lambda$CDM model and peculiar expansion in the
timescape model, and which drives the Hubble tension (Di Valentino et al.,
2021).

Peculiar velocity modelling based on FLRW assumptions is routinely included in
many SNe Ia data releases. We needed to remove such modelling for our non-FLRW
analysis of Pantheon+. While our Bayesian comparison based on SHSα appears
sufficient to resolve the largest contributions to the Hubble tension, Fig. 6
shows that new selection criteria based on the SHSα should already be
sufficient to isolate smaller higher order effects of cosmic inhomogeneity in
the Pantheon+ catalogue. This is supported by Fig. 11: for P+1690 we find $\ln
B\sim 0$ beyond the SHSα instead of $\ln B\sim 1$.

The most significant result is that the overall Bayesian evidence from
bootstrapped resampling of Pantheon$+$ takes a value $\ln B=1.02\pm 0.04$ in
favour of timescape, contrary to the naïve inference one might make by visual
inspection of the bottom panel of Fig. 11. While this is only at the
marginally positive boundary of the Jeffrey’s scale, systematic effects such
as Malmquist biases have to still to be accounted for using the broadened
definition of model independence. Quite simply: by removing the source of the
Hubble tension, the Bayesian evidence in favour of timescape has increased.

It should be stressed that this is completely independent of the goodness of
fit of the timescape cosmology.

##  5 Discussion

The updated analysis of the timescape and spatially flat $\Lambda$CDM models,
utilising new data of over twice the size, has uncovered significant
discoveries regarding both the structure of the catalogue and the evaluation
of the models.

The Pantheon$+$ results use a covariance matrix of distance moduli, assuming a
fiducial $w$CDM cosmology. In order to compare the fundamental differences
between cosmological models, or to identify any inherent biases in the current
standard of cosmology, it is crucial to perform cosmology-dependent steps as
late as possible in the data-reduction pipeline. Due to the cosmology-
dependent bias-correction simulations that were introduced by Pantheon$+$ and
applied to the SALT2mu algorithm, along with the inclusion of peculiar
velocity modelling in their final Hubble diagram, it was not possible to
incorporate these calibrations in the calculation of the covariance matrix.
Consequently, there remains the possibility of uncorrected selection biases in
these results. However, we presume that these biases are of a similar
magnitude for both models. This approach allows a fair comparison between
cosmological models not built upon the FLRW geometry. In order to mitigate
potential issues and biases, we cross-reference our results with the
observations reported in the JLA catalogue by Dam et al. (2017).

Prior to the marginalisation process, we recover the differences in the SALT2
surfaces reported by Taylor et al. (2021), when comparing the JLA and
Pantheon$+$ catalogues. It is important to account for the inherent variations
within the published datasets when making comparisons between different SNe Ia
catalogues.

As there are significant differences in the results of the Pantheon$+$/JLA
common subsample and the full Pantheon$+$ sample, we conduct the analysis for
random subsamples (see Section 4.3 for more details) for both the JLA and the
Pantheon$+$ datasets to investigate the selection effects. The dispersion of
the randomised results for the Pantheon$+$ sample is greater than JLA, having
significant scatter at low-redshift cuts. Notably, $42\%$ of unique SNe Ia in
the Pantheon$+$ catalogue have a redshift of $z_{\lower
2.0pt\hbox{$\scriptstyle\rm CMB$}}\leq 0.075$. Given the large number of low-
redshift observations, there is considerable scatter in the $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ and $c$ parameters for $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}<0.0225$. Hence, light-curve parameters
deduced using these values must consider the uncertainties stemming from the
variations in the observed parameters at low-redshift ranges.

One reason for the scatter is the fact that the supernova observations within
the Pantheon$+$ catalogue are derived from diverse samples and thus the value
adopted for $H_{\lower 2.0pt\hbox{$\scriptstyle 0$}}$ differs. To account for
this, we include marginalisation over the zero-point offset in the definition
of the likelihood function. Such an approach effectively controls the level of
scatter.

As the redshift cut increases, the scatter becomes negligible. However, the
light-curve parameters only converge beyond the scale of statistical
homogeneity. To evaluate this quantitatively, we rely on a power–law method
fitted to $\alpha x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ presented in
Section 4.4 and quantify its approach to convergence.

For this data driven approach we take a percentage range of the $\alpha
x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ limiting value for high $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}$ which is motivated by a few key reasons.
One rationale is the need for this method to be applicable to future datasets,
and we cannot assume uniformity in the level of noise across diverse datasets.
Our method is also more robust to datapoints that have not yet reached the
SHSα or display numerical instability.

Different levels of statistical noise may lead to earlier or later
convergence. We chose to use 5% as it provides more constraints on the SHSα
for the convergence while having a higher tolerance for noise and fitting
constraints than a 3% uncertainty does. In Pantheon$+$, the 5% bound for the
chosen range corresponds to 2.6$\sigma$ for $\alpha x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ using a spatially flat $\Lambda$CDM model.
However, because of the larger values of $\beta$ a $1$–$3$% bound may be more
appropriate for SHSβ which was not possible to define here due to the lack of
convergence.

Nonetheless, this approach does exhibit several limitations. Firstly,
power–law functions tend to approach negative infinity as $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}\to 0$ and it is obvious that the $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ parameter remains finite and bounded across all
redshifts. This highlights the important need for a well-defined lower bound
when applying this model-fitting method. Modifying the power–law by
introducing a fourth parameter, such as $-a(z-z_{\lower
2.0pt\hbox{$\scriptstyle 0$}})^{-n}+C$, does marginally improve the $\chi^{2}$
statistic and, for negative values of $z_{\lower 2.0pt\hbox{$\scriptstyle
0$}}$ prevents the power–law from diverging for $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}\to 0$. However, this modification does not
significantly alter the resulting SHSα and instead increases the Akaike
information criterion (AIC)181818The Akaike information criterion (AIC) is a
criterion for the evaluation of statistical models based on the goodness of
fit and the number of estimated parameters. For more information see Akaike
(1974) and Dam et al. (2017, Appendix C)..

The second limitation pertains to the fitting process itself. While fitting a
power–law, any numerical spikes present in the data can significantly
influence the resulting best-fit curve.

Lastly, the power–law method relies on the convergence of the $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ parameter, which was not the case for JLA. The
$x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ convergence issue introduces a
challenge when comparing results across different catalogues that may not
control the redshift dependence of $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$
as effectively as Brout et al. (2022b) did in Pantheon$+$. Due to the
degeneracy of $\alpha$ and $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$,
convergence issues within one of the parameters can be counterbalanced by the
other, making the treatment more robust when considering the product.

As discussed previously in Section 4.4 and shown in Fig. 7, the parameters
$\beta$ and $c$ do not converge and thus a SHSβ based on these is not
considered as an alternative to the SHSα. This is likely due to the Malmquist
biases present at redshifts greater than $z_{\lower
2.0pt\hbox{$\scriptstyle\rm min$}}\sim 0.05$ which we could not yet account
for in the timescape model. Therefore, we also exclude the correction from
spatially flat $\Lambda$CDM for consistency. Applying corrections for this
bias to timescape would be an important step to improve our sample and confirm
the SHSα by a second set of parameters.

The result we obtained for the SHSα with the power–law method significantly
differs from the one found by Dam et al. (2017). These differences stem from
inherent biases within the distinct datasets, emphasising the need for robust
controls to assess their convergence behaviours effectively. Such calibration
was done effectively by Pantheon$+$ for the $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ parameter. Further cosmology-independent
treatment of such issues is needed for a more comprehensive examination of the
SHSα.

Moreover, random sampling has a tangible impact on the parameter values. For
$\alpha$, $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$ and $c$, these changes are
not significant and the results converge towards that of the full sample when
increasing the sample size. However, the most significant change is in the
$\beta$ parameter for varying subsample sizes. For this parameter, we do not
observe a convergence, as the final value differs significantly from the ones
found for lower $N$. Considering marginalised likelihoods does not resolve
this feature as it changes the absolute values but not the behaviour for
varying sample sizes. Consequently, the results for $\beta$ cannot be
attributed to the choice of $H_{\lower 2.0pt\hbox{$\scriptstyle 0$}}$. As
shown in Fig. 5, selecting subsamples with different redshift distributions
resolves this issue. Therefore, for a thorough investigation into the
convergence of the light-curve parameters, a sample with uniform redshift
distribution would be needed. Constructing such a sample will require more SNe
Ia particularly at $z_{CMB}\gtrsim 0.04$.

A noteworthy observation regarding the light-curve parameters of the
Pantheon$+$ dataset, specifically $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}$,
is that they converge for $z_{\lower 2.0pt\hbox{$\scriptstyle\rm
min$}}=0.075^{+0.007}_{-0.009}$. In contrast, the dependency on redshift for
the light-curve parameters in the JLA catalogue indicates the presence of an
accelerated expansion. This feature has been in heavy debate for the JLA
sample. Nielsen et al. (2016) enforced redshift independence in $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ and $c$ in their analysis, which was criticised
by Rubin & Hayden (2016). However, Rubin & Hayden (2016) added twelve
additional parameters to the Bayesian Hierarchical Model, and their inclusion
has been challenged by several papers, notably, Dam et al. (2017) and Colin et
al. (2020). Dam et al. (2017) found that the extra parameters did not increase
$\ln B$, while Colin et al. (2020) criticises their addition by analysing the
Bayesian information criterion (BIC)191919The Bayesian information criterion
(BIC) is closely related to the Akaike information criterion (AIC) but
introduces a larger penalty term for $k\geq 8$. The BIC is given by: ${\rm
BIC}=k\ln{N}-2\ln{\mathcal{(L)}}$. For more information see Kass & Raftery
(1995) and Dam et al. (2017, Appendix C).. When the same supernova analysis as
in Dam et al. (2017) is conducted on the Pantheon$+$ data, the lack of
redshift dependence indicates that these features from JLA were due to
different calibration systematics and incomplete datasets. In particular,
since $x_{\lower 2.0pt\hbox{$\scriptstyle 1$}}\to\mathop{\hbox{\,const}}$ for
$z_{\lower 2.0pt\hbox{$\scriptstyle\rm min$}}>\mathop{\hbox{\,SHS}}_{\alpha}$,
the systematic error in the JLA catalogue that was at the centre of the debate
has been controlled in Pantheon$+$.

The Bayesian analysis of both the JLA data and the Pantheon$+$/JLA common
subsample agree with the results by Dam et al. (2017). When we extend the
analysis to encompass the complete P+1690 dataset, certain artefacts emerge,
stemming from distribution-related effects. In particular, the Pantheon$+$
catalogue incorporates a substantial number of low-redshift observations that
fall below the SHSα. In this low-redshift regime, the fundamental assumption
of statistical homogeneity does not hold, and consequently, the timescape
model demonstrates a more favourable fit to the data. The P+1690 global
minimum $\ln B_{\lower 2.0pt\hbox{$\scriptstyle\rm P+1690$}}=-2.82$, (c.f.,
Section 4.7.2, Item (iii)) which naïvely favours $\Lambda$CDM was found to be
a distribution effect. This preference arises due to specific structural
features present in the data and the partial reliance on the spatially flat
$\Lambda$CDM cosmological framework to account for these features.

Future SNe Ia surveys such as those from Euclid, the Rubin Observatory, and
the Roman Space Telescope will vastly expand on the SNe Ia sample size and
allow for more detailed examination of large scale structures (Euclid
Collaboration et al., 2020, 2022; Ivezić et al., 2019; Rose et al., 2021). For
these detailed analyses to be possible, the data products from the surveys
must be presented free from underlying cosmological models, as was achieved
here with Pantheon$+$. For example, subtle biases that favour spatially flat
$\Lambda$CDM can be introduced when correcting for sample related biases such
as the Malmquist bias.

##  6 Conclusions

We have performed a from–first–principles reanalysis of the best available SNe
Ia catalogue, Pantheon+ (Scolnic et al., 2022). For the first time in 25
years, we find evidence that this catalogue already contains such a wealth of
data that with further reanalysis, a genuine paradigm shift may soon emerge.
Supernovae Ia were the events that first convinced a majority of the community
of the need for an observational paradigm shift (Riess et al., 1998;
Perlmutter et al., 1999). That paradigm shift invoked a mysterious “dark
energy”, which has always been recognised as a place-holder for new physics.
All major paradigm shifts in physics have involved the complex interplay of
concepts, experimental phenomenology and mathematical logic. While such
questions are likely to still play out for decades, in our view some
theoretical new physics features of the “dark energy” place-holder are already
emerging.

The observed spatial isotropy and homogeneity of the most distant sources has
been a deep theoretical conundrum ever since relativistic cosmology was first
constructed a century ago. By necessity, the first mathematical cosmologists
put the answer into their models as an assumption, and elevated the assumption
to the Cosmological Principle. Rather than assuming spatial isotropy and
homogeneity, any appeal to a deeper physical principle should seek to explain
such puzzling observations. The timescape cosmology does this by returning to
fundamental questions left unanswered in the first formulations of
relativistic cosmology, via its Cosmological Equivalence Principle (CEP)
(Wiltshire, 2008), an extension of the suite of Einstein’s equivalence
principles. Rather than seeking new physics in direct violation of the known
equivalence principles, in the timescape model each of Einstein’s principles
remains intact and embedded within the others according to its scale of
validity.

In directly addressing the question “what is the origin of scale?”, the CEP
provides a deeper principle for operationally determining the relevant scales
SHSα, SHSβ – as required by the new selection criteria that emerged from our
data-driven bootstrapping in Sections 4.7.2 and 4.7.3. As a direct consequence
of our analysis, timescape’s proposed resolution of the Hubble tension202020A
potential resolution of the Hubble tension is a generic mathematical feature
of inhomomogeneous GR cosmologies (Bolejko, 2018). The important question for
the timescape scenario has been: Is it compelling? has uncovered deeper
understanding of many issues that may signal the beginning of a paradigm
shift. The term paradigm shift was first applied in relation to fundamental
physics (Kuhn, 1962), making our results a prime test case for the other
important criterion of the philosophy of science: falsifiability (Popper,
1935).

In the conceptual framework of Kuhn (1962), major paradigm shifts occur over
timescales of decades to centuries when multiple failures of “normal science”
(incremental research) are followed by the resolution of anomalies through the
novel predictions of a new paradigm. Since the preconditions of a genuine
paradigm shift in cosmology are already well documented (Aluri et al., 2023;
Peebles, 2022), let us propose the

null hypothesis that:

the timescape resolution of the Hubble tension (TRHT) is not accompanied by
any novel selection criteria in fundamental cosmology.

Following Popper’s (1935) criterion, we should seek just one compelling
falsification of the null hypothesis to rule it out. We may verify the TRHT in
this manner, but not ever prove it.

As yet, we are not able to rule out the timescape. Rather the following
quantitative observations of the TRHT resolve puzzles noted over the past 25
years with smaller sparser datasets:

  * •

The results of Smale & Wiltshire (2011) showed slight Bayesian preference for
timescape over $\Lambda$CDM using MLCS2k2 and slight Bayesian preference for
$\Lambda$CDM over timescape using SALT/SALT2 fitters as then understood. The
maximum evidence in “favour” of either model had $|\ln B|\lesssim 1.68$,
depending on precise $z_{\rm min}$ cuts, which outliers were excluded etc.
Detailed discussion of two decades of advances in light curve fitting, and the
historical treatment of systematic uncertainties, is beyond the scope of this
paper. Nonetheless, looking into details one finds, e.g., that with further
analysis a sequence of orderings of $|\ln B|$ analogous to Section 4.7.1 (iii)
is to be found in Smale & Wiltshire (2011, Tables 6–9, Fig. 8).

  * •

Smale & Wiltshire (2011) also discussed the fact that for timescape to be
preferred in Bayesian analysis, the extinction by dust in the Milky Way is
needed to be typical compared to other galaxies, which is in agreement with
independent direct measurement of Finkelman et al. (2008, 2010). In two
studies of nearby galaxies these authors found values of the reddening
parameter $R_{\lower 2.0pt\hbox{$\scriptstyle\rm V$}}=2.82\pm 0.38$,
$R_{\lower 2.0pt\hbox{$\scriptstyle\rm V$}}=2.71\pm 0.43$ for their two
samples of 7 and 8 galaxies respectively, indicating that the Milky Way value
$R_{\lower 2.0pt\hbox{$\scriptstyle\rm V$}}=3.1$ was within a standard
deviation of expectation. The Milky Way value had been used in the MLCS based
values, but was not included as a parameter in the SALT/SALT2 methods leading
to potential degeneracies with free parameters.

While direct comparison of the numerical value for the $\beta$ parameter that
was favoured over a decade ago is complicated by the advances in light curve
fitting, what remains consistent is the fact that higher relative values of
the $\beta$ parameter were favoured by timescape as compared to $\Lambda$CDM
(Smale & Wiltshire, 2011). In the present analysis, the resulting $\beta$
values are much higher than considered previously. However, the
$\Lambda$CDM/timescape differences are negligible compared to the deviations
caused by choosing different statistical approaches (Table 1). In particular,
using a likelihood that marginalises over the zero-point offset degenerate
with the Hubble constant, significantly reduces the value of $\beta$. This is
thus related to the Hubble tension.

  * •

Furthermore, for the timescape model to be a consistent description of the
sample, an apparent “Hubble bubble” – as mooted by Zehavi et al. (1998) –
should be present in the data, as discussed by Smale & Wiltshire (2011).
Conceptually this strikes at the heart of the peculiar motion versus expansion
dichotomy that the timescape cosmology resolves. Once again, this relates to
the $\beta$ values obtained by the different models and statistical analyses,
which in turn may resolve the Hubble tension, as described before.

Our in depth comparison between the spatially flat $\Lambda$CDM and timescape
cosmological models was only possible with Pantheon$+$ as we were able to
remove all cosmology-dependent calibration and correction steps. Without the
standard corrections, biases such as the Malmquist bias are present in our
analysis, however, we expect the resulting differences to be small when
comparing models. From this extensive dataset of 1690 SNe Ia observations we
follow the analysis of Dam et al. (2017) and construct a covariance matrix
from the SALT2 parameters augmented with relevant perturbations and the
statistical covariance from Pantheon$+$.

In our analysis for the parameters $\alpha$, $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ and $c$, we find significant differences
between JLA and Pantheon$+$. We find a convergent behaviour of $x_{\lower
2.0pt\hbox{$\scriptstyle 1$}}$ with $z_{\lower 2.0pt\hbox{$\scriptstyle\rm
min$}}$ due to resolved systematic errors which resolves a previous debate
(Nielsen et al., 2016; Rubin & Hayden, 2016; Dam et al., 2017; Colin et al.,
2020). In contrast, the $\beta$ parameter does not converge for varying sample
size. To fully understand the nature of the $\beta$ convergence further tests
are required with SNe Ia datasets that have a uniform redshift distribution.

Furthermore, we consider the SHS in terms of a data-driven definition of
SHS${}_{\alpha}=0.075^{+0.007}_{-0.009}$ which is substantially larger than
the conventional SHS${}_{{\rm g}2}\mathop{\sim}\limits 0.033$. By relying
solely on observations, our approach gives a true representation of the
underlying structure in the Universe. In particular, while the conventional
2-point galaxy correlation function measure SHSg2 picks out gross features of
a transition to average homogeneity (when refined) which does coincide with
known large scale structures such as the Great Attractor, the SHSα statistic
reveals further “fine structure” in the transition to homogeneity at scales
corresponding to more distant structures. For future analyses of the SHS, an
additional consideration of the SHSβ would enable a more thorough
investigation. For this to be possible, a treatment of the Malmquist bias’,
which was unaccounted here, is needed.

Considering the Bayesian comparison of the two cosmological models, we find
that timescape is positively favoured in the JLA740, JLA580 and P+580
datasets. However, the extension of the full Pantheon$+$ catalogue, reveals
important artefacts that highlight the need for analysing the underlying
cosmic structures through which we observe SNe Ia. Instead of assuming uniform
homogeneity throughout the entire Universe, the indications favouring the
timescape model at lower redshifts highlight the necessity of accounting for
inhomogeneities. The significant changes in the value for the $\ln B$ can be
explained from the cosmic structures. To further investigate this and to
determine which model is to be definitively preferred, more SNe Ia at high
redshifts are required to be calibrated in a cosmology-independent manner.

As final remarks we observe that one of the hallmarks of a genuine paradigm
shift is that at tipping points, many puzzles tend to fall into place in rapid
succession. This article is the first of several anticipated in the coming
year. Some of the research directly related to the new characterisation of
statistical homogeneity which is underway includes:

  * •

A deeper phenomenological understanding of the Newtonian limit in cosmology
without a background FLRW geometry promises to be found. This may explain why,
despite violation of the Buchert-Ehlers theorem (Buchert & Ehlers, 1997), the
heuristic computational scheme of Rácz et al. (2017) nonetheless produces
results which are close to those of the timescape cosmology.

  * •

Szekeres models (Szekeres, 1975) have a rich structure that can be used as
effective models on $\mathop{\sim}\limits$SHS scales which can be directly
constrained by large scale peculiar velocity / peculiar expansion datasets
(Hills, 2022; O’Keeffe, 2023). Future precision in such work relies on the
definition of SHSβ which was not yet possible due to the lack of cosmology-
independent treatment of Malmquist bias. This is left for future work.

  * •

The BAO scale has been extracted directly from BOSS data (without FLRW
assumptions) for the timescape cosmology (Heinesen et al., 2019). This is
important for future tests at scales of interest. However, calibration of this
scale relative to the BAO at the CMB last-scattering surface requires
understanding the interplay general relativistic nonlinearities with matter at
the galactic level. First steps to achieve this require understanding why
particular simple classes of previously studied models and novel
generalisations of those solutions fail (Galoppo, 2023), leaving open the
possibilities that may work in future.

General relativity is a deep and beautiful theory, whose principles remain to
be fully explored – and maybe extended – to better understand the Universe.
The hypotheses tested and further advanced here, if correct, indicate that we
will likely be exploring those consequences for decades to come.

## Acknowledgements

We are deeply grateful to the Pantheon+ team for supplying the systematic
perturbation files used in their analysis and for providing extensive feedback
on its implementation. In particular, we thank Dillon Brout for his perceptive
insights and hands–on assistance with our revisit of the Pantheon+ data
reduction pipeline beyond conventional perturbed FLRW models. RRH and ZGL are
supported by the Marsden Fund administered by the Royal Society of New
Zealand, Te Apārangi under Fast Start Grant MFP-UOC2204. RRH is also supported
by the Rutherford Foundation Postdoctoral Fellowship under the grant RFT-
UOC2203-PD. AS thanks the German Academic Exchange Service for supporting her
stay at the University of Canterbury. DLW thanks Subir Sarkar for hosting a
Visiting Oxford Erskine Fellowship, and Linacre College and the Rudolf Peierls
Centre for Theoretical Physics for their generous support. We are grateful to
all members of the University of Canterbury Gravity and Cosmology group for
stimulating discussions: in particular Marco Galoppo, Christopher Harvey-
Hawes, Morag Hills, Shreyas Tiruvaskar and Michael Williams. Finally, we thank
Lawrence Dam, Tamara Davis, Asta Heinesen, Roya Mohayee, Eoin Ó Colgáin,
Mohamed Rameez, David Rubin, Subir Sarkar, Dan Scolnic, Nathan Secrest and
Shahin Sheikh-Jabbari for correspondence, discussions, debates and the
dispassionate scientific consensus–building that made this analysis possible.

## Data Availability

A complete set of the codes and details used for our analysis and how to use
them can be found at https://github.com/antosft/SNe-PantheonPlus-Analysis.

The JLA catalogue and their covariance can be found at
http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/568/A22 and
https://supernovae.in2p3.fr/sdss_snls_jla/ReadMe.html respectively. The
Pantheon$+$ calibration/systematic perturbation files and data that were used
in this paper were obtained via private communication with Dillon Brout and
Dan Scolnic.

## References

  * Abbott et al. (2017) Abbott B. P., et al., 2017, Astrophysical Journal Letters, 848, L12
  * Akaike (1974) Akaike H., 1974, IEEE Transactions on Automatic Control, 19, 716
  * Alam et al. (2017) Alam S., et al., 2017, Monthly Notices of the Royal Astronomical Society, 470, 2617
  * Alcock & Paczynski (1979) Alcock C., Paczynski B., 1979, Nature, 281, 358
  * Aluri et al. (2023) Aluri P. K., et al., 2023, Classical and Quantum Gravity, 40, 094001
  * Bentivegna & Bruni (2016) Bentivegna E., Bruni M., 2016, Physical Review Letters, 116, 251302
  * Betoule et al. (2014) Betoule M., et al., 2014, Astronomy & Astrophysics, 568, A22
  * Blas et al. (2011) Blas D., Lesgourgues J., Tram T., 2011, J. Cosmology Astropart. Phys., 2011, 034
  * Bolejko (2018) Bolejko K., 2018, Physical Review D, 97, 103529
  * Bolejko et al. (2016) Bolejko K., Nazer M. A., Wiltshire D. L., 2016, Journal of Cosmology and Astroparticle Physics, 06, 035
  * Brout et al. (2022a) Brout D., et al., 2022a, Astrophysical Journal, 938, 110
  * Brout et al. (2022b) Brout D., et al., 2022b, Astrophysical Journal, 938, 111
  * Buchert (2000) Buchert T., 2000, General Relativity and Gravitation, 32, 105
  * Buchert (2001) Buchert T., 2001, General Relativity and Gravitation, 33, 1381
  * Buchert & Ehlers (1997) Buchert T., Ehlers J., 1997, Astronomy & Astrophysics, 320, 1
  * Buchert & Räsänen (2012) Buchert T., Räsänen S., 2012, Annual Review of Nuclear and Particle Science, 62, 57
  * Buchert et al. (2020) Buchert T., Mourier P., Roy X., 2020, General Relativity and Gravitation, 52, 27
  * Buchner et al. (2014) Buchner J., et al., 2014, Astronomy & Astrophysics, 564, A125
  * Cheng et al. (1979) Cheng E. S., Saulson P. R., Wilkinson D. T., Corey B. E., 1979, Astrophysical Journal, 232, L139
  * Clarkson et al. (2008) Clarkson C., Bassett B., Lu T. H.-C., 2008, Physical Review Letters, 101, 011301
  * Colin et al. (2020) Colin J., Mohayaee R., Rameez M., Sarkar S., 2020, A response to Rubin & Heitlauf: "Is the expansion of the universe accelerating? All signs still point to yes" (arXiv:1912.04257)
  * Condon et al. (1998) Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., Perley R. A., Taylor G. B., Broderick J. J., 1998, Astronomical Journal, 115, 1693