# New Constraints on Anisotropic Expansion from Supernovae Type Ia

W. Rahman1, R. Trotta1,2,3, S. S. Boruah4, M. J. Hudson5,6,7, and D. A. van
Dyk2,8
1Astrophysics Group, Physics Department, Blackett Laboratory, Imperial College
London, Prince Consort Rd, London SW7 2AZ
2Data Science Institute, William Penney Laboratory, Imperial College London,
London SW7 2AZ
3International School for Advanced Studies (SISSA), Theoretical and Scientific
Data Science Group, Physics Department, Via Bonomea 265, 34136 Trieste, Italy
4Department of Astronomy and Steward Observatory, University of Arizona, 933 N
Cherry Ave, Tucson, AZ 85719, USA
5Department of Physics and Astronomy, University of Waterloo, Waterloo, ON,
N2L 3G1, Canada
6Waterloo Centre for Astrophysics, Waterloo, ON, N2L 3G1, Canada
7Perimeter Institute for Theoretical Physics, 31 Caroline St. N., Waterloo,
ON, N2L 2Y5, Canada
8Statistics Section, Mathematics Department, Huxley Building, Imperial College
London, Queen’s Gate, London SW7 2AZ E-mail<EMAIL_ADDRESS>
(Accepted 2022 April 22. Received Received 2022 April 22; in original form
2021 August 16)

###### Abstract

We re-examine the contentious question of constraints on anisotropic expansion
from Type Ia supernovae (SNIa) in the light of a novel determination of
peculiar velocities, which are crucial to test isotropy with SNe out to
distances $\sol 200/h$ Mpc. We re-analyze the Joint Light-Curve Analysis (JLA)
Supernovae (SNe) data, improving on previous treatments of peculiar velocity
corrections and their uncertainties (both statistical and systematic) by
adopting state-of-the-art flow models constrained independently via the 2M$++$
galaxy redshift compilation. We also introduce a novel procedure to account
for colour-based selection effects, and adjust the redshift of low-$z$ SNe
self-consistently in the light of our improved peculiar velocity model.

We adopt the Bayesian hierarchical model BAHAMAS to constrain a dipole in the
distance modulus in the context of the $\Lambda$CDM model and the deceleration
parameter in a phenomenological Cosmographic expansion. We do not find any
evidence for anisotropic expansion, and place a tight upper bound on the
amplitude of a dipole, $|D_{\mu}|<5.93\times 10^{-4}$ (95% credible interval)
in a $\Lambda$CDM setting, and $|D_{q_{0}}|<6.29\times 10^{-2}$ in the
Cosmographic expansion approach. Using Bayesian model comparison, we obtain
posterior odds in excess of 900:1 (640:1) against a constant-in-redshift
dipole for $\Lambda$CDM (the Cosmographic expansion). In the isotropic case,
an accelerating universe is favoured with odds of $\sim 1100:1$ with respect
to a decelerating one.

###### keywords:

supernovae: general, methods: statistical, (cosmology:) dark energy,
cosmology: observations, (cosmology:) cosmological parameters

††pubyear: 2021††pagerange: New Constraints on Anisotropic Expansion from
Supernovae Type Ia–A

##  1 Introduction

A fundamental assumption underpinning the cosmological concordance model is
the cosmological principle, namely that the universe is homogeneous and
isotropic on sufficiently large scales. Given the ubiquity of the cosmological
principle, an observational test of this assumption is an important step
towards validating our best description of the large scale universe. Testing
of homogeneity is hampered by the need of surveying extremely large scales
(see Maartens (2011)), although recent studies have found the transition to
homogeneity at high ($z\sim 2$) redshift consistent with expectations from the
$\Lambda$CDM cosmological concordance model (Gonçalves et al., 2018, 2020).

The assumption of isotropy has been tested over a range of redshifts and with
many different probes, from the relatively local universe out to the redshift
of recombination. Analyses of Cosmic Microwave Background (CMB) anisotropies
data obtained by the Wilkinson Microwave Anisotropy Probe (WMAP) and the
Planck satellite found up to $\sim 3\sigma$ evidence of breaches of
statistical isotropy in the form of power asymmetry between hemispheres,
multipole alignments, anomalous clustering of directions, although the
significance of these results is difficult to assess, partially because of
issues of a posteriori testing (Bennett et al., 2013; Akrami et al., 2014;
Schwarz et al., 2016; Planck Collaboration et al., 2016). Quasar polarization
directions also appear to be aligned along anomalous directions in the CMB
(Hutsemékers et al., 2005) and with coherence scales in excess of 500 Mpc
(Hutsemékers, D. et al., 2014), in potential disagreement with the
cosmological principle. Investigating the distribution of galaxies on large
scales, Sarkar et al. (2019) found however good agreement between the
predictions of $\Lambda$CDM and the Sloan Digital Sky Survey data, with a
transition to isotropy observed beyond a length scale of $200h^{-1}$ Mpc
(where $h$ is the dimensionless Hubble-Lemaître parameter). More recently,
Secrest et al. (2020) reported a one-sided 4.9-$\sigma$ rejection of the
hypothesis that the dipole in a sample of 1.3 million quasars is purely due to
our motion with respect to the CMB.

Supernovae Type Ia (SNIa) can be used to test the second expression of the
cosmological principle, namely that the expansion of the universe is
isotropic. SNIa are a sub-class of supernovae (SNe), resulting from the
thermonuclear explosion of CO white dwarfs accreting mass near the
Chandrasekhar limit, whose spectra exhibit no hydrogen lines but strong
silicon lines. A series of corrections can be applied to account for
correlations of absolute peak magnitude with their light curve decline rate
and colour. The correlation with light curve decline rate were first observed
by Rust (1974) and further corroborated later by Pskovskii (1977, 1984) and
Phillips (1993) using larger samples of SNe type Ia. Later a correlation with
color was also observed and it was noticed that brighter SNIa had a bluer
colour (Hamuy et al., 1995; Riess et al., 1996; Perlmutter et al., 1997).
Empirical corrections that account for correlations of absolute peak magnitude
with their decline time and colour can then be applied to standardise the SNIa
data (Tripp, 1998; Phillips et al., 1999). Light curves of SNe type Ia can be
standardised so that the residual scatter of their peak $B$-band magnitude is
sufficiently small ($\sim 0.1$ mag) to infer cosmological parameters, as was
first demonstrated by Riess et al. (1998) and Perlmutter et al. (1999). SNe
type Ia observations can also be used to test the hypothesis of isotropy in
the expansion of the universe underpinning the Friedemann-Lemaître-Robertson-
Walker (FLRW) metric of the concordance cosmological model, which exhibits an
isotropic scale factor $a(t)$. To this end, various authors have analyzed
increasingly large SNIa compilations with different statistical approaches,
often with sharply discordant results.

After early works (Kolatt & Lahav, 2001; Schwarz & Weinhorst, 2007; Gupta et
al., 2008), Cooke & Lynden-Bell (2010) analysed a subset of 250 SNIa from the
Union compilation (Kowalski et al., 2008) with $z>0.2$ with a maximum
likelihood approach to constrain a dipolar modulation to the luminosity
distance, finding no significant deviation from isotropy. However, Cai & Tuo
(2012) claimed that the deceleration parameter shows a preferred direction in
the Union2 (Amanullah et al., 2010) compilation of 557 SNIa, a result
corroborated by the analysis of Antoniou & Perivolaropoulos (2010), who
combined SNIa data with other cosmological probes. Jiménez et al. (2015)
analysed the same Union2 data, additionally including the SNLS3 data, and
showed that previous claims of anisotropy disappear if one accounts for
correlations among the observations by including the full data covariance
matrix in the analysis. Other null results of anisotropic expansion include
Heneka et al. (2014) and Lin et al. (2016a), who investigated the Joint Light-
curve Analysis (JLA) compilation (Betoule et al., 2014) of 740 SNIa. (See also
Sun & Wang (2019), who obtain discrepant results from three different
compilations of SNIa namely Union 2.1, JLA and Constitution (Hicken et al.,
2009). Similarly Andrade et al. (2018a), find the JLA data prefers isotropy,
with the results being inconclusive on the Union2.1 data. Javanmardi et al.
(2015) however, find that the null hypothesis of isotropy cannot be rejected
unless ones specifically takes into account its alignment with the dipole of
the Cosmic Microwave Background (CMB), in which case the null hypothesis can
be rejected at a level ranging from 95% to 99% confidence. An important
distinction between these data sets is that Union2/2.1 and Constitution have
no corrections for the peculiar velocities of the SNe host galaxies, whereas
JLA does. Bernal et al. (2017) found from Union2 and LOSS data potential
differences between hemispheres in the isotropy of the deceleration parameter.

The situation becomes more confused when considering the largest SNIa
compilation to date, Pantheon, encompassing 1048 objects in the redshift range
$0.01<z<2.3$ (Scolnic et al., 2018b). A major hurdle to any re-analysis that
uses Pantheon is the lack of a publicly available full correlation matrix for
its SNIa light curve standardisation coefficients, which hampers a principled
statistical approach. Nevertheless, several papers have attempted to use
Pantheon-derived measurements of the distance modulus as a function of
redshift to investigate potential deviations from an isotropic expansion,
finding isotropy is still favoured (Sun & Wang, 2018b; Soltis et al., 2019;
Zhao et al., 2019; Andrade et al., 2018b).

Despite the existence of substantially larger compilations like Pantheon, the
JLA remains a useful data set for analyses of this kind, because all the
necessary statistical and systematics covariance matrixes are publicly
available, unlike e.g. the Pantheon set. Recently, Colin et al. (2019b)
(henceforth C19) claimed 3.9-$\sigma$ evidence for a dipole in the
deceleration parameter from a maximum likelihood analysis of JLA data, leading
to a lack of statistical evidence for acceleration in the expansion. The claim
was disputed by Rubin & Heitlauf (2020) (henceforth, RH20), who pointed out
the incorrect use of heliocentric redshifts in C19 and other technical
assumptions about selection effects which, when corrected, remove the
preference for a dipole and restore the high significance for an accelerated
expansion. (See also the discussion in Rubin & Hayden (2016), itself a
rebuttal of Nielsen et al. (2016).) A further reply by Colin et al. (2019a)
appears to concede some technical points, but not the overall conclusion on
the actual lack of statistical significance for an accelerated expansion.

The aim of this work is to clarify the status of claims for a statistically
significant dipole in the accelerated expansion of the universe, especially in
light of the ongoing controversy. In so doing, we also revisit the important
question of the level of statistical evidence in favour of the accelerated
expansion in an isotropic universe from SNIa data alone. We address the
criticisms of the published JLA data made by C19, who claimed that the
peculiar velocity corrections made to the JLA SNIa data based on local bulk
flows are incorrect. In this paper, we introduce a state-of-the-art treatment
of peculiar velocities, which are independently constrained using the 2M++
galaxy catalogue (Lavaux & Hudson, 2011; Boruah et al., 2020a, b), re-derive
correlated peculiar velocity uncertainties (both statistical and systematic)
from a fully consistent flow model, and upgrade the Bayesian hierarchical
model BAHAMAS Shariff et al. (2016) to include an new treatment of residual
colour-based selection effects in SNIa data.

The remaining of this paper is structured as follows: section 2 introduces the
cosmological model, the anisotropy model, our Bayesian framework, the data
used, our new peculiar velocities treatment and our new colour-based selection
effects correction. Section 3 demonstrates the performance of our method on
simulated data. Our results from the JLA data, both in terms of parameter
inference and Bayesian model comparison, are presented in section 4. Our
conclusions are given in section 5.

##  2 Methodology and Data

###  2.1 Cosmological model and dipole modulation

We investigate the isotropy of the expansion in both a model-specific and a
model-independent way: firstly, we consider the $\Lambda$CDM model for the
underlying cosmology; secondly, we use the so-called “Cosmographic expansion”
(i.e., a Taylor expansion of the scale factor as a function in time) as a
model-independent description of the underlying matter-energy density of the
universe.

The $\Lambda$CDM model has cosmological parameters
${\mathscr{P}}_{1}=\\{\Omega_{m},\Omega_{\Lambda},H_{0}\\}$, where
$\Omega_{m}$ and $\Omega_{\Lambda}$ are the density parameters of matter (both
baryonic and dark) and cosmological constant, $\Lambda$, respectively, in
units of the critical energy density; $H_{0}$ is the Hubble-Lemaître constant
(which we fix to $H_{0}=72$ km/s/Mpc, as it is exactly degenerate with the
SNIa intrinsic magnitude). The curvature parameter $\Omega_{\kappa}$ is given
by

$\Omega_{\kappa}=1-\Omega_{m}-\Omega_{\Lambda}$ (1)

and we assume a universe with constant dark energy equation of state,
$w(z)=-1$. We denote by $\bar{z}$ the redshift of a comoving galaxy seen by an
observer who is also at rest w.r.t. to the CMB restframe (i.e., the
“cosmological” redshift, with no peculiar velocities from either the source or
the observer) and by $z_{\rm hel}$ the redshift for an observer in the Sun’s
frame of reference111We neglect the distinction between geocentric and
heliocentric frames of references (the difference due to the $\sim 30$ km/s
orbital speed of the Earth is of order $\Delta z\sim 10^{-5}$), since redshift
measurements are routinely reported in the heliocentric frame and also already
corrected for atmospheric refraction.. The measured redshift in our
heliocentric frame of reference is given by $\hat{z}_{\rm hel}$, and it
differs from $z_{\rm hel}$ by measurement noise. The redshift in our
heliocentric frame of reference, $z_{\rm hel}$, differs from the redshift of a
comoving observer, $\bar{z}$, by virtue of peculiar velocities of the source
and the observer, and gravitational red/blueshifts due to the local
gravitational potential at the location of the source and observer. In the
following, we neglect gravitational effects, which are subdominant (see
however Calcino & Davis (2017)) and focus instead on the impact of peculiar
velocities.

The relationship between heliocentric redshift, $z_{\rm hel}$, and the
redshift of a comoving galaxy as seen by an observer at rest w.r.t. the CMB,
$\bar{z}$ is given by:

$\displaystyle(1+z_{\rm hel})$ $\displaystyle=(1+z_{\rm\tiny{CMB}})(1+z_{\rm
pec}^{\odot})$ (2) $\displaystyle(1+z_{\rm\tiny{CMB}})$
$\displaystyle=(1+\bar{z})(1+z_{\rm pec}^{\rm SN})$ (3)

where $z_{\rm pec}^{\odot}$ is the redshift induced by the peculiar velocity
of the Solar System w.r.t. the CMB restframe, while $z_{\rm pec}^{\rm SN}$ is
the redshift caused by the peculiar velocity of the SNIa w.r.t. the CMB frame.
The second equality introduces the redshift in the CMB restframe,
$z_{\rm\tiny{CMB}}$, i.e., the frame in which our motion w.r.t. the CMB has
been removed222A source of confusion in the literature is the widespread use
of the term “CMB restframe” to denote what we call $\bar{z}$ (i.e., the
cosmological redshift, with no peculiar motions from either source nor
observer). This misleading nomenclature is for example used by Betoule et al.
(2014), as well as in the data products of the JLA data release.. With the
above definitions, we can write the luminosity distance to redshift $\bar{z}$,
as (Davis et al., 2011)

$\begin{split}d_{L}(\bar{z},z_{\rm pec}^{\odot},z_{\rm pec}^{\rm
SN},{\mathscr{P}}_{1})&=\frac{c}{H_{0}}\frac{(1+\bar{z})(1+z_{\rm
pec}^{\odot})(1+z_{\rm pec}^{\rm SN})^{2}}{\sqrt{|\Omega_{\kappa}|}}\times\\\
&\mathrm{sinn}\left\\{\sqrt{|\Omega_{\kappa}|}\int_{0}^{\bar{z}}\frac{dz}{E(z)}\right\\}.\end{split}$
(4)

$E(z)$ depends on our choice of cosmology and for the $\Lambda$CDM universe,
is given by

$E^{2}(z)=\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}+\Omega_{\kappa}(1+z)^{2}.$ (5)

The sinn$(x)$ function is defined as

$\text{sinn}(x)=\begin{cases}x&\quad\text{if}\quad\Omega_{\kappa}=0\\\
\text{sin}(x)&\quad\text{if}\quad\Omega_{\kappa}<0\\\
\text{sinh}(x)&\quad\text{if}\quad\Omega_{\kappa}>0\\\ \end{cases}.$ (6)

In our model-independent approach, we follow Visser (2004) and Taylor-expand
the scale factor of the FLRW metric up to third order in time around $t_{0}$
(today), as:

$a(t)=a_{0}\\{1+H_{0}(t-t_{0})-\frac{1}{2}q_{0}H_{0}^{2}(t-t_{0})^{2}+\frac{1}{3!}j_{0}H_{0}^{3}(t-t_{0})^{3}+O([t-t_{0}]^{4})\\}$
(7)

where $q_{0}$ is the dimensionless deceleration parameter, defined as

$q_{0}=-\frac{1}{a}\frac{\mathrm{d}^{2}a}{\mathrm{d}t^{2}}\left[\frac{1}{a}\frac{\mathrm{d}a}{\mathrm{d}t}\right]^{-2}_{t=t_{0}}$
(8)

and $j_{0}$ is the so-called “jerk”,

$j_{0}=+\frac{1}{a}\frac{\mathrm{d}^{3}a}{\mathrm{d}t^{3}}\left[\frac{1}{a}\frac{\mathrm{d}a}{\mathrm{d}t}\right]^{-3}_{t=t_{0}},$
(9)

which is also dimensionless. This model-independent expansion only relies on
the FLRW metric but makes no assumption about the underlying matter-energy
density, and leads to the following form of the luminosity distance at
redshift $\bar{z}$:

$\begin{split}&d_{L}(\bar{z},z_{\rm pec}^{\odot},z_{\rm pec}^{\rm
SN},{\mathscr{P}}_{2})=\frac{(1+z_{\rm pec}^{\odot})(1+z_{\rm pec}^{\rm
SN})^{2}}{1+\bar{z}}\frac{c\bar{z}}{H_{0}}\times\\\
&\qquad\left[1+\frac{1}{2}(1-q_{0})\bar{z}-\frac{1}{6}(1-q_{0}-3q_{0}^{2}+j_{0}-\Omega_{\kappa})\bar{z}^{2}+O(\bar{z}^{3})\right],\end{split}$
(10)

where $c$ is the speed of light and the model-independent parameters are
${\mathscr{P}}_{2}=\\{H_{0},q_{0},j_{0},\Omega_{\kappa}\\}$. From expressions
(4) or (10), we obtain the isotropic distance modulus, $\mu_{I}$, using the
standard formula

$\mu_{I}(\bar{z},z_{\rm pec}^{\odot},z_{\rm pec}^{\rm
SN},{\mathscr{P}}_{a})=25+5\log_{10}\frac{d_{L}(\bar{z},z_{\rm
pec}^{\odot},z_{\rm pec}^{\rm SN},{\mathscr{P}}_{a})}{1\text{ Mpc}},$ (11)

where $a=1,2$ depending on the chosen parameterization.

There are several different ways one can parameterize the possibility of
anisotropic expansion, depending on the underlying physical origin for the
effect. A spherical harmonics expansion introduces, to lowest order, a dipolar
modulation in the direction of a SNIa situated at redshift $\bar{z}$ and in
direction ${\mathbf{n}}_{\rm\tiny{SN}}$ in the sky, with
${\mathbf{n}}_{\rm\tiny{SN}}$ a unit vector pointing from the centre of the
coordinate system (the Earth) to the location of the SNIa on the celestial
sphere. Different authors have taken different approaches in the literature,
with no consensus as to which quantity should be modulated: one could expand
the scale factor $a(t)$, the luminosity distance, the comoving distance, the
Hubble parameter, the matter density, the cosmological constant density, or
the distance modulus. Each of these possibilities leads to a different
anisotropic imprint onto the Hubble-Lemaître law. A dipole moment that is
constant with distance, $r$, in the peculiar velocity field (i.e., a bulk
flow) leads to

$c\bar{z}\approx H_{0}r+D_{v}({\mathbf{n}}_{\rm
dip}\cdot{\mathbf{n}}_{\rm\tiny{SN}}),$ (12)

By contrast, a constant dipole in either $H_{0}$ or $r$ leads to

$c\bar{z}\approx H_{0}r+H_{0}D_{H}r({\mathbf{n}}_{\rm
dip}\cdot{\mathbf{n}}_{\rm\tiny{SN}}),$ (13)

which increases linearly with distance. Another possibility is to modulate the
distance modulus directly:

$\mu=\mu_{I}(\bar{z},z_{\rm
hel},{\mathscr{P}}_{a})\left(1+D_{\mu}F(\bar{z})({\mathbf{n}}_{\rm
dip}\cdot{\mathbf{n}}_{\rm\tiny{SN}})\right),$ (14)

where $F(z)$ is a function of redshift alone which can be used to localize the
dipole at a given length scale. Yet another approach, adopted by C19, is to
model the dipole on the deceleration parameter, $q_{0}$, in a Cosmographic
expansion:

$q_{0}(z)=q_{m}+D_{q_{0}}F(\bar{z})({\mathbf{n}}_{\rm
dip}\cdot{\mathbf{n}}_{\rm\tiny{SN}}).$ (15)

In this paper, we add the dipole to either the distance modulus, Eq. (14), or
to the deceleration parameter, Eq. (15), and consider both $F(z)=1$ and,
following C19, an exponentially decaying function of redshift with
characteristic scale given by the free parameter $S$, namely
$F(z)=\exp{(-z/S)}$. These two forms have the effect of either creating a
dipole that is constant in redshift or constrained to a local scale which
could arise for reasons such as existing within a cosmic void.

As noted in previous works, a phenomenological approach as the one taken here
that perturbs an underlying FLRW metric may not be entirely consistent. An
alternative route would require specifying a physical model for the
anisotropy, and then derive the ensuing predictions for the distance modulus
and compare those with observations, as done for example in the context of an
ellipsoidal universe from Bianchi type I models (Campanelli et al., 2011).
However, the advantage of a purely phenomenological approach is that it
remains agnostic about the underlying cause of any anisotropy, and it provides
constraints on the level of anisotropy that can then be applied to other
models.

###  2.2 Bayesian Hierarchical Model

In this paper, we improve on previous works constraining anisotropy from SNe
data by adopting a fully Bayesian hierarchical model (BHM) for the statistical
analysis of SNIa data, called BAHAMAS. “Hierarchical” here refers to the
existence, within the model, of a layer of unobserved (so-called “latent”)
variables for each SNIa, corresponding to the true value of their light-curve-
derived properties (as opposed to the noisy measured value). The latent
variables are marginalised over in the inference, and are constrained in
virtue of the fact that they are all generated from the same underlying
population distribution, which is modelled with a set of hyperparameters,
themselves determined from the data. This approach has been shown to have
better coverage statistics overall than the $\chi^{2}$ method traditionally
employed, and leads to a reduction of mean squared errors for the recovered
cosmological parameters by a factor of $\sim$ 2 - 3 when deployed on simulated
data (March et al., 2011; Shariff et al., 2016). Furthermore, the BAHAMAS
framework can be used for principled Bayesian model comparison of the kind we
perform in section 4.2, as it correctly marginalizes out all nuisance
parameters. By contrast, the heuristic $\chi^{2}$ approach traditionally
adopted is an approximation to the BAHAMAS likelihood, and cannot be used to
compute Bayes factors or for formal Bayesian model comparison.

The Bayesian methodology pioneered in March et al. (2011) has been adopted and
extended in several other papers, including e.g., UNITY (Rubin et al., 2015),
Steve Hinton et al. (2018) and Simple-BayeSN (Mandel et al., 2017). (See also
Nielsen et al. (2016) for a profile likelihood interpretation.) Here we
briefly summarize our Bayesian Hierarchical model, which builds on BAHAMAS,
referring the reader to (March et al., 2011; Shariff et al., 2016) for fuller
details.

We denote with a hat symbol observed quantities, in order to distinguish them
from the latent (i.e., unobserved) variables in our model. For each SNIa $i$,
the data $d_{i}$ can be summarised by a vector

$d_{i}=\\{\hat{z}_{{\rm
hel},i},\hat{c}_{i},\hat{x}_{1_{i}},\hat{m}_{B_{i}},\hat{\Sigma}_{C,i}\\},$
(16)

where $\hat{z}_{{\rm hel},i}$ is the observed heliocentric redshift,
$\hat{m}_{B_{i}}$ is the observed peak B-band apparent magnitude,
$\hat{x}_{1_{i}}$ and $\hat{c}_{i}$ are observed ‘stretch’ and ‘colour’
corrections, which are summary statistics of the lightcurve of the SNIa
obtained with the lightcurve fitter SALT2 (Guy et al., 2005, 2007) during the
standardization procedure. Furthermore,
$\hat{\Sigma}_{C,i}=\text{Cov}(\hat{c}_{i},\hat{x}_{1_{i}},\hat{m}_{B_{i}})$
is a $3\times 3$ variance-covariance matrix that describes the measurement
error on the observables. On the standard deviation scale, the measurement
error for redshift for SNIa with spectroscopic follow-up is $\sigma_{z}^{\rm
sp_{SN}}\sim 5\times 10^{-3}$ when the redshift is determined from the SNIa
spectrum alone, and $\sigma_{z}^{\rm sp_{host}}\sim 5\times 10^{-4}$ when it
is obtained from host-galaxy spectra (Zheng et al., 2008). The redshift
measurements are independent from each other and from all other observables.
We discuss this uncertainty further in section 2.6. (See also Calcino & Davis
(2017) for the potentially important impact of systematic redshift errors as
small as $\Delta z\sim 10^{-4}$.)

In BAHAMAS, we introduce latent variables for each SNIa in order to model each
source of uncertainty according to its origin: measurement error, population
scatter and intrinsic (residual) variability. A probabilistic hierarchical
model is built as follows: each SNIa has latent variables $M_{i}^{\epsilon}$,
$x_{1_{i}}$ and $c_{i}$, representing the objects’ “true” (i.e., noiseless)
absolute magnitude, stretch correction and colour correction, respectively.
These latent variables follow normal distributions, representing population
variability of the SNe and parameterised by their means and variances:

$\displaystyle x_{1_{i}}$
$\displaystyle\sim\mathcal{N}\left(x_{1\star},R_{x_{1}}^{2}\right),$ (17)
$\displaystyle c_{i}$
$\displaystyle\sim\mathcal{N}\left(c_{\star},R_{c}^{2}\right),$ (18)
$\displaystyle M_{i}^{\epsilon}$
$\displaystyle\sim\mathcal{N}\left(M_{0}^{\epsilon},\sigma_{\rm
res}^{2}\right),$ (19)

where $x_{1\star},c_{\star}$ and $M_{0}^{\epsilon}$ are the population means
and $R_{x_{1}}^{2},R_{c}^{2}$ and $\sigma_{\rm res}^{2}$ are the population
variances, all of which are also estimated from the data. We collect the
population-level parameters in a vector of variables
$\vartheta\equiv\\{x_{1\star},c_{\star},M_{0}^{\epsilon},R_{x_{1}}^{2},R_{c}^{2},\sigma_{\rm
res}^{2}\\}$. The intrinsic magnitude of each SNIa $M_{i}$, is modified by
applying the linear ‘Tripp relation’ (Tripp, 1998), so that $M_{i}\rightarrow
M_{i}^{\epsilon}\equiv M_{i}+\alpha{x}_{1_{i}}-\beta{c}_{i}$, where $\alpha$
and $\beta$ are nuisance parameters that control the slope of the stretch and
colour correction, respectively. Therefore, $M_{i}^{\epsilon}$ is a linear
function of $M_{i}$ that features a lower population variance, represented by
$\sigma_{\rm res}^{2}$. In astrophysical parlance, $M_{i}^{\epsilon}$ is
referred to as “the corrected intrinsic magnitude” of the SNIa. Thanks to the
standardization procedure, the residual standard deviation of the SNe type
Ia’s corrected intrinsic magnitude can be sufficiently reduced so that they
can be used as luminosity distance indicators. Shariff et al. (2016) used the
BAHAMAS model to determine the residual standard deviation of the corrected
intrinsic magnitude, finding (from JLA data) a value $\sigma_{\rm res}\sim
0.104\pm 0.005$, similar to (if somewhat smaller than) the value obtained with
a $\chi^{2}$ method by Betoule et al. (2014). At the latent level, the
apparent peak magnitude $m_{Bi}$ is related to the standardised intrinsic
magnitude $M_{i}^{\epsilon}$ via the isotropic distance modulus of Eq. (11):

$m_{Bi}=\mu_{I}(z_{{\rm
hel},i},{\bar{z}},_{i},{\mathscr{P}}_{a})-\alpha{x}_{1_{i}}+\beta{c}_{i}+M_{i}^{\epsilon}.$
(20)

Finally, the observed values of
$\\{\hat{m}_{B_{i}},\hat{x}_{1_{i}},\hat{c}_{i}\\}$ are modelled as normally
distributed around their latent values, with variance-covariance matrix given
by $\hat{\Sigma}_{C,i}$. Additionally, a systematic errors covariance matrix
that correlates different SNIe (for example, because of calibration
uncertainties common between SNIe within the same survey) is included when
available. Inference is based on the marginal distribution of the quantities
of interest, ${\mathscr{P}}_{a}$, which includes uncertainty at all levels of
the hierarchy.

It has become common practice to split the SNe into two groups, based on their
host-galaxy stellar mass (Sullivan et al., 2010), with the mass threshold
between the two groups being around value of $\log_{10}(M_{g})=10$, where
$M_{g}$ is the host-galaxy mass measured in solar masses. Smith et al. (2020)
found a difference of up to 0.04 mag in the average intrinsic magnitude of the
two groups, and application of the Tripp relation to the two groups separately
may further reduce the residual dispersion in SNIa absolute magnitudes. (See
also Thorp et al. (2021) for a similar result obtained using a Bayesian
hierarchical approach.) However, Brout & Scolnic (2021) have recently cast
doubt on the robustness of this mass-step correction, which they ascribed
instead to incorrect dust modeling. The ultimate origin of the mass step
remains unclear, and it might relate to stellar population age (Childress et
al., 2014) and metallicity (Sullivan et al., 2010). In any case, Shariff et
al. (2016) showed that adding a mass-step or a more general linear covariate
as a function of host-galaxy mass has little impact on the ensuing
cosmological parameters inference. Therefore, in this paper we do not adopt a
mass-step correction. We refer the reader to March et al. (2011) and Shariff
et al. (2016) for the full mathematical detail of the BAHAMAS model,
marginalization procedure and algorithms used for sampling the resulting
posterior distribution.

###  2.3 Data

The largest SNIa compilation to date is the ‘Pantheon sample’ (Scolnic et al.,
2018a), which contains 1048 spectroscopically confirmed SNIa. This compilation
includes 279 new SNIa discovered by the Pan-STARRS1 (PS1) Medium Deep Survey
(Chambers et al., 2016) in addition to the previous SNe discovered by previous
catalogues to create the total.

Unfortunately, the Pantheon sample only provides estimates (and associated
uncertainties) of the distance modulus for each SNIa but does not include the
covariance matrices of both the measurement error for each SNIa and the
systematic covariance matrix across the whole data set. Because our Bayesian
hierarchical model also requires these covariances over the light curve fit
parameters as opposed to the covariance over distance modulus provided by the
Pantheon data, we instead use the smaller ‘Joint Light-Curve Analysis’ (JLA)
compilation (Betoule et al., 2014). Very recently, a ‘Pantheon+’ dataset has
been presented in Brout et al. (2022), increasing the set of spectroscopically
confirmed SNe to 1550 in the redshift span $0.001$ to $2.26$. The data
products of this larger set have not yet been released, so unfortunately we
cannot use these data in our framework. The JLA data contains 740 SNIa
including 374 SNIa from the SDSS-II survey (Frieman et al., 2008; Sako et al.,
2018), 239 from SNLS (Conley et al., 2011; Sullivan et al., 2011), a low-$z$
sample of 118 SNIa at $z<0.1$ which is comprised of numerous smaller surveys
and nine Hubble Space Telescope SNIa. These have been fit and standardised
using the SALT2 Light-Curve Fitter (Guy et al., 2007).

An overview of the distribution of JLA objects in the sky is provided in Fig.
1. The long stripe in the lower left hemisphere is from the SDSS objects. It
is clear to see that the distribution of the JLA objects in the sky is highly
anisotropic.

Figure 1: A plot of the JLA SNe showing their location in the sky in galactic
coordinates, as well as the direction of the CMB dipole (left-pointing
triangle) and bulk flow (right-pointing triangle), as determined by Boruah et
al. (2020a).

###  2.4 Accounting for Colour-dependent Selection Effects

Towards the high end of the redshift range of a given survey, SNe that are
intrinsically brighter (smaller $M_{i}$) or bluer (smaller $c_{i}$) are more
likely to be observed and to be followed up spectroscopically to confirm their
type. This selection bias must be accounted for to avoid bias in the estimates
of the cosmological parameters. Ignoring magnitude-based selection effects
leads to estimates of the distance modulus that are biased low: at a given
redshift, the average observed magnitude is smaller (i.e., the observed peak
flux brighter) than the population mean, which leads to an estimated distance
modulus that is biased low. This effect reduces or even obliterates the
preference for a non-zero cosmological constant. Traditionally, this has been
addressed by “correcting back” the estimates of the distance modulus by the
average bias in each redshift bin, established with forward simulations of
data subject to selection effects. This approach is adopted e.g., by Betoule
et al. (2014). More recently, this method has been extended and refined with
the so-called “BEAMS with Bias Corrections” (BBC) method (Kessler & Scolnic,
2017). The JLA analysis only corrected for magnitude-based selection effects,
concluding that no additional correction were necessary for colour (Betoule et
al. (2014), Fig. 11), despite observing a downward trend in observed colour
with redshift for SNSL and SDSS.

In the context of a Bayesian analysis, however, selection effects are treated
differently: the posterior is conditional on the observed data (Kelly, 2007),
which leads to a re-weighing factor increasing the statistical weight of SNe
that are less likely to be observed (see Eq. (23) below). Rubin et al. (2015)
introduced a general formalism for the selection function that was further
developed by Hinton et al. (2018). In practice, however, this formalism
typically requires several simplifying assumptions that may be difficult to
justify (e.g., a well-sampled SNIa redshift distribution, a selection function
that is described by a normal cumulative distribution function (CDF),
independence of the selection probability from the underlying cosmology). We
will present an improved selection effects treatment in an upcoming, dedicated
work.

Rubin & Hayden (2016) argue that uncorrected-for colour-dependent selection
effects remaining in the JLA data (after bias correction of the data) should
be addressed by introducing a population colour mean that is both redshift-
and survey-dependent. However, at the top level of the BHM we would like a
_physically-meaningful_ population mean that describes the underlying,
population-level latent mean colour, itself a reflection of the physical
properties of the SNe. Such a colour mean can be a function of redshift, to
reflect evolution in the physical properties of SNe with lookback time, but it
cannot be survey-dependent, for clearly the latter dependency is caused by
survey-specific selection effects and, in virtue of being survey-specific,
cannot be the consequence of changing underlying physical properties of the
SNe being observed. Thus, by using population-level variables to address a
survey-induced selection effect, the method advocated by Rubin & Hayden (2016)
goes against the physical interpretability of the BHM; we prefer the
population-level variables to be tied to the physics of the SNIa explosion
mechanism rather than the survey-induced selection effects. Furthermore, as
pointed out by Dam et al. (2017), modeling residual colour drift with redshift
as advocated by Rubin & Hayden (2016) introduces undesirable degeneracies with
cosmological parameters in the cosmographic expansion.

Here, we account for residual colour-based selection effects with an
approximate method that captures the spirit of the correct Bayesian procedure.
(We will present a more complete exposition in an upcoming work.) In general,
we aim to base the likelihood function on the distribution of the data $d_{i}$
for SNIa $i$ conditional on it having been observed, which is denoted by an
indicator variable $I_{i}=1$. ($I_{i}=0$ would indicate that SNIa $i$ is not
observed.) Denote by $\hat{D}=\\{\hat{d}_{1},\dots,\hat{d}_{n}\\}$ a random
sample of SNe where all SNe have equal probability of being observed (i.e.,
there are no selection effects in $\hat{D}$). In the presence of selection
effects, we observe a non-representative sample of $n_{\rm obs}<n$ SNe from
$\hat{D}$, whereby a SNIa $i$ is observed with probability given by the
selection function $p(I_{i}=1|\hat{d}_{i},\Psi)$, where $\hat{d}_{i}$ are the
observed (noisy) data, and $\Psi$ are parameters describing the selection
function (assumed known). The distribution of the observed $\hat{d}_{i}$,
conditional on SNIa $i$ being observed, is given by

$p(\hat{d}_{i}|I_{i}=1,\Psi,\Theta)=\frac{p(I_{i}=1|\hat{d}_{i},\Psi)p(\hat{d}_{i}|\Theta)}{p(I_{i}=1|\Psi,\Theta)},$
(21)

where $\Theta=\\{{\mathscr{P}}_{a},\vartheta\\}$ are the parameters of the
hierarchical model (including the parameterization of the distance modulus,
${\mathscr{P}}_{a}$, and the population-level distribution parameters,
$\vartheta$). We assume that the selection probability conditional on the
observed data, $\hat{d}_{i}$, i.e., $p(I_{i}=1|\hat{d}_{i},\Psi)$ in the
numerator of (21), does not depend on $\Theta$ and thus is an ignorable
constant if (as we assume here) the selection function and its parameter
$\Psi$ are known. The quantity $p(\hat{d}_{i}|\Theta)$ is the likelihood in
the absence of selection effects, and the denominator gives the probability of
observing a SNIa, irrespective of the value of the data:

$p(I_{i}=1|\Psi,\Theta)=\int{\rm
d}\hat{d}_{i}p(I_{i}=1|\hat{d}_{i},\Psi)p(\hat{d}_{i}|\Theta).$ (22)

Omitting the ignorable multiplicative constants in the numerator of (21)
yields the likelihood function of $\Theta$ including selection effects for a
sample of $n_{\rm obs}$ observed SNe, $\hat{d}_{\rm
obs}\equiv\\{\hat{d}_{1},\dots,\hat{d}_{n}{}_{\rm obs}\\}$,

$p(\hat{d}_{\rm obs}|\\{I_{i}=1\\}_{i=1}^{n_{\rm
obs}},\Psi,\Theta)\propto\frac{p(\hat{d}_{\rm
obs}|\Theta)}{p(I=1|\Psi,\Theta)^{n}_{\rm obs}},$ (23)

where we have dropped the dummy index $i$ in the denominator. Therefore,
selection effects are accounted for by dividing the likelihood function of the
observed data, $p(\hat{d}_{\rm obs}|\Theta)$, by a “correction factor” that
gives the probability of making $n_{\rm obs}$ observations.

Thus far, we have been entirely general. Next we specify the form of the
selection function, $p(I_{i}=1|\hat{d}_{i},\Psi)$, entering in Eq. (22). The
probability of a SNIa being selected, and spectroscopically followed-up to
determine its type, depends primarily on its magnitude and colour. (SNe with
larger stretch parameter $x_{1}$ are slower declining and thus remain visible
and potentially detectable for longer, but this effect is subdominant.) The
data correction procedure in Betoule et al. (2014) should in principle account
for both magnitude- and colour-based selection, but their discussion makes it
clear that there are large uncertainties in the determination of the selection
probability that enters their forward simulation of the data. For example,
Betoule et al. (2014) mention that the SDSS spectroscopic follow-up target
selection favours intrinsically bluer SNe, introducing complex colour-
dependency in the selection function. In light of such difficult-to-simulate
selection effects, we advocate a method that estimates any residual selection
effect (after the data correction procedure of Betoule et al. (2014)) directly
from the observed data.

We wish to account for residual colour-based selection effects that may remain
in the data. Therefore, we assume that $p(I_{i}=1|\hat{d}_{i},\Psi)$ depends
only on $\hat{c}_{i}$ and $\hat{z}_{{\rm hel},i}$, with the $N_{s}$ redshift
bins for survey $s$ chosen as discussed in section 2.3, and factorize both the
selection function and the likelihood in a product over redshift and survey
bins, assumed independent of each other. Within each redshift, and survey bin,
we allow a different selection function, which is derived below. With these
assumptions, the probability of observing $n_{\rm obs}$ SNe is (with the
shorthand notation $I_{n}{}_{\rm obs}=1$ denoting
$\\{I_{1}=1,\dots,I_{n}{}_{\rm obs}=1\\}$):

$p(I_{n}{}_{\rm
obs}=1|\Psi,\Theta)=\prod_{s=1}^{4}\prod_{j=1}^{N_{s}}\left(\int{\rm
d}\hat{d}_{i}p_{sj}(I_{i}=1|\hat{d}_{i},\Psi)p_{sj}(\hat{d}_{i}|\Theta)\right)^{N_{sj}},$
(24)

where $N_{s}$ is the number of bins for survey $s$ and $N_{sj}$ the number of
observed SNe in bin $sj$. Within each redshift bin for survey $s$, we
parameterize the selection function as a normal cumulative distribution
function (CDF), and assume that we observe a SNIa $i$ with colour
$\hat{c}_{i}$ with probability:

$p_{sj}(I_{i}=1|\hat{d}_{i},\Psi)=\Phi\left(\frac{c^{\rm
obs}_{sj}-\hat{c}_{i}}{\sigma^{\rm obs}_{sj}}\right),$ (25)

where

$\Phi(x)=\int_{-\infty}^{x}{\mathcal{N}}_{y}(0,1)dy$ (26)

is the CDF of a standard normal, and ${\mathcal{N}}_{y}(\mu,\sigma^{2})$ is a
Gaussian distribution in $y$ with mean $\mu$ and variance $\sigma^{2}$. In Eq.
(25), $c^{\rm obs}_{sj}$ is the colour value at which there is a 50%
probability of observing a SNIa in redshift bin $j$ and for survey $s$;
$\sigma^{\rm obs}_{sj}$ denotes the width of the transition from the regime
where all objects are observed, i.e. for $(c^{\rm
obs}_{sj}-\hat{c}_{i})/\sigma^{\rm obs}_{sj}\ll 0$, to the regime where no
objects are observed, where $(c^{\rm obs}_{sj}-\hat{c}_{i})/\sigma^{\rm
obs}_{sj}\gg 0$, for the bin $sj$ being considered. In Eqs. (31)-(32) below we
show how to estimate the selection function parameters $\Psi=\\{(c^{\rm
obs}_{sj},\sigma^{\rm obs}_{sj})\\}$ ($s=1,\dots,4,j=1,\dots.N_{s}$). As an
approximation, we ignore uncertainty in the resulting estimates and assume the
parameters are known exactly333We ignore uncertainty in our estimates of
$\Psi$, and we defer the evaluation of the impact of this approximation to an
upcoming, dedicated work..

In principle, we would like to use the likelihood function of March et al.
(2011) and Shariff et al. (2016) as the second term of the integrand in Eq.
(22). Since we only wish to account for residual colour-based selection
effects, however, we ignore the part of the BAHAMAS likelihood that relates
colour to magnitude via the Tripp linear relation (Eq. (C2) in March et al.
(2011)), and instead only consider the distribution of colour values that one
would obtain when integrating out the latent colour variables conditional on
all other variables in the BHM, leading to the simple expression for the
likelihood entering into Eq. (24):

$p_{sj}(\hat{d}_{i}|\Theta)={\mathcal{N}}_{\hat{c}}{}_{i}(c_{\star},R_{c}^{2}+\bar{\sigma}_{c,sj}^{2}),$
(27)

where $\bar{\sigma}_{c,sj}$ is the average colour measurement error for the
$n_{sj}$ SNe in bin $sj$ (for simplicity, we assume all $n_{sj}$ SNe in bin
$sj$ have the same colour measurement error, given by $\bar{\sigma}_{c}$; we
also ignore correlation between colour and stretch and magnitude)

Eq. (27) features a redshift-independent conditional expectation of colour,
described by $c_{\star}$. This formalizes the assumption of Rubin & Hayden
(2016) within BAHAMAS that the observed drift to bluer SNe with redshift
within a survey is a consequence of selection effects and not of a change in
the underlying population colour distribution with redshift.

With the above elements, we can compute the probability of observing SNIa $i$
in redshift bin $j$ for survey $s$ by integrating over its colour,
$\hat{c}_{i}$, obtaining 444It is useful to recall that
$\int_{-\infty}^{\infty}\Phi\left(\frac{\mu-x}{\sigma}\right){\mathcal{N}}_{x}(\nu,\tau^{2})dx=\Phi\left(\frac{\mu-\nu}{\sqrt{\sigma^{2}+\tau^{2}}}\right)$.

$p_{sj}(I_{i}=1|\Psi,\Theta)=\Phi\left(\frac{c^{\rm
obs}_{sj}-c_{\star}}{\sqrt{(\sigma^{\rm
obs}_{sj})^{2}+R_{c}^{2}+\bar{\sigma}_{c,sj}^{2}}}\right).$ (28)

For a given $sj$ bin, when $c_{\star}\ll c^{\rm obs}_{sj}$ selection effects
are irrelevant, because the survey is seeing the entire colour population, and
accordingly $p(I_{i}=1|\Psi,\Theta)\rightarrow 1$ from Eq. (28). However,
$p(I_{i}=1|\Psi,\Theta)$ becomes smaller for values of $c_{\star}>c^{\rm
obs}_{sj}$, with the difference measured in units of the total standard
deviation, i.e., when the survey is preferentially seeing the bluer part of
the population because of colour-based selection bias. In this case, the
observed distribution of the $N_{sj}$ objects in the $sj$ bin deviates from
the latent distribution, and the likelihood that ignores selection effects
would incorrectly penalize this value of $c_{\star}$. According to Eq. (23),
the correction factor in the denominator of Eq. (23) increases the weight
given to observed SNe with $c_{\star}>c^{\rm obs}_{sj}$.

Finally, there remains the issue of determining the value of $\Psi$, the
selection function parameters. Ideally, one would do so from forward
simulation of surveys, but this is unpractical for our purposes, and
unfeasible for the low-$z$ sample, which is obtained from a collection of
telescopes with poorly understood selection functions. Furthermore, as noted
earlier, colour-based selection effects might escape ab initio modeling of
this kind, as indicated in e.g. Betoule et al. (2014). As an alternative, we
estimate the value of $\Psi$ in each redshift- and survey-bin $sj$, by
matching the first and second moment of the empirical colour distribution
within the bin to the marginal distribution of colour based on the right-hand-
side of (23), understood as a distribution over observed colour values
$\hat{c}_{i}$ within each bin. Let $m_{C}(t)$ denote the moment generating
function of the random variable $C$ (the observed colour), defined as

$m_{C}(t)\equiv\int_{-\infty}^{\infty}e^{ty}f_{C}(y)dy,$ (29)

where $f_{C}(c)$ in our case is given by (in the bin under consideration)

$f_{C}(c)=\frac{{\mathcal{N}}_{c}(c_{\star},R_{c}^{2}+\bar{\sigma}_{c,sj}^{2})\Phi\left(\frac{c^{\rm
obs}_{sj}-c}{\sigma^{\rm obs}_{sj}}\right)}{p_{sj}(I_{i}=1|\Psi,\Theta)},$
(30)

with denominator given by Eq. (28). Our strategy, known as the _method of
moments_ , is to analytically compute the first two moments of the
distribution, set them equal to the empirical moments, and solve the resulting
system of equations to obtain estimates of the selection function parameters
in each survey and redshift bin, $\\{c^{\rm obs}_{sj},\sigma^{\rm
obs}_{sj}\\}$. Details of the calculation are presented in Appendix A, where
we show that the moment generating function is given by Eq. (54), its first
moment, $dM_{C}(t)/{dt}|_{t=0}$, by Eq. (67), and its second moment,
$d^{2}M_{C}(t)/{dt^{2}}|_{t=0}$, by Eq. (68). Specifically, in each bin, we
set the first moment equal to the empirical mean, and set the second central
moment equal to the empirical variance:

$\displaystyle\frac{dM_{C}(t)}{dt}{\Big{|}}_{t=0}=\langle\hat{c}_{i}\rangle,$
(31)
$\displaystyle\frac{d^{2}M_{C}(t)}{dt^{2}}{\Big{|}}_{t=0}-\frac{dM_{C}(t)}{dt}{\Big{|}}_{t=0}^{2}=\frac{1}{1-N_{sj}}\sum_{i=1}^{N_{sj}}(\hat{c}_{i}-\langle\hat{c}_{i}\rangle)^{2},$
(32)

where
$\langle\hat{c}_{i}\rangle=\frac{1}{N_{sj}}\sum_{i=1}^{N_{sj}}\hat{c}_{i}$ is
the empirical mean. We set the population mean and standard deviation to
$c_{\star}=-0.0022$ and $R_{c}=0.0758$, the empirical mean and standard
deviation from the lowest two bins of the SDSS, SNLS and Low-$z$ surveys,
where colour-based selection effects are expected to be negligible. We then
solve the resulting coupled equations to determine $\\{c^{\rm
obs}_{sj},\sigma^{\rm obs}_{sj}\\}$.

In cases where ${N_{sj}}$ is small $(\sim 10$) and the sample variance is
small, the method of moments yields an estimate of $\sigma^{\rm obs}_{sj}$
near or equal to zero. This is similar to a well-known pathology in the
estimation of the shape parameter in the skew-normal distribution (see e.g.
Azzalini & Arellano-Valle (2012)). A simple solution is to impose a lower
cutoff to the value of $\sigma^{\rm obs}_{sj}$. We chose a cutoff value
$\sigma^{\rm obs}_{sj}>0.01$, since any smaller value leads to a selection
function indistinguishable from a step function. Tests of our method of
moments on simulated data are provided in Appendix A.

Because our selection effect model assumes in Eq. (27) that the colour
observations are independent, we set the corresponding covariances in the
variance-covariance matrix for the systematic effects to zero. This has a
minor effect on our estimates of the cosmological parameters, as those terms
are subdominant with respect to other off-diagonal terms.

We apply the method of moments to estimate $\Psi$ to each the four sub-surveys
separately, and we verify its accuracy by simulating colour observations from
the model of Eq. (27), assuming a Normal constant-in-redshift latent colour
distribution with mean $c_{\star}=0$ and standard deviation $R_{c}=0.1$. We
then apply the reconstructed selection function in each redshift bin, and
compare the resulting distribution of simulated SNe with the observed sample
within that bin. This comparison is shown in Fig. 2. For the simulation study
and the real JLA data in Section 3.2, we divided each sub-survey in 5
approximately equal-spaced redshift bins, with the following exceptions: the
HST data is treated in one single bin owing to the the small number of SNIa in
HST. For the same reason, in the Low-z survey data, the two highest-redshift
bins are combined into a single bin. For the SDSS survey, the right edge of
the $4^{\rm th}$ bin has its boundary shifted 0.015 in redshift space to the
right, in order to account for a discontinuity in the data, where the
population of objects have a lower colour and are separated by a gap in
redshift. The agreement between real data and simulation is excellent, thus
validating our approach.

Figure 2: Colour-based selection function in binned JLA data, with redshift
bins boundaries indicated by the vertical black lines.The blue circles are the
individual SNe, the blue errorbars represent the data mean and standard
deviation within each top-hat bin, while the orange errorbars give the mean
and standard deviation of simulated data from the model using the
reconstructed selection function in that bin (shifted horizontally for ease of
comparison).

###  2.5 A New Derivation of Peculiar Velocity Corrections

Our motion w.r.t. the CMB frame is measured precisely by the temperature
dipole observed in CMB anisotropies. The most accurate result is from Planck
Collaboration et al. (2020), giving a velocity $v_{\odot-{\rm CMB}}=369.82\pm
0.11$ km/s in the direction $l=264.021^{\circ}\pm
0.011^{\circ},b=48.253^{\circ}\pm 0.005^{\circ}$. This induces in the non-
relativistic limit a redshift correction (Davis et al., 2011)

$z_{\rm pec}^{\odot}\approx-\frac{v_{\odot-{\rm
CMB}}}{c}({\mathbf{n}}_{\rm\tiny{CMB}}\cdot{\mathbf{n}}_{\rm\tiny{SN}})\lesssim
10^{-3},$ (33)

where ${\mathbf{n}}_{\rm\tiny{CMB}}$ is a unit vector in the direction of the
CMB dipole and ${\mathbf{n}}_{\rm\tiny{SN}}$ is a unit vector in the direction
of the SNIa. Given the small uncertainties in the measurement for
$v_{\odot-{\rm CMB}}$ and the CMB dipole direction, we can consider $z_{\rm
pec}^{\odot}$ as known exactly, and thus neglect measurement error on this
quantity (as it is $\sim 3\times 10^{-7}$). Estimating $z_{\rm pec}^{\rm
SN}=v_{\rm pec}^{{\rm SN}}/c$ requires knowledge of the peculiar velocity of
the SNIa in the CMB frame, $v_{\rm pec}^{{\rm SN}}$. This can be measured
either from a peculiar velocity field survey (for example, by using the
Fundamental Plane (FP) relation or the Tully-Fisher (TF) relation to measure
the distance to a galaxy, and then subtracting from the observed velocity the
expansion component obtained from the Hubble-Lemaître law) or derived from
linear perturbation theory applied to a smoothed density field. The latter
approach has a long history, originally having been used to predict the
peculiar velocities of FP and TF samples (Hudson, 1993; Strauss & Willick,
1995; Davis et al., 1996) but more recently applied to peculiar velocity data
that include SNe type Ia. (Riess et al., 1997; Radburn-Smith et al., 2004;
Pike & Hudson, 2005; Neill et al., 2007; Turnbull et al., 2012; Carrick et
al., 2015; Boruah et al., 2020b; Lilow & Nusser, 2021; Stahl et al., 2021).

While the peculiar velocity of the SNIa becomes rapidly negligible for
$z\gtrsim 0.1$, it is important for local objects ($z\ll 0.1$), where it can
be significant w.r.t. the expansion velocity (up to $\sim 30\%$) and where it
leads to much larger changes in the apparent magnitude, due to the steeper
gradient of the distance modulus at low redshift. For example, at $z=0.01$ an
uncorrected SNIa peculiar velocity $v_{\rm pec}^{{\rm SN}}$ induces a redshift
systematic error $\delta z=v_{\rm pec}^{{\rm SN}}/c$, which corresponds to a
significant change in the theoretical distance modulus
$\delta\mu\approx\frac{d\mu}{dz}\delta z\approx 5/\ln(10)(\delta v/(cz))=0.14$
mag for $v_{\rm pec}^{{\rm SN}}=200$ km/s. To avoid difficulties with peculiar
velocities, earlier SNIa cosmological analyses routinely adopted a lower
redshift cutoff $z_{\rm cut}$, removing SNe below $z_{\rm cut}$; for example,
Kessler et al. (2009) used $z_{\rm cut}=0.02$; Riess et al. (2007) used
$z_{\rm cut}=0.023$. Recently, Huterer (2020) estimated the impact of
uncorrected peculiar velocities on the Pantheon sample from numerical N-body
simulations, and recommended a cutoff $z_{\rm cut}=0.02$ to protect against
significant bias to cosmological parameters. However, a better way that does
not discard useful data at low redshift is to assign uncertainties that scale
with distance, as we do here.

The JLA sample contains 37 SNe with $z_{\rm hel}<0.02$, and 110 with $z_{\rm
hel}<0.05$, for which an appropriate treatment of peculiar velocities is
required if they are to be used in the cosmological analysis – particularly in
our case, where we wish to use them to constrain a local dipole in the
expansion. To first order in redshift Eq. (2) gives

$\bar{z}=z_{\rm hel}-z_{\rm pec}^{\odot}-z_{\rm pec}^{\rm SN},$ (34)

meaning that the redshift of a comoving observer, $\bar{z}$, is obtained from
the measured heliocentric redshift by subtracting our local dipole ($z_{\rm
pec}^{\odot}$) and the redshift due to the SNIa peculiar velocity, $z_{\rm
pec}^{\rm SN}$.

The model used in Betoule et al. (2014) to estimate $z_{\rm pec}^{\rm SN}$ has
been criticised by C19, who highlighted potential bulk flow velocity
discontinuities at $z=0.04$, pointed out that peculiar velocity corrections
arbitrarily disappear beyond $200/h$ Mpc ($z\sim 0.067$, the limit of the
galaxy density field measurements from which the peculiar velocities were
derived) and that the residual uncorrelated velocity dispersion of
$\sigma_{v}=150$ km/s might be underestimated. While RH20 pointed out
technical flaws with the analysis of C19, it is important in the light of this
valid criticism to revisit the issue of low-redshift peculiar velocity
corrections here.

To this end, in this work we replace the peculiar velocity corrections used by
Betoule et al. (2014) – which rested on the IRAS PSCz catalogue from Branchini
et al. (1999) – with the more recent ones obtained by Carrick et al. (2015).
We follow Boruah et al. (2020a, b), who carried out a thorough comparison
between density reconstruction from galaxy redshift surveys and kernel
smoothing of peculiar velocity data methods. We adopt here their peculiar
velocity field inferred from 69,160 galaxies from the 2M++ galaxy redshift
catalogue (Lavaux & Hudson, 2011). The catalogue covers almost the entire sky
(with the notable exception of the plane of the galaxy), is highly complete
out to 200$/h$ Mpc ($z\sim 0.067$) in the region covered by 6dF and SDSS, and
out to 125$/h$ Mpc ($z\sim 0.041$) in the region covered by 2MRS. We thus
remove555Differently from C19, we do not remove the magnitude bias corrections
made to the JLA SNe, as they are important to account for selection effects,
nor do we neglect the contribution of peculiar velocities uncertainty to the
covariance matrix, which we re-derive for our case. the SNIa peculiar velocity
corrections for the low-$z$ JLA sample that are in common with the A2 sample
of Boruah et al. (2020a) (107 objects), and replace them with new values
obtained as follows.

The radial peculiar velocity for a SNIa at comoving distance $r$ and direction
${\mathbf{n}}_{\rm\tiny{SN}}$ is obtained from the luminosity-weighted density
field ${\bf v}_{\rm gal}$ as

$v_{\rm pec}^{{\rm
SN}}(r,{\mathbf{n}}_{\rm\tiny{SN}},\theta)={\mathbf{n}}_{\rm\tiny{SN}}\cdot\left(\beta_{v}{\bf
v}_{\rm gal}(r,{\mathbf{n}}_{\rm\tiny{SN}})+{\bf V}_{\rm ext}\right)$ (35)

where $\theta=\\{\beta_{v},{\bf V}_{\rm ext}\\}$, with a rescaling parameter
$\beta_{v}=0.411\pm 0.020$ and external residual bulk flow velocity (in
Galactic Cartesian coordinates) ${\bf V}_{\rm ext}=[52\pm 20,-163\pm 21,49\pm
16]$ km/s (how we treat and propagate the uncertainties in these values is
addressed in section 2.6). We follow the methodology of Boruah et al. (2020a),
with the difference that we only use the SFI$++$ peculiar velocity sample
(therefore not including A2 SNe data) in order to avoid circularity (i.e.
using SNe data to predict the peculiar velocity correction for the same SNe
data).

We do not wish to use the distance modulus information from a SNIa at this
stage of the analysis, only its observed redshift in the CMB restframe,
$\hat{z}_{\rm\tiny{CMB}}\equiv\hat{z}_{\text{hel}}-z_{\rm pec}^{\odot}$.
Firstly, the observed redshift in the CMB restframe is corrected to the
average redshift of the group to which the host galaxy belongs. Differently
from JLA, we correct the CMB restframe redshift for all host galaxies,
including those in clusters and poorer groups. This is necessary to suppress
the highly non-linear velocity contribution to the observed redshift, and it
leads to deviations of a few percent in $\hat{z}_{\rm\tiny{CMB}}$ in most
low-$z$ SNe, compared with the value used by JLA, see Fig. 3. However, there
are 6 SNe that show much larger changes in their CMB frame redshift, up to
$\sim 30\%$ (highlighted in green in Fig. 3); two of them (sn2007ci and
sn2001cz) are in common with the outliers in peculiar velocity, identified in
Fig. 4. It is worth noting that Carr et al. (2021) reviewed Type Ia SNe
literature and discovered that some SNIa had been reported with incorrect
redshifts and/or positions due to misprints. Two of the outliers in Fig. 4,
sn2008bf and sn1996c, are among the SNIa identified by Carr et al. (2021) as
requiring updates: sn2008bf suffers from an uncertain host galaxy
identification, while sn1996c had an incorrect sky location. Many other SNIa
had their redshifts and positions updated by Carr et al. (2021), but their
dataset is not publicly available yet and therefore we cannot make use of
their findings.

We compute the expected peculiar velocity by marginalizing over the unknown
comoving distance of the SNIa, $r$:

$\langle v_{\rm pec}^{{\rm SN}}\rangle=\int
drp(r|c\hat{z}_{\rm\tiny{CMB}})v_{\rm pec}^{{\rm
SN}}(r,{\mathbf{n}}_{\rm\tiny{SN}}),$ (36)

where $v_{\rm pec}^{{\rm SN}}$ is computed self-consistently from the flow
model, Eq. (35), and $p(r|\hat{z}_{\rm\tiny{CMB}})$ is the probability density
function (pdf) for $r$ given the observed redshift. This can be linked via a
variable transformation to the pdf for the true (latent) CMB redshift of the
SNIa, $z_{\rm\tiny{CMB}}(r)$, via:

$p(r|c\hat{z}_{\rm\tiny{CMB}})=p(cz_{\rm\tiny{CMB}}(r)|c\hat{z}_{\rm\tiny{CMB}})\bigg{|}\frac{\partial
cz_{\rm\tiny{CMB}}(r)}{\partial r}\bigg{|},$ (37)

where from Eq. (3)

$z_{\rm\tiny{CMB}}(r)=\bar{z}(r)+(1+\bar{z}(r))v_{\rm pec}^{{\rm
SN}}(r,{\mathbf{n}}_{\rm\tiny{SN}})/c,$ (38)

and the cosmological redshift at comoving distance $r$, $\bar{z}(r)$, is
computed for the $\Lambda$CDM concordance model, with parameters as in Table
LABEL:table:priors. Note that the transformation between $r$ and
$z_{\rm\tiny{CMB}}$ may not be unique due to the existence of triple-valued
regions. However, in practice we checked that the reconstructed velocity field
indeed leads to unique transformation between $r$ and $z_{\rm\tiny{CMB}}$ for
the parameters under consideration. Nevertheless, it can lead to ‘flat’
regions in redshift space (i.e, where $\frac{\partial
z_{\rm\tiny{CMB}}}{\partial r}\approx 0$), which result in large uncertainties
in the expected peculiar velocity.

The 2M++ reconstruction employs linear perturbation theory to predict the
velocities. As shown in Carrick et al. (2015), this leads to an uncertainty
due to non-linearities of $\sigma_{v}^{\text{\tiny{NL}}}=150$ km/s. Therefore,
assuming a Gaussian uncertainty, we can write the probability of the predicted
redshift in the CMB frame for a SNIa at comoving distance $r$ given its
observed redshift transformed in the CMB frame in Eq. (37) as:

$p(cz_{\rm\tiny{CMB}}(r)|c\hat{z}_{\rm\tiny{CMB}})={\mathcal{N}}(c\hat{z}_{\rm\tiny{CMB}},(\sigma_{v}^{\text{\tiny{NL}}})^{2})$
(39)

As a check for the robustness of our method, we also estimated the predicted
peculiar velocity using an iterative prescription. In this alternative method,
we start from the observed CMB restframe redshift as an initial, rough
approximation for the comoving distance (under the fiducial $\Lambda$CDM
assumption), taking $v_{\rm pec}^{{\rm SN}}=0$. In subsequent iterations, an
updated estimate for the comoving distance is obtained using Eq. (38) and the
velocity estimate we get from the reconstruction at the given comoving
distance. This step is repeated until convergence. The iterative method gives
broadly consistent results as the marginalization-based method described
above. However, the iterative prescription may underestimate the uncertainty
in the vicinity of triple-valued regions, and therefore we elected to use the
marginalization method instead.

Our new peculiar velocity corrections obtained from Eq. (36) are compared
against those used in Betoule et al. (2014) in Fig. 4, which are obtained from
the JLA data products via666In Eq. (40) we use the notation introduced in this
paper but we notice that our $\bar{z}$ is (confusingly) denoted
“$z_{\rm\tiny{CMB}}$” in the Betoule et al. (2014) and associated data
products.

$v_{\rm corr,JLA}=c\left(\frac{1+z_{\rm hel}}{(1+\bar{z})(1+z_{\rm
pec}^{\odot})}-1\right).$ (40)

The correlation coefficient between our peculiar velocity corrections and that
used in Betoule et al. (2014) is only $\sim 0.60$. There are also 6 SNe
(labeled in Fig. 4) that are more than 3 standard deviations of the sample
away from the identity line. Several reasons can explain the differences
between our peculiar velocities and those used in JLA: firstly, the density
field used to predict peculiar velocities in JLA was based on the IRAS PSCz
survey, which is likely to be noisier than 2M++ everywhere. Secondly, in the
JLA analysis only galaxies in rich clusters are corrected to the mean redshift
of the cluster, whereas here we correct the CMB restframe redshift for all
host galaxies, including those in clusters and poorer groups.

Figure 3: Left panel: comparison of the CMB restframe redshift used in this
work for the low-redshift sample (horizontal axis, 107 SNe) with that of the
JLA analysis (vertical axis). There are 6 SNe (highlighted in green) with a
difference exceeding 3 standard deviations around the identity line (dashed).
Right panel: the same comparison but showing the fractional differences
between redshifts on the vertical axis. Figure 4: Comparison of the peculiar
velocity corrections used in this work (horizontal axis) with those adopted in
the JLA analysis (vertical axis) for the 107 SNe in the A2 low-redshift sample
common to the JLA data. The colour coding gives the CMB restframe redshift
used in this work. Outliers (defined as more than 3 standard deviations away
from the dashed identity line) are labeled.

For SNe beyond $z=0.067$ (in the SNLS, SDSS and high-$z$ samples), the
relative importance of the peculiar velocity corrections diminishes as the
recession velocity of the Hubble flow increases, so the detailed treatment
becomes less important. In the JLA analysis, peculiar corrections have been
set to 0 beyond redshift $z\sim 0.067$, the limit of the galaxy survey from
which said corrections were derived. We do the same here, noting that beyond
redshift $z\sim 0.067$ the exact treatment of peculiar velocity correction
becomes less crucial since their relative importance diminished quickly.

###  2.6 Peculiar Velocity Correction Uncertainties

As we have replaced the peculiar velocities of JLA with our own, new
estimates, we also update their error analysis associated with the peculiar
velocities treatment. First, we remove the diagonal term from the JLA
statistical covariance matrix representing the uncertainty from peculiar
velocity corrections. We also remove the “pecvel” contribution to the
systematic covariance matrix. We replace them with the following terms,
flowing from our updated treatment of peculiar velocity corrections.

From Eq. (34), the covariance of the cosmological redshift value for SNe $i,j$
is given by (without including the negligible error in $z_{\rm pec}^{\odot}$)
is given by:

$\displaystyle\Xi^{2}_{z,ij}$ $\displaystyle=\delta_{ij}^{\rm
K}\sigma^{2}_{z,i}+\Sigma^{2}_{z,ij},$ (41) $\displaystyle\Sigma^{2}_{z,ij}$
$\displaystyle=\delta_{ij}^{\rm K}\left(\sigma^{2}_{\rm NL}+\sigma^{2}_{\rm
2M++}(z_{\rm\tiny{CMB}})+\sigma^{2}_{v,i}\right)/c^{2}+C^{\rm
flow}_{ij}/c^{2},$ (42)

where $\sigma_{z,i}$ is the spectroscopic redshift measurement uncertainty,
$\sigma_{{\rm NL}}=150$ km/s is the uncertainty in the peculiar velocity due
to non-linearities, $\sigma_{\rm 2M++}$ is the redshift-dependent uncertainty
due to survey incompleteness, $\sigma_{v,i}$ is the standard deviation of the
average peculiar velocity prediction, Eq. (36), and $C^{\rm flow}_{ij}$ is the
correlated covariance coming from uncertainty in the flow model
($\delta_{ij}^{\rm K}$ is the Kronecker delta). We address each term in turn.

For the 107 low-$z$ SNe in our analysis, the largest reported statistical
uncertainty in their redshift measurement in the JLA data release is
$\max_{i}\sigma_{z,i}^{\rm JLA}=0.0014$, corresponding to a velocity
uncertainty of 420 km/s. The JLA data release also has $\sigma_{z,i}^{\rm
JLA}=0$ for 10 of the 107 low-$z$ SNe and $\sigma_{z,i}^{\rm JLA}<0$ for 42
SNIa. In order to resolve the issue of 0 or negative redshift uncertainties
and to be conservative, we adopt the following prescription for the standard
deviation of the spectroscopic uncertainty:

$\sigma_{z,i}=\max(\sigma_{z,i}^{\rm JLA},5\cdot 10^{-4}).$ (43)

where the floor value of $5\cdot 10^{-4}$ represents the typical uncertainty
of spectroscopic redshift determination from host galaxies spectra.

The term $\sigma^{2}_{{\rm NL}}=150$ km/s represents uncertainty in the linear
velocity prediction due to unaccounted-for non-linearities, which we fix at
the value recommended in Carrick et al. (2015). However, the uncertainty in
our reconstructed peculiar velocity increases with redshift, an effect that
was ignored in previous work: firstly, the predicted peculiar velocities for
tracers near the outer edge of the 2M$++$ catalogue (200 Mpc/$h$) have larger
uncertainty because of unaccounted-for structures outside of the survey
limits, as well as because of lack of coverage beyond 125 Mpc/$h$ for part of
the sky (Hollinger & Hudson, 2021); secondly, the noise increases at larger
distances due to the smaller number of galaxies with larger weights that are
used to represent the density field (Lilow & Nusser, 2021). We capture these
effects via the redshift-dependent term:

$\sigma_{\rm
2M++}(z_{\rm\tiny{CMB}})=\begin{cases}\sigma_{1}z_{\rm\tiny{CMB}}&\mbox{ for
}z_{\rm\tiny{CMB}}<z_{400},\\\ \sigma_{1}z_{400}&\mbox{ for
}z_{\rm\tiny{CMB}}\geq z_{400},\end{cases}$ (44)

where $z_{400}=0.138$ is the redshift corresponding to a radial comoving
distance of 400$/h$ Mpc, and $\sigma_{1}$ is chosen so that the total peculiar
velocity rms beyond $z_{400}$, i.e. $(\sigma^{2}_{{\rm NL}}+\sigma_{\rm
2M++}^{2}(z_{400}))^{1/2}$, equals $380$ km/s. This prescription also
approximately matches the $\Lambda$CDM prediction at the 2M++ boundary,
$z=0.067$, where $(\sigma^{2}_{\rm NL}+\sigma^{2}_{\rm
2M++}(0.067))^{1/2}=227$ (km/s). This is in contrast with the original JLA
analysis which uses a redshift-independent 150 km/s uncertainty throughout the
redshift range.

The term $\sigma^{2}_{v,i}$ is the variance of $v_{\rm pec}^{{\rm SN}}$ under
the distribution given by Eq. (37), i.e.

$\sigma^{2}_{v,i}=\langle(v_{\rm pec}^{{\rm SN}})^{2}\rangle-\langle v_{\rm
pec}^{{\rm SN}}\rangle^{2}.$ (45)

Finally, we translate the redshift covariance of Eq. (41) into a magnitude
covariance via linear propagation of errors using the isotropic distance
modulus of Eq. (11), i.e.

$\sigma^{2}_{m,ij}=\Sigma^{2}_{z,ij}\frac{\partial\mu_{I}}{\partial\bar{z}}{\bigg{\rvert}}_{\begin{subarray}{l}z_{\rm
hel}=\hat{z}_{{\rm hel},i}\\\
\bar{z}=\bar{z}_{i}\end{subarray}}\frac{\partial\mu_{I}}{\partial\bar{z}}{\bigg{\rvert}}_{\begin{subarray}{l}z_{\rm
hel}=\hat{z}_{{\rm hel},j}\\\
\bar{z}=\bar{z}_{j}\end{subarray}}+\delta_{ij}^{\rm
K}\sigma^{2}_{z,i}\left(\frac{\partial\mu_{I}}{\partial z_{\rm
hel}}{\bigg{\rvert}}_{\begin{subarray}{l}z_{\rm hel}=\hat{z}_{{\rm hel},i}\\\
\bar{z}=\bar{z}_{i}\end{subarray}}\right)^{2}$ (46)

where $\bar{z}_{i}$ is computed from Eq. (34), and the distance modulus
derivatives are evaluated at the fiducial cosmological parameter values given
in Table LABEL:table:priors.

So far, we have considered a fixed value for the flow parameters,
$\theta=\\{\beta_{v},{\bf V}_{\rm ext}\\}$ entering in Eq. (35). The
uncertainties in the inferred flow parameters lead to correlated uncertainties
in the peculiar velocities which needs to be accounted for, and that in
previous work are usually considered a source of systematic error. Our
parameterised flow model allows us to translate them into a statistical error,
as follows. In order to estimate the covariance coming from uncertainty in the
flow model parameters, we draw $10^{4}$ posterior samples of the flow
parameters $\theta_{k}$ ($k=1,\dots,10^{4}$) from the fitted flow model, using
the method of Boruah et al. (2020b), and we calculate the average peculiar
velocity, $\langle{v_{\rm pec}^{{\rm SN}}},_{i}\rangle$ for all the 107 SNe in
our low-$z$ sample from those samples. We then estimate the covariance of the
average peculiar velocity between SNIa $i$ and $j$ as
$C_{ij}^{\rm{flow}}=\text{Cov}(\langle{v_{\rm pec}^{{\rm
SN}},}_{i}\rangle,\langle{v_{\rm pec}^{{\rm SN}},}_{j}\rangle)$, where the
covariance matrix is computed from the $k$ samples. Since the value of
$\sigma_{v,i}$ above varies among the $k$ samples (although the variation is
small, $\sol 10\%$), we use the average of $\sigma_{v,i}$ from the $10^{4}$
posterior samples. In accord with terminology used in the literature, we call
this term the “systematic uncertainty”, although as noted above we have
actually translated it into a statistical uncertainty. There are no changes to
the $x_{1}$ and $c$ terms of the covariance and these are left unchanged from
the original JLA analysis.

We show in the left panel of Fig. 5 the square root of the diagonal entries of
the peculiar velocities covariance matrix, translated into magnitude
covariance,i.e. $\sigma_{m,ii}$ in Eq. (46)). The right panel shows the square
root of the diagonal entries of the total magnitude covariance matrix
(including all other magnitude uncertainties). Our values are compared with
the original JLA values on the same figure. The largest difference in the
total value of $\sigma_{m}$ appears in the low redshift range, where our re-
analysis modifies the associated peculiar velocities, which are dominant in
this redshift range. In general, the net effect is to increase the statistical
uncertainty while decreasing the systematic uncertainty with respect to the
JLA analysis: at the median redsfhit of the low-$z$ sample, $z=0.0243$, the
average diagonal $\sigma_{m}$ due to statistical uncertainty in the peculiar
velocities is $0.076$ in our analysis (vs $0.045$ JLA), while the average
systematic diagonal error is $0.020$ in our work (vs $0.039$ JLA). Overall,
the total magnitude uncertainty due to peculiar velocities is increased by
$\sim 30\%$ in our analysis (at the median redshift) compared to JLA. In the
right panel of Fig. 5, we compare the total uncertainties on the apparent
magnitude (including all statistical and systematic uncertainties) between
this work and the JLA analysis, showing that our magnitude uncertainties are
generally larger, especially at low redshifts where the new peculiar velocity
uncertainties dominate the error budget.

Figure 5: Left panel: statistical and systematic (diagonal entries of the
covariance matrix only) uncertainties in the apparent magnitude induced by
peculiar velocity corrections for our re-analysis compared with JLA. Right
panel: total diagonal $\sigma_{m}$ for our reanalysed and original JLA data.
These components include the statistical error from the SALT2 fits and
peculiar velocities as well systematics from several other components that are
outlined in Betoule et al. (2014).

###  2.7 Choice of Priors

As always in a Bayesian analysis, particular attention must be paid to priors,
especially in the present case where we are interested in performing not only
parameter inference but also model comparison (for an overview of the issue,
see e.g. Trotta (2008)).

The priors for the cosmological parameters ${\mathscr{P}}_{1}$ and for the
other parameters in the hierarchical model are chosen as in Shariff et al.
(2016), to which we refer for fuller details. A summary is provided in Table
LABEL:table:priors. For the cosmographic expansion parameters
${\mathscr{P}}_{2}$, the $\mathrm{Uniform}(0,2)$ priors in
$\Omega_{m},\Omega_{\Lambda}$ translate into the following non-uniform prior
for the deceleration parameter $q_{0}=\Omega_{m}/2-\Omega_{\Lambda}$:

$p(q_{0})=\begin{cases}\frac{2}{3}(q+2)&\mbox{ for }-2\leq q_{0}<-1,\\\
2/3&\mbox{ for }-1\leq q_{0}<0,\\\ \frac{2}{3}(1-q)&\mbox{ for }0\leq
q_{0}\leq 1,\\\ 0&\mbox{ otherwise}.\end{cases}$ (47)

Given that the likelihood’s support is almost entirely within the region
$-1.0\lesssim q_{0}\lesssim 0$, we choose the slightly simplified, uniform
prior $q_{0}\sim\mathrm{Uniform}(-2,1)$. Since the jerk $j_{0}$ and spatial
curvature $\Omega_{\kappa}$ appear in the degenerate combination
$j_{0}-\Omega_{\kappa}$ in the second-order term in the cosmographic
expansion, Eq. (10), we adopt a uniform prior on this combination, with ranges
motivated by the range of physically plausible curvature values.

The dipole amplitude parameter, $D_{\mu}$, is a positive quantity for which
the obvious choices of priors are a uniform or log-uniform prior. The lower
boundary of the former is naturally 0, while for the latter a lower cutoff
must be imposed in order for the prior to be proper (i.e., normalizable). This
however is problematic for the model selection outcome: as the likelihood
becomes flat (i.e., insensitive to the value of $D_{\mu}$) once $D_{\mu}$
drops below a certain threshold, the marginal posterior becomes equal to the
marginal prior for arbitrarily lower values of $D_{\mu}$ with a log-uniform
prior. This has an influence on the Bayes factor, as a larger prior range
(i.e., a lower $D_{\mu}$ cut-off in a log-uniform prior) leads to a less
favourable model selection outcome for the anisotropic model. In order to
avoid this difficulty, we choose a uniform prior on $D_{\mu}$ itself. The
upper prior cutoff is set by the characteristic scale expected by a dipole
signal. This could be gleaned from a theoretical model, or, in a
phenomenological approach such as ours, guided by the order of magnitude of
previous upper limits on the parameter, which is of order $\sim 10^{-3}$ (see
Table LABEL:table:compiled_author_results). Such upper limits can however be
considerably relaxed in the case of a dipole that is decaying with redshift,
leading to upper limits of order $\sim 10^{-1}$ even when no dipole is present
(see our simulated case of Fig. 8). In order to accommodate such a scenario,
we choose a uniform prior $D_{\mu}\sim\mathrm{Uniform(0,0.2)}$.

We choose to sample the area of the sky in a uninformative manner since we do
not have any prior belief of the directions a dipole might be pointing to.
Requiring rotational invariance on the surface of the 2-sphere leads to a
uniform distribution on the Galactic longitude of the dipole vector,
$l_{d}\sim\mathrm{Uniform}(0,2\pi)$ (in radians), and a uniform distribution
on the cosine of the latitude of the dipole vector,
$\cos(b_{d})\sim\mathrm{Uniform}[0,1]$, with $b_{d}\in[0,\pi/2]$. Flipping the
sign of $b_{d}$ is equivalent to the transformation $l_{d}\rightarrow
l_{d}+\pi\,\text{mod}\,2\pi$ and $D_{\mu}\rightarrow-D_{\mu}$. Hence in order
to cover the possibility of a dipole pointing in a direction in the southern
Galactic hemisphere, we extend the dipole amplitude to negative values, and
therefore our prior is modified to $D_{\mu}\sim\mathrm{Uniform(-0.2,0.2)}$.
Similar considerations lead to a prior for the dipole amplitude on the
deceleration parameter $D_{q_{0}}\sim\mathrm{Uniform(-30,30)}$.

For the prior on the exponential scale parameter, $S$, we need to select a
lower boundary (lest $D_{\mu}$ becomes unidentifiable and to stop pathologies
associated with $S=0$), which we take to be the scale of the lowest redshift
SNe in our data, namely $S=0.01$; for the upper boundary, we take $S=0.1$ as
it is known that the bulk flow does not disappear at least out to $z\sim
0.067$. In summary, our prior is thus $S\sim\textrm{Uniform}[0.01,0.1]$.

##  3 Simulations and Tests of Methodology

###  3.1 Simulated Data

Table 1: Supernova population parameters, cosmological parameters and dipole model parameters adopted in this work, together with the prior choices and fiducial values for simulation studies. ‘SD’ stands for standard deviation SNIa population distributions and covariates | Symbol | Fiducial value | Prior distribution
---|---|---|---
Mean absolute magnitude of SNIa | $M_{0}$ | -19.3 | $M_{0}\sim\mathrm{Normal}(-19.3,2^{2})$
Residual scatter of SNIa magnitude after corrections | $\sigma_{\rm res}$ | 0.1 | $\sigma_{\rm res}\sim\mathrm{InvGamma(0.003,0.003)}$
Coefficient of stretch covariate | $\alpha$ | 0.14 | $\alpha\sim\mathrm{Uniform(0,1)}$
Coefficient of colour covariate | $\beta$ | 3.2 | $\beta\sim\mathrm{Uniform(0,4)}$
Mean of stretch | $x_{*}$ | 0.0 | $x^{*}\sim\mathrm{Normal}(0,10^{2})$
Mean of colour | $c_{*}$ | 0.0 | $c^{*}\sim\mathrm{Normal}(0,1^{2})$
SD of stretch distribution | $R_{x}$ | 1.0 | $R_{x}\sim\mathrm{LogUniform}(-5,2)$
SD of colour distribution | $R_{c}$ | 0.1 | $R_{c}\sim\mathrm{LogUniform}(-5,2)$
Parameters controlling the expansion history |  |  |
Matter energy density | $\Omega_{M}$ | 0.3 | $\Omega_{M}\sim\mathrm{Uniform}(0,2)$
Dark energy density | $\Omega_{\Lambda}$ | 0.7 | $\Omega_{\Lambda}\sim\mathrm{Uniform}(0,2)$
Deceleration Parameter | $q_{0}$ | -0.55 | $q_{0}\sim\mathrm{Uniform}(-2,1)$
Jerk and spatial curvature | $j_{0}-\Omega_{K}$ | 1 | $j_{0}-\Omega_{K}\sim\mathrm{Uniform}(-2,2)$
Hubble-Lemaître constant [km/s/Mpc] | $H_{0}$ | 72 | Fixed
Anisotropic expansion parameters |  |  |
Galactic longitude of dipole [rad] | $l_{d}$ | $\\{-,4.60\\}$ | $l_{d}\sim\mathrm{Uniform}(0,2\pi)$
Galactic latitude of dipole [rad] | $b_{d}$ | $\\{-,0.84\\}$ |  $\cos(b_{d})\sim\mathrm{Uniform}(0,1)$, $b_{d}\in[0,\pi/2]$
Amplitude of dipole on $\mu$ | $D_{\mu}$ | { 0, 0.02} | $D_{\mu}\sim\mathrm{Uniform}(-0.2,0.2)$
Amplitude of dipole on $q_{0}$ | $D_{q_{0}}$ | { 0, 10} | $D_{q_{0}}\sim\mathrm{Uniform}(-30,30)$
Exponential scale of dipole | $S$ | { –, 0.026} | $S\sim\mathrm{Uniform}(0.01,0.10)$

To test the ability of our setup to recover the anisotropy parameters within
the BAHAMAS framework, we forward simulate data from our model in a similar
manner to March et al. (2011), adding however JLA-like Galactic coordinates to
our simulated SNe in order to carry out inference on a potential dipole. The
fiducial values for the parameters are listed in Table LABEL:table:priors. We
consider a case with no dipole ($D_{\mu}=0$) and a case with a large dipole
($D_{\mu}=0.02$), pointing in the approximate direction of the CMB dipole and
with an exponential scale, $S=0.026$, matching the value preferred by the
results of C19. The steps to generate the simulated data are as follows:

  1. 1.

Draw a value for the latent CMB restframe redshift, $z_{\text{CMB},i}$, for
each SNIa. $z_{\rm hel}$ is assumed to equal this CMB redshift. We draw 740
objects at the same redshifts as the real JLA data to ensure the binning of
the data is the same when applying selection effects correction on simulated
data as it is on real data. Notice that our simulations do not include the
issue of peculiar velocity corrections, which are assumed to have already been
performed.

  2. 2.

Compute $\mu_{i}(\hat{z},\Theta)$ using the fiducial values for our chosen
cosmology, whether $\Lambda$CDM or the Cosmographic expansion.

  3. 3.

Apply the dipolar modulation to the distance modulus using (14).

  4. 4.

Draw the latent parameters $x_{1i}$, $c_{i}$ and $M_{i}$ from the normal
distributions, $x_{1i}\sim N(x_{*},R_{x}^{2})$, $c_{i}\sim N(c_{*},R_{c}^{2})$
and $M_{i}\sim N(M_{0},\sigma_{\rm res}{}^{2})$ respectively.

  5. 5.

Compute $m_{Bi}$ using $x_{1_{i}}$, $c_{i}$, $M_{i}$ and the Phillip’s
relation equation in (11).

  6. 6.

Draw the value of the standard deviations $\sigma_{x_{1i}}$, $\sigma_{c_{i}}$,
and $\sigma_{m_{i}}$, from the appropriate normal distributions fitted to the
errors in the JLA data and use them to construct the $3\times 3$ covariance
matrix for each SNe as
$C=\mathrm{diag}(\sigma_{c_{i}}^{2},\sigma_{x_{1i}}^{2},\sigma_{m_{i}}^{2})$.

  7. 7.

Draw the observed SALT2 parameters from $\hat{x}_{1i}\sim
N(x_{1i},\sigma_{x_{1i}}^{2})$, $\hat{c}_{i}\sim N(c_{i},\sigma_{ci}^{2})$ and
$\hat{m}_{Bi}\sim N(m_{Bi},\sigma_{mi}^{2})$.

  8. 8.

Apply the selection function on the colour values drawn in the previous step.
We use the values of $\sigma^{\rm obs}_{sj}$ and $c^{\rm obs}_{sj}$ inferred
from the real JLA data in section 2.4 for the redshift bins the SNIa falls in.
If a given SNIa is not selected, we cycle back to step (iv) and redraw that
SNIa. The process ends when all 740 objects are selected.

  9. 9.

Generate positions for SNIa in the sky which match the positions of the JLA
data. The non-isotropic distribution of the data has an effect on our ability
to constrain a dipole, so it is important to match the real JLA SNIa positions
for a realistic simulation.

Figure 6: Simulated SNIa data generated using the JLA dataset as a reference
point. The simulations assume a $\Lambda$CDM model and no anisotropy in the
distance modulus for the left most model. A dipole of value $D_{\mu}=0.02$ is
present for the second plot. The third plot also has this value of the dipole,
but restricted to a local scale ($z\sim 0.1$) by multiplying the dipole term
by the function $F(z)=\exp{(-z/0.026)}$. For this third plot the redshift has
been truncated to only show the redshift range where the dipole is noticeable.
Figure 7: Simulated SNIa data generated using the JLA dataset as a reference
point, assuming $\Lambda$CDM as a fiducial model: distribution of colour
($c$), stretch ($x_{1}$) and the respective standard deviation entering the
observational error covariance matrix, as well as the SNIa’s sky coordinates.

An example of the simulated data, compared with the JLA dataset, is shown in
Fig. 6 and 7, for a $\Lambda$CDM realisation. Some banding of the data can be
seen of the plots for $\sigma_{x_{1}}$, $\sigma_{c}$, and $\sigma_{m_{B}}$.
This results from the different survey components comprising the JLA data. We
can see that our simulated data does not capture this banding in the errors
because we draw our errors from Gaussian distributions, but this approximation
has no quantitative effect on our simulations. The well defined ‘stripe’ in
the plot of the galactic coordinates of the data is from the SDSS component of
the JLA data.

We test the ability of our method to retrieve the true fiducial values listed
in Table LABEL:table:priors from the simulated data under the 7 different
scenarios below:

  1. 1.

$\Lambda$CDM-Isotropic simulates an isotropic universe ($D_{\mu}=0$) from
$\Lambda$CDM and the inference also assumes a $\Lambda$CDM cosmology.

  2. 2.

$\Lambda$CDM-D simulates a dipole ($D_{\mu}=0.02$) with $F(z)=1$ from
$\Lambda$CDM and the inference also assumes a $\Lambda$CDM cosmology.

  3. 3.

$\Lambda$CDM-D-exp simulates a dipole ($D_{\mu}=0.02$) with
$F(z)=\exp{(-z/S)}$ from $\Lambda$CDM and inference also assumes a
$\Lambda$CDM cosmology, with the additional free parameter $S$.

  4. 4.

Cosmographic-Isotropic simulates an isotropic universe ($D_{\mu}=0$) from a
cosmographic expansion and inference also assumes a cosmographic expansion.

  5. 5.

Cosmographic-D simulates a dipole with $F(z)=1$ from a cosmographic expansions
and inference also assumes a cosmographic expansion.

  6. 6.

Cosmographic-D-exp simulates a dipole ($D_{\mu}=0.02$) with
$F(z)=\exp{(-z/S)}$ and inference also assumes a cosmographic expansion, with
the additional free parameter $S$.

  7. 7.

Cosmographic-D-exp* simulates data as in $\Lambda$CDM-D-exp but the
reconstruction adopts the cosmographic expansion instead. This serves to
assess the bias in parameter reconstruction from the cosmographic expansion
when the reconstruction model is miss-specified.

For each parameterisation we generate 10 realisations of the data; each
realisation contains 740 SNIa objects to be similar in size to the JLA data.
The posterior results are averaged over the 10 realisations for each
parameterisation when we reconstruct the parameters in section 3.2. We do not
use the Gibbs sampler of Shariff et al. (2016) for posterior sampling, but
rather adopt PyMultiNest (Buchner et al., 2014), an implementation of the
Nested Sampling algorithm Multinest (Feroz & Hobson, 2008; Feroz et al., 2009;
Feroz et al., 2013). The benefit is that we can also compute the Bayesian
evidence which we will use for Bayesian model comparison in section 4.2.

###  3.2 Parameter Reconstruction from Simulations

Table 2: Difference between the 1D marginal posterior mean and the true value
used to simulated JLA-data, averaged over 10 data realizations each, $\Delta
x\equiv x_{\rm true}-\bar{x}$ for each parameter $x$. In parenthesis, the
difference is expressed in units of the 1D posterior standard deviation.The
top (middle) section shows the case where $\Lambda$CDM (Cosmographic
expansion) is assumed both for the simulation and the inference; the bottom
row assumes $\Lambda$CDM for the simulation and the the Cosmographic expansion
for the inference – a case of model miss-specification.

Model | $\Delta q_{0}$ | $\Delta(j_{0}-\Omega_{\kappa})$ | $\Delta\Omega_{M}$ | $\Delta\Omega_{\Lambda}$ | $\Delta l(\mathrm{rad)}$ | $\Delta b(\mathrm{rad)}$ | $|\Delta D_{\mu}|(\times 10^{-2})$ | $\Delta S(\times 10^{-2})$
---|---|---|---|---|---|---|---|---
Simulation: $\Lambda$CDM; reconstruction: $\Lambda$CDM
(ia) $\Lambda$CDM-Isotropic, $F(z)=0$ | - | - | 0.03(0.51$\sigma$) | 0.05(0.36$\sigma$) | - | - | < $1.00\times 10^{-1}$(2$\sigma$) | -
(ib) $\Lambda$CDM-Isotropic, $F(z)=\exp(-z/S)$ | - | - | 0.03(0.51$\sigma$) | 0.05(0.36$\sigma$) | - | - | < $4.96\times 10^{-1}$ (2$\sigma$) | < 9.29 (2$\sigma$)
(ii) $\Lambda$CDM-D | - | - | -0.02(0.19$\sigma$) | 0.02(0.10$\sigma$) | -0.003(0.20$\sigma$) | -0.003(0.25$\sigma$) | 0.036(0.13$\sigma$) | -
(iii) $\Lambda$CDM-D-exp | - | - | 0.03(0.91$\sigma$) | 0.07(0.72$\sigma$) | 0.02(0.20$\sigma$) | -0.02(0.17$\sigma$) | -0.16(0.28$\sigma$) | -0.10(0.21$\sigma$)
Simulation: Cosmographic; reconstruction: Cosmographic
(iv) Cosmographic-Isotropic, $F(z)=0$ | -0.02(0.20$\sigma$) | 0.03(0.01$\sigma$) | - | - | - | - | - | -
(v) Cosmographic-D | -0.06(0.56$\sigma$) | 0.27(0.58$\sigma$) | - | - | -0.0116(0.37$\sigma$) | -0.003(0.26$\sigma$) | -0.01(0.35$\sigma$) | -
(vi) Cosmographic-D-exp | -0.01(0.10$\sigma$) | -0.05(0.14$\sigma$) | - | - | -0.04(0.20$\sigma$) | 0.01(0.08$\sigma$) | -0.13(0.30$\sigma$) | 0.16(0.37$\sigma$)
Simulation: $\Lambda$CDM; reconstruction: Cosmographic
(vii) Cosmographic-D-exp | -0.11(1.05$\sigma$) | -0.63(1.37$\sigma$) | - | - | 0.11(0.60$\sigma$) | 0.038(0.37$\sigma$) | -0.128(0.44$\sigma$) | -0.165(0.24$\sigma$)

We use BAHAMAS to construct posterior distributions (averaged over 10 data
realizations) for the set of cosmological parameters ${\mathscr{P}}_{1}$
(i.e., $\Lambda$CDM) or ${\mathscr{P}}_{2}$ (i.e., cosmographic expansion) and
anisotropy parameters $\\{l_{d},b_{d},D_{\mu},S\\}$. Although they are sampled
over during reconstruction, SNe population parameters and SALT2 coefficients
are numerically marginalised over in corner plots and not visualised as they
are not the focus of this paper. The difference between the 1D marginal
posterior mean (averaged over realizations) and the true value of each
parameter is displayed in Table LABEL:table:simulated_results. We observe that
in all cases except for the scenario Cosmographic-D-exp* the difference is a
fraction of a standard deviation, hence entirely within realization and
sampling noise. The model misspecification of Cosmographic-D-exp*, however,
does lead to shifts of up to $\sim 1.4\sigma$ in the reconstructed
cosmographic parameters, a reflection of the fact that the data have been
generated under a different model, namely $\Lambda$CDM, than has been assumed
in the reconstruction. However, the difference for the anisotropy parameters
remains below $0.5\sigma$.

It is instructive to investigate the expected constraints on the dipole
amplitude when the simulated data are from a isotropic universe (scenario
$\Lambda$CDM-Isotropic), shown in Fig. 8. The 1-sided $95\%$ upper limit we
can place on $D_{\mu}$ for this simulated data is $D_{\mu}\leq 8.08\times
10^{-4}$ (left panel). When however we introduce the additional freedom of an
exponential scale parameter $S$, the constraints in $D_{\mu}$ degrade by two
orders of magnitude, as a small value of $S$ confines any anisotropy to very
small redshifts where the statistical power of our data is small and therefore
degenerate with many values of $D_{\mu}$ which leads to the entire prior space
on $D_{\mu}$ being well explored right up to the prior edge (right panels of
Fig. 8.). The 95% upper bound becomes $D_{\mu}<4.36\times 10^{-2}$ and
$S<4.56\times 10^{-2}$. Qualitatively similar results apply for the
Cosmographic expansion.

The above simulations do not include colour-based selection effects for
simplicity. A simulation study including colour-based selection appears in
Appendix A (see Fig. 14) and shows that averaging the posterior distributions
across replicates yields 2-dimensional and 1-dimensional marginal
distributions that are centered on the true values of the parameters.

Figure 8: Posterior distributions from simulated JLA-like data, averaged over
10 data realizations. The simulations assume a $\Lambda$CDM isotropic universe
($D_{\mu}=0$); the reconstructions allow for a dipole in the distance modulus
with $F(z)=1$ (left) and $F(z)=\exp{(-z/S)}$ (right). Contours in 2D plots
enclose joint 1, 2 and 3 $\sigma$ highest posterior density (HPD) credible
regions; the vertical dashed line in the 1D marginals for $D_{\mu}$ and $S$
delimits the 2$\sigma$ upper limit. Notice the different scales in the axis of
$D_{\mu}$ between the two cases.

##  4 Results

###  4.1 Parameter Inference

We begin by presenting the impact of our new treatment of peculiar velocities
and colour selection effects correction on the constraints on
$\Omega_{m},\Omega_{\Lambda}$ for the case of an isotropic universe, as shown
in Fig. 9. In the left panel, we compare the constraints using the old JLA
peculiar velocities and $z_{\rm\tiny{CMB}}$ as in Betoule et al. (2014)
(green) to the ones obtained using our new, group-corrected values of the CMB
restframe redshifts (blue), and additionally replacing the JLA peculiar
velocities with our newly derived ones (orange). The constraints on the
parameters of our models from JLA data are summarised in Table
LABEL:table:final_jla_results, for both the $\Lambda$CDM model and the
Cosmographic expansion. In the top section we also investigate the impact of
our newly derived values for $z_{\rm\tiny{CMB}}$, peculiar velocity
corrections and colour-based selection effects on $\Omega_{m}$,
$\Omega_{\Lambda}$ constraints in an isotropic universe. Starting from the
same treatment as Betoule et al. (2014), we find $\Omega_{m}=0.306\pm 0.087$,
in good agreement with the value in Betoule et al. (2014),
$\Omega_{m}=0.295\pm 0.034$ but with significantly larger uncertainty, perhaps
on account of the different statistical approach. When replacing the CMB
restframe redshifts used in Betoule et al. (2014) with the ones presented
here, we find $\Omega_{m}=0.253\pm 0.089$, a shift of about half a standard
deviation according to our uncertainty, but of $1.5\sigma$ in units of the
standard deviation quoted by Betoule et al. (2014). The effect of using the
new peculiar velocity corrections (with their newly derived associated
covariance matrix) while maintaining the value of CMB restframe redshift from
Betoule et al. (2014) results in a more modest shift, $\Omega_{m}=0.297\pm
0.089$ (case 2M++, old $z_{\rm\tiny{CMB}}$). When using both the new redshift
values and our newly derived peculiar velocity corrections in combination, we
obtain $\Omega_{m}=0.273\pm 0.090$. All these results do not use our new
colour-based selection effects corrections; once those are included, the
constraint on the matter density shifts back to a value close to the original
JLA analysis, namely $\Omega_{m}=0.290\pm 0.091$ (but notice the larger
uncertainty on our result).

In the right panel of Fig. 9, we observe a shift in the posterior towards
lower $\Omega_{\Lambda}$ and larger $\Omega_{m}$ when adding the systematic
covariance matrix (including our new peculiar velocity covariances) to the
statistical covariance matrix, as already noticed by Shariff et al. (2016).
Further adding the correction for colour selection effects shifts the
posterior only slightly.

Figure 9: Impact of our new treatment of peculiar velocities and colour
selection effects correction on the constraints on
$\Omega_{m},\Omega_{\Lambda}$ for the case of an isotropic universe. Left
panel: comparison of marginal posteriors $\Omega_{m},\Omega_{\Lambda}$ using
the old JLA peculiar velocity uncertainties and $z_{\rm\tiny{CMB}}$ (green) to
the ones obtained using our new, group-corrected values of the CMB restframe
redshifts (blue), and additionally replacing the JLA peculiar velocity
uncertainties with our newly derived ones (orange). Right panel: posterior
using only the statistical covariance matrix (green), adding the systematic
covariance matrix (including our new peculiar velocity covariances, blue) and
further adding colour selection effects corrections (orange). In this panel,
we adopt our newly derived, group-corrects CMB restframe redshifts.  Table 3:
Parameter constraints from 1D marginal posterior distributions from JLA data,
and outcome of Bayesian model comparison (last column) for the models and data
considered. The second column indicates the peculiar velocity corrections
adopted: ‘PSCz’ is the same treatment as in Betoule et al 2014; ‘2M++’
indicates the new treatment in this paper; ‘PSCz, new $z_{\rm\tiny{CMB}}$’
uses the peculiar velocity corrections of Betoule et al 2014 but our newly
derived CMB restframe redshift, while ‘2M++, old $z_{\rm\tiny{CMB}}$’ uses our
newly derived peculiar velocity corrections but the same $z_{\rm\tiny{CMB}}$
as Betoule et al 2014. The last column is the difference between the model’s
log-evidence and log-evidence of the $\Lambda$CDM-Isotropic or Cosmographic-
Isotropic case (with or without colour correction, as appropriate), i.e., the
log of the Bayes factor. A negative value indicates that the model under
consideration is disfavoured w.r.t. the isotropic universe. The $\star$ on
some of the data for $\Delta\ln{(Z)}$ are to indicate that the calculation
involved a comparison for a Cosmographic expansion with a restricted prior on
the deceleration parameter of $q_{0}\sim[-2,0)$ which limits it to an
accelerating universe only.

Model | $v_{\rm pec}^{{\rm SN}}$ | Selection | $\Omega_{m}$ | $\Omega_{\Lambda}$ | $l_{d}~{}\mathrm{[rad]}$ | $b_{d}~{}\mathrm{[rad]}$ | $|D_{\mu}|$ | Dipole | Quantity | $\Delta\ln{(Z)}$ |
---|---|---|---|---|---|---|---|---|---|---|---
correction | correction | scale, $S$ | modulated |
$\Lambda$CDM-Isotropic | PSCz | No | $0.306\pm 0.087$ | $0.562\pm 0.135$ | - | - | - | - | - |  |
$\Lambda$CDM-Isotropic | PSCz, new $z_{\rm\tiny{CMB}}$ | No | $0.253\pm 0.089$ | $0.460\pm 0.139$ | - | - | - | - | - |  |
$\Lambda$CDM-Isotropic | 2M++, old $z_{\rm\tiny{CMB}}$ | No | $0.297\pm 0.089$ | $0.547\pm 0.137$ | - | - | - | - | - |  |
$\Lambda$CDM-Isotropic | 2M++ | No | $0.273\pm 0.090$ | $0.508\pm 0.140$ | - | - | - | - | - | $0.0$ |
$\Lambda$CDM-Isotropic | 2M++ | Yes | $0.290\pm 0.091$ | $0.500\pm 0.141$ | - | - | - | - | - | $0.0$ |
CDM only-Isotropic | 2M++ | No | $0.037\pm 0.028$ | $0$ | - | - | - | - | - | $-5.65$ |
CDM only-Isotropic | 2M++ | Yes | $0.037\pm 0.028$ | $0$ | - | - | - | - | - | $-4.85$ |
$\Lambda$CDM-D | 2M++ | No | $0.263\pm 0.090$ | $0.493\pm 0.140$ | Unconstrained | Unconstrained | $<5.75\times 10^{-4}\,(95\%)$ | - | $\mu$ | $-8.68\pm 0.01$ |
$\Lambda$CDM-D | 2M++ | Yes | $0.285\pm 0.094$ | $0.490\pm 0.143$ | Unconstrained | Unconstrained | $<5.93\times 10^{-4}\,(95\%)$ | - | $\mu$ | $-6.81\pm 0.01$ |
$\Lambda$CDM-D-Exp | 2M++ | No | $0.269\pm 0.087$ | $0.499\pm 0.137$ | Unconstrained | Unconstrained | $<9.81\times 10^{-3}\,(95\%)$ | $<8.75\times 10^{-2}\,(95\%)$ | $\mu$ | $-5.27\pm 0.08$ |
$\Lambda$CDM-D-Exp | 2M++ | Yes | $0.286\pm 0.091$ | $0.487\pm 0.142$ | Unconstrained | Unconstrained | $<1.05\times 10^{-2}\,(95\%)$ | $<8.71\times 10^{-2}\,(95\%)$ | $\mu$ | $-3.49\pm 0.01$ |
$\downarrow$ | $\downarrow$ | $\downarrow$ | $q_{0}$ | $j_{0}-\Omega_{\kappa}$ | $\downarrow$ | $\downarrow$ | $\downarrow$ | $\downarrow$ | $\downarrow$ | $\downarrow$ |
Cosmographic-Isotropic | PSCz | No | $-0.352\pm 0.092$ | $0.114\pm 0.391$ | - | - | - | - | - |  |
Cosmographic-Isotropic | PSCz, new $z_{\rm\tiny{CMB}}$ | No | $-0.288\pm 0.095$ | $-0.132\pm 0.370$ | - | - | - | - | - |  |
Cosmographic-Isotropic | 2M++, old $z_{\rm\tiny{CMB}}$ | No | $-0.342\pm 0.093$ | $0.071\pm 0.390$ | - | - | - | - | - |  |
Cosmographic-Isotropic | 2M++ | No | $-0.320\pm 0.094$ | $-0.254\pm 0.383$ | - | - | - | - | - | $0.0$ |
Cosmographic-Isotropic | 2M++ | Yes | $-0.302\pm 0.090$ | $-0.023\pm 0.357$ | - | - | - | - | - | $0.0$ |
Cosmographic-($q_{0}=0$) | 2M++ | No | $0$ | $-0.953\pm 0.140$ | - | - | - | - | - | $-3.90\pm 0.06^{\star}$ |
Cosmographic-($q_{0}=0$) | 2M++ | Yes | $0$ | $-0.881\pm 0.134$ | - | - | - | - | - | $-3.28\pm 0.09^{\star}$ |
Cosmographic-($q_{0}\geq 0$) | 2M++ | No | $<0.067\times 10^{-2}\,(95\%)$ | $-0.922\pm 0.141$ | - | - | - | - | - | $-7.39\pm 0.05^{\star}$ |
Cosmographic-($q_{0}\geq 0$) | 2M++ | Yes | $<0.074\times 10^{-2}\,(95\%)$ | $-0.922\pm 0.144$ | - | - | - | - | - | $-7.00\pm 0.09^{\star}$ |
Cosmographic-D | 2M++ | No | $-0.307\pm 0.093$ | $-0.085\pm 0.375$ | Unconstrained | Unconstrained | $<6.43\times 10^{-4}\,(95\%)$ | - | $\mu$ | $-6.36\pm 0.02$ |
Cosmographic-D | 2M++ | Yes | $-0.291\pm 0.095$ | $-0.070\pm 0.371$ | Unconstrained | Unconstrained | $<6.32\times 10^{-4}\,(95\%)$ | - | $\mu$ | $-6.47\pm 0.02$ |
Cosmographic-D-Exp | 2M++ | No | $-0.307\pm 0.089$ | $-0.069\pm 0.353$ | Unconstrained | Unconstrained | $<1.00\times 10^{-2}\,(95\%)$ | $<8.68\times 10^{-2}\,(95\%)$ | $\mu$ | $-2.96\pm 0.05$ |
Cosmographic-D-Exp | 2M++ | Yes | $-0.293\pm 0.094$ | $-0.048\pm 0.362$ | Unconstrained | Unconstrained | $<9.63\times 10^{-3}\,(95\%)$ | $<8.94\times 10^{-2}\,(95\%)$ | $\mu$ | $-3.26\pm 0.01$ |
$\downarrow$ | $\downarrow$ | $\downarrow$ | $q_{0}$ | $j_{0}-\Omega_{\kappa}$ | $\downarrow$ | $\downarrow$ | $|D_{q}|$ | $\downarrow$ | $\downarrow$ | $\downarrow$ |
Cosmographic-D | 2M++ | No | $-0.305\pm 0.094$ | $-0.083\pm 0.376$ | Unconstrained | Unconstrained | $<6.06\times 10^{-2}\,(95\%)$ | - | $q_{0}$ | $-7.12\pm 0.02$ |
Cosmographic-D | 2M++ | Yes | $-0.286\pm 0.094$ | $-0.084\pm 0.367$ | Unconstrained | Unconstrained | $<6.29\times 10^{-2}\,(95\%)$ | - | $q_{0}$ | $-7.45\pm 0.08$ |
Cosmographic-D-Exp | 2M++ | No | $-0.316\pm 0.097$ | $-0.040\pm 0.370$ | Unconstrained | Unconstrained | $<22.46\,(95\%)$ | $<5.59\times 10^{-2}\,(95\%)$ | $q_{0}$ | $-1.22\pm 0.02$ |
Cosmographic-D-Exp | 2M++ | Yes | $-0.296\pm 0.093$ | $-0.043\pm 0.359$ | Unconstrained | Unconstrained | $<23.45\,95\%$ | $<6.00\times 10^{-2}\,(95\%)$ | $q_{0}$ | $-1.42\pm 0.02$ |

The second section of Table LABEL:table:final_jla_results presents our
constraints on the distance modulus dipole parameters, also comparing the
impact of using the colour selection effect correction (as indicated in the
third column), which is found to be quite minor on all constraints. Posterior
1D and 2D distributions are shown in Fig. 10 for the $\Lambda$CDM model and in
Fig. 11 for the Cosmographic expansion. The posterior distribution for the
dipole $|D_{\mu}|$ in the left panel of Fig. 10 for the $\Lambda$CDM model
with $F(z)=1$ peaks at 0, and we set a 1-tailed $95\%$ upper bound of
$|D_{\mu}|<5.93\times 10^{-4}$ (95.45% probability) – a factor of $\sim 2$
tighter than the limits derived by Lin et al. (2016a) from the same data,
namely $|D_{\mu}|<1.98\times 10^{-3}$. The dipole direction is correspondingly
unconstrained.

In the second case, shown in the right panel of Fig. 10, a scale function of
the form $F(z)=\exp{(-z/S)}$ is used to constrain the dipole to local region.
As expected from our simulations, we find a degeneracy between the dipole
parameter and its scale, whose effect is to degrade the upper limits on the
dipole amplitude to $|D_{\mu}|<1.05\times 10^{-2}$ for the $\Lambda$CDM case.
We can see $|D_{\mu}|$ is well explored again, right up to its prior edge
because of this degeneracy, which contributes to these weaker constraints. The
limits on the dipole scale $S$ are also weak, with the 1D marginal
distribution stretching all the way to the prior boundary ($S=0.10$), but
peaking near the lower prior boundary. The very weak preference for a non-zero
dipole (seen in the peak away from 0 in the 1D marginal distribution) could be
an indication of a residual effects of the bulk flow, which points broadly in
the same direction as the more prominent peak in the $l_{d}$ posterior
distribution. Such departures from perfect isotropy are weak, and not
dissimilar from what we observed in our isotropic universe simulations (Fig.
8, right panel) – hence they can be easily ascribed to the result of random
noise. With regards to constraints on cosmological parameters, with respect to
the $\Lambda$CDM-Isotropic case, we observe only a very mild shift in their
value as a consequence of the introduction of a potential dipole in the model,
$|\Delta\Omega_{M}|=0.005$ and $|\Delta\Omega_{\Lambda}|=0.006$. Both are
shifts of less than $0.1$ standard deviations of the posterior, and are
similar in scale for the case with $F(z)=1$. A similar result is seen in the
Cosmographic case with $|\Delta q_{0}|=0.009$ which is also quite minor
relative to the posterior standard deviation.

Figure 10: Posterior inference allowing for a distance modulus dipole from JLA
data, assuming the $\Lambda$CDM model with $F(z)=1$ (left) and
$F(z)=\exp{(-z/S)}$ (right), including colour-based selection effects
correction. Contours in 2D plots enclose joint 1, 2 and 3 $\sigma$ HPD
credible regions; the vertical dashed line in the 1D marginals for $D_{\mu}$
and $S$ delimits the 2$\sigma$ upper limit.

For the isotropic Cosmographic expansion, the posterior mean for the
deceleration parameter increases from $q_{0}=-0.352\pm 0.092$ when using the
old PSCz velocity corrections and $z_{\rm\tiny{CMB}}$ to $q_{0}=-0.302\pm
0.090$ when updating both velocity corrections and CMB restframe redshifts to
the new values we derived here (the error indicates the standard deviation of
the posterior, not the uncertainty on the mean). When introducing the
possibility of a dipole, the posterior mean for $q_{0}$ hovers around $-0.30$,
depending on the dipole model and whether we adopt the colour selection
effects corrections. These values are quite a bit larger than the expectation
under $\Lambda$CDM, namely $q_{0}=-0.55$, but not as large as the results
reported in Table 2 of Colin et al. (2019b) ($q_{0}=-0.157$ in our notation),
who ascribed the shift of the deceleration parameter towards 0 to evidence for
an anisotropic universe. Our analysis (with its improved treatments of
peculiar velocities and colour selection effects) disagrees with those
conclusions: the marginal posterior probability (obtained from our posterior
samples, not a Gaussian approximation to the posterior) for $q_{0}\geq 0$
(i.e., no acceleration) for the $F(z)=1$ case with selection effects
corrections and a distance modulus dipole is $9.3\times 10^{-4}$. We return on
this question from the point of view of Bayesian model comparison in section
4.2.

The constraints on the dipole parameters for the Cosmographic expansion model
are qualitatively similar to those presented for the $\Lambda$CDM model, as
shown in Fig. 11 and detailed in the central two sections of Table
LABEL:table:final_jla_results. For the $F(z)=1$ case, the posterior dipole
amplitude peaks at 0 and we set a $95\%$ upper limit $|D_{q_{0}}|<6.32\times
10^{-4}$. There is no significant evidence for a dipole moment in the
Cosmographic expansion case under our data and models. A qualitatively similar
picture holds for the $F(z)=\exp{(-z/S)}$ case, albeit with weaker limits on
the dipole amplitude owing to the degeneracy with the scale parameter $S$.

For a more direct comparison with the results of C19, we have also
investigated the same model as C19, where the dipole modulation is applied to
the deceleration parameter $q_{0}$, as in Eq. (15), rather than to the
distance modulus $\mu$. Differently from C19, we did not remove the bias
corrections to the magnitude, we kept the direction of the dipole free (as
opposed to being fixed in the CMB dipole direction), used our new peculiar
velocity corrections and CMB restframe redshifts (as opposed to heliocentric
redshifts) and applied our new colour selection effects. The resulting
posteriors are shown in Fig. 12 and constraints presented in the bottom
section of Table LABEL:table:final_jla_results.

We see qualitatively similar results to the case where we model the dipole on
$q_{0}$ in the Cosmographic expansion. For the $F(z)=1$ case, the posterior
dipole amplitude again peaks at $0$ and we set a $95\%$ upper limit
$|D_{q_{0}}|<6.29\times 10^{-2}$. A qualitatively similar picture holds for
the $F(z)=\exp{(-z/S)}$ case with a preference for a small $S$ value. This is
in contrast with Colin et al. (2019b), who saw a preference of a scale value
$S=0.0262$, a likely consequence of their removal of peculiar velocity
correction and the use of heliocentric redshifts.

Figure 11: Posterior inference allowing for a distance modulus dipole from JLA
data, assuming the Cosmographic expansion model with $F(z)=1$ (left) and
$F(z)=\exp{(-z/S)}$ (right), including colour-based selection effects
correction. Contours in 2D plots enclose joint 1, 2 and 3 $\sigma$ HPD
credible regions; the vertical dashed line in the 1D marginals for $D_{\mu}$
and $S$ delimits the 2$\sigma$ upper limit.

Figure 12: Posterior inference allowing for a distance modulus dipole from JLA
data, assuming the Cosmographic expansion model with $F(z)=1$ (left) and
$F(z)=\exp{(-z/S)}$ (right), including colour-based selection effects
correction. Contours in 2D plots enclose joint 1, 2 and 3 $\sigma$ HPD
credible regions; the vertical dashed line in the 1D marginals for $D_{q_{0}}$
and $S$ delimits the 2$\sigma$ upper limit. The dipole here is modelled
directly on the deceleration parameter rather than the distance modulus,
similar to C19. Table 4: Comparison of our results with previous searches for
a dipolar modulation. We only include comparable results that utilize the
“dipole fitting” approach (as opposed to the “hemisphere comparison” method).
Reported detections (at higher than $95\%$ significance) are highlighted in
boldface. The third column ‘$v_{\rm pec}$ corr’states whether peculiar
velocity corrections were adopted, and if so from which data set they were
derived. It should be noted here, previous significant detections involved
data that had no peculiar velocity corrections.

Reference | Data |  $v_{\rm pec}$ corr |  $l_{d}$ (deg) |  $b_{d}$ (deg) | Dipole amplitude | Quantity modulated
---|---|---|---|---|---|---
Cooke & Lynden-Bell (2010) | Union | None | $309$ | $43$ | $0.14\pm 0.12$ | $d_{L}$
Lin et al. (2016b) | Union2.1 | None | $171.8_{-42.0}^{+42.0}$ | $9.9_{-20.3}^{+20.3}$ | $0.160\pm 0.115$ | $\Omega_{m}$
Mariano & Perivolaropoulos (2012) | Union2 | None | $309_{-18}^{+18}$ | $-15.1_{-11.5}^{+11.5}$ | $(1.3\pm 0.6)\times 10^{-3}$ | $\mu$
Yang et al. (2014) | Union2.1 | None | $307.1_{-16.2}^{+16.2}$ | $-14.3_{-10.1}^{+10.1}$ | $(1.2\pm 0.5)\times 10^{-3}$ | $\mu$
Wang & Wang (2014) | Union2.1+GRB | None | $309.2_{-15.8}^{+15.8}$ | $-8.6_{-10.5}^{+10.5}$ | $(1.37\pm 0.57)\times 10^{-3}$ | $\mu$
Lin et al. (2016a) | JLA | PSCz | $316_{-110}^{+107}$ | $-5_{-60}^{+41}$ | $<1.98\times 10^{-3}\,(95\%)$ | $\mu$
Sun & Wang (2018a) | Union2.1 | None | $309_{-15.7}^{+15.5}$ | $-8.9_{-9.8}^{+11.2}$ | $(1.46\pm 0.56)\times 10^{-3}$ | $\mu$
Sun & Wang (2018a) | Constitution | None | $67.0_{-66.2}^{+66.5}$ | $-0.6_{-26.3}^{+25.2}$ | $(4.4\pm 5.0)\times 10^{-4}$ | $\mu$
Sun & Wang (2018a) | JLA | PSCz | Unconstrained | Unconstrained | Unconstrained | $\mu$
Sun & Wang (2018c) | Pantheon | 2M++ | $329_{-28}^{+101}$ | $37_{-21}^{+52}$ | $(3.7_{-3.7}^{+2,5})\times 10^{-4}\,(95\%)$ | $\mu$
Zhao et al. (2019) | Pantheon | 2M++ | $306.00_{-125.01}^{+82.95}$ | $-34.20_{-54.93}^{+16.82}$ | $<1.16\times 10^{-4}\,(95\%)$ | $\mu$
This work | JLA | 2M++ | Unconstrained | Unconstrained | $<5.93\times 10^{-4}(95\%)$ |  $\mu$, $F(z)=1$
This work | JLA | 2M++ | Unconstrained | Unconstrained | $<1.05\times 10^{-2}(95\%)$ |  $\mu$, $F(z)=\exp{(-z/S)}$
Colin et al. (2019b) | JLA | None |  $264.021$ (fixed) |  $48.253$ (fixed) | $-8.03_{-2.05}^{+2.05}$ |  $q_{0}$, $F(z)=\exp(-z/S)$
Rubin & Heitlauf (2020) | JLA | None |  $264.021$ (fixed) |  $48.253$ (fixed) | $-8.65^{+2.2}_{-2.6}$ |  $q_{0}$, $F(z)=\exp(-z/S)$
Rubin & Heitlauf (2020) | JLA | PCSz |  $264.021$ (fixed) |  $48.253$ (fixed) | $-1.1^{+3.2}_{-3.4}$ |  $q_{0}$, $F(z)=\exp(-z/S)$
This work | JLA | 2M++ | Unconstrained | Unconstrained |  $<6.32\times 10^{-2}$ ($95\%$) |  $q_{0}$, $F(z)=1$
This work | JLA | 2M++ | Unconstrained | Unconstrained |  $<13.68$ ($95\%$) |  $q_{0}$, $F(z)=\exp(-z/S)$

###  4.2 Bayesian Model Comparison

We compare the isotropic expansion model to the alternatives featuring a
dipole via Bayesian model comparison, and report the difference in the log of
the Bayesian evidence (i.e., the log of the Bayes factor) in Table
LABEL:table:final_jla_results:

$\Delta\ln(Z)=\ln B_{D}-\ln B_{I},$ (48)

where $B_{I}$ is the Bayesian evidence for the isotropic model (either
$\Lambda$CDM or Cosmographic expansion) and $B_{D}$ is the evidence for a
model featuring a dipole, with priors as in Table LABEL:table:priors. A value
of $\Delta\ln(Z)<0$ indicates a preference for the isotropic model. According
to the Jaynes’ scale for the strength of evidence, values of
$|\Delta\ln(Z)|=2.5(5.0)$ correspond to moderate (strong) evidence for one of
the models being compared. For equal prior probability for the models, the
quantity $\exp(\Delta\ln(Z))=B_{D}/B_{I}$ gives the posterior odds between the
isotropic and the dipole expansion models, which are approximately 12:1
(150:1) for moderate (strong) evidence (see Trotta (2008) for details on
Bayesian model comparison). When considering the numerical odds derived from
the Bayesian evidence, it is important to bear in mind that these can be very
sensitive to the choice of prior distribution, particularly the prior on the
parameters that only appears in the more complicated model (i.e., the model
featuring the dipole). While we believe that our choices of prior
distributions are well justified in Section 2.7, researchers that make other
choices for their prior distributions may compute odds ratios that differ from
those that we report.

In Table LABEL:table:priors we only carry out pair-wise model comparisons
between models that use the same data and same treatment of colour selection
effects, for comparing models with different data and/or assumptions about the
data generating process would be meaningless. The evidences (and associated
uncertainties) are estimated with PyMultinest, which was run with 400 live
points and an evidence tolerance of 0.5. We observe that the isotropic model
is favoured over all others: in the $\Lambda$CDM scenario, the constant dipole
model is disfavoured with odds ranging from 900:1 to almost 6000:1, depending
on the adoption of colour selection effects. The model with a dipole falling
off with redshift is also disfavoured, albeit with smaller odds ranging
between 32:1 and 194:1.

A similar pictures holds for comparison between models in the Cosmographic
expansion case. Here, the odds against the anistropic models are generally
smaller than in $\Lambda$CDM, owing to the smaller parameter space volume
ratio between posterior and prior, which control the strength of the Occam’s
razor effect in favour of the isotropic model. We also note that when
introducing a dipole scale parameter $S$, despite the larger number of free
parameters in this model w.r.t. the case where the dipole is constant in
redshift, the Bayes factor against it is actually smaller than the constant-
in-redshift dipole case. This can be explained by noting that the introduction
of $S$ as a free parameter leads to much less stringent limits on the dipole
amplitude because of the degeneracy explained above. Since $S$ itself is only
weakly constrained, the Occam’s razor effect for these two parameters is
weakened, leading to a weaker preference for the isotropic model. As a
consequence, the Cosmographic-D-Exp model where the dipole is on the
deceleration parameter is only very mildly disfavoured (odds smaller than 5:1)
w.r.t. the isotropic model – a consequence of the fact that the additional
parameters for this model cannot be strongly constrained, and hence the
posterior odds remain approximately equal to the prior odds.

As another case, we reproduce the setup used in C19, namely removing peculiar
velocity corrections entirely, using heliocentric redshifts (i.e., $z_{\rm
hel}$ instead of $\bar{z}$ in Eq. (34), a choice that imprints the dipole due
to the solar system’s motion onto the data, as pointed out by Rubin & Heitlauf
(2020)) and removing from the covariance matrix all uncertainties associated
with peculiar velocity corrections. In this setup, we compare the evidence for
an isotropic Cosmographic expansion with that of a dipolar modulation of the
form $F(z)=\exp{(-z/S)}$ either on the deceleration parameter, Eq. (15) (as in
C19), or on the distance modulus, Eq. (14). We adopt a Gaussian prior with a
standard deviation of $10^{\circ}$ on the dipole direction, centred on the
bulk flow direction from Boruah et al. (2020b), namely $l_{bf}=301^{\circ}\pm
4^{\circ},b_{bf}=0^{\circ}\pm 3^{\circ}$ (in excellent agreement with the
results of Said et al. (2020), obtained using the Fundamental Plane relation,
namely $l_{bf}=304^{\circ}\pm 4^{\circ},b_{bf}=1^{\circ}\pm 4^{\circ}$). When
the dipole is modelled on the deceleration parameter, as in C19, the
Cosmographic anisotropic model is still disfavoured with respect to the
isotropic one, with odds of approximately 17:1 ($\Delta\ln(Z)=-2.84\pm 0.08$).
Although the anisotropic model achieves a better quality of fit by absorbing
the dipole in the data, from an Occam’s razor persepective it remains
disfavoured due to its additional, unwarranted model complexity. We observe a
similar effect (if stronger) when the dipole is modelled instead on distance
modulus $\mu$, with odds of approximately 150:1 ($\Delta\ln(Z)=5.03\pm 0.03$)
in favour of the isotropic model. We can repeat this comparison in the
isotropic $\Lambda$CDM case and compare that with an anisotropic $\Lambda$CDM
model with $F(z)=\exp{(-z/S)}$, finding $\Delta\ln(Z)=-4.92\pm 0.15$, which
again favours the isotropic model with odds of 136:1.

Another interesting question is the strength of evidence in favour of an
accelerating universe in the isotropic expansion case. The Bayes factor
between the isotropic $\Lambda$CDM model and an isotropic model with no dark
energy ($\Omega_{\Lambda}=0$) disfavours the latter with odds in excess of
120:1 (including selection effects corrections). For completeness (and to
compare the above Bayesian model comparison results with a hypothesis testing
approach), we have also computed the log-likelihood difference for the best-
fit parameter values:

$\Delta\ln(L)=\ln\hat{L}_{\Omega_{\Lambda}=0}-\ln\hat{L}_{\Lambda\text{CDM}}$
(49)

where $\hat{L}_{\Lambda\text{CDM}}$ is the maximum likelihood value for the
$\Lambda$CDM model and $\hat{L}_{\Omega_{\Lambda}=0}$ is the maximum
likelihood value for a universe with no dark energy. We find
$\Delta{\ln(L)}=-6.35$ (with selection effects correction), which we translate
into a $p$-value using Chernoff theorem (as the hypothesis being tested,
$\Omega_{\Lambda}=0$, lies on the boundary of the allowed parameter space),
obtaining a $p$-value of $5.9\times 10^{-3}$ for the hypothesis that
$\Omega_{\Lambda}=0$.

In the phenomenological Cosmographic expansion setting, the accelerating
isotropic model, with a uniform prior on $q_{0}\in[-2,0)$ is favoured with
odds of approximately 26:1 when compared to a coasting universe, i.e.,
$q_{0}=0$. The accelerating universe is preferred with odds of almost $\sim$
1100:1 when compared with a decelerating model, i.e. one with a uniform prior
$q_{0}\in(0,1]$.

##  5 Summary and Conclusions

We have revisited the question of a dipolar anisotropy in the expansion of the
universe, and derived new constraints on a possible dipolar modulation from
SNIa data. Our approach builds on the Bayesian hierarchical model BAHAMAS,
which has been extended to include a new approximate correction for residual
colour-based selection effects. We have also upgraded the treatment of
peculiar velocities and host galaxy redshifts from the original JLA paper, by
adopting state-of-the-art flow models constrained using the 2M++ galaxy
catalogue. Finally, we have improved the treatment of both statistical and
systematic uncertainties pertaining to peculiar velocity corrections – the
dominant source of error for $z\lesssim 0.1$ SNIa, which are all-important for
a robust, accurate and precise measurement of anisotropy in the local
expansion.

We did not find any evidence for a deviation from isotropy, either in the
framework of $\Lambda$CDM or in phenomenological Cosmographic expansion. We
placed tight constraints on the amplitude of a possible dipole both in the
distance modulus and on the deceleration parameter. Our upper bounds are more
stringent by a factor of $\sim 2$ than the results previously obtained from
the same data sets with a comparable approach. We note that all previous
searches that have claimed a significant detection of anisotropy have
neglected peculiar velocity corrections.

We have used the framework of Bayesian model comparison to evaluate the Bayes
factor between models featuring a dipole and an isotropically expanding
universe (both in $\Lambda$CDM and in the Cosmographic expansion). We found
moderate to strong Bayesian evidence against an anisotropic expansion. We have
also evaluated the evidence in favour of acceleration, finding that a non-zero
cosmological constant is preferred, using JLA SNe alone, by odds of 120:1, a
result corroborated by a more traditional $p$-value approach based on a
frequentist hypothesis test, which rejects $\Omega_{\Lambda}=0$ with a
$p$-value of $5.9\times 10^{-3}$. In the Cosmographic expansion, a
decelerating universe is disfavoured with odds of almost $1100:1$ w.r.t. an
accelerating one. We conclude that the preferred model remains the
$\Lambda$CDM isotropically expanding model.

Our work does not address the important question of the so-called “$H_{0}$
tension” – the fact that the value of the Hubble-Lemaître constant obtained
from local distance ladder measurements shows a highly statistical significant
discrepancy with respect to the the smaller value obtained via CMB
anisotropies at high redshift (see, e.g., Di Valentino et al. (2021) for a
review). SNe type Ia, on their own, only provide a relative distance
measurement, given the exact degeneracy between $H_{0}$ and $M_{0}$. However,
in combination with calibrated distance indicators, such as Cepheid variables,
they deliver an absolute distance scale that is used to obtain tight
constraints on $H_{0}$, see Riess et al. (2021) for the latest results of the
SH0ES collaboration, and Freedman (2021), who finds that calibrating SNIa via
the tip of the Red Giant Branch (TRGB) instead, the tension with the CMB
anisotropies-determined value is no longer significant. The novel peculiar
velocity corrections presented here have potentially a bearing on this
important problem, in that a change in the estimated peculiar velocity implies
a horizontal shift of the SNIa in the Hubble diagram – potentially even
correlated across several SNe. This in turn could affect the linear fit that
gives, via the slope, the value of $H_{0}$. While we leave the investigation
of this effect to future work, we do provide our new corrections to the JLA
data publicly (see Data Availability section at the end), so that others
wishing to use it for studies of $H_{0}$ can do so.

In this work we adopted the JLA compilation since we were unable to use the
most recent Pantheon sample, owing to the full covariance matrix (including
systematics) not being publicly available. We plan to apply our new framework,
including up-to-date peculiar velocity corrections, to the Foundation sample
(Foley et al., 2018; Jones et al., 2019), together with the recent data
release of the Dark Energy Survey (DES Collaboration et al., 2018). In the
near future, we expect to be able to obtain even tighter constraints on
possible anisotropy in the expansion from upcoming, larger SNIa surveys like
the one that will be delivered by the Vera Rubin Observatory Legacy Survey of
Space and Time (LSST) (LSST Science Collaboration et al., 2009).

## Acknowledgements

RT was partially supported by STFC grants ST/P000762/1 and ST/T000791/1. WR
was supported by PhD scholarship by the Science and Technology Facilities
Council (STFC), UK and partially supported by a Quantitative Research grant
provided by G Research, which is gratefully acknowledged. MJH acknowledges
support from an NSERC Discovery Grant.

## Data Availability

A Python implementation of BAHAMAS for public used in this article is being
worked on and will be made available soon. The data products generated from
our analysis from the article is available from
https://zenodo.org/record/5854639

## References

  * Akrami et al. (2014) Akrami Y., Fantaye Y., Shafieloo A., Eriksen H. K., Hansen F. K., Banday A. J., Górski K. M., 2014, ApJ, 784, L42
  * Amanullah et al. (2010) Amanullah R., et al., 2010, ApJ, 716, 712
  * Andrade et al. (2018a) Andrade U., Bengaly C. A. P., Alcaniz J. S., Santos B., 2018a, Phys. Rev. D, 97, 083518
  * Andrade et al. (2018b) Andrade U., Bengaly C. A. P., Santos B., Alcaniz J. S., 2018b, Astrophys. J., 865, 119
  * Antoniou & Perivolaropoulos (2010) Antoniou I., Perivolaropoulos L., 2010, Journal of Cosmology and Astroparticle Physics, 2010, 012
  * Azzalini & Arellano-Valle (2012) Azzalini A., Arellano-Valle R. B., 2012, arXiv e-prints, p. arXiv:1203.2376
  * Bennett et al. (2013) Bennett C. L., et al., 2013, The Astrophysical Journal Supplement Series, 208, 20
  * Bernal et al. (2017) Bernal C., Cárdenas V. H., Motta V., 2017, Physics Letters B, 765, 163
  * Betoule et al. (2014) Betoule M., et al., 2014, A&A, 568, A22
  * Boruah et al. (2020a) Boruah S. S., Hudson M. J., Lavaux G., 2020a, arXiv e-prints, p. arXiv:2010.01119
  * Boruah et al. (2020b) Boruah S. S., Hudson M. J., Lavaux G., 2020b, Monthly Notices of the Royal Astronomical Society, 498, 2703
  * Branchini et al. (1999) Branchini E., et al., 1999, MNRAS, 308, 1
  * Brout & Scolnic (2021) Brout D., Scolnic D., 2021, ApJ, 909, 26
  * Brout et al. (2022) Brout D., et al., 2022, arXiv e-prints, p. arXiv:2202.04077
  * Buchner et al. (2014) Buchner J., et al., 2014, A&A, 564, A125
  * Cai & Tuo (2012) Cai R.-G., Tuo Z.-L., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 004
  * Calcino & Davis (2017) Calcino J., Davis T., 2017, Journal of Cosmology and Astroparticle Physics, 2017, 038
  * Campanelli et al. (2011) Campanelli L., Cea P., Fogli G. L., Marrone A., 2011, Phys. Rev. D, 83, 103503
  * Carr et al. (2021) Carr A., Davis T. M., Scolnic D., Said K., Brout D., Peterson E. R., Kessler R., 2021, The Pantheon+ Analysis: Improving the Redshifts and Peculiar Velocities of Type Ia Supernovae Used in Cosmological Analyses, doi:10.48550/ARXIV.2112.01471, https://arxiv.org/abs/2112.01471
  * Carrick et al. (2015) Carrick J., Turnbull S. J., Lavaux G., Hudson M. J., 2015, Monthly Notices of the Royal Astronomical Society, 450, 317
  * Chambers et al. (2016) Chambers K. C., et al., 2016, arXiv e-prints, p. arXiv:1612.05560
  * Childress et al. (2014) Childress M. J., Wolf C., Zahid H. J., 2014, Monthly Notices of the Royal Astronomical Society, 445, 1898
  * Colin et al. (2019a) Colin J., Mohayaee R., Rameez M., Sarkar S., 2019a, arXiv e-prints, p. arXiv:1912.04257
  * Colin et al. (2019b) Colin J., Mohayaee R., Rameez M., Sarkar S., 2019b, Astronomy & Astrophysics, 631, L13
  * Conley et al. (2011) Conley A., et al., 2011, ApJS, 192, 1
  * Cooke & Lynden-Bell (2010) Cooke R., Lynden-Bell D., 2010, MNRAS, 401, 1409
  * DES Collaboration et al. (2018) DES Collaboration et al., 2018, arXiv e-prints,
  * Dam et al. (2017) Dam L. H., Heinesen A., Wiltshire D. L., 2017, MNRAS, 472, 835
  * Davis et al. (1996) Davis M., Nusser A., Willick J. A., 1996, ApJ, 473, 22
  * Davis et al. (2011) Davis T. M., et al., 2011, The Astrophysical Journal, 741
  * Di Valentino et al. (2021) Di Valentino E., et al., 2021, Class. Quant. Grav., 38, 153001
  * Feroz & Hobson (2008) Feroz F., Hobson M. P., 2008, Monthly Notices of the Royal Astronomical Society, 384, 449
  * Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, Monthly Notices of the Royal Astronomical Society, 398, 1601
  * Feroz et al. (2013) Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2013, arXiv e-prints,
  * Foley et al. (2018) Foley R. J., et al., 2018, MNRAS, 475, 193
  * Freedman (2021) Freedman W. L., 2021, The Astrophysical Journal, 919, 16
  * Frieman et al. (2008) Frieman J. A., et al., 2008, AJ, 135, 338