Probing the independence within the dark sector in the fluid approximation

TL;DR

This paper investigates potential interactions between dark matter and dark energy in the universe using a non-parametric approach, finding no current evidence but forecasting future surveys could detect or rule out such interactions.

Contribution

It introduces a non-parametric reconstruction method for dark sector interactions and assesses the detectability with upcoming large-scale surveys.

Findings

No significant evidence for dark sector interaction with current data.

Future surveys like LSST and DESI could constrain interaction to 20% precision.

Detection of interaction would challenge the standard non-interacting dark sector model.

Abstract

The standard model of cosmology is based on two unknown dark components that are uncoupled from each other. In this paper we investigate whether there is evidence for an interaction between these components of cold dark matter (CDM) and dark energy (DE). In particular, we focus on a minimal extension and reconstruct the interaction history at low-redshifts non-parametrically using a variation of the commonly used principal component analysis. Although we focus on the interaction in the dark sector, any significant deviation from the standard model that changes the expansion history of the Universe, should leave imprints detectable by our analysis. Thus, detecting signatures of interaction could also be indicative of other non-standard phenomena even if they are not the results of the interaction. It is thus interesting to note that the results presented in this paper do not provide support for the interaction in the dark sector, although the uncertainty is still quite large. In so far as interaction is present but undetectable using current data, we show from a Fisher forecast that forthcoming LSST and DESI surveys will be able to constrain a DM-DE coupling at $20\%$ precision --- enough to falsify the non-interacting scenario, assuming the presence of a modest amount of interaction.

Tables

Table 1. Table 1 : Summary statistics of the compressed CMB likelihood [ 61 ] , describing the marginalised mean values and their 68 % percent 68 68\% confidence limits. The last three columns give the pairwise correlations.

			Correlation coefficients
	$\hat{X}$	$\sigma$	$ℛ$	${\ell }_A$	${\omega }_b$
$ℛ$	1.7382	0.0088	1.0	0.64	-0.75
${\ell }_A$	301.63	0.15	-	1.0	-0.55
${\omega }_b$	0.02262	0.00029	-	-	1.0

Equations

$${\nabla }_{\mu }{T}_{(A)}^{\mu \nu }=Q_{(A)}^{\nu },$$

$$Q_{(A)}^{\nu }=Q_{(A)}{u}_{(A)}^{\nu }+f_{(A)}^{\nu },g_{\mu \nu }{f}_{(A)}^{\nu }{u}_{(A)}^{\nu }=0,$$

$$Q^{\mu }=Q{u}_c^{\mu },$$

$$Q(a)=q(a)H(a){\rho }_X(a),$$

$${\dot{\rho }}_c+3H{\rho }_c=-Q,{\dot{\rho }}_X=Q,$$

$$q(a)=\sum _{i=1}^n{q}_i{T}_i(a),T_i(a)=\{\begin{array}{ll}1,& a_{i-1}\le a<a_i,\\ 0,& \text{otherwise}.\end{array}$$

$$H^2(a)/H_0^2={\Omega }_m(a)+{\Omega }_X(a)+{\Omega }_K(a),$$

$${\Omega }_X(a)$$

$${\Omega }_m(a)$$

$${\dot{\delta }}_m+\frac{1}{2}\dot{h}=(Q/{\rho }_m){\delta }_m,$$

$$\ddot{h}+2H\dot{h}=-8\pi G(\delta \rho +3\delta p),$$

$${\ddot{\delta }}_m+\left(2H-Q/{\rho }_m\right){\dot{\delta }}_m-\left[2HQ/{\rho }_m+\frac{d}{dt}(Q/{\rho }_m)\right]{\delta }_m=4\pi G(\delta \rho +3\delta p).$$

$$D^{\prime\prime}+\frac{1}{a}\left(3+\frac{d\ln H}{d\ln a}-\Gamma\right)D^{\prime}=\frac{1}{a^{2}}\left[\frac{3}{2}\Omega_{m}(a)+\frac{1}{H}\frac{d}{d\ln a}\big{(}\Gamma H\big{)}+2\Gamma\right]D,$$

$$d_{z}(z)\equiv r_{s}(z_{\mathrm{d}})/D_{V}(z),\qquad D_{V}(z)\equiv\big{[}(1+z)^{2}d_{A}^{2}(z)\times cz/H(z)\big{]}^{1/3},$$

$$\mu (z)=m-M=5{\log }_{10}\left[d_L(z)/\text{10 pc}\right],$$

$$\mu =m_B^{\ast }-M+\alpha x_1-\beta c,$$

$$\pi (\mathbf{q})=\frac{1}{(2\pi {)}^{n/2}\sqrt{\det {\mathbf{C}}_{\pi }}}\exp \left[-\frac{1}{2}(\mathbf{q}-{\mathbf{q}}^{\mathrm{fid}}{)}^{\mathrm{T}}{\mathbf{C}}_{\pi }^{-1}(\mathbf{q}-{\mathbf{q}}^{\mathrm{fid}})\right],$$

$$q(a)\to q_i={\int }_0^{1}da{W}_i(a)q(a),$$

$$({\mathbf{C}}_{\pi }{)}_{ij}$$

$$={\int }_0^{1}da{W}_i(a){\int }_0^{1}d{a}^{\prime }{W}_j(a^{\prime })⟨(q(a)-q^{\mathrm{fid}}(a))(q(a^{\prime })-q^{\mathrm{fid}}(a^{\prime }))⟩.$$

$$\xi(a,a^{\prime})\equiv\left\langle\big{(}q(a)-q^{\text{fid}}(a)\big{)}\big{(}q(a^{\prime})-q^{\text{fid}}(a^{\prime})\big{)}\right\rangle.$$

$$\xi (|a-a^{\prime }|)=\frac{{\xi }_0}{1+(|a-a^{\prime }|/a_c{)}^2}.$$

$${\mathbf{F}}_{\pi }{\mathbf{v}}_i={\lambda }_i{\mathbf{F}}_p{\mathbf{v}}_i,$$

$$q(a)=\sum _{i=1}^n{\alpha }_i{e}_i(a),e_i(a)=\sum _{j=1}^n{A}_{ij}{T}_j(a),$$

$$C\equiv-2\big{(}D_{\mathrm{KL}}[p,\pi]-\widehat{D}_{\mathrm{KL}}\big{)},$$

$$D_{\mathrm{KL}}[p,\pi ]=\int d\theta p(\theta |D)\ln \left(\frac{p(\theta |D)}{\pi (\theta )}\right),$$

$$C=n-\mathrm{tr}({\mathbf{F}}_p^{-1}{\mathbf{F}}_{\pi })=n-\sum _{i=1}^n{\lambda }_i=n-\sum _{i=1}^n\frac{1}{(S/N{)}_i}.$$

$$\theta =\{{\Omega }_{m0},{\Omega }_{b0},{\sigma }_{80},q_1,q_2,\dots ,q_{20},\alpha ,\beta ,x_{\ast },c_{\ast },M_0,R_x,R_c,R_M\}.$$

$$\mathrm{MSE}=\sum _{i=1}^n{\left(q_i^{(m)}-{\bar{q}}_i\right)}^2+\sum _{i=1}^n{\left({\sigma }_i^{(m)}\right)}^2,$$

$$\sigma^{(m)}_{i}=\bigg{[}\sum_{j=1}^{m}e^{2}_{j}(a_{i})\,\sigma^{2}_{\alpha_{j}}\bigg{]}^{1/2}$$

$$p[q(a)]\equiv p(q_1,q_2,\dots ,q_n)=p({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n).$$

Full text

1 Introduction

The $\Lambda$ Cold Dark Matter ($\Lambda$CDM) model is now firmly established as the standard paradigm of cosmology, having fitted a wide range of observations [1]. It is, nevertheless, a phenomenological model and thereby provides no explanation for non-baryonic CDM, and the cosmological constant that appears to drive cosmic acceleration [2, 3]. Although these two ingredients make up $95\%$ of the matter-energy content of the Universe, little is known about their non-gravitational nature.

In the standard model it is assumed that (i) CDM is pressureless and dilutes with the cosmic expansion as $\rho_{c}\propto a^{-3}$ and (ii) DE has negative pressure, and is undiluting with an equation of state $w=-1$, i.e. a cosmological constant. In the development of the concordance cosmology the need for these two components arose at different times, from different lines of evidence, both astrophysical and cosmological. They are both considered independent of each other with different evolutions, and thus it is interesting to note that the coincidence problem — i.e. the observation that DM and DE have comparable densities only recently when for most of the lifetime of the Universe they were different — occurs roughly when the validity of the fluid approximation might be questioned.

The fluid approximation, i.e. the assumption that the cosmic density field behaves as an ideal fluid, is used at both early and late times of cosmological evolution. While this is a reasonable assumption in the former regime, when gravitational structures had yet to form and matter and energy existed as a thermal bath of particles, it might be questioned whether the fluid approximation still holds in the late Universe [4]. From the primordial epoch to the epoch preceding the first gravitational structures the evolution of the Universe could simply be described by thermal physics. By contrast, the late Universe we observe today is significantly more complicated, being composed of a complex hierarchy of nonlinear gravitational structures. Tracers of the cosmic density are no longer simply parcels of fluid particles but are instead galaxies that follow the underlying density field in nontrivial ways [5].

Models that involve interacting dark energy have long been studied as alternatives to $\Lambda$ (see the reviews [8, 6, 7] and references therein). Early work centred on quintessence models of DE, in which the scalar field is coupled to either matter [9, 10] or only dark matter [11, 12, 13]. If the interaction is contained in the dark sector, however, a particularly appealing feature is that it can provide an explanation of the coincidence problem [11, 10]. Recent work though has largely shifted to the question of whether DE is time-varying (i.e. with a possible time-dependent DE equation of state) and models now tested against data are typically based on the fluid picture of DE [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. In all of these models, however, the fundamental mechanism giving rise to the coupling is generic with only the effects of interaction on linear scales usually studied (though $N$-body simulations have recently begun investigating nonlinear scales [35, 36]). While the detection of interaction may be physical in nature it could also be argued to signal a breakdown of the fluid approximation. Even if the fluid approximation breaks down we might still expect a fluid-like evolution, just one that is different from the usual scalings.

These models can also be motivated by the fact that gravitational probes are sensitive only to the total energy-momentum tensor $T_{\mu\nu}$ [37], with the splitting of $T_{\mu\nu}$ into different constituents typically based on physical considerations and the strength of gravity insensitive to any coupling. At late times it is reasonable to think ordinary matter and radiation (baryons, neutrinos, photons etc) in the cosmic fluid do not couple. Moreover, interactions between Standard Model particles and the dark sector are strongly constrained by experimental data. However, without a fundamental theory behind the dark sector there is a priori no reason to split into non-interacting DM and DE.

Recently, such models have seen renewed interest as a possible solution to tensions in the measured values of $H_{0}$ and $\sigma_{8}$ [27, 30, 31, 32, 38, 34]. Deviations from the $\Lambda$CDM scenario hint at a breakdown of the usual assumptions of the dark sector, which could be interpreted in several different ways. In particular, if we take DM and DE to be as yet undetected particles or fields then detection of interaction is to be understood at face value, i.e. physical in nature. Alternatively, if we take the dark sector to be phenomenological artefacts required for concordance with observations, then it might call into question the reality of the dark components. It is also timely to revisit some of the assumptions of the $\Lambda$CDM model given that observational cosmology is poised to see an influx of data from next generation experiments.

Without guidance from theoretical arguments, all parametrisations of the interaction investigated so far are necessarily phenomenological. Naturally, simple forms of the interaction have largely been pursued [6, 16, 19, 29, 31, 33]. Alternatively, we can take a less rigid approach by directly reconstructing the quantity of interest from, e.g. a set of basis functions or Padé approximants. However, this approach of introducing a large number of degrees of freedom generally suffers from parameter degeneracies leading to slow convergence using standard parameter estimation methods or weak constraints. In this work we will take a more data-driven approach. We use a generalised principal component analysis, allowing us to reconstruct the interaction history from a statistically decorrelated basis of eigenfunctions. The main advantage is that we are able to constrain only those features that are actually being probed, and thus let the data decide on the best functional form for the interaction.

The plan of this paper is as follows. In Section 2 we review aspects of interacting models, discuss the relevant theory, and set out the model to be analysed. In Section 3 we describe the data, statistical methods and tools used in the analysis. Section 4 presents the results; in Section 5 we study the detectability of interaction in upcoming surveys; finally, in Section 6 we summarise our main findings.

2 Preliminaries

We consider an energy-momentum tensor consisting of multiple fluids, labelled $A$. Typically it is assumed that each fluid species $A$ satisfies its own energy-momentum conservation equation, $\nabla_{\mu}T^{\mu\nu}_{(A)}=0$. In general, if we allow the transfer of energy-momentum between species then

[TABLE]

where $Q^{\nu}_{(A)}$ is the covariant interaction of species $A$. The conservation of the total energy-momentum $T^{\mu\nu}=\sum_{A}T^{\mu\nu}_{(A)}$ implies the balance condition $\sum_{A}Q^{\nu}_{(A)}=0$. Separating the energy and momentum part by decomposing $Q^{\nu}_{(A)}$ relative to the fluid 4-velocity $u^{\mu}_{(A)}$ of each fluid component we can write

[TABLE]

where $Q_{(A)}$ is the rate of energy transfer, $f^{\nu}_{(A)}$ the rate of momentum transfer, $u^{\nu}_{(A)}$ the 4-velocity, and $g_{\mu\nu}$ is the metric tensor. It can be observed that $Q^{\mu}$ does not appear in Einstein’s equations as they depend on $T_{\mu\nu}$ and not its derivative. The interaction term enters through the fluid equations only, which are in general modified from there usual forms.

As we consider only interaction between CDM ($c$) and DE ($X$) we have $Q^{\nu}:=Q_{c}^{\nu}=-Q_{X}^{\nu}$. Following [18] the form of the covariant interaction we assume to be

[TABLE]

where $Q$ is time-dependent only and $u^{\mu}_{c}$ is the 4-velocity of CDM. In this simple model there is no net momentum transfer in the rest frame of DM; any transfer that takes place is along the geodesic flow of DM. Consequently, in the synchronous gauge the peculiar velocity of DM (and also baryons) vanishes and we have $u_{c}^{\mu}=(1,0,0,0)$. Moreover, as there is no momentum transfer, no spatial gradients arise in the density of DE and we have $\delta\rho_{X}=0$, i.e. DE is spatially homogeneous. This model is known as the geodesic interaction model and is among the simplest interacting extension to $\Lambda$CDM. In this model the fluid equations retain their usual non-interacting forms and it is the interaction model we consider in this work.

2.1 Model specifications

The form of the function $Q$ in (2.2) we assume to be of the following form

[TABLE]

where $q(a)$ is the dimensionless interaction history, $H(a)$ is the Hubble constant and $\rho_{X}(a)$ is the DE density. The chosen form for $Q(a)$ is for convenience only; any arbitrariness of $Q(a)$ is absorbed into $q(a)$. Note however that it does not necessarily lead to solutions that remain physical into the future. For instance, if $q>0$ then DM decays unbounded at a rate proportional to the DE density, eventually going negative. Therefore, we consider (2.3) an ansatz valid for the late-epoch that we focus on in this work.

Given that the matter density and DE density are approximately equal in the recent past we expect that any interaction will be greatest at low-redshifts. This is generally realised by assuming a logistic-like interaction parametrisation in which $Q$ only becomes appreciable at late-times when $\rho_{X}\simeq\rho_{c}$.

The convention we use here is that positive values of $Q$ (or $q$) gives a universe in which DM decays to DE, while negative values gives the reverse behaviour. For a fixed $\Omega_{m0}$, models with $q>0$ will have a greater fraction of matter through all epochs. In such a case, the growth of structure is enhanced relative to $\Lambda$CDM, as the universe is more matter dominated than in the non-interacting case. For $q<0$ the growth of structure is suppressed relative to $\Lambda$CDM as the onset of DE occurs earlier.

Assuming a DE equation of state $w_{X}=-1$, the coupled continuity equations becomes

[TABLE]

where overdots denote differentiation with respect to cosmic time. The cosmological constant is of course recovered when $Q=0$ so $\rho_{X}\propto\Lambda$. All other fluid components evolve in the standard, non-interacting way.

The issue of how to parametrise $q$ can be likened to that of determining the DE equation of state. In the absence of any plausible $q$ from theory we will take a model-independent approach and reconstruct it directly from data. We divide $q$ up into $n$ bins and constrain the amplitudes $q_{i}$ of each bin. The bins are chosen to be uniformly spaced in scale factor $a$, with edges $a_{0}<a_{1}<a_{2}<\ldots<a_{n}$. The $i^{\mathrm{th}}$ bin spans the interval $[a_{i-1},a_{i})$ and we set $a_{n}=a_{\mathrm{max}}=1$. Since $q(a)$ has a piecewise constant amplitude in each bin we represent it as

[TABLE]

The amplitudes $q_{1},q_{2},\ldots,q_{n}$ are dimensionless parameters characterising the interaction strength. We set $q(a)=0$ outside the binning range $[a_{0},a_{n}]$. The low-redshift window chosen is motivated by the fact that the onset of cosmic acceleration occurs in the recent past.

The Friedmann equation is given by

[TABLE]

where $H_{0}$ is the present-day value of the Hubble constant, $\Omega_{K}(a)=\Omega_{K0}a^{-2}$ with $\Omega_{K0}$ the spatial curvature parameter, and $\Omega_{m}(a)=\Omega_{b}(a)+\Omega_{c}(a)$ together with $\Omega_{X}(a)$ have modified time-dependence given by solving (2.4). With the specific form given by (2.3) we find

[TABLE]

where $a\in[a_{j-1},a_{j})$, $\Omega_{m0}$ is the present matter density parameter, and $\Omega_{X0}$ is the present DE density parameter. If $a<a_{\mathrm{min}}=a_{0}$ then the densities recover their usual forms, $\rho_{c}\propto a^{-3}$ and $\rho_{X}=\mathrm{const}$. The binning strategy is chosen to effectively impose standard $\Lambda$CDM evolution, up until the onset of cosmic acceleration at late-times where we are most interested.

At the level of perturbations the presence of interaction modifies the continuity equation to allow an exchange of energy between fluid species. For the total matter fluctuation $\delta_{m}$ we have, in the synchronous gauge,

[TABLE]

where $\rho_{m}=\rho_{b}+\rho_{c}$ is the total matter density in the background. From Einstein’s equations, the metric perturbation $h$ satisfies

[TABLE]

where $G$ is the gravitational constant, $\delta\rho=\sum_{A}\delta\rho_{(A)}$ and $\delta p=\sum_{A}\delta p_{(A)}$ are the total density and pressure fluctuations, respectively. An ordinary, second-order differential equation can be obtained from (2.7) and (2.8) that is closed in the perturbation variable $\delta_{m}$:

[TABLE]

At late-times the effect of radiation is negligible so we neglect it and treat the total density perturbations as composed of baryons, CDM and DE. In comoving synchronous gauge it can be shown that $\delta\rho_{X}=0$, i.e. dark energy is spatially homogeneous and nonclustering [24]. Equation (2.7) implies the amplitudes of fluctuations evolve uniformly. Thus we can write $\delta_{m}(\bm{x},a)=D(a)\delta_{m}(\bm{x},a=1)$, where $D(a)$ is the growing mode. With a change of independent variable $t\to a$ the above equation becomes

[TABLE]

where primes denotes differentiation with respect to scale factor, and we have defined the dimensionless function $\Gamma\equiv Q/(\rho_{m}H)$ and $\Omega_{m}(a)\equiv 8\pi G\rho_{m}(a)/3H^{2}(a)$. Note we have suppressed all dependence on the scale factor, except for $\Omega_{m}(a)$, to avoid confusion with the matter density parameter $\Omega_{m0}$. We solve this equation numerically with the initial conditions $D(a_{\mathrm{init}})=a_{\mathrm{init}}$ and $D^{\prime}(a_{\mathrm{init}})=1$, taking $a_{\mathrm{init}}=0.03$ (or $z\simeq 30$), that is we begin integration at a time deep in the matter-dominated epoch when $D\sim a$. Having then solved (2.9) we normalise the growth factor to unity today. The quantity of interest is not $D(a)$, which is not observable, but the growth rate $f\equiv d\ln D/d\ln a$. Without solving this equation we can understand the effect of interaction qualitatively. If $Q>0$ there are two competing effects: on the one hand a conversion of DM to DE produces a faster cosmic expansion and a lower matter fraction but on the other we have $\Gamma>0$ driving the dissipation down and reinforcing the source term. In effect this means we can always compensate a high $\Omega_{m0}$ with a low $q$. Clearly, $D(a)$ no longer just depends on the cosmic expansion $H(a)$, but also $Q(a)$. This implies that a detection of non-zero $Q$ can be translated to mean a violation of a consistency relation if the assumptions of $\Lambda$CDM do not hold [40, 41, 42].

3 Methodology

3.1 Data sets

In this section we describe the data used and emphasize the physics they probe.

Baryon Acoustic Oscillations.

Galaxy surveys exhibit enhancements at a certain length scale in the clustering of matter, due to Baryon Acoustic Oscillations (BAO). From the anisotropic 2-point correlation function BAO surveys commonly report the distilled quantity

[TABLE]

where $z$ is an effective redshift found by taking a weighted average of galaxies in a given slice, $z_{\mathrm{d}}$ is the redshift at the baryon drag epoch, $d_{A}(z)=(1+z)^{-1}\int^{z}_{0}dz^{\prime}/H(z^{\prime})$ is the angular diameter distance for a spatially flat FLRW metric, and for a sound speed $c_{s}(z)$ in the baryon-photon fluid $r_{s}(z)=\int^{\infty}_{z}dz^{\prime}\,c_{s}(z^{\prime})/H(z^{\prime})$ is the comoving sound horizon. For the particular interacting model being considered the sound speed remains unchanged from $\Lambda$CDM (see Appendix A).

We use data from the Six-degree Field Galaxy Survey (6dFGS) [43], the Sloan Digital Sky Survey (SDSS) DR7 Main Galaxy Sample [44], the LOWZ and CMASS galaxy samples of the Baryon Oscillation Spectroscopic Survey (BOSS) DR12 [45], and the SDSS Luminous Red Galaxies (LRG) [46]. We also use the three correlated measurements reported by WiggleZ with the provided covariances [47]. The data comprise eight measurements at different effective redshifts of the distilled parameter given variously as $d_{z}$, its reciprocal $D_{V}(z)/r_{s}(z_{\mathrm{d}})$, or sometimes normalised to a fiducial cosmology.

Redshift Space Distortions.

The growth of structure in the Universe depends on its energy contents through its effect on background expansion (and importantly also interaction). The competition between cosmic expansion and the tendency for inhomogeneous regions to be further enhanced by gravity gives a useful dynamical probe at linear scales, particularly of DE and modified gravity. Galaxy surveys exploit the anisotropies induced on the power spectrum from redshift space distortions (RSD) to measure the normalised growth rate $f\sigma_{8}$, where $\sigma_{8}(z)$ is the root-mean-square of the amplitude of matter fluctuation averaged in a spherical volume of radius $8h^{-1}\,\textrm{Mpc}$. In the linear regime $\sigma_{8}(z)$ simply scales with $D(a)$ so $\sigma_{8}(z)=D(a)\sigma_{80}$, where $\sigma_{80}$ is its present-day value and can be treated as a free parameter.

In this analysis we use the same data as compiled by Planck (2018) [1], which consists of measurements from 6dFGS [70, 71], SDSS MGS [72], SDSS DR7 LRG [73], GAMA [74], BOSS DR12 [75], WiggleZ [76], VIPERS [77], FastSound [78], and BOSS DR14 quasars [79].

Type Ia Supernovae.

Type Ia supernovae (SNe Ia) are standardisable candles that can be used to probe the expansion history. The distance to SNe Ia is given by the distance modulus defined as

[TABLE]

where $m$ is the apparent magnitude, $M$ is the absolute magnitude, and $d_{L}(z)=(1+z)^{2}d_{A}(z)$ is the luminosity distance.

The cosmological observable is the distance moduli $\mu(z)$, which is to be compared to the measured value given by the Tripp relation

[TABLE]

with the peak apparent magnitude $m_{B}^{*}$ and intrinsic magnitude $M$ are given in $B$-band. The two additional terms — known as Phillips corrections — are the time stretch parameter ${x}_{1}$ and the colour correction parameter at maximum brightness ${c}$ [50, 51]. The stretch and colour coefficients, $\alpha$ and $\beta$, are nuisance parameters, i.e. parameters to be fitted simultaneously with the cosmological parameters.

The data for each SN Ia are the parameters $\hat{m}^{*}_{B}$, $\hat{x}_{1}$ and $\hat{c}$ as well as the heliocentric redshift $\hat{z}$ produced using the SALT2 light-curve fitting procedure [52]. We use the Joint Light-curve Analysis (JLA) catalogue of $740$ spectroscopically confirmed SNe Ia with redshifts ranging from $z=0.01$ to 1.3 [48].111http://supernovae.in2p3.fr/sdss_snls_jla/ Though the more recent Pantheon sample [49] is larger, with over 1000 SNe Ia, it does not yet include ancillary data products we need for the statistical method used here (described in Section 3.2).

Cosmic Chronometers.

The basic idea of the cosmic chronometers (CC) as a probe is to directly measure the cosmic expansion history by using the kinematic form $H(z)=-(1+z)^{-1}dz/dt$ (assuming a Friedmann-Lemaître-Robertson-Walker (FLRW) background redshift $1+z=1/a$). In principle, because it is not an integrated quantity like $d_{L}(z)$ it provides greater sensitivity to cosmological parameters. The difficulty however is in obtaining accurate estimates of the differential changes in redshift. At present, measurements are at about $6\%$ precision and limited by systematics. Observations of the age differences of passively evolving old, elliptical galaxies are typically serve as the “standard clocks” of this method. This analysis makes use of 26 measurements of $H(z)$ with redshifts ranging from $z\simeq 0.01$ to $z\simeq 2$ [62, 63, 64, 65].

Cosmic Microwave Background.

To complement the low-redshift probes we also include data from Planck. Although the CMB is mainly a probe of the early Universe (when dark energy was negligible), it does contain some valuable geometric information through the angular scale of the sound horizon [56]. At the background level, the CMB data provides a precise determination of the distance to last scattering.

We include Planck data using the compressed CMB likelihood method that considers a handful of parameters (thought of as observables) that summarise key features of the CMB power spectrum [57, 58, 59]. The data consists of the following: (i) the CMB shift parameter $\mathcal{R}\equiv\sqrt{\Omega_{m}H_{0}^{2}}D_{A}(z_{*})$, where $D_{A}(z)\equiv(1+z)d_{A}(z)$ is the comoving angular diameter distance evaluated at the redshift of last scattering $z_{*}$; (ii) the angular scale of the sound horizon at last scattering, $\ell_{A}\equiv\pi D_{A}(z_{*})/r_{s}(z_{*})=\pi/\theta_{*}$, where $r_{s}(z_{*})$ is the comoving sound horizon and $\theta_{*}=r_{s}(z_{*})/D_{A}(z_{*})$ is the angular size of the sound horizon; (iii) the physical baryon density $\omega_{b}$. With a precision of $\lesssim 1\%$, these three quantities are among the most precisely determined by Planck. As they summarise key geometric features of the CMB angular power spectrum they are sometimes referred to as the CMB distance prior. The observable $\mathcal{R}$ determines the distance to the last scattering surface independent of $H_{0}$, $\ell_{A}$ is closely related to the position of the first acoustic peak, and $\omega_{b}$ sets the relative heights of odd to even peaks. When combined with other data sets there is no significant loss of information using the compressed likelihood versus the full likelihood [59]. Although they are not strictly cosmology-free measurements, but rather constrained quantities obtained from a CMB analysis assuming a given model, they can be considered early Universe observables independent of the (late-)DE model assumed [61]. We use the Planck 2015 data release [61] of the lensing amplitude marginalised-compressed likelihood.222It is well-known that combinations of $\Omega_{\Lambda}h^{2}$ and $\Omega_{K0}h^{2}$ that give the same $\mathcal{R}$ will produce a near identical CMB spectrum at high multipoles [60]. However, with constraints from CMB lensing this degeneracy is broken. This likelihood is summarised in Table 1.

3.2 Statistical modelling

The parameters are sampled from the posterior formed from the joint likelihood of BAO, CMB, CC, RSD, and SNe Ia. For BAO, CC, and RSD we take the likelihoods to be Gaussian distributed in the data. As the Planck joint posterior distribution of the data $\hat{\bm{X}}=(\mathcal{R},\ell_{A},\omega_{b})^{\mathrm{T}}$ is near-Gaussian we also take the CMB likelihood to be Gaussian with mean given by the second column of Table 1 and covariance matrix $\mathbf{\Sigma}_{\mathrm{CMB}}=\mathbf{\Sigma}^{1/2}{\mathbf{D}}\mathbf{\Sigma}^{1/2}$, with $\mathbf{\Sigma}^{1/2}\equiv\mathrm{diag}(\sigma_{\mathcal{R}},\sigma_{\ell_{A}},\sigma_{\omega_{b}})$ and ${\mathbf{D}}$ the correlation matrix.

The SN data consists of standardisation outputs from the SALT2 template and as a result the regression model demands a more principled approach. We adopt a recent Bayesian hierarchical approach in which the dependencies are constructed within a probabilistic framework and has been shown to also deliver tighter constraints on parameters [54, 53].

In the hierarchical approach we introduce the latent variables $\bm{M}$, $\bm{x}_{1}$, $\bm{c}$ and $\bm{z}$ (represented as vectors) as the true variables that we do not observe. The SN likelihood involves additional hyperparameters that describe the distributions of latent or unobserved variables relating to the colour, stretch, and absolute magnitude. These distributions are taken to be Gaussian with means and standard deviations included in the set of free parameters $\theta$. The likelihood of the SALT2 outputs has the probabilistic form $\mathcal{L}_{\mathrm{SN}}({\theta};D)\equiv p(\hat{\bm{z}},\hat{\bm{m}}_{B}^{*},\hat{\bm{x}}_{1},\hat{\bm{c}}\,|\,{\theta})$. The set of parameters ${\theta}$ includes the nuisance parameters $\alpha$ and $\beta$. The data $D$ includes $\hat{\bm{m}}_{B}^{*}$, $\hat{\bm{x}}_{1}$, $\hat{\bm{c}}$ and $\hat{\bm{z}}$. The details of the marginalised likelihood used in this analysis are given in Appendix B.

To summarise, our regression model consists of the base cosmological parameters $\Omega_{m0}$, $\Omega_{b0}$, $\sigma_{80}$ the reconstruction parameters $q_{1},\ldots,q_{n}$, the SN nuisance parameters $\alpha$, $\beta$ and additional SN parameters $M_{0},x_{*},c_{*},R_{M}^{2},R_{x}^{2},R_{c}^{2}$ that describe the Gaussian priors of the latent variables.

3.3 Priors

3.3.1 Smoothing priors on reconstruction

Given the emphasis of this analysis on the parameters $q_{1},q_{2},\ldots,q_{n}$ it is important that appropriate priors are chosen. Typically flat priors are used to allow the inference to be driven by the data and provided the likelihood is informative and supported by the prior this is usually a reasonable choice. However in reconstruction there are necessarily a large number of degrees of freedom and these flat priors become informative relative to the likelihood. For sparse data it is unlikely that all parameters can be constrained and, because neighbouring bins are uncorrelated, often leads to a noisy reconstruction. Increasing the number of bins only introduces more unconstrained degrees of freedom and as a result the posterior is multimodal and convergence of MCMC methods is slow. This is one of the basic problems that reconstruction methods face. One popular method is to add a regularisation or penalty term to the $\chi^{2}$, designed in a way to favour smooth reconstructions. The method we adopt here instead incorporates this information into the prior [66]. The prior is assumed to be of Gaussian form

[TABLE]

where $\bm{q}=(q_{1},q_{2},\ldots,q_{n})^{\mathrm{T}}$ and similarly for the fiducial model $\bm{q}^{\mathrm{fid}}$. The off-diagonal entries of $\mathbf{C}_{\pi}$ can be specified so as to favour smoothness, ensuring that the $q_{i}$s do not change abruptly between bins. Since the discretisation of the unknown $q(a)$ is achieved by a simple average over each bin we can write

[TABLE]

where $W_{i}(a)=T_{i}(a)/(a_{i}-a_{i-1})$ is the normalised top-hat function. With this the covariance matrix can be written in component form as

[TABLE]

Thus covariances between bins are encoded in temporal correlations of $q(a)$ given by the two-point correlation function:

[TABLE]

Here $\xi(a,a^{\prime})$ is a function that we are free to specify. As a matter of convenience we assert that correlations are invariant under time translations and reversals, which implies that the arguments of $\xi$ depends only on the magnitude of the difference between $a$ and $a^{\prime}$, i.e. $\xi(a,a^{\prime})=\xi(|a-a^{\prime}|)$. A physically plausible $\xi$ should take into account the fact that correlations should be strongest for small separations then fall off with “distance”. In this analysis we use the CPZ correlation function [66], which effects a $\sim 1/r^{2}$ fall-off:

[TABLE]

The tuning parameter $a_{c}$ sets the characteristic correlation length while $\xi_{0}$ determines the overall strength of correlations. Larger values of $a_{c}$ correspond to stronger correlations between bins, vice-versa, and in the limit $a_{c}\to 0$ there are no cross-correlations between bins, which implies ${\xi(|a-a^{\prime}|)\to\delta_{D}(a-a^{\prime})}\xi_{0}$. On the other hand, in the limit $a_{c}\to\infty$ the denominator approaches unity and ${\xi(|a-a^{\prime}|)\to\xi_{0}}$ for all $a$ and $a^{\prime}$, and we effectively recover a flat prior.

In this analysis we consider two fiducial models: (i) $\bm{q}^{\mathrm{fid}}$ determined by a five bin running average of $\bm{q}$, which we will call Prior I and (ii) a $\Lambda$CDM-biased prior with $\bm{q}^{\mathrm{fid}}=0$, which we call Prior II. We remark that in using Prior II we can write $\bm{q}^{\mathrm{fid}}={\mathbf{R}}\bm{q}$, for some constant matrix ${\mathbf{R}}$, i.e. $\bm{q}^{\mathrm{fid}}$ is now a function of $\bm{q}$. Technically, the PDF (3.4) is not normalised to unity upon inserting $\bm{q}^{\mathrm{fid}}={\mathbf{R}}\bm{q}$, but instead must be rescaled by multiplying it by $\det[({\mathbf{I}}-{\mathbf{R}})({\mathbf{I}}-{\mathbf{R}})^{\mathrm{T}}]$. The prior remains Gaussian but is now centred about $\bm{q}=0$ with a new covariance matrix that more easily allows the recovery of low-frequency features.

3.3.2 Other priors

For the rest of the parameters, the prior distributions and ranges used in this analysis are as follows: the matter density parameter $\Omega_{m0}\in[0,1]$, the DE density parameter $\Omega_{X0}\in[0,1]$, the baryon density parameter $\Omega_{b0}\in[0,0.4]$, the log of the clustering amplitude $\log_{10}\sigma_{80}\in[-5,2]$ (a log-uniform prior of $\sigma_{80}$). For the SN-specific parameters we choose the intervals $\alpha\in[0,1]$, $\beta\in[0,4]$, $\log_{10}R_{M}^{2}\in[-10,4]$, $\log_{10}R_{x}^{2}\in[-10,4]$, $\log_{10}R_{c}^{2}\in[-10,4]$; for the means we choose $M_{0}\sim\mathcal{N}(19.3,2^{2})$, $x_{*}\sim\mathcal{N}(0,10^{2})$ and $c_{*}\sim\mathcal{N}(0,1^{2})$.

3.4 Identifying a data-oriented basis

The binning of $q(a)$ introduces a large number of extra parameters and it is often the case that the data can support only a few parameters, with the rest being unconstrained. In such cases a common approach to reconstruction is to perform a principal component analysis (PCA) to tease out the features most sensitive to data (usually at low redshifts where data are more abundant). The PCA approach has been widely used to reconstruct the DE equation of state using real data and mock data in forecasting exercises [80, 81, 66, 82, 83, 67, 68, 84, 69]. The principal components (PCs) represent orthogonal directions in parameter space and therefore provide a statistically decorrelated basis. These directions are obtained by performing an eigendecomposition of the Fisher matrix, either of the prior distribution for constraint forecasting or the posterior when using real data. When considered as an expansion in the top-hat basis the PCs are eigenfunctions constructed from an average of $q(a)$ weighted according to those features that are actually being probed by the data. In this case the number of bins $n$ is usually taken to be large so as to allow the recovery of smooth PCs. Using the correlated prior (4.7) allows us to employ a greater number of bins than would otherwise be possible, and effectively smooths out the discrete artefacts that tend to arise when using real data. In this approach the $n$ top-hat basis functions are considered as an intermediate basis from which we construct an uncorrelated basis.

Typically a basis is obtained by finding eigenvectors of the posterior parameter covariance $\mathbf{C}_{p}$. When using the correlated prior (3.4) this approach suffers from an inherent difficulty in finding a natural way to order the modes, as they do not cleanly separate according to their eigenvalues. The modes with the largest eigenvalues generally are unconstrained high-frequency modes, and the number of modes that are actually probed by the data are sensitive to the chosen $\xi_{0}$.

In order to identify a new basis that captures the features the data are probing we find it useful to expand $q(a)$ in terms of a basis expressing the signal-to-noise. This is found by solving the Karhunen-Loève (KL) eigenvalue problem [87]

[TABLE]

with the Fisher information matrices ${\mathbf{F}}_{\pi}=\mathbf{C}^{-1}_{\pi}$ and ${\mathbf{F}}_{p}=\mathbf{C}^{-1}_{p}$ of the of the prior ($\pi$) and posterior ($p$), respectively. ($\mathbf{C}_{\pi}$ can be computed directly from (3.6) but $\mathbf{C}_{p}$ is estimated from the MCMC samples.) The set of generalised eigenvectors $\{\bm{v}_{i}\}$ form the new basis and we order them from highest to lowest signal-to-noise, which we define here as $(S/N)_{i}=\lambda_{i}^{-1}$, for noise given by ${\mathbf{F}}_{\pi}$ and signal by ${\mathbf{F}}_{p}$. Because low-frequency modes have the highest $S/N$, we can consider (3.9) a low-pass filter. In the space of functions spanned by $T_{i}(a)$ we can write (2.5) in this basis as

[TABLE]

where $A_{ij}=(\mathbf{v}_{i})_{j}$, i.e. the rows of ${\mathbf{A}}=(A_{ij})$ are given by the generalised eigenvectors $\mathbf{v}_{i}$. Note $\{\mathbf{v}_{i}\}$ form a complete spanning set of the original basis but are not mutually orthogonal. The basis coefficients $\alpha_{i}$ however are uncorrelated by virtue that its covariance matrix is $({\mathbf{A}}{\mathbf{F}}_{p}{\mathbf{A}}^{\mathrm{T}})^{-1}={\mathbf{I}}$, where ${\mathbf{I}}$ is the identity matrix. It can also be seen that all $\alpha_{i}$ have unit variance and their values can be recovered from $\bm{\alpha}={\mathbf{A}}^{-\mathrm{T}}\bm{q}$.

It is important to quantify how many $\alpha_{i}$ are really being constrained by the data to avoid fitting for the noise of the reconstruction. To do this we compute the Bayesian complexity [85]

[TABLE]

where, for a given posterior $p(\theta|D)$ and prior $\pi(\theta)$,

[TABLE]

is the Kullback-Leibler divergence, which quantifies the amount of information gained from the data, and $\widehat{D}_{\mathrm{KL}}$ is a point estimator of $D_{\mathrm{KL}}[p,\pi]$. The Bayesian complexity effectively tells us how many parameters are being constrained by the data. Restricting to the parameters $\theta=\{q_{1},q_{2},\ldots,q_{n}\}$ by marginalising over all others, since $\pi(\theta)$ is a Gaussian we can write the effective number of $\alpha_{i}$ being constrained as

[TABLE]

If the $S/N$ is high for all modes then the sum is approximately zero and the effective number of parameters is equal to the total number of model parameters. It should be noted however that the formula defined by the first equality of (3.13) only holds if both the prior and posterior distributions are Gaussians [88]. Note that in the case of Prior I, the prior Fisher matrix ${\mathbf{F}}_{\pi}$ is not the same as the inverse of (3.6) but is instead given by ${\mathbf{F}}_{\pi}=({\mathbf{I}}-{\mathbf{R}})^{\mathrm{T}}\mathbf{C}_{\pi}^{-1}({\mathbf{I}}-{\mathbf{R}})$.

4 Analysis

In this section we analyse an interacting model with 20 uniformly spaced bins between $a_{\mathrm{min}}=0.4$ and $a_{\mathrm{max}}=1$. The reconstruction of the interaction is thus restricted to this range, which corresponds to non-uniformly spaced bins from a redshift of $z_{\mathrm{min}}=0$ to $1.5$. Bins containing no data provide little information on the corresponding $q_{i}$ so to ensure that each bin contains at least one data point we have set $a_{\mathrm{min}}=0.4$. For the correlated $\bm{q}$ prior we choose a smoothing length of $a_{c}=0.12$, corresponding to a characteristic correlation with the nearest $\sim 4$ bins. We will consider $\xi_{0}=0.2$, which corresponds to a standard deviation of $q(a)$ when averaged between $a=0.4$ and $a=1$ of $\sigma_{\bar{q}}\approx 0.4$.

The parameters are estimated from the joint posterior given by Bayes’s theorem $p({\theta}|D)\propto\mathcal{L}({\theta};D)\pi({\theta})$, with $\pi({\theta})$ being the joint prior and $\mathcal{L}({\theta};D)$ the joint likelihood formed from the BAO, CMB, CC, RSD, and SN likelihoods. For all data except SN (see Appendix B) we take the data to be Gaussian distributed. The parameter set consists of the following parameters:

[TABLE]

We fix the Hubble constant to $H_{0}=67.3\,\mathrm{km}\,\mathrm{s}^{-1}\mathrm{Mpc}^{-1}$ and specialise to a spatially flat geometry, $\Omega_{X0}=1-\Omega_{m0}$. The posterior is sampled using the affine-invariant MCMC sampler emcee [89].

4.1 Reconstruction

The individual constraints for the 20 bin interaction model are summarised in Figure 1, which are plotted together with the prior probability distribution to show the improvement in each bin. It is clear the first 10 amplitudes $q_{1},q_{2},\ldots,q_{10}$ covering a redshift range from $z=0$ to $z\approx 0.4$ furnish the best constraints. The next few bins show mild improvements, while the last five indicate very little constraining power at those redshifts. As we have noted before this is due to these bins containing more, and better quality, data than those at higher redshifts (and this is even after we used logarithmically spaced bins in redshift space to account for the unevenly distributed data). It should be noted that the CPZ prior requires bins of uniform size (in redshift, scale factor, etc) to ensure $\mathbf{C}_{\pi}$ is positive definite. As Figure 1 shows, Prior I is uninformative as to the amplitudes size of the amplitudes and is essentially flat over an interval of moderate values of $q_{i}$s.

In Figure 2 we reconstructed $q(a)$ (though shown as a function of redshift) and show the probability density function (PDF) of $q(z)$. Although the constraints are not strong it is clear that only for the first few $q_{i}$ are the limits tightest, and in the case of Prior I the variance rapidly grows with redshift. We recall that in the absence of any new information provided by the data the joint posterior distribution is identical to the joint prior of $\bm{q}$.

We have checked that our reconstruction is robust to the number of bins by changing to 10 and 30 bins — both reproducing the basic features seen in Figure 2.

The reconstructed models will of course give a better fit to data than flat $\Lambda$CDM. To see just how much the fit improves we can compare the $\chi^{2}\equiv-2\ln\mathcal{L}(\theta;D)$, with lower values being preferred. When evaluated at the mean parameter values we find a $\Delta\chi^{2}\equiv\chi^{2}-\chi^{2}_{\mathrm{\Lambda CDM}}=-2.2$ when using Prior I and $\Delta\chi^{2}=-3.6$ when using Prior II. These represent very modest improvements in the quality of fit, given we have 20 more parameters than $\Lambda$CDM. We present a model comparison in Section 4.2.1.

4.1.1 How many modes?

Having reconstructed $q(a)$ the question arises as to just how many KL eigenmodes $m$ should be retained in the expansion (3.10): too many and we risk fitting the noise, while too few we fail to capture the physical features being suggested by the data. To assess this trade-off we compute the mean squared error given by $\mathrm{MSE}=\mathrm{bias}^{2}+\mathrm{variance}$ (see, e.g. [80]). Generally, the bias decreases with $m$, while the variance increases with $m$. In our decorrelated basis it is given by

[TABLE]

where

[TABLE]

is the error of the reconstructed $q(a)$ in the $i^{\mathrm{th}}$ bin keeping only the first $m$ modes, each having a variance $\sigma^{2}_{\alpha_{j}}=1$ in the KL basis. Moreover, $q_{i}^{(m)}$ is the associated mean of $\sigma^{(m)}_{i}$, whereas $\bar{q}_{i}$ is the mean in the $i^{\mathrm{th}}$ bin in the original parametrisation (i.e. keeping all modes). A caveat to using (4.2) however is that we assume the true $q(a)$ is given by the mean of the full reconstruction. The optimal number of modes $m$ to keep is thus determined by minimising (4.2). For both Priors I and II we find the MSE is minimised when $m=1$.

We next compute the Bayesian complexity $C$, using the Gaussian approximation given by (3.13). Out of the 20 parameters we find $C=5.2$ for Prior I and $C=3.2$ for Prior II. We show in Figure 4 the first four modes; above $i=5$ for Prior I and $i=3$ for Prior II the modes generally become less smooth as we expect. The third PC in fact provides little information above $z=0.2$. It is clear that the most interesting features of $q(z)$ in our reconstruction is for low-redshifts. This could be for two reasons. The first is that this could be a data effect caused by the specific cosmological probes being used. Unsurprisingly, these two describe the low-redshift features we see in Figure 2. For reference we list the three best constrained eigenmodes: (Prior I) $\alpha_{1}=1.5\pm 1.0$, $\alpha_{2}=0.23\pm 1.00$ and $\alpha_{3}=0.48\pm 1.00$, and (Prior II) $\alpha_{1}=1.2\pm 1.0$, $\alpha_{2}=0.49\pm 1.00$ and $\alpha_{3}=0.89\pm 1.00$. It is clear all except the first mode is consistent with zero (and the constraints only get worse for larger $i$). However, the deviation of $\alpha_{1}$ from zero is not statistically significant (at the $\sim 1\sigma$ level).

From the posterior $p(q_{1},q_{2},\ldots,q_{20})$ estimated from all MCMC samples we can obtain the joint posterior of $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$ from

[TABLE]

To do this we condition on the $q_{i}$s and marginalise over $\alpha_{i}$, for $i=m+1,m+2,\ldots,n$, and use that the joint PDF of $\alpha_{i}$s are separable:

[TABLE]

where $\delta$ is the Dirac delta function. In effect, we project $q(a)$ onto a subspace spanned by a subset of the KL modes, which we achieve in practice by discarding modes $m+1$ and higher. Indeed, in the case $m=20$ we have $q^{(m)}(a)=q(a)$ so that we recover the results of Figure 2.

Figure 5 shows the reconstruction $q^{(m)}(a)$ for the leading $m$ KL modes, the leading $m-1$ KL modes, and also the PDF of the reconstruction (4.5) for both priors. This is the reconstruction being constrained by the data, i.e. after eliminating the extraneous prior modes. The number of modes $m$ are determined by the Bayesian complexity (3.13), which tells us heuristically how many KL modes are being probed by data. It can be seen that without the higher $m$ modes the variance at redshift is significantly reduced. The fifth mode of Prior I, however, dominates the variance above a redshift of $z=0.6$. As can be seen in Figure 4, this mode has $S/N\gtrsim 1$ indicating this mode ought to be discarded, too. Generally we find that the Bayesian complexity provides a useful way to identify the relevant modes.

As expected, the features that the data are actually constraining (i.e. the low-redshift features) in Figure 5 are preserved after discarding the higher $m$ modes. As can be further seen, the fifth mode of Prior I does not possess significant explanatory power as dropping this still reproduces the basic form of the full reconstruction. By constrast, the third mode of Prior I is crucial in reproducing the low-$z$ feature ($z\lesssim 0.4$) of the full reconstruction.

As mentioned in Section 4.1, the improvement in $\chi^{2}\equiv-2\ln\mathcal{L}$ for the full reconstruction gives $\Delta\chi^{2}\approx-3$. This is because any deviations of the best-fit ${q}_{i}$ from zero are small, leading to changes in the DM and DE density of less than $5\%$ at redshifts $z\lesssim 0.6$. It is evident that subspaces of the full reconstruction of $q(a)$ will only yield improvements to the $\chi^{2}$ of $|\Delta\chi^{2}|\lesssim 3$. From the Akaike Information Criterion $\mathrm{AIC}=\chi^{2}+2k$, for $k$ parameters, we see that the introduction of an additional parameter must decrease the $\chi^{2}$ by more than four to be competitive. We thus see that the case for any non-zero $q_{i}$ is not strong.

4.1.2 Sensitivity analysis

In reconstructing $q(a)$ using (3.8) we have specified the tuning parameter $\xi_{0}$ that sets the strength of the correlations. We have chosen $\xi_{0}=0.2$ for both Priors I and II, but different choices are of course possible and it is worth exploring how our inference depends on $\xi_{0}$. Rather than perform the analysis for a range of different $\xi_{0}$ we can also marginalise over it. Thus we expand the hypothesis space to include $\xi_{0}$ as a free parameter. This parameter is however unconstrained by data, as the joint likelihood does not depend on it. Nevertheless we can assign it a prior $\pi(\xi_{0})$ and fold it into the original prior (3.4):

[TABLE]

where $\pi(\bm{q}|\xi_{0})$ will be given by (3.4). As $\xi_{0}$ is a scale parameter we assign it a logarithmically uniform distribution $\pi(\xi_{0})\propto 1/\xi_{0}$. For convenience we define the $\xi_{0}$-independent covariance $\bar{\mathbf{C}}_{\pi}\equiv\mathbf{C}_{\pi}/\xi_{0}$, where $\mathbf{C}_{\pi}$ can be either the covariance of Prior I or II. In general, for $n$ interaction parameters, by integrating $\xi_{0}$ over the range $[\xi_{a},\xi_{b}]$ we obtain

[TABLE]

where $x^{2}(\bm{q})\equiv(\bm{q}-\bm{q}^{\mathrm{fid}})^{\mathrm{T}}\bar{\mathbf{C}}^{-1}_{\pi}(\bm{q}-\bm{q}^{\mathrm{fid}})/2$ and $\Gamma(s,x)$ is the incomplete gamma function. We consider a wide interval with limits $\xi_{a}$ and $\xi_{b}$ that enclose $\xi_{0}=0.2$, as used in (3.4). Compared with the previous priors used, (4.7) has a heavier tail so that regions in parameter space far from $\bm{q}^{\mathrm{fid}}$ are more easily explored. Moreover, since $\xi_{a}<0.2$ we now have a more narrowly peaked mode at $\bm{q}=\bm{q}^{\mathrm{fid}}$, reflecting more confidence in the fiducial model. (In the case of Prior II, (4.7) can be centred on $\bm{q}=0$ by defining a new $\mathbf{C}_{\pi}$ that absorbs the shift.) If $\xi_{b}\gg\xi_{a}$, then (4.7) has faster than Gaussian dropoff since $\pi(\bm{q})\sim(\xi_{b}^{-1}x^{2})^{-1}e^{-\xi_{b}^{-1}x^{2}}$, as $x^{2}\to\infty$ (ignoring multiplicative constants), so that $\pi(\bm{q})$ will be more strongly peaked than the Gaussian prior (3.4).

Figure 6 shows the full reconstruction using the $\xi_{0}$ marginalised prior with bounds $\xi_{a}=0.02$ and $\xi_{b}=2$, i.e. enclosed by an interval with limits an order of magnitude smaller or larger than $\xi_{0}=0.2$ used previously. Using $\bm{q}^{\mathrm{fid}}=0$ with (4.7) we see this prior is more constraining at high redshifts as the lower limit $\xi_{a}$ can be seen to have the effect of causing $q(z)$ at redshifts $z\gtrsim 0.2$ to hew more closely to $q=q^{\mathrm{fid}}$ and suppress the variance; c.f. Figure 2. The slight deviation for $z\lesssim 0.2$ is robust to this prior, and also a more conservative choice of $\xi_{0}=2$ (i.e. $\xi_{0}$ ten times larger than before). Compared with Figure 2, Prior II shows larger deviations of the mean values of $q$ from $q=0$, particularly around $z\simeq 0.8$; this is due to $q$ being less constrained to explore regions away from $\bm{q}=\bm{q}^{\mathrm{fid}}$, but results in a larger variance around the mean. Where there are large deviations (at higher redshifts) they are always accompanied by substantial uncertainties related to the choice of prior.

4.2 Simple one parameter extensions

The previous sections have shown that while we are able to extract a handful of modes with moderate to high $S/N$, only the best mode provides a constraint on the associated $\alpha_{1}$ that is not totally overwhelmed by its uncertainty. This indicates the data are able to support at most one parameter. The improvement in chi-square is marginal and we expect robust model comparison with $\Lambda$CDM to strongly disfavour the full reconstruction model.

Instead of reanalysing the data using a subset of modes (which could be considered using the data twice) we instead consider a one parameter model covering a single wide bin spanning $z_{\mathrm{min}}=0$ to $z_{\mathrm{max}}=1.5$. We will focus in particular on constraints on the associated parameter $q$ and also ask whether the simplest extension can be competitive with flat $\Lambda$CDM.

For comparison with the ansatz (2.3) we will also consider constraints from a physically motivated interaction given by

[TABLE]

In contrast to (2.3), which has an abrupt transition to non-interacting $\Lambda$CDM dynamics, this model smoothly interpolates in a logistic manner between $\Lambda$CDM at early times ($Q\to 0$), to one with DM-DE interaction ($Q\sim qH\rho_{X}$; c.f. (2.3)) at late times.333The interacting model specified by (4.8) can alternatively be viewed as one in which the interacting DM-DE is described by a single fluid with an equation of state $p=-A\rho^{-\alpha}$, where $A$ is a positive constant. Such an exotic fluid is known as a generalised Chaplygin gas [90, 91], and is notable for having asymptotic behaviour that mimics CDM at early times and a cosmological constant at late times [15, 92, 25]. We will call this model $q$XCDM and the one bin model with (2.3) $q_{1}$XCDM. For both models we adopt flat priors on $q$, and again assume a spatially flat background.

Combining all data we find $q=0.039\pm 0.031$ for $q_{1}$XCDM and a slightly stronger constraint of $q=0.041\pm 0.027$ for $q$XCDM, though both are consistent with no interaction at $95\%$ confidence level. The datum $\omega_{b}$ of the compressed CMB data set represents a highly informative prior. Leaving this out of the data set we find a much weaker constraint of $q=0.021\pm 0.029$, a $\approx 85\%$ shift in the mean value.

Figure 7 shows the constraints from individual data sets obtained from their respective likelihoods. It is clear that individual data sets are not sensitive probes to $q$, except for RSD data, which however shows a strong degeneracy with $\Omega_{m0}$. It is therefore necessary to exploit the complementarity of data sets to break the degeneracy. A positive $q$ yields relatively more matter and less DE relative to $\Lambda$CDM. The matter-radiation equality, on which the sound horizon depends, will occur earlier. As the window between the epoch of matter-radiation equality and the time of decoupling is wider, density fluctuations have more time to decay in the radiation-dominated epoch and the overall effect being a suppression of the baryon acoustic peaks [22].

Since the CMB data depends on an assumed cosmological model, it is interesting to consider constraints on $q$ from local cosmological probes only (BAO, CC, RSD, SN), which have the virtue of being model independent. Though CMB data comprises three precise measurements, we find by leaving it out we recover significantly worse constraints on $q$. In addition cosmic chronometers data, while a direct measurement of the expansion history, provide only slight improvements in constraints owing to the large systematic errors and the fact that SN probes roughly the same physics (see Figure 7).

The origin of the dark degeneracy was discussed in [37]. The expansion history as probed by BAO and SN is sensitive only to the evolution of the total energy content, which determines $H$ through Friedmann’s equation, and not to couplings between individual components. The same evolution can be produced in more than one way. For instance, a coupled dark sector is observationally indistinguishable from one that is uncoupled with DE having a certain equation of state. The interaction parameters only enter through the derivative of $T_{\mu\nu}$ and not through Einstein’s equation. As a result the degeneracy persists regardless of whether we consider background cosmological observables or perturbation observables [37]. Although RSD places narrow constraints it comes with an almost exact degeneracy with $\Omega_{m0}$ as seen in Figure 7. We can see this quantitatively as follows. Noting that $q$ is small (of order $10^{-2}$) we linearise (2.9) about $q=0$ and solve perturbatively. To do this we first reformulate (2.9) in terms of $f=d\ln D/d\ln a$. We can thus recast (2.9) as

[TABLE]

where $p(a)$ and $r(a)$ are given by the coefficients of the second and third term of (2.9), respectively. We write the interacting solution as the sum of non-interacting $\Lambda$CDM component $f^{(0)}(a)$ and the interacting component $f^{(1)}(a)$: $f(a)\approx f^{(0)}(a)+qf^{(1)}(a)$. The zeroth-order solution is $f^{(0)}(a)=\Omega_{m}(a)^{\gamma}$, with the matter density obeying the usual scaling and the growth index $\gamma$ equal to 0.55 for flat $\Lambda$CDM [94]. Inserting this into (4.9) then formally integrating we find

[TABLE]

where $u(a)=a^{4}H(a)$, $r^{(1)}(a)$ is the linear term of $r\approx r^{(0)}+qr^{(1)}$ and likewise for $p^{(1)}$. As in [93] we have discarded the quadratic term $(f^{(0)})^{2}$, and here additionally $f^{(0)}f^{(1)}$. In the case of $Q$ given by (4.8) we have $f^{(1)}\sim a^{4}$, as $a\to 0$, i.e. like $\Lambda$CDM  $G$ is small in the early matter-dominated era. Thus we can understand this degeneracy from Figure 8: less matter implies a greater abundance of DE which suppresses the growth of structure. However, this can be offset by a positive interaction rate ($q>0$), which raises the matter abundance across all epochs compared to $\Lambda$CDM (for the same $\Omega_{m0}$) — this explains the anticorrelation between $q$ and $\Omega_{m0}$ as seen in Figure 7.

We can also compare how the growth function changes in the presence of interaction. Recall in flat $\Lambda$CDM the growth function $f(a)$ is completely specified by the expansion history $H(a)$ (assuming also $\rho_{m}\propto a^{-3}$). With interaction, $Q(a)$ alters $H(a)$ and $\rho_{m}(a)$, but also modifies the growth equation through additional terms. We can connect (4.10) to the growth index formalism [94, 93] by integrating $f=d\ln D/d\ln a$ to obtain

[TABLE]

where we defined $A(a;q)=e^{q\int^{a}_{0}d\ln a^{\prime}f^{(1)}(a^{\prime})}$ (since $f^{(1)}\sim a^{4}$ as $a\to 0$ this prefactor tends to unity as $a\to 0$; c.f. Figure 8). It is clear that the growth factor depends on $q$ in addition to the growth index parameter $\gamma$. Thus if we attempt to fit $\gamma$ assuming $f=\Omega_{m}(a)^{\gamma}$ we are liable to find a systematic bias, resulting in a shift of $\gamma$ from its $\Lambda$CDM value of $0.55$. The growth index formalism is commonly used to probe modified gravity but care must be taken when interpreting $\gamma$ as this analysis shows it is possible to have both a scale-indepedent growth factor $D(a)$ and a value of $\gamma\neq 0.55$, and still have the theory of gravity be given by general relativity.

4.2.1 Model comparison

Any model that is an extension of $\Lambda$CDM is guaranteed to give an equally good or better fit to data. The question that must be asked of the model is whether the improvement in fit justifies the additional flexibility. To address this we compute the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). For a model with $k$ parameters and maximum likelihood $\mathcal{L}_{\mathrm{max}}$ the AIC is given by $\mathrm{AIC}=-2\ln\mathcal{L}_{\mathrm{max}}+2k$, and the BIC is given by $\mathrm{BIC}=-2\ln\mathcal{L}_{\mathrm{max}}+k\ln N$, with $N$ the size of the data set. It should be noted that while the AIC and BIC are very similar, differing only in their penalty term, they are obtained from different starting points: the AIC has its origins in information theory and is based on an approximation of the Kullback-Leibler divergence (c.f. (3.12)), while the BIC estimates the (logarithm) Bayesian evidence under the assumption that the likelihood is Gaussian and $N$ is large. The quantity of interest however is the difference, and we will take positive values of $\Delta\mathrm{AIC}$ and $\Delta\mathrm{BIC}$ to indicate preference for the interacting model.

We find $\Delta\mathrm{AIC}=-0.69$ when comparing spatially flat $q_{1}$XCDM and $\Lambda$CDM with $k=12$ and 11 parameters, respectively. This indicates a slight preference for $\Lambda$CDM. We find also $\Delta\mathrm{BIC}=-5.4$, indicating strong preference for $\Lambda$CDM when assessed on Jeffreys’s scale. In the case of spatially flat $q$XCDM and $\Lambda$CDM, we find $\Delta\mathrm{AIC}=0.22$ and $\Delta\mathrm{BIC}=-4.5$. Again the BIC is decisively in favour of the non-interacting scenario. That the BIC is more penalising than the AIC is not surprising: the AIC is generally more accommodating of additional parameters when $N$ is large, but unlike the BIC it is dimensionally inconsistent, in that the tendency of AIC to select the more complex model does not fall as the size of the data set grows [95].

As a further comparison, we will also compare $\Lambda$CDM with the previous reconstruction of $q(a)$. While this may be seen as post-hoc tuning of the interacting model, it is nevertheless interesting to see whether this model can be competitive under favourable circumstances. We will thus focus on a minimal model that keeps just the best measured mode, parametrised by $\alpha_{1}$. This model is a one parameter extension of $\Lambda$CDM upon marginalising over all other $\alpha_{i}$. The model selection measure of choice we take to be the Bayesian evidence $E=\int d\theta\,\mathcal{L}(\theta;D)\pi(\theta)$, which readily embodies Occam’s razor. Given two models, model selection is decided by computing the Bayes factor $B$ given by the ratio of the evidences. We will take $\ln B<0$ to mean preference for the reconstructed model and $\ln B>0$ to mean preference for $\Lambda$CDM. Since non-interacting models are nested within interacting models, we can compute the Bayes factor using the Savage-Dickey density ratio formula $B={p(\alpha_{1}=0|D)}/{\pi(\alpha_{1}=0)}$. We recall in the KL basis $\alpha_{1},\alpha_{2},\ldots,\alpha_{20}$ are uncorrelated variables and the joint prior is also a Gaussian that can therefore be separated: $\pi(\bm{\alpha})=\pi(\alpha_{1})\pi(\alpha_{2})\ldots\pi(\alpha_{20})$. Marginalising over $\alpha_{2},\alpha_{3},\ldots,\alpha_{20}$ is trivial for Gaussian distributions and we have that in this decorrelated KL basis $\pi(\alpha_{1})$ is given by a 1-dimensional Gaussian with a mean of zero and a variance of $1/\lambda_{1}$ (for both Priors I and II). Parameters shared by both models will also have the same priors. We find $\ln B=-1.4$ for Prior I and $\ln B=-1.9$ for Prior II, indicating weak to moderate evidence for $\Lambda$CDM as judged on Jeffreys’s scale, although not as decisively as in the one parameter models considered earlier. The slightly less favourable evidence in the case of Prior II is to be expected as this prior takes $\bm{q}^{\mathrm{fid}}=0$, which necessarily drives the regression towards reconstructions consistent with null interaction.

5 Future prospects: the Fisher forecast

As is clear from the reconstruction we are not yet able to obtain tight constraints on a possible DM-DE coupling. A more detailed reconstruction of the finer features of $q(a)$ may be possible in the future with upcoming stage-IV surveys, such as with the Large Synoptic Survey Telescope (LSST) [96] and the Dark Energy Spectroscopic Instrument (DESI) [97]. In anticipation of this we use the Fisher framework to forecast improvements in the constraints obtained in this analysis.

We assume a SN redshift distribution of the form

[TABLE]

For the LSST survey we take $z_{0}=0.04$ and $\alpha=0.7$. As LSST is expected to deliver $>10^{4}$ SNe Ia per year for ten years we consider what constraints we might obtain with a one year sample with a total number of $N_{\mathrm{tot}}=5\times 10^{4}$ SNe Ia and the full ten year sample with $N_{\mathrm{tot}}=5\times 10^{5}$ SNe Ia. We thus draw $N_{\mathrm{tot}}$ samples distributed according to the probability density function $p(z)=(1/N_{\mathrm{tot}})dN/dz$. For simplicity we adopt the conventional SN analysis in which the data are Gaussian distributed and the Fisher matrix is given by

[TABLE]

where the indices $i$ and $j$ runs over all parameters. We take the total error of the apparent magnitude $\sigma_{m_{k}}$ to be made up of an intrinsic scatter $\sigma_{\mathrm{int}}$ and assume a Gaussian redshift uncertainty $\sigma_{z}$. Propagating the redshift error by (3.2) the magnitudes have a total squared error given by

[TABLE]

Here $\sigma_{\mathrm{int}}=0.12$ and we assume photometric redshift errors modelled by a linear drift, $\sigma_{z_{k}}=0.05(1+z_{k})$.

In addition to SN constraints from LSST we also include constraints from RSD data from DESI and the DESI Bright Galaxy survey [97]. We use the projected RSD constraints based on a sky coverage of $14,000\,\text{deg}^{2}$ and the pessimistic wavenumber cutoff of $k_{\mathrm{max}}=0.1h\,\text{Mpc}^{-1}$. The DESI estimates are of $f\sigma_{8}$ at 18 redshifts between $z=0.05$ and $z=1.85$ (see Tables 2.3 and 2.5 in [97]). We further supplement this with DESI projected errors on the transverse and radial BAO scales, $d_{A}/r_{s}$ and $Hr_{s}$, respectively, and which are correlated measurements with a correlation coefficient of 0.4. We assume negligible correlations between current and future galaxy surveys, i.e. we assume the data are not being double counted.

We forecast constraints using the Fisher information matrix

[TABLE]

where $F_{ij}^{\mathrm{data}}$ is computed from the joint likelihood and represents hypothetical future prior information obtained in this work from current data. To obtain a sense of the improvement in constraints possible we focus on a single wide bin between $z=0$ and $z=1.5$. We evaluate the Fisher matrix at the maximum likelihood estimate using $q$XCDM, which we find to be $\Omega_{m0}=0.3116$, $\Omega_{b0}=0.04977$, $q=0.041$, and $\sigma_{80}=0.767$. As $q$ is most strongly correlated with $\Omega_{m0}$ we show in Figure 9 the joint constraints in the $(q,\Omega_{m0}$)-plane having marginalised over all other parameters. We find a constraint of $\Delta q=\pm 0.019$ and $\pm 0.0078$ ($68\%$ C.L.) for the one and ten year SN sample, respectively. For the full ten year LSST observing run we expect the constraints on $q$ to improve by a factor of $\approx 3.5$. Although this represents only a $\approx 20\%$ determination of $q$ it is still sufficient to rule out the no-interaction hypothesis at $99\%$ confidence level.

6 Conclusions

In this work we have investigated in a model-independent way a range of minimal extensions to $\Lambda$CDM that relax the usual assumption that DM and DE are non-interacting species. We adopted a non-parametric approach and reconstructed directly from data the interaction history within a low-redshift window of $z<1.5$. Using a range of primarily low-redshift cosmological probes we do not see statistically significant evidence for coupling between DM and DE, although we note that the tightest constraints obtained in the reconstruction are for $z\lesssim 0.4$, which also happen to show the strongest signs of a departure from $\Lambda$CDM. While these departures are $\lesssim 2\sigma$ in statistical significance, they suggest a very mild late-time breakdown of the assumptions of $\Lambda$CDM.

We find that although several Karhunen-Loève modes of the reconstruction can be cleanly extracted from the data, only one mode can be satisfactorily constrained. That most of the interaction parameters have substantial uncertainties is not surprising given the weak constraints on $q(a)$ can be likened to that of the dark energy equation of state $w_{X}(a)$. Current data are not yet able to rule out the cosmological constant scenario using the CPL parametrisation [98, 99] $w_{X}(a)=w_{0}+w_{a}(1-a)$, which is the simplest dynamical dark energy model. Like $w_{X}$ the interaction parameter $q$ directly modify the scaling relations of the energy density and so enter $\rho_{X}$ and $\rho_{c}$ through an integral. On the other hand geometric probes are built from distance measures that are integrals over the expansion and thus sensitive to $q$ and $w_{X}$ through a double integral. In both cases this limits the effectiveness of such probes and emphasizes the need for dynamical probes, such as RSDs. However, in the case of interaction it is also necessary to combine probes to break the ‘dark’ degeneracy to obtain tight constraints [37]. Therefore the situation with interaction should not be expected to yield strong constraints even for the minimal models considered here.

The added flexibility of an extra degree of freedom should always be weighed against the increased model complexity, and in this work we have also assessed the viability of interacting models. We find that $\Lambda$CDM remains the preferred model over all interacting scenarios considered here. Even under favourable conditions, the most competitive interacting model of the reconstructions (using the best constrained KL mode only) we find that $\Lambda$CDM is still favoured but with weak Bayesian evidence ($|\ln B|=1.4$).

Our constraints on $q$ for the minimal models (Section 4.2) are consistent with those in the recent literature using ansatz forms of $Q$. In [100] four parametrisations of $Q$ were investigated, finding interaction to be statistically insignificant in each. More recently, in [32], a Bayesian comparison is performed using two novel parametrisations of $Q$ proportional to the dark matter density that specially overcome the early-time instability [18]. For various data combinations, the authors find that $\Lambda$CDM is always favoured with Bayesian evidence ranging from “positive” to “very strong”. Similarly in [31], interaction is investigated in the context of relieving the $H_{0}$ tension. The authors find moderate Bayesian evidence $\ln B=1.8$ for interaction using Planck 2015 and a prior on $H_{0}$ given by the SH0ES 2016 determination [101]. However, by also adding BAO and SN Ia data the Bayes factor reduces to $\ln B=0.53$, i.e. inconclusive evidence for interaction. In [102] a more elaborate quintessence-type interaction model is studied, in which a spin-$\frac{1}{2}$ fermionic field (DM) is coupled to a canonical scalar field (DE) by the Yukawa interaction, following the earlier work of [103]. This model is described by two additional parameters for the Yukawa potential ($\lambda$) and the coupling ($r$). Though KiDS weak lensing and CMB lensing mesaurements are considered, the tightest constraints are found using only the Planck 2015 data set (TT, TE, and EE): $r=-0.01074^{+0.0424}_{-0.0426}$ at $68\%$ confidence level, with $\lambda$ found to be poorly constrained. (Interestingly, using only the Planck data set, the authors find this model is moderately favoured over $\Lambda$CDM using the Deviance Information Criterion.)

The notable exception to these results reaffirming $\Lambda$CDM is that of Salvatelli et al. [27], who reported a late-time interaction that excluded $\Lambda$CDM at $99\%$ confidence level. While these results generated much interest in the phenomenology of interacting models as viable extensions to $\Lambda$CDM, but to date the results of [27] have not been borne out by recent data. In particular, the recent work of [104] reanalysed the same interacting models using recent data sets and found no significant deviation from $\Lambda$CDM. They suggest that the differences in results with the earlier analysis are due simply to differences in the Planck 2013 and 2015 data sets, with the latter no longer favouring interaction. However, we also note that the constraints of [104] on the same model are more stringent, coming from a combined analysis of CMB and the low redshift probes of RSD, BAO, and SN. By contrast the results of [27] are based on CMB (Planck 2013 temperature fluctuations and WMAP polarization) and RSD data. While interacting scenarios can resolve tensions in local analyses [31, 32], $\Lambda$CDM is still robust against a wide variety of probes.

Though we have aimed to be as model-independent as possible, we have nevertheless had to make assumptions about the type of interaction. We have thus assumed the covariant interaction 4-vector is directed along the geodesic flow of CDM (i.e. no momentum transfer) and also that DE has an equation of state $w=-1$. With these assumptions we introduce no additional dynamical degrees of freedom and maintain the picture of DE as a vacuum energy, only now with the possibility of an energy exchange with CDM. Further, the inclusion of an equation of state $w\neq-1$ typically degrades the parameter constraints as extra parameters generally increase parameter uncertainties overall. (It is not surprising that an extended parameter space is able to partially ease cosmological tensions.) Though this gives a more general description of the dark sector we note that the freedom enabled by a variable DE equation of state will to some extent be captured by the model-independent form of $Q$ considered here. We leave this more comprehensive analysis to future work when better data become available.

As a caveat to using the CPZ correlated prior (3.4) we are forced to assume a fiducial model of $q(a)$, which necessarily introduces some bias into the reconstruction. Of course, some structure must be introduced into the problem if we are to make progress. We have obtained results against two fiducial models: one with a running average of $\bm{q}$ as in [67], and a more conservative choice that favours null interaction with $\bm{q}^{\mathrm{fid}}=0$. Reassuringly, features of the reconstruction where data are abundant ($z\lesssim 0.4$) are largely the same meaning that the choice of $\bm{q}^{\mathrm{fid}}$ is mostly irrelevant. By contrast, going out to higher redshifts we find the reconstruction is prior dominated, being highly dependent on the choice of $\bm{q}^{\mathrm{fid}}$.

The next generation of surveys will certainly allow a more elaborate, fine-scale reconstruction of the interaction history than considered here. A Fisher forecast shows that the constraints from low-redshift data can be expected to improve by a factor of two, for a minimal one parameter interaction model. In summary, current data do not reveal any statistically significant deviations from $\Lambda$CDM. If however there is interaction to be found at the level $q\gtrsim 0.04$, we anticipate that a one year sample of LSST SN data combined with DESI BAO and RSD data will be capable of distinguishing from $\Lambda$CDM at a statistical significance of $\simeq 3\sigma$.

Acknowledgements

We thank Jan Hamann and David Parkinson for comments and discussions during the early stages of this work, and Florian List for his reading and comments on an early manuscript. LD is supported by the Australian government Research Training Program. KB acknowledges the support of the Australian Research Council through the Future Fellowship FT140101270. This work has made use of the publicly available codes CAMB [106], emcee [89], and getdist.444https://github.com/cmbant/getdist/ The authors acknowledge the use of Artemis at The University of Sydney for providing HPC resources that have contributed to the research results reported within this paper.

Appendix A On using standard fitting formulae

In general no analytical form exists for $z_{*}$ and $z_{\mathrm{d}}$ — though widely used fitting formulae have been given [105]. We justify the use of these fitting formulae for our interacting model by noting the following. The sound horizon depends on physics before the time of recombination physics ($z\sim 1000$) going back to the big bang. While the standard formula given usually assume $\Lambda$CDM the fact that the interaction considered here is a late-time effect meaning that the evolution tracks $\Lambda$CDM up until $z\lesssim 2$. Thus, up to the time of last scattering DM and DE evolve in their usual way. In non-interacting models $z_{*}$ has a weak dependence on $\omega_{b}$ and $\omega_{c}$; we find in our interacting model that it also has a weak dependence on $q$, as well. Moreover, in standard recombination we have

[TABLE]

where $R(z)\equiv\dot{\rho}_{b}/\dot{\rho}_{\gamma}$ and $R_{\mathrm{eq}}=R(z_{\mathrm{eq}})$, i.e. evaluated at matter-radiation equality $z_{\mathrm{eq}}$ (which is in general different from $\Lambda$CDM). The logarithm of (A.1) is order one and insensitive to cosmological parameters, because the sound horizon is determined largely by the physical densities $\omega_{b}$ and $\omega_{c}$.

In Figure 10 we have confirmed using CAMB [106] that the error for $z_{\mathrm{d}}$ and $z_{*}$ is less than $0.1\%$ for a wide range of $q$. Although this is for a single interaction parameter, we expect this approximation to also hold for models with multiple $q$, provided they are all of similar magnitude.

Appendix B Details of the SNe Ia likelihood

The commonly used chi-squared method has been shown to be inadequate, leading to a systematic shift in the cosmological parameters [53]. This shift is due to the uncertainties of the colour parameter $c$ being roughly of the same size as its scatter [55] and the issue originates from the implicit flat priors on the stretch and colour parameters. The solution is to adopt informative priors and below we sketch the construction of the likelihood using a recent Bayesian hierarchical approach.

The dependence of the observables on the latent variables can be made explicit by writing the likelihood as

[TABLE]

Note that we have ignored measurement errors in the redshift as it has little effect on parameter estimation [53]. (The redshift errors, however, are still propagated to the apparent magnitude using (3.2).) Thus we put $\hat{\bm{z}}=\bm{z}$ and have suppressed the dependence in redshift. Further conditioning on the latent stretch and colour parameters we have

[TABLE]

Since $\bm{x}_{1}$ and $\bm{c}$ can be considered independent variables the joint prior is separable, i.e. $p(\bm{x}_{1},\bm{c}|{\theta})=p(\bm{x}_{1}|{\theta})p(\bm{c}|{\theta})$. As $p(\bm{m}_{B}^{*}|\bm{M},\bm{x}_{1},\bm{c},{\theta})$ expresses the deterministic relation (3.3) we have

[TABLE]

where $\delta$ is the Dirac delta function and $\bm{m}_{B}^{*}(\bm{M},\bm{x}_{1},\bm{c};{\theta})=\bm{\mu}(\bm{z};{\theta})+\bm{M}-\alpha\bm{x}_{1}+\beta\bm{c}$. Marginalising over $\bm{m}_{B}^{*}$ (B.1) becomes

[TABLE]

The SALT2 outputs for the $i^{\mathrm{th}}$ SN Ia $\{\hat{m}_{B,i}^{*}$, $\hat{x}_{1,i}$, $\hat{c}_{i}\}$ are correlated variables and so we construct $p(\hat{\bm{m}}_{B}^{*},\hat{\bm{x}}_{1},\hat{\bm{c}}|\bm{m}_{B}^{*},\bm{x}_{1},\bm{c},{\theta})$ to be a $3N$-dimensional multivariate Gaussian with mean given by the corresponding latent variables ($\bm{m}_{B}^{*}$, $\bm{x}_{1}$ and $\bm{c}$) with a $3N\times 3N$ covariance block diagonal matrix $\mathbf{\Sigma}_{\mathrm{stat}}$. Each SN Ia, having latent variables $M_{i}$, $x_{1,i}$ and $c_{i}$, are plausibly assumed to be drawn from independent and identical Gaussian distributions, i.e. $M_{i}\sim\mathcal{N}(M_{0},R_{M}^{2})$, $x_{1,i}\sim\mathcal{N}(x_{*},R_{x}^{2})$ and $c_{i}\sim\mathcal{N}(c_{*},R_{c}^{2})$. With these considerations the integral (B.1) reduces to a convolution of Gaussians, which can be analytically resolved to give the marginalised negative log-likelihood

[TABLE]

where

[TABLE]

are vectors of length $3N$, $\mathbf{\Sigma}_{\mathrm{stat}}$ is the covariance matrix of statistical uncertainties from the light-curve fit, and $\mathbf{\Sigma}_{\mathrm{sys}}$ is the covariance matrix of systematic uncertainties, including from calibration, the light-curve model, dust extinction, and bias uncertainty [48]. Finally, ${\mathbf{B}}=\mathrm{diag}({\mathbf{J}},{\mathbf{J}},\ldots,{\mathbf{J}})$ is a $3N\times 3N$ block diagonal matrix where each block is identical with

[TABLE]

Note $\mathbf{\Sigma}_{\mathrm{stat}}$ has parameter dependence so the normalisation term of (B.2) cannot be neglected in MCMC parameter estimation. Moreover despite the Gaussian form of (B.2) the data $\{\hat{\bm{m}}_{B}^{*},\hat{\bm{x}}_{1},\hat{\bm{c}}\}$ are not Gaussian distributed. For a more detailed derivation, including analytic marginalisation of hyperparameters, we refer the reader to [53].