# Joint reconstructions of growth and expansion histories from Stage-IV
surveys with minimal assumptions I: Dark Energy beyond $\Lambda$

Rodrigo Calderón<EMAIL_ADDRESS>Korea Astronomy and Space Science
Institute, Daejeon 34055, Korea  Laboratoire Charles Coulomb, Université de
Montpellier & CNRS, 34095 Montpellier, France  Benjamin L’Huillier
<EMAIL_ADDRESS>Department of Physics and Astronomy, Sejong University,
Seoul 05006, Korea  Department of Astronomy & Space Science, Yonsei
University, Seoul 03722, Korea  David Polarski<EMAIL_ADDRESS>Laboratoire Charles Coulomb, Université de Montpellier & CNRS, 34095
Montpellier, France  Arman Shafieloo<EMAIL_ADDRESS>Korea Astronomy and
Space Science Institute, Daejeon 34055, Korea  University of Science and
Technology, Daejeon 34113, Korea  Alexei A. Starobinsky<EMAIL_ADDRESS>L. D. Landau Institute for Theoretical Physics RAS, Moscow 119334, Russia

###### Abstract

Combining Supernovae, Baryon Acoustic Oscillations and Redshift-Space
Distortions data from the next generation of (Stage-IV) cosmological surveys,
we aim to reconstruct the expansion history up to large redshifts using
forward-modeling of
$f_{\mathrm{DE}}(z)=\rho_{\mathrm{DE}}(z)/\rho_{\mathrm{DE,0}}$ with Gaussian
processes (GP). In order to reconstruct cosmological quantities at high
redshifts where few or no data are available, we adopt a new approach to GP
which enforces the following minimal assumptions: a) Our cosmology corresponds
to a flat Friedman-Lemaître-Robertson-Walker (FLRW) universe; b) An Einstein
de Sitter (EdS) universe is obtained on large redshifts. This allows us to
reconstruct the perturbations growth history from the _reconstructed_
background expansion history. Assuming various DE models, we show the ability
of our reconstruction method to differentiate them from $\Lambda$CDM at
$\gtrsim 2\sigma$.

##  I Introduction

Despite its success, the concordance model of cosmology, the $\Lambda$CDM
model, is facing numerous challenges [1, 2, 3, 4]. On the theory side, the
nature of dark energy (DE) and dark matter remain a mystery though their
phenomenology is constrained from the data with increasing accuracy. In
addition, while general relativity (GR) has been very successful on solar
system scales, alternative gravity models with observable signatures on large
cosmic scales remains a challenging possibility. Observationally, despite all
the successful predictions of the concordance model, the arrival of data of
ever increasing quantity and quality has led to the rise of tensions, most
notably, the $H_{0}$ and the $S_{8}$ tensions [5, 6, 7, 8]. This is another
incentive to question the nature of gravity itself [9]. As it is well-known,
combining background and perturbations can be a decisive tool to constrain
gravity and hence the nature of dark energy [10]. Very generally, extraction
of physical information from the data in a way that is as reliable as possible
remains an important challenge.

To face these issues, two main approaches have been adopted by the community.
A common approach is to confront new theories to the data, hence constraining
the model parameters. On the other hand, one can perform a bottom-up approach
from the data to the phenomenology [_e.g._ 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22]. This latter approach aims to retrieve the “correct” phenomenology
from the data as independently as possible of theoretical assumptions. In this
work, we take a hybrid approach, in-between traditional model-fitting and
model-independent approaches. At the core of the work presented here lies the
reconstruction of the expansion history from the data with Gaussian processes
(GP). Indeed Gaussian processes are a very effective tool when one follows an
approach where no or little theoretical assumptions are made on the possible
DE phenomenologies. However this method can lead to large and undesirable
uncertainties at large redshifts, from $z\gtrsim 1-2$ on [19]. This is a
consequence of the method itself when only a small amount of data are
available. For example one should remember that no SNIa data are available
when one comes close to $z\sim 2$. Nevertheless sticking to this
reconstruction method in view of its advantages and elegance, we first devise
a way to get reasonable uncertainties on very large redshifts: we simply
impose that our universe should be close to an Einstein de Sitter (EdS)
universe on large redshifts (we have in mind here redshifts still obeying
$z\ll z_{\rm eq}$). We assume a flat FLRW universe, where for $z\ll
z_{\mathrm{eq}}$, the dimensionless expansion rate $h(z)\equiv H(z)/H_{0}$ is
given by

$h^{2}(z)=\Omega_{\mathrm{m,0}}\;(1+z)^{3}+(1-\Omega_{\mathrm{m,0}})\;f_{\mathrm{DE}}\,(z)\;,$
(1)

where
$f_{\mathrm{DE}}(z)\equiv\frac{\rho_{\mathrm{DE}}(z)}{\rho_{\mathrm{DE,0}}}$
encodes the evolution of the dark energy density. Equality (1) covers many
physical situations. As we have in mind a joint reconstruction of both the
background expansion and the perturbations, we postpone a discussion of the
scope and limitations of (1), and generally of the method presented here,
after we address the evolution of the perturbations (see Eqs. (2) and (3)
below). Hereafter, EdS universe will refer to the limit when the second term
of (1) becomes negligibly small. While this is obvious theoretically, its
implementation can be problematic when using Gaussian Processes. We will
reconstruct $f_{\mathrm{DE}}(z)$ using GP with reasonable minimal assumptions.
In this way only these expansion histories are reconstructed for which the
second term in (1) is negligible compared to the dust-like matter term at
large-$z$. While this can be seen, and indeed is, a theoretical prior, we
believe it is altogether reasonable and leaves the possibility to consider
many dark energy models which is the key point. In other words, this prior can
still fit a very large number of dark energy models and these are the models
we are interested in here. Another advantage of this approach is to render
this assumption explicit rather than tacit. So (1) fixes the theoretical
framework in which we want to get the dark energy phenomenology and the
reconstruction is such that we exclude behaviours at large redshifts which are
in our opinion less interesting. In this work we will also assume a spatially
flat universe as it is obvious from (1) and we will comment in the conclusion
on relaxing this assumption. In the next sections we present the method and
check its successful use. Last but not least, we stick to General Relativity
(GR) in this work, leaving modified gravity models for future work. The method
and data are described in Section II, Section III shows the validity of the
method, and Section IV is devoted to the reconstructions of other fiducial
cosmologies. We discuss our results in Section V.

##  II Method and Data

###  II.1 Method

In this paper, we aim to reconstruct the expansion and growth histories
without assumptions concerning the dark energy component and assuming GR. To
do so, as we have said earlier, we assume a flat FLRW universe, where for
$z\ll z_{\mathrm{eq}}$, the dimensionless expansion rate $h(z)\equiv
H(z)/H_{0}$ is given by (1). We note that the $f_{\mathrm{DE}}$ term is
essentially free to capture a large variety of DE behaviours.

In General Relativity, the equation for the growth of perturbations in the
Newtonian approximation valid on observable (sub-Hubble radius) cosmic scales
is given by

$\displaystyle\ddot{\delta}+2H~{}\dot{\delta}$ $\displaystyle=4\pi
G~{}\rho~{}\delta,$ (2)

where a dot stands for a derivative with respect to time, $\delta$ is the
relative energy density fluctuation of dust-like matter and $\rho$ its
background energy density. The above equation can also be re-written in terms
of the growth factor $f\equiv\delta^{\prime}/\delta$, to give

$\displaystyle
f^{\prime}+\left(f+2+\frac{h^{\prime}}{h}\right)f-\frac{3}{2}\Omega_{\text{m}}=0~{},$
(3)

where a prime stands for a derivative with respect to $\ln a$. Let us discuss
now the variety of situations for which our theoretical set-up holds. In
models where DE perturbations are negligible, we can use Eqs. (2) and (3) as
they stand. From this point of view, the two essentially dark components
appearing in (1) are divided in the following way: only the first term
includes baryonic matter and cold dark matter, which undergo gravitational
clustering at the comoving scale corresponding to $8\,h^{-1}\rm Mpc$ today.
The second term can even include hot non-relativistic matter, for example
massive neutrinos, for which the free-streaming scale is much larger than this
scale. Of course, strictly speaking, while $\Omega_{\mathrm{m,0}}$ appearing
in the first term of (1) corresponds to clustered matter – crucially, it is
this quantity which would appear in the driving term of the equation for the
growth of matter perturbations (2) and (3) – the parameter
$\Omega_{\mathrm{DE,0}}$ on the other hand includes all remaining unclustered
components besides DE, and the same applies to $f_{\mathrm{DE}}$. But these
effects remain tiny because $\Omega_{\nu}\ll 1$. In the same way we do not
include here the effect of massive neutrinos on the growth of perturbations
which translates into a very tiny damping during the matter era with
$\delta_{m}\propto a^{1-\frac{3}{5}\Omega_{\nu}}$ for $\Omega_{\nu}\ll 1$.
This cumulative effect becomes however relevant if we integrate back from very
high redshifts of order $z_{\mathrm{eq}}$, which is not the case here. In this
context, we note that the presence of any hot matter can be pinpointed by
considering matter perturbations on larger comoving scales, but there are
unfortunately few such accurate data available at the present time.

We note further that equations (1–3) can in principle apply to non-interacting
tracking (early) dark energy. In that case DE tends to a dust-like behaviour
with $\Omega_{\mathrm{DE}}\to\Omega_{\mathrm{DE},\infty}\neq 0$
($\Omega_{m}\to\Omega_{m,\infty}$), with
$\Omega_{\mathrm{DE},\infty}\ll\Omega_{m,\infty}$ from observations (see e.g.
[23]). We stress again that when tracking DE or massive neutrinos behaves like
dust, their contribution to (1) is very small compared to the first term.

Interestingly, the second term in (1) can in principle also include the energy
density of dark matter-dark energy interactions in models where these exist.
However this term is known to be small because observations severely restrict
deviations of $w_{\mathrm{DE}}$ from $w_{\Lambda}=-1$ on one hand, but also
because these interactions, whether physical or purely gravitational as in
dark energy models based on scalar-tensor gravity [24], would result in a
small gravitational clustering of the components represented by the second
term due to clustering of the dust-like matter components, and we do not
consider this clustering here.

Observationally, redshift-space distortion (RSD) provides us with growth rate
measurements of the (unbiased) quantity111We note that the pipelines used to
extract $f\sigma_{8}$ themselves assume a fiducial cosmology, which may
introduce a bias.

$f\sigma_{8}\equiv\frac{\sigma_{8,0}}{\delta_{0}}f\delta=\frac{\sigma_{8,0}}{\delta_{0}}\delta^{\prime},\;\;\;\text{
with }\;\;\;\delta_{0}=\delta(z=0).$ (4)

To validate the method, we proceed as follows.

Figure 1: Diagramatic representation of our method. Given a set of
hyperparameters $(\sigma_{f},\ell_{f})$ we obtain a behaviour of
$f_{\mathrm{DE}}$, which together with a value of $\Omega_{\text{m},0}$,
yields an expansion history $h(z)$. We then solve the growth equation assuming
a $G_{\mathrm{eff}}$ (though in the present work we stick to GR, where $G_{\rm
eff}={\rm const}=G$) and a given $\sigma_{8,0}$ to obtain the growth history
$f\sigma_{8}$. The values for the parameters (in light-blue) are obtained
through MCMC sampling.

  1. 1.

Generate the mock data:

    1. (a)

Choose a fiducial cosmology

    2. (b)

Generate the background (SN, BAO) and perturbations (RSD) mock data according
to this cosmology.

  2. 2.

Reconstruct the expansion and growth histories using MCMC: Choose the
cosmological parameters $(\Omega_{\text{m},0},\sigma_{8,0})$ and the hyper
parameters $(\sigma_{f},\ell_{f})$ via MCMC.

For each point in the parameter space:

    1. (a)

Generate a single realization of
$f_{\mathrm{DE}}\sim\mathcal{GP}(-1,K(\sigma_{f},\ell_{f}))$—imposing $f_{\rm
DE,0}\equiv 1$.

    2. (b)

Obtain the expansion history $h(z)$ via Eq. (1) for
$f_{\mathrm{DE}},\Omega_{\text{m},0}$

    3. (c)

Solve the growth equation (3) and obtain the growth history $f\sigma_{8}(z)$

    4. (d)

Confront the expansion and the growth histories to the mock data and obtain a
likelihood for this choice of (hyper)parameters.

A diagrammatic representation of our method is displayed on Figure 1. In §
III, we will show the effect of different choices of mock datasets on the
reconstructions of the expansion history.

###  II.2 Gaussian Process modelling

Parameter | $\Omega_{\rm m,0}$ | $\sigma_{8,0}$ | $\sigma_{f}$ | $\ell_{f}$
---|---|---|---|---
Prior | $[0.1,0.6]$ | $[0.6,1.2]$ | $[10^{-2},2.5]$ | $[10^{-2},7]$
Table 1: Uniform priors for the parameters used in the MCMC analyses.

In this study, we need to know the expansion history up to high redshift in
order to solve the growth equation (2). Gaussian Processes (GP) have been
extensively used in the literature [25, 26, 27, 15] to reconstruct the
expansion history directly from the data.

A clear downside of the usual application of GP is that one can only
reconstruct a quantity where there are data, as in the absence of data, the GP
will fall back to its mean function. Therefore, modeling $h$ as a GP can (i)
lead to an incorrect asymptotic behaviour of $h$ (because of the lack of data
at high-$z$) and (ii) does not disentangle the matter term
$\Omega_{\text{m},0}$ from the dark energy term , which is necessary to solve
the growth equation (2) in a consistent way. Hence, we must find a way to
impose $h^{2}(z\gg 1)\sim\Omega_{\text{m},0}(1+z)^{3}$: we enforce this EdS
limit regardless of the value of $\Omega_{\text{m},0}$, by modelling
$f_{\mathrm{DE}}$ as a GP with mean $\bar{f}=1$ and imposing
$f_{\mathrm{DE}}(z=0)=1$.

This ensures the correspondence between the $\Omega_{\text{m},0}$ part of $h$
and $\Omega_{\text{m},0}$ appearing in the right-hand side of the growth
equation (2).

We choose an exponential squared kernel as covariance function [28] such that

$f_{\mathrm{DE}}(z)\sim\mathcal{GP}\left(\bar{f}=1,k=K(\sigma_{f},\ell_{f})\right)$
(5)

where $\sigma_{f}$ and $\ell_{f}$ are the usual hyperparameters controlling
the deviations from the mean and the correlation length in dataspace,
respectively. We can then draw samples of $f_{\mathrm{DE}}(a)$ from this
multivariate (Gaussian) distribution. Substituting $f_{\mathrm{DE}}$ in (1)
provides us with a family of samples of $h(z)$ that we can use to solve for
the perturbations $\delta_{m}(z)$ using (2) (or equivalently (3)) without
assuming a particular DE model. This approach uses _untrained_ GP rather than
_trained_ GP, and calculate the goodness of fit of each sample. It is similar
to the approach in [29] and [19]. The advantage of using _untrained_ GP over
trained GP is that one can perform forward modeling, and therefore is not
restricted to the redshift range where the data are. This is particularly
interesting in our case: since we need to solve the growth equation (2) with
the correct initial conditions. By modeling $f_{\mathrm{DE}}$ as a GP with
mean function one, we enforce that at high-redshift, the expansion history is
EdS, and we can solve the growth equation. In this approach, rather than
finding the best-fit hyperparameters that maximize the marginal likelihood as
is commonly done[28], we use a more Bayesian approach and marginalize over the
hyperparameters [29, 19, 30]. For a detailed comparison between trained and
untrained GP and on the effect of marginalizing over the hyperparameters, see
[31].

We perform our analysis in increasing levels of complexity. We first fix the
cosmological parameters ($\Omega_{\text{m},0}$, $\sigma_{8,0}$) to their
fiducial value, and sample the (hyper)parameter space using Markov Chain Monte
Carlo (MCMC) methods as implemented in emcee [32]. We then include the
cosmological parameters one by one in the MCMC analysis to asses their
individual effect on the reconstructions. Table 1 shows the prior ranges used
throughout this work. The hyperparameters $(\sigma_{f},\,\ell_{f})$ will
determine the shape of $f_{\mathrm{DE}}$, which together with a given value of
$\Omega_{\text{m},0}$, determines the expansion history $h(z)$ through Eq.(1).
Solving Eq. (2) and translating into $f\sigma_{8}(z)$ further requires the
knowledge of $\sigma_{8,0}$ and a given shape of $G_{\mathrm{eff}}$$(k,z)$.
Again, in this work we restrict ourselves to the case of GR.

Figure 2: Left: reconstructed $f_{\mathrm{DE}}$ from various combinations of
the data, assuming perfect knowledge of $(\Omega_{\text{m},0},\sigma_{8,0})$.
Right: reconstructed $h(z)/(1+z)$. Solid lines and shaded regions correspond
to the median and 68% confidence intervals around it, respectively. The
dashed-black line correspond to the true (fiducial) cosmology used to generate
the mock data.

###  II.3 Data

####  II.3.1 Mocks data

To validate our method, we proceed as follows. We first generate mock data
from a fiducial CPL cosmology, with
$\theta^{\mathrm{fid}}=\\{\Omega_{\mathrm{m,0}}^{\mathrm{fid}}=0.28,w_{0}^{\mathrm{fid}}=-0.95,w_{a}^{\mathrm{fid}}=0.3,\sigma_{8,0}^{\mathrm{fid}}=0.81\\}$.
In this fiducial model one has [33, 34]

$\frac{\Omega_{\rm DE}}{\Omega_{m}}=\frac{\Omega_{\rm
DE,0}}{\Omega_{m,0}}(1+z)^{3(w_{0}+w_{a})}\exp{\left(-3w_{a}\frac{z}{1+z}\right)}~{}.$
(6)

Hence for the interesting range $-1<w_{0}+w_{a}<0$, while $\Omega_{\rm DE}\to
0$ in the past, it does so more slowly than if DE was in the form of a
cosmological constant $\Lambda$. As an illustration, we get for any given
values $\Omega_{\rm DE,0}$ and $\Omega_{m,0}$ at $z=99$ (chosen for
convenience)

$\frac{\Omega^{\rm fid}_{\rm DE}}{\Omega_{m}}(z=99)=\frac{\Omega_{\rm
DE,0}}{\Omega_{m,0}}~{}f^{\rm fid}_{99}\times 10^{-6}~{},$ (7)

where we have defined

$f_{99}\equiv f_{\rm DE}(z=99)~{},$ (8)

Using $f_{99}^{\rm fid}=51.64$, it is seen that the ratio $\frac{\Omega^{\rm
fid}_{\rm DE}}{\Omega_{m}}(z=99)$, while very tiny, is still larger by a
factor $51.64$ than it would be for a cosmological constant $\Lambda$. As we
shall later discuss in Section III.4, for the particular fiducial values
chosen here, this can lead to a biased determination of $\Omega_{\text{m},0}$
and $\sigma_{8,0}$. We generate Mock SNIa, BAO, and RSD data as follows:

  * •

SNIa: we generate Roman-like (formerly “WFIRST) data following Tab.7 in [35].
Due to the high-computational nature of the problem, we use this compressed
likelihood made of 9 points in $h(z)$. We note that these assume a flat
Universe. We will refer to this dataset as “SN”.

  * •

BAO: we simulate a DESI-like survey, following Tabs. 2.3-2.7 in [36]. More
specifically, we consider the $14,000\deg^{2}$ survey as baseline. The BAO
measurements are given in terms of the combinations $H(z)r_{\text{d}}$ and
$d_{A}(z)/r_{\text{d}}$. Since we do not include any early universe
information in our analysis, we express our $\chi^{2}$ in terms of the product
$H_{0}r_{\text{d}}$. For the method validation with mock data, and for
computational reasons, we fixed $H_{0}r_{d}$ to its fiducial value222We tested
including the combination $H_{0}r_{d}$ as a free parameter, finding that a
(Stage-IV) SN+BAO+RSD analysis is able to constrain the product to be within
the $1\sigma$ region of its true (fiducial) value.. In the case of (real)
current data, it is left as a free parameter, as we ignore its true value.

  * •

RSD: we generate $f\sigma_{8}$ measurements for a DESI-like survey, always
following Tabs. 2.3-2.7 in [36] from the fiducial cosmology. Note that these
require assuming a given $G_{\mathrm{eff}}$. We refer to these measurements as
“RSD”.

We note at this stage that the acoustic peaks in the CMB allow for a precise
determination of the quantity
$100\,\theta_{s}(z_{\star})=100\,r_{s}(z_{\star})/r(z_{\star})=1.04110\pm
0.00031$ where $r_{s}(z_{\star})$ is the comoving sound horizon and
$r(z_{\star})$ is the comoving distance at recombination. We use this
additional constraint at recombination, hence at redshifts much larger than
our numerical calculations, to check that all our models have an acceptable
behaviour at very large $z$. Using the Planck values $z_{\star}=1089.92$,
$r_{s}(z_{\star})=$144.4\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$$ and
$H_{0}=$70\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}\,\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$,
a calculation of the quantity $100\,\theta_{s}$ for the fiducial CPL model
considered here yields 1.0427 close to the one reported by Planck [37]. Since
we deal with mock data and we do not include any early-universe information in
the analysis, we expect that inclusion of data at high-$z$ will only confirm
and refine our results. Similar values are obtained for the rest of the
fiducial cosmologies considered in Section IV.

We performed a series of separate and joint analysis, using SN, BAO and RSD
data. In the next section, we illustrate our reconstructions for a fiducial
cosmology with a substantial departure from $\Lambda$CDM —depicted by the
dashed line in Fig. 2— to test a “worst-case” scenario. This enables us to
asses whether our framework allows for an accurate reconstruction of
$f_{\mathrm{DE}}$ and thus, of the expansion history $h(z)$.

##  III Validation of the Method with Mock Data

In this section, we apply our formalism to the mock data generated above in
order to validate the method in the case of GR. We first assume perfect
knowledge of $(\Omega_{\text{m},0},\sigma_{8,0})$, then we will only assume
$\sigma_{8,0}$, and finally we do not assume any parameter to reconstruct the
growth and expansion histories.

Figure 3: Posteriors of the kernel hyperparameter $\sigma_{f}$ vs.
$\Omega_{\text{m},0}$ from various chains—assuming perfect knowledge of
$\sigma_{8,0}$. It is seen that in all cases, those samples of $f_{\rm DE}$
with a value of $\Omega_{\rm m,0}=0.28$ (corresponding to the truth) require a
deviation from the mean (_i.e._ $\sigma_{f}\neq 0$). Furthermore, “background-
only” probes (SN & BAO) are consistent with no deviations from the mean
(_i.e._ $\sigma_{f}=0$) provided $\Omega_{\rm m,0}\sim 0.31$ \- close to
$\Lambda$CDM’s best-fit value. However, RSD are more sensitive to the
evolution of $f_{\mathrm{DE}}$ and do not allow for $f_{\mathrm{DE}}=1$
(_i.e._ $\sigma_{f}\neq 0$), though the value of $\Omega_{\text{m},0}$ is
essentially unconstrained.

###  III.1 ($\Omega_{\rm m,0}$, $\sigma_{8,0}$) known

The left panel of Fig. 2 shows the reconstructed $f_{\mathrm{DE}}$ assuming
perfect knowledge of $\Omega_{\text{m},0}$ and $\sigma_{8,0}$ using various
combinations of SNIa, BAO and RSD data. The lower panels throughout this paper
show the percentage error with respect to the fiducial cosmology used to
generate the data—for a given quantity X, we define $\frac{\Delta
X}{X}\equiv\left(\frac{X-X^{\rm fid}}{X^{\rm fid}}\right)\times 100$. The DE
evolution $f_{\mathrm{DE}}$ is well reconstructed at low-$z$ where the data
are present. At high-$z$, our choice of mean function $\bar{f}_{\rm DE}=1$
ensures a EdS behaviour as expected. However, at intermediate $1.5\lesssim
z\ll 20$, $f_{\mathrm{DE}}$ is not reconstructed well. This is understandable
because this redshift range is in the matter-dominated era, where the
constraining power on $f_{\mathrm{DE}}$ becomes negligible. However, because
of the knowledge of the fiducial $\Omega_{\text{m},0}$, the expansion history
is reconstructed relatively well as shown on the right-hand panel.

Figure 4: Upper panels: Percent errors in the reconstructed background
quantities $f_{\mathrm{DE}}$ (left) and $h(z)$ (right) assuming perfect
knowledge of $\sigma_{8,0}$, marginalized over $\Omega_{\text{m},0}$. Lower
panels: Reconstructed (perturbed) quantities
$\gamma\equiv\ln{f}/\,\ln{\Omega_{\rm m}}$ and $f\sigma_{8}$ from joint
analyses. Solid lines and shaded regions show the median and $1\sigma$
confidence regions around it, respectively. As shown in the upper panels, when
$\Omega_{\text{m},0}$ is free to vary, $f_{\mathrm{DE}}$ is not well
reconstructed with SN+BAO, as background-only probes can compensate for the
“lack of DE” at large-$z$ with larger values of $\Omega_{\text{m},0}$. This
bias in $\Omega_{\text{m},0}$ is reflected in the $f\sigma_{8}$
reconstructions shown in the bottom-right panel, being marginally consistent
with the fiducial solution. Adding data from RSD (in red) yields a better
reconstruction of the quantities of interest.

In this particular choice of fiducial cosmology, $f_{\mathrm{DE}}$ grows with
redshift, and thus contributes more to the expansion rate of the universe at
large-$z$, compared to $\Lambda$CDM. Our choice of mean function $\bar{f}_{\rm
DE}=1$ induces a slight bias towards $\Lambda$CDM, as in the absence of data
$f_{\mathrm{DE}}$ will go back to its mean, which explains why $h(z)$ is
underestimated at intermediate redshifts. At higher-redshifts, this bias
becomes negligible and an accurate estimation of $\Omega_{\text{m},0}$ is more
relevant, since $h(z\gg
1)\sim\left[\Omega_{\text{m},0}(1+z)^{3}\right]^{1/2}$. However, this bias
induced in $h(1.5<z<25)$ is below $2\%$.

###  III.2 $\sigma_{8,0}$ known

We now study the effect of releasing the assumption of perfect knowledge of
$\Omega_{\text{m},0}$. Instead, we now vary $\Omega_{\text{m},0}$ together
with the hyperparameters in the MCMC runs. The reconstructions are shown in
Fig. 4. Without perfect knowledge of $\Omega_{\text{m},0}$, the uncertainties
on the reconstructions increase, as expected. The SN+BAO reconstructions
underestimate $f_{\mathrm{DE}}$ while the RSD-only reconstruction overestimate
it for $z\leq 1$. Combining the three probes results in a more accurate and
precise reconstruction of $f_{\mathrm{DE}}$ on this redshift range.
Interestingly, while the SN+BAO reconstructions of $h$ are overestimated at
high-redshift, and the RSD-only reconstructions overestimate $h$ over the
whole redshift range, the reconstruction of $h$ with the three probes is
accurate and precise over the whole range of redshifts. We also note that the
relative errors for $f\sigma_{8}(z)$ measurements are larger than those of BAO
and SNIa at low-$z$, which further explains the larger uncertainty at low-$z$.

Figure 5: Residuals from a joint analysis, marginalized over
$\Omega_{\text{m},0}$ and $\sigma_{8,0}$, assuming GR. We also compare the
reconstructions with and without $f_{99}\equiv f_{\rm DE}(z=99)$— _c.f._
Section III.4. Although the uncertainties increase significantly, the
reconstructions including $f_{99}$ are in $1\sigma$ agreement with the
fiducial solutions throughout the redshift range $0<z\lesssim 10^{2}$.

In the lower panels of Fig. 4, we show the reconstructions of the growth index
$\gamma(z)={\ln f}/{\ln\Omega_{\text{m}}}$ and $f\sigma_{8}(z)$ when
marginalizing over $\Omega_{\text{m},0}$, assuming GR. The dashed line shows
the fiducial solution, for which $\gamma\simeq 0.55$, with a slight redshift
evolution—see [38, 39, 40]. As explained before, the background probes
(SN+BAO) are not able to constrain the value of $\Omega_{\text{m},0}$ (and
thus to accurately reconstruct $f_{\mathrm{DE}}$—see upper panels of Fig.4)
which in turn is reflected in the reconstructed growth history $f\sigma_{8}$,
particularly at low-$z$. When including RSD, and assuming the fiducial value
of $\sigma_{8,0}$, the reconstructions (in orange) are much more accurate,
both in $\gamma$ and $f\sigma_{8}$. From the well-known relation
$f\simeq\Omega_{\text{m}}^{\gamma}$, we see that larger values of
$\Omega_{\text{m},0}$ will require lower values of $\gamma$ in order to keep
the agreement with observations, as reflected in the bottom-left panel of
Fig.4 for the case of SN+BAO.

The bias in the background-only (SN+BAO) reconstructions is understood by
looking at the posterior distributions of $(\sigma_{f},\Omega_{\text{m},0})$,
shown in Fig. 3. Each point in parameter space correspond to a given
realization of $f_{\mathrm{DE}}$ and expansion history characterized by
certain combinations of $(\sigma_{f},\ell_{f})$ and
$\Omega_{\text{m},0}$—while $\sigma_{8,0}$ is fixed to its fiducial value. For
the background-only probes (SN+BAO), the data are consistent with no
deviations from the mean ($\sigma_{f}=0$), provided a high value of
$\Omega_{\text{m},0}$$\sim 0.31$, while those samples with the true (fiducial)
value of $\Omega_{\text{m},0}$$=0.28$ require $\sigma_{f}\neq 0$. However,
redshift-space distortions are more sensitive to the nature of DE, and thus do
not allow for a $\Lambda$CDM-like evolution ($f_{\mathrm{DE}}=1$) and require
$\sigma_{f}\neq 0$ \- at the price of losing the constraining power on
$\Omega_{\text{m},0}$. As mentioned before, when combining the three different
datasets, these effects compensate and we get an accurate reconstruction of
$f_{\mathrm{DE}}$, and a clearer detection of $\sigma_{f}\neq 0$ (_i.e._ a
deviation from the mean—$\Lambda$CDM), by narrowing down the allowed range of
$\Omega_{\text{m},0}$.

###  III.3 ($\Omega_{\rm m,0}$, $\sigma_{8,0}$) unknown

We now release the assumption of perfect knowledge of $\Omega_{\text{m},0}$
and $\sigma_{8,0}$ and marginalize over them. Fig. 5 shows the relative errors
in the reconstructions of $f_{\mathrm{DE}}$ and $h$ after marginalizing over
the two cosmological parameters and the two hyperparameters. As previously,
$f_{\mathrm{DE}}$ cannot be reconstructed accurately at $z\geq 1$. The RSD-
only reconstruction underestimate $f_{\mathrm{DE}}$ from low redshift, and the
uncertainties are quite large. The smaller uncertainties in the combined
probes case yields a reconstruction within $2\sigma$ of the expansion history,
with a slight overestimation at high-redshifts. This stems from the fact that
we are biased towards $\Lambda$CDM’s best-fit value of
$\Omega_{\text{m},0}$$\sim 0.31$. Furthermore, since $\sigma_{8,0}$ is now a
free parameter, and is degenerate with $\Omega_{\text{m},0}$, we loose the
ability to “narrow down” $\Omega_{\text{m},0}$ to its “true” (fiducial) value
and thus, to accurately reconstruct $f_{\mathrm{DE}}$. This bias will be
reflected on the growth evolution $f\sigma_{8}$, which depends not only on
$h(z)$ but also on $\Omega_{\text{m}}(z)$ via the source term for the
perturbations $\delta$\- see lower panels of Fig. 5, particularly the left
panel showing $\gamma(z)$.

The lower panels of Fig. 5 also show the reconstructed perturbed quantities
when marginalizing over both $\Omega_{\text{m},0}$ and $\sigma_{8,0}$. In this
case, as understood from the background reconstructions in the upper panel,
the degeneracy between $\Omega_{\text{m},0}$ and $\sigma_{8,0}$ causes the
reconstructions of $f_{\mathrm{DE}}$ to be less accurate, even in the joint
SN+BAO+RSD case. From the $h(z)$-data standpoint, having “less” Dark Energy
($f_{\mathrm{DE}}<f_{\mathrm{DE}}^{\rm fid}$) at large-$z$ can be compensated
by having more matter $\Omega_{\text{m},0}>\Omega_{\text{m},0}^{\rm fid}$.
Which enforces $\sigma_{8,0}<\sigma_{8,0}^{\rm fid}$, while keeping the
agreement with RSD observations. Thus, when $\sigma_{8,0}$ is also free to
vary, we lose the constraining power on $\Omega_{\text{m},0}$ provided by RSD
measurements and therefore we do not reconstruct $f_{\mathrm{DE}}(z\geq 0.6)$.
The $2\%$ percent bias induced in $h(z\geq 2)$ translates into a bias in
$f\sigma_{8}(z\geq 2)$ as can be seen in the lower-right panel of Fig. 5.
Despite the $2\sigma$ agreement with both $h(z)$ and $f\sigma_{8}(z)$
measurements, the bias in the values of $\Omega_{\text{m},0}$ and
$\sigma_{8,0}$ are clearly seen in the evolution of $\gamma(z)$ — shown on the
left panels of Fig. 5.

Figure 6: Posterior distributions for the cosmological parameters in a
combined analysis. By allowing for a non-negligible contribution from DE at
large-$z$, encoded in $f_{99}$, we remove the bias seen before. Dashed lines
indicate the fiducial values.

###  III.4 Removing the Biases

As discussed before, assuming $\Lambda$CDM (_i.e. $f_{\mathrm{DE}}=1$_) at
high-$z$ is a reasonable assumption but can induce a bias in the
reconstructions if the contribution from DE is non-negligible at such
redshifts. As can be seen from Fig. 6, our contours are biased and shifted
towards $\Lambda$CDM’s best-fit value - _i.e._ a higher $\Omega_{\text{m},0}$
and lower $\sigma_{8,0}$ to compensate for the lack of DE. To deal with such a
bias, we introduce an additional degree of freedom in the MCMC runs. Namely,
$f_{99}\equiv f_{\mathrm{DE}}(z=99)$ characterizing the DE density at $z\sim
10^{2}$, when we start our integration. More specifically, for each step in
the MCMC, we impose a value for $f_{\mathrm{DE}}$ at both ends of the redshift
range, namely $f_{\mathrm{DE}}=1$ at $z=0$ and $f_{\mathrm{DE}}=f_{99}$ at
$z=99$—while $f_{99}$ is free to vary in the range $f_{99}\in[0.01,300]$. The
behaviour of $f_{\mathrm{DE}}$ between $0<z<99$ is solely determined by the
data at low-redshift, without biasing our parameter inferences by assuming
$\Lambda$ in the past—particularly those of $\Omega_{\text{m},0}$ and
$\sigma_{8,0}$. This enables us to cover a wide variety of DE models,
including those where the DE density grows with redshift and contributes as
much as $\sim 0.1\%$ of the matter density at $z=10^{2}$. The upper bound in
our prior is motivated by the fact that observations require an epoch of
matter domination for structures to efficiently form at early times
($\delta\propto a(t)$), and the latest Pantheon+ analysis [41] on CPL
cosmologies constrain $-1.1\lesssim w_{0}\lesssim-0.6$ and $-1.3\lesssim
w_{a}\lesssim 0.2$. We note that for our choice of fiducial cosmology (_i.e.
$w_{0}=-0.95,w_{a}=0.3$_), the normalized energy density of DE is
$f_{\mathrm{DE}}(z=99)\simeq 52$.

Figure 7: Reconstructions of various (different) fiducial cosmologies. Solid lines and shaded regions correspond to median and $95\%$ C.L. around it, respectively. Dashed-colored lines in the upper panels depict the “true” behaviour of $f_{\mathrm{DE}}(z)$ for the different cosmologies. In the lower panels, the corresponding observables $H(z)$ and $f\sigma_{8}(z)$ normalized to the correct solution in the fiducial cosmology. Datasets | $\Omega_{m,0}$ | $\sigma_{8,0}$ | $f_{99}$
---|---|---|---
SN+BAO+RSD | $0.2920^{+0.0073}_{-0.0077}$ | $0.7972^{+0.008}_{-0.007}$ | —
SN+BAO+$f_{99}$ | $0.283^{+0.019}_{-0.023}$ | — | $21^{+83}_{-17}$
SN+BAO+RSD+$f_{99}$ | $0.282^{+0.011}_{-0.012}$ | $0.806^{+0.013}_{-0.011}$ | $15^{+45}_{-11}$
Table 2: Forecasted constraints on cosmological parameters from joint analyses
using data from Stage-IV surveys. We remind the reader that the fiducial
values used to generate the data are $\Omega_{\text{m},0}^{\rm fid}=0.28$,
$\sigma_{8,0}^{\rm fid}=0.81$ and $f_{99}^{\rm fid}=51.64$.

In Fig. 5 we compare the reconstructions with and without $f_{99}$. The biased
case (in dark-red) is in $1\sigma$ agreement with the background observable
$h(z)$ and within $2\sigma$ agreement with the growth observable
$f\sigma_{8}(z)$, despite the clear discrepancy with the true DE evolution at
$z\gtrsim 1$, depicted in the left panels. The dark gray line in Fig. 5 shows
the reconstructions when including the $f_{99}$ parameter in the runs,
increasing the uncertainty on the reconstructions but being consistent with
$f_{\mathrm{DE}}^{\rm fid}$ throughout the redshift range.

##  IV Reconstructing various Fiducial Cosmologies

In this section, we apply our methodology to reconstruct other (different)
cosmologies. In particular, we focus on cosmologies where the behaviour of
$f_{\mathrm{DE}}$ is different than the one presented in the previous section
of this paper. First, we reconstruct a CPL model with the best-fit parameters
as obtained in the latest joint Pantheon+ analysis [41] (including BAO and
Planck data). We also consider two additional examples, namely the
Phenomenological Emergent Dark Energy (PEDE) [42] and the Transitional Dark
Energy (TDE) [43] models.

  * •

$\rm CPL^{\rm\star}$: We take the best-fit values for the CPL cosmology in the
latest Pantheon+ analysis [41]. These correspond to $\theta^{\rm
fid}=\\{\Omega_{\text{m},0}^{\rm fid}=0.31,w_{0}^{\rm fid}=-0.85,w_{a}^{\rm
fid}=-0.6\\}$ and we fix $\sigma_{8,0}^{\rm fid}=0.81$. We denote this $\rm
CPL^{\rm\star}$ to distinguish it from the fiducial cosmology presented in
§III.

  * •

PEDE: Similarly, for the PEDE model we fix $\theta^{\rm
fid}=\\{\Omega_{\text{m},0}^{\rm fid}=0.31,\sigma_{8,0}^{\rm fid}=0.81\\}$
with $w(z)$ given by Eq. (6) in [42]

  * •

TDE: For the TDE model [43], the cosmological parameters are $\theta^{\rm
fid}=\\{\Omega_{\text{m},0}^{\rm fid}=0.31,w_{0}^{\rm fid}=-0.8,w_{1}^{\rm
fid}=-2.0,z_{t}^{\rm fid}=1.0,\Delta_{z}^{\rm fid}=0.2,\sigma_{8,0}^{\rm
fid}=0.81\\}$ with the expression for $w(z)$ given by Eq. (4.2) in [44].

The $f_{\mathrm{DE}}(z)$ and $h(z)$ reconstructions for the corresponding
cosmologies are shown in Fig.7, when marginalizing over the cosmological+hyper
parameters and using SN+BAO+RSD data.

The $f_{\mathrm{DE}}$ (and thus $h(z)$) reconstructions in the CPL and PEDE
cosmologies are accurate within $1\sigma$ across the entire redshift range
$0<z<4$. In such cosmologies, the energy density of DE smoothly decreases with
redshift, and as such, contribute even less than the $\Lambda$-term to the
energy budget of the universe at $z\sim 10^{2}$, and thus our assumption of an
EdS universe at such redshift is well-justified. In the case of larger
departures from $\Lambda$CDM, as in the case of TDE shown in red, the
reconstructions are able to capture the “increasing” trend in
$f_{\mathrm{DE}}$ at low-$z$ ($z\lesssim 1$), where the data are abundant and
DE dominates, but tend to go back to the mean ($f_{\mathrm{DE}}\sim 1$) where
the constraining power on DE is lost—both because the quality of the data
decreases, and because matter becomes the dominant component. This translates
into $\lesssim 2\%$ deviations in terms of $h(z)$, while being consistent
within $2\sigma$ with the fiducial expansion and growth histories.
Interestingly, in all cases the $f_{\mathrm{DE}}$ reconstructions are in
tension—at $\gtrsim 2\sigma$—with $\Lambda$CDM ($f_{\mathrm{DE}}\equiv 1$) for
$z\lesssim 1$.

##  V Discussions and Conclusions

In this work a method is presented for the joint reconstruction of the cosmic
expansion and growth histories up to large redshifts $z\lesssim 100$ using
Gaussian processes. When using Gaussian processes, the problem to be solved is
the scarcity or even the total absence of data for $z\gtrsim 2$. On the other
hand as we expect most DE models to behave like a quasi Einstein de Sitter
universe on these intermediate redshifts, we chose to model
$f_{\mathrm{DE}}\sim\mathcal{GP}$ with mean function 1 while only assuming
flat FLRW+GR. This is a sufficient condition to impose an EdS behaviour at
high-redshifts. These reasonable assumptions leave the DE evolution to be
solely determined by the data at low-$z$. Hence we emphasize the flexibility
of our method with respect to the possible DE phenomenologies at low-$z$ where
DE dominates and where specific models leave their imprint in observations. We
tested the efficiency of our method in jointly reconstructing the cosmic
expansion and growth histories by simulating mock data for the upcoming
generation of cosmological surveys and for various dark energy models. In §III
we extensively discuss the reconstruction of a fiducial CPL cosmology with
$w_{0}=-0.95$ and $w_{a}=0.3$ (hence with a rather strong evolution of $w$
upwards for growing $z$) which is at the limit of being excluded by current
data (_e.g._ Pantheon+ [41]), in other words we choose to test our method with
a “worst-case” scenario which departs substantially from a cosmological
constant at large redshifts while still obeying an EdS like behaviour at
high-$z$. The successful reconstruction of this fiducial cosmology implemented
through the introduction of the parameter $f_{99}$ means that it will work as
well, perhaps even without making use of $f_{99}$, for most DE models
including Phenomenological Emergent Dark Energy (PEDE) or Transitional Dark
Energy (TDE) models (see Fig. 7) which generically do not contribute much to
the energy budget at high-$z$. We have shown that $f_{\rm DE}$ reconstruction
from stage-IV surveys can potentially detect deviations from $\Lambda$CDM at
more than $2\sigma$. We believe the joint reconstruction method presented
here, which was refined in order to capture a small amount of DE at high
redshifts when $\Omega_{\mathrm{DE}}\to 0$ though not as quickly as in
$\Lambda$CDM, can address such tracking (non-interacting smooth) DE as well.
We conjecture that by including information at higher-$z$ (_e.g._ CMB distance
priors), we could potentially “constrain” the contribution of (early) DE at
$z\sim z_{\rm rec}$ by further constraining $f_{99}$. We also expect that our
approach can be extended to more cosmologies, and that a non-vanishing spatial
curvature or the effect of massive neutrinos can be constrained accurately. We
finally expect interesting results can be obtained for modified gravity DE
models and we leave all these developments for future work.

## Acknowledgements

BL acknowledges the support of the National Research Foundation of Korea
(NRF-2019R1I1A1A01063740) and the support of the Korea Institute for Advanced
Study (KIAS) grant funded by the government of Korea. AS would like to
acknowledge the support by National Research Foundation of Korea
NRF2021M3F7A1082053, and the support of the Korea Institute for Advanced Study
(KIAS) grant funded by the government of Korea. AAS was partly supported by
the project number 0033-2019-0005 of the Russian Ministry of Science and
Higher Education.

## References

  * Weinberg [1989] S. Weinberg, The Cosmological Constant Problem, Rev. Mod. Phys. 61, 1 (1989).
  * Sahni and Starobinsky [2000] V. Sahni and A. A. Starobinsky, The Case for a positive cosmological Lambda term, Int. J. Mod. Phys. D 9, 373 (2000), arXiv:astro-ph/9904398 .
  * Bullock and Boylan-Kolchin [2017] J. S. Bullock and M. Boylan-Kolchin, Small-scale challenges to the $\Lambda$CDM paradigm, Annual Review of Astronomy and Astrophysics 55, 343 (2017), https://doi.org/10.1146/annurev-astro-091916-055313 .
  * Perivolaropoulos and Skara [2021] L. Perivolaropoulos and F. Skara, Challenges for $\Lambda$CDM: An update (2021), arXiv:2105.05208 [astro-ph.CO] .
  * Riess  _et al._ [2022] A. G. Riess, W. Yuan, L. M. Macri, D. Scolnic, D. Brout, S. Casertano, D. O. Jones, Y. Murakami, L. Breuval, T. G. Brink, A. V. Filippenko, S. Hoffmann, S. W. Jha, W. D. Kenworthy, J. Mackenty, B. E. Stahl, and W. Zheng, A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km/s/Mpc Uncertainty from the Hubble Space Telescope and the SH0ES Team (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Di Valentino  _et al._ [2021a] E. Di Valentino _et al._ , Snowmass2021 - letter of interest cosmology intertwined II: The Hubble constant tension, Astroparticle Physics 131, 102605 (2021a).
  * Di Valentino  _et al._ [2021b] E. Di Valentino _et al._ , Cosmology intertwined III: $f\sigma_{8}$ and $S_{8}$, Astroparticle Physics 131, 102604 (2021b).
  * Di Valentino  _et al._ [2021c] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, In the realm of the Hubble tension—a review of solutions, Classical and Quantum Gravity 38, 153001 (2021c).
  * Clifton  _et al._ [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Modified gravity and cosmology, Physics Reports 513, 1–189 (2012).
  * Starobinsky [1998] A. A. Starobinsky, How to determine an effective potential for a variable cosmological term, JETP Lett. 68, 757 (1998), arXiv:astro-ph/9810431 .
  * Shafieloo  _et al._ [2006] A. Shafieloo, U. Alam, V. Sahni, and A. A. Starobinsky, Smoothing supernova data to reconstruct the expansion history of the Universe and its age, MNRAS 366, 1081 (2006), astro-ph/0505329 .
  * Shafieloo [2007] A. Shafieloo, Model-independent reconstruction of the expansion history of the Universe and the properties of dark energy, MNRAS 380, 1573 (2007), astro-ph/0703034 .
  * Alam  _et al._ [2009] U. Alam, V. Sahni, and A. A. Starobinsky, Reconstructing Cosmological Matter Perturbations using Standard Candles and Rulers, Astrophys. J. 704, 1086 (2009), arXiv:0812.2846 [astro-ph] .
  * Shafieloo and Clarkson [2010] A. Shafieloo and C. Clarkson, Model independent tests of the standard cosmological model, Phys. Rev. D 81, 083537 (2010), arXiv:0911.4858 [astro-ph.CO] .
  * Shafieloo  _et al._ [2012] A. Shafieloo, A. G. Kim, and E. V. Linder, Gaussian process cosmography, Phys. Rev. D 85, 123530 (2012), arXiv:1204.2272 [astro-ph.CO] .
  * Shafieloo  _et al._ [2013] A. Shafieloo, A. G. Kim, and E. V. Linder, Model independent tests of cosmic growth versus expansion, Phys. Rev. D 87, 023520 (2013), arXiv:1211.6128 [astro-ph.CO] .
  * L’Huillier and Shafieloo [2017] B. L’Huillier and A. Shafieloo, Model-independent test of the FLRW metric, the flatness of the Universe, and non-local estimation of $H_{0}r_{\mathrm{d}}$, J. Cosmology Astropart. Phys 1, 015 (2017), arXiv:1606.06832 .
  * L’Huillier  _et al._ [2018] B. L’Huillier, A. Shafieloo, and H. Kim, Model-independent cosmological constraints from growth and expansion, MNRAS 476, 3263 (2018), arXiv:1712.04865 .
  * L’Huillier  _et al._ [2020] B. L’Huillier, A. Shafieloo, D. Polarski, and A. A. Starobinsky, Defying the laws of gravity I: model-independent reconstruction of the Universe expansion from growth data, MNRAS 494, 819 (2020), arXiv:1906.05991 [astro-ph.CO] .
  * Raveri  _et al._ [2021] M. Raveri, L. Pogosian, K. Koyama, M. Martinelli, A. Silvestri, G.-B. Zhao, J. Li, S. Peirone, and A. Zucca, A joint reconstruction of dark energy and modified growth evolution (2021), arXiv:2107.12990 [astro-ph.CO] .
  * Arjona  _et al._ [2021] R. Arjona, A. Melchiorri, and S. Nesseris, Testing the $\Lambda$CDM paradigm with growth rate data and machine learning, ArXiv e-prints (2021), arXiv:2107.04343 [astro-ph.CO] .
  * Nesseris  _et al._ [2021] S. Nesseris _et al._ , Euclid: Forecast constraints on consistency tests of the $\Lambda$CDM model, ArXiv e-prints (2021), arXiv:2110.11421 [astro-ph.CO] .
  * Ade  _et al._ [2016] P. A. R. Ade _et al._ (Planck), Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Boisseau  _et al._ [2000] B. Boisseau, G. Esposito-Farèse, D. Polarski, and A. A. Starobinsky, Reconstruction of a Scalar-Tensor Theory of Gravity in an Accelerating Universe, Phys. Rev. Lett. 85, 2236 (2000), arXiv:gr-qc/0001066 [gr-qc] .
  * Holsclaw  _et al._ [2010a] T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib, and D. Higdon, Nonparametric reconstruction of the dark energy equation of state, Phys. Rev. D 82, 103502 (2010a), arXiv:1009.5443 [astro-ph.CO] .
  * Holsclaw  _et al._ [2010b] T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib, and D. Higdon, Nonparametric Dark Energy Reconstruction from Supernova Data, Phys. Rev. Lett. 105, 241302 (2010b), arXiv:1011.3079 [astro-ph.CO] .
  * Holsclaw  _et al._ [2011] T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib, and D. Higdon, Nonparametric reconstruction of the dark energy equation of state from diverse data sets, Phys. Rev. D 84, 083501 (2011), arXiv:1104.2041 [astro-ph.CO] .
  * Rasmussen and Williams [2006] C. Rasmussen and C. Williams,  _Gaussian Processes for Machine Learning_, Adaptative computation and machine learning series (University Press Group Limited, 2006).
  * Joudaki  _et al._ [2018] S. Joudaki, M. Kaplinghat, R. Keeley, and D. Kirkby, Model independent inference of the expansion history and implications for the growth of structure, Phys. Rev. D 97, 123501 (2018), arXiv:1710.04236 [astro-ph.CO] .
  * Ruiz-Zapatero  _et al._ [2022] J. Ruiz-Zapatero, C. García-García, D. Alonso, P. G. Ferreira, and R. D. P. Grumitt, Model-independent constraints on $\Omega_{m}$ and $H(z)$ from the link between geometry and growth, ArXiv e-prints (2022), arXiv:2201.07025 [astro-ph.CO] .
  * Hwang  _et al._ [2022] S.-g. Hwang, B. L’Huillier, R. Keeley, A. Shafieloo, and M. J. Jee (2022).
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, Publications of the Astronomical Society of the Pacific 125, 306–312 (2013).
  * Chevallier and Polarski [2001] M. Chevallier and D. Polarski, Accelerating Universes with Scaling Dark Matter, Int. J. Mod. Phys. D 10, 213 (2001), gr-qc/0009008 .
  * Linder [2003] E. V. Linder, Exploring the Expansion History of the Universe, Physical Review Letters 90, 091301 (2003), astro-ph/0208512 .
  * Riess  _et al._ [2018] A. G. Riess _et al._ , Type Ia Supernova Distances at Redshift $>1.5$ from the Hubble Space Telescope Multi-cycle Treasury Programs: The Early Expansion Rate, Astrophys. J. 853, 126 (2018), arXiv:1710.00844 [astro-ph.CO] .
  * Aghamousa  _et al._ [2016] A. Aghamousa _et al._ (DESI), The DESI Experiment Part I: Science,Targeting, and Survey Design, ArXiv e-prints (2016), arXiv:1611.00036 [astro-ph.IM] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Polarski  _et al._ [2016] D. Polarski, A. A. Starobinsky, and H. Giacomini, When is the growth index constant? (2016), arXiv:1610.00363 [astro-ph.CO] .
  * Calderon  _et al._ [2019] R. Calderon, D. Felbacq, R. Gannouji, D. Polarski, and A. A. Starobinsky, Global properties of the growth index of matter inhomogeneities in the universe, Phys. Rev. D 100, 083503 (2019), arXiv:1908.00117 [astro-ph.CO] .
  * Calderon  _et al._ [2020] R. Calderon, D. Felbacq, R. Gannouji, D. Polarski, and A. A. Starobinsky, Global properties of the growth index: mathematical aspects and physical relevance, Phys. Rev. D 101, 103501 (2020), arXiv:1912.06958 [astro-ph.CO] .
  * Brout  _et al._ [2022] D. Brout _et al._ , The Pantheon$+$ Analysis: Cosmological Constraints, ArXiv e-prints (2022), arXiv:2202.04077 [astro-ph.CO] .
  * Li and Shafieloo [2019] X. Li and A. Shafieloo, A simple phenomenological emergent dark energy model can resolve the hubble tension, The Astrophysical Journal 883, L3 (2019).
  * Keeley  _et al._ [2019] R. E. Keeley, S. Joudaki, M. Kaplinghat, and D. Kirkby, Implications of a transition in the dark energy equation of state for the $h_{0}$ and $\sigma_{8}$ tensions, Journal of Cosmology and Astroparticle Physics 2019 (12), 035–035.
  * Koo  _et al._ [2022] H. Koo, R. E. Keeley, A. Shafieloo, and B. L’Huillier, Bayesian vs frequentist: Comparing bayesian model selection with a frequentist approach using the iterative smoothing method (2022), arXiv:2110.10977 [astro-ph.CO] .


SIS/R model on Bi-Uniform hypergraphs
Rohan Nuckchady

Abstract

This report is based on the work in (1). We first review definitions and
notation developped there and provide derivations for the exact mathematical
description of an SIS epidemic on a hypergraph. We then generalise the work in
(1) to a new class of models that en-compass SIS and SIR models. The exact
differential equations are derived for the expected values of the population
of each state. Focusing on Bi-uniform hypergraphs, we make suitable
approximations obtain numerical solutions to those equations. These are
compared with stochastic simulations of the model for various systems.

Introduction

Models for epidemic propagation have usually been made on standard graphs.
However these assume pair-wise relationships between nodes, where as in
reality nodes interact in groups of more than 2. In order to account for this,
network based models often contain complete subgraphs to model those larger
groups, however these lead to a loss of information in the system. Recently,
work was done in (1), (2) to develop epidemic models on hypergraphs. Primarily
SIS models were derived however these do not always provide a good model of
epidemic propagation. In particular the type of hypergraph and relative sizes
and number of hyper edges should heavily affect the reliability of the model.
Often, epidemics propagate in a none-binary way-i.e there are more than two
possible states per node. Such a model on graphs has been heavily studied in
literature and is known as the SIR model. However, these assume that I-type
nodes cannot transition back to S-type nodes. Hence these do not properly
model epidemics where immunity is only gradually gained, or where the epidemic
is deadly, infectious and where immunity is not guaranteed. This motivates the
need to develop a model for these classes of epidemics.

We begin by introducing basic ideas and definitions behind the model and
reviewing the SIS model in (1). We provide an exact mathematical description
of the SIS model and derivations of the equations for the means of size of the
population of each node type. Similar methods are then applied to derive
equations for our new model (SISR). Our model will depend on two parameters
and by adjusting their values both the SIS and SIR models can be recovered.
Next, we describe the simulation methods used to produce stochastic
simulations of the model on Bi-uniform hypergraphs. We also make a mean-field
approximation valid for this class of hypergraphs. Finally, we compare the
resultsof the stochastic simulations to numerical solutions to our system of
equations.

Preliminaries

In this section we provide some basic definitions and develop the the SIS and
SISR model on hypergraphs.

Definition 1. $V$ is a set and its elements are called nodes. $E$ is another
set whose elements consists of subsets of $V$ . The pair $(V,E)$ is called a
hypergraph with vertices $v\in V$ and hyperedges $e\in E$.

Thus in the case where $|e|=2$ for all $e\in E$, we recover a graph. In the
context of epidemic propagation, nodes on the hypergraph represent individuals
and the hyperedges they are contained in represent the communities they are a
member of. In the SIS model, each node can be in 2 possible states: i)Infected
ii)Susceptible. There are two processes that are occuring in parallel on the
hypergraph: i)Recovery ii) Infection. Recovery is the process by which an
infected node becomes susceptible again and infection is the converse.

Definition 2. $S$ denotes a susceptible node and $I$ an infected node.
Infection and Recovery are modelled by poisson processes with parameters
$r_{S}$ (node dependent) and $\gamma$ respectively such that the probability
infection and recovery in a small amount of time $\delta t$ are given by
$1-\exp\left(-r_{S}\delta t\right)$ and $1-\exp(-\gamma\delta t)$
respectively.

Note with the above definition, $\gamma$ is taken to be independent of the
nodes in the environment. It remains however to define how the environment
(i.e. the hyperedges) affect the infection rate of a particular node. One
could approximate a hypergraph by replacing each hyperedge with a complete
subgraph. In doing so, the infection rate of a susceptible node would be the
some $f(n)\tau$ , where $\tau$ is the infection rate per node, $n$ is the
number of infected neighbours of the susceptible node and $f$ is an arbitrary
function. In the case of a hypergraph though, each hyperedge containing the
node provides its own infection pressure and thus we provide a different
definition.

Definition 3. Let $S$ be a susceptible node and $H_{S}=\\{h\in$ $E\mid S\in
h\\}$. Let also $f:\mathbb{Z}\rightarrow\mathbb{R}$ and $\tau$ be the
infection rate per node such that if a susceptible node has $n$ infected
neighbours in a particular hyperedge, its infection rate is given by
$f(n)\tau$. For $S$, the total infection rate $r_{S}$ is given by

$r_{S}=\sum_{h\in H_{S}}f\left(n_{h}\right)$

Where $n_{h}$ is the number of infectious nodes in $h\in E$.

Exact description of SIS model on hypergraphs

The state space for hypergraph with $|V|=N$ is $\\{S,I\\}^{N}$ containing
$2^{N}$ elements. We partition this state space into $N+1$ classes.

Definition 4. A class $S^{k}$ corresponds to all $\left(\begin{array}[]{l}N\\\
k\end{array}\right)$ states containing $k$ infected nodes. $\left\\{S^{k}\mid
0\leq k\leq N\right\\}$ is a partition of $\\{S,I\\}^{N}\cdot S_{i}^{k}\in
S^{k}$ are indered elements of $S^{k}$ with $1\leq i\leq$
$\left(\begin{array}[]{l}N\\\
k\end{array}\right)=c_{k}.S_{i}^{k}(j)\in\\{S,I\\}$ denotes the state of the
$j^{\text{th }}$ node.

In this language we describe the processes of infection and recovery
transitions:
i) Infection: An $S_{j}^{k}\rightarrow s_{i}^{k+1}$ transition, such that for
some $p\in\\{1,\ldots,N\\}$, we have $S_{j}^{k}(p)=S,S_{i}^{k}(p)=I$ and
otherwise $S_{j}^{k}(m)=S_{i}^{k-1}(m)$ for all $m\neq
p,m\in\\{1,\ldots,N\\}.$ Also, there is a $q\neq p$ such that $S_{j}^{k}(q)=I$
and $p,q$ are both contained in the same hyperedge. ii)Recovery: An
$S_{j}^{k}\rightarrow S_{i}^{k-1}$ transition, such that there is some
$p\in\\{1,\ldots,N\\}$ with $S_{j}^{k}(p)=I,S_{i}^{k-1}(p)=S$ and for all
$m\neq p,m\in\\{1,\ldots,N\\}$ we have $S_{j}^{k}(m)=s_{i}^{k-1}(m)$

These give rise to a linear system of differential equations for the vectors
of probabilities $X^{k}(t)=\left(X_{1}^{k},\ldots,X_{r_{k}}^{k}\right)$. Each
$X_{i}^{k}$ denotes the probability that the system is in the state
$S_{i}^{k}$ at time $t$. In particular this differential equation can be
written as:

$\frac{dX^{k}}{dt}=A^{k}X^{k-1}+B^{k}X^{k}+C^{k}X^{k+1}$

$A^{k}$ describes recovery, $B^{k}$ describes the behaviour if nothing happens
and $C^{k}$ describes infection. These can be written in the form of a
Kolomogrov equation as $\frac{dX}{dt}=PX$, where $P$ is in block tri-diagonal
form:

$P=\left(\begin{array}[]{ccccccc}B^{0}&C^{0}&0&0&0&\cdots&0\\\
A^{1}&B^{1}&C^{1}&0&0&\cdots&\vdots\\\ 0&A^{2}&B^{2}&C^{2}&0&\cdots&\vdots\\\
0&0&A^{3}&B^{3}&C^{3}&\vdots&\vdots\\\
\vdots&\vdots&\cdots&\cdots&\cdots&\vdots&\vdots\\\
0&0\ldots&\cdots&\cdots&\cdots&A^{n}&B^{N}\end{array}\right)$

The matrix element $A_{ij}^{k}$ is the rate of transition from $S_{j}^{k-1}$
to $S_{i}^{k}$ and is therefore given by:

$A_{ij}^{k}=\begin{cases}\tau\sum_{h_{i\in
h}h}f\left(N_{h}\left(S_{j}^{k-1}\right)\right)&\text{ if an
}S_{j}^{k-1}\rightarrow S_{1}^{k}\\\ &\text{ transition exists and }l\text{ is
the }\\\ 0&\text{ newly infected node }\\\ &\text{ otherwise }\end{cases}$

With $N_{h}$ denoting the number of infected nodes in $h$.

Similarly $C_{ij}^{k}$ is the transition rate from $S_{j}^{k}$ to
$S_{j-1}^{k-1}$ and

$C_{ij}^{k}=\left\\{\begin{array}[]{l}\gamma\quad\text{ if an
}S_{j}^{k-1}\rightarrow S_{i}^{k}\\\ \quad\text{ transition exints and
}l\text{ is the }\\\ \quad\text{ previously infected node }\\\ 0\quad\text{
otherwise }\end{array}\right.$

$B^{k}$ describes transitions between $S_{j}^{k}$ and $S_{i}^{k}$, which is
only possible in our model if $i=j$. Since the sum of elements in each column
of $P$ is zero (3), then it follows that the diagonal elements of $B^{k}$ are
given by

$B_{u}^{k}=-\sum_{\lambda IH}A_{mu}^{k+1}-\sum_{AH}C_{mu}^{*-1}$

We denote the first term in this equation by $N_{St}$, and note that it
represents the sum of
$\sum_{\mathrm{A}=\mathrm{n}}f\left(N_{\mathrm{A}}\left(S_{j}^{k-1}\right)\right.$
over all susceptible nodes in $S_{j}^{k-1}$. The second term can be evaluated
as $\gamma(k+1)$ since the number of non zero terms equals to the number of
infected nodes in $S_{\mathrm{m}}^{\text{k }+1}$.

$B_{ij}^{k}=\begin{cases}N_{Su}-\gamma(k+1)&\text{ if }i=j\\\ 0&\text{
otherwise }\end{cases}$

While these provide the exact evolution of the system, it is very difficult to
solve exactly and numerically for large hypergraphs. We therefore focus on
developing approximations for the expected values of the number of infected
and susceptible nodes.

Theorem 5 Let $[X]=\mathbb{E}(X)$ where $X=l$ ar $S$ and
$[SI]=\sum_{k=0}^{N}\sum_{j=1}^{c}N_{St}X_{j}^{k}(t)$ with
$\alpha=\left(\begin{array}[]{l}N\\\ k\end{array}\right)$. Then the following
equations hold for arbitrary hypergraphs:

$\displaystyle\frac{d[S]}{dt}=\gamma[I]-\tau[SI]$
$\displaystyle\frac{d[I]}{dt}=-\gamma[I]+\tau[SI]$

Proof. Introduce $c_{k}=(1,1,\ldots,1)$ as vector with $\alpha$ ones. We can
then write $[l]=\sum_{k=0}^{N}ke_{k}X^{k}$ and
$|S|=\sum_{k=0}^{N}(N-k)e_{k}X^{k}$. We can also write 4 as
$B_{4}^{k}=-\left(e_{k+1}A^{k+1}\right)_{i}-\left(e_{k-1}C^{k-1}\right)_{i}=$
$\left(c_{i}B^{k}\right)_{i}$. So it follows that

$a_{k+1}A^{k+1}+c_{k-1}C^{k-1}+c_{k}B^{k}=0$

Next we compute

$\displaystyle\frac{d[S]}{dt}$
$\displaystyle=\sum_{k=0}^{N}(N-k)e_{k}\frac{dX^{k}}{dt}$
$\displaystyle=\sum_{k=0}^{N}(N-k)\left(e_{k}A^{k}X^{k-1}+e_{k}B^{k}X^{k}+e_{k}C^{k}X^{k+1}\right)$
$\displaystyle=\sum_{k=0}^{N}\left((N-k-1)e_{k+1}A^{k+1}+(N-k)e_{k}B^{k}\right.$
$\displaystyle\left.+(N-k+1)C^{k-1}\right)X^{k}$
$\displaystyle=-\sum_{k=0}^{N}\left(e_{k+1}A^{k+1}-e_{k-1}C^{k-1}\right)X^{k}$

In the first line of the above simplification we used 1, in the second line we
re-labled the dummy variables and in the third line we used 8 twice to
simplify. To continue we note as in 5 that
$\left(e_{k-1}C^{k-1}\right)_{4}=\gamma_{k}=\left(\gamma k_{u}\right)_{2}$.
Also note that
$\left(e_{k+1}A^{k+1}\right)_{t}=\sum_{m=0}^{k+1}A_{mi}^{k+1}=FN_{Si}$ so that
$\sum_{k=0}^{N}e_{k+1}A^{k+1}X^{k}=$ $r\sum_{k=0}^{N}N_{sZ}X^{*}=r|ST|$.
Combining them we obtain

$\frac{d[S]}{dt}=\gamma[l]-r[SI]$

The proof of the second equation is similar and is done in (1). Alternatively
one could note that $\mathbb{E}(I+S)=N$ and thus that
$\frac{d[I]}{dt}=-\frac{d[S]}{dt}$.

It is important to note at this stage that these equations do not provide a
closed set of differential equations as $[SI]$ is not known. Various
approximations can be made to estimate [SI] depending on the structure of the
hypergraph and the form of $f$.

Exact description of SISR model on hypergraphs.

We introduce a new possible state for each node denoted as D. Our state space
now has dimension $3^{N}$. We keep the same notation as in the $SIS$ model to
denote our state space and introduce a new transition ’deactivation’:

Definition 6. The state space in the SISR model on a hypergraph with $|V|=N$
is $\\{S,I,D\\}^{N}$. We keep the notation of $0.4$ and introduce a class
$R^{k}$ which correspond to all $2^{N-k}\left(\begin{array}[]{c}N\\\
k\end{array}\right)=b_{k}$ states containing $k$ susceptible nodes.
$R_{i}^{k}\in R^{k}$ are indexed elements of $R^{k}$ with $1\leq i\leq
c_{k}.R_{i}^{k}(j)\in\\{S,I,D\\}$ denotes the state of the $j^{\text{th }}$
node.

The transitions in the model are identical to the SIS model but we add another
a deactivation transition: An $S_{j}^{k}\rightarrow S_{i}^{k-1}$ transition,
such that there exists a $p\in\\{1,\ldots,N\\}$ satisfying $S_{j}^{k}(p)=I$
and $S_{i}^{k-1}(p)=D$ and $S_{j}^{k}(m)=S_{m}^{k-1}$ for all $m\neq
p,m\in\\{1,\ldots,N\\}$.

Definition 7. A deactivation is a poisson process between $S$ and $D$ with
parameter $\beta$ such that the probability of deactivation in a small time
$\delta t$ is $1-\exp(-\beta\delta t)$

Next we introduce $Y^{k}=\left(Y_{1}^{k},\ldots,Y_{c_{k}}^{k}\right)$ to
denote the vector of probabilities with $Y_{i}^{k}$ being the probability that
the system is in state $R_{i}^{k}$. We now have a pair of differential
equations:

$\displaystyle\frac{dX^{k}}{dt}=A_{1}^{k}X^{k-1}+B_{1}^{k}X^{k}+C_{1}^{k}X^{k+1}$
$\displaystyle\frac{dY^{k}}{dt}=M^{k}Y^{k-1}+N^{k}Y^{k}+L^{k}Y^{k+1}$

$A_{1}$ again describes the rate of infection, so $A_{1}^{k}=A^{k}$ with some
abuse of notation. $C_{1}$ describes recovery and deactivation processes, both
of which correspond to $S^{k+1}\rightarrow S^{k}$ processes. We have

$C_{ij}^{k}=\begin{cases}\gamma&\text{ if an }S_{j}^{k-1}\rightarrow
S_{i}^{k}\text{ infection }\\\ &\text{ transition exists. }\\\ \beta&\text{ if
an }S_{j}^{k-1}\rightarrow S_{i}^{k}\text{ deactivation }\\\ &\text{
transition exists. }\\\ 0&\text{ otherwise }\end{cases}$

The matrix $B_{1}^{k}$ is determined identically as in the case of the SIS
model. $M_{ij}^{k}$ is the recovery rate from $R_{j}^{k-1}$ to $R_{i}^{k}$, so
we have

$M_{ij}^{k}=\begin{cases}\gamma&\text{ If a recovery transition
}R_{j}^{k-1}toR_{i}^{k}\\\ &\text{ exists }\\\ 0&\text{ otherwise
}\end{cases}$

$L_{ij}^{k}$ describes either infection of a single node. Hence we have;

$L_{ij}^{k}=\begin{cases}\tau\sum_{h:l\in
h}f\left(N_{h}\left(R_{j}^{k+1}\right)\right)&\text{ If an infection
transition }\\\ &\text{ is possible and where }l\text{ is the }\\\ &\text{
newly infected node }\\\ 0&\text{ otherwise }\end{cases}$

The matrix $N^{k}$ is determined as in the case of the SIS model. This system
has higher dimension than the SIS model and we therefore focus on expected
valued. We denote $[SI]_{t}=\mathbb{E}(SI)$

Theorem 6. With all the definitions of $0.4$ and $[D]=$ $\mathbb{E}(D)$ the
following hold for an arbitrary hypergraph:

$\frac{d[I]}{dt}=\tau[SI]-(\gamma+\beta)[I]$
$\begin{gathered}\frac{d[S]}{dt}=\gamma[I]-\tau[SI]\\\
\frac{d[D]}{dt}=\beta[I]\end{gathered}$

Note that the evolution equation for $[S]$ is as expected, the same as in the
SIS model.

Proof. The proof is similar to that of the proof of $Theorem\,5$ or as done in
(1). We will only show how the additional $\beta$ factor is obtained. This
difference arises in the evaluation of
$\left(e_{k-1}C^{k-1}\right)_{j}=\sum_{i}C_{ij}$. for fixed $j,C_{ij}^{k}$ is
none zero $k$ times as $S_{j}^{k}$ has $k$ infected nodes. Each one of these
can undergo the same transitions and each correspond to a unique term in
$S^{k-1}$. So
$\left(e_{k-1}C^{k-1}\right)_{j}=\left(k(\gamma+\beta)e_{k}\right)_{j}$.

The proof 16 is similar, but we reproduce it here due to the different
notation. We introduce $e_{k}=(1,\ldots,1)$ with $b_{k}$ ones and similarly to
the SIS model we have: $N_{ii}^{k}=\left(e_{k}N^{k}\right)_{i}=$
$-\left(e_{k+1}M^{k+1}\right)_{i}-\left(e_{k-1}L^{k-1}\right)_{i}$. Computing
the derivative:

$\displaystyle\frac{d[S]}{dt}$
$\displaystyle=\sum_{k=0}^{N}ke_{k}\frac{dY^{k}}{dt}$
$\displaystyle=\sum_{k=0}^{N}\left(e_{k+1}M^{k+1}-e_{k-1}L^{k-1}\right)Y^{k}$

The manipulations required are identical to the proof of $0.4.1$. Now we look
at each term separately. Let $N_{I}(G)$ denote the number of infected nodes in
state $G$.

$\displaystyle\sum_{k=0}^{N}e_{k+1}M^{k+1}Y^{k}$
$\displaystyle=\sum_{k=0}^{N}\sum_{allj}\sum_{alli}M_{ij}^{k+1}Y_{j}^{k}$
$\displaystyle=\sum_{k=0}^{N}\sum_{allj}\gamma
N_{I}\left(R_{j}^{k}\right)Y_{j}^{k}$ $\displaystyle=\gamma[I]$

With $N_{I}$ denoting the number of infected nodes in the state. The last
equality follows from the definition of expectation values in two state system
since the classes $Y^{k}$ form a partition of the state space.

$\displaystyle\sum_{k=0}^{N}e_{k-1}L^{k-1}Y^{k}$
$\displaystyle=\sum_{k=0}^{N}\sum_{allj}\sum_{alli}L_{ij}^{k-1}Y_{j}^{k}$
$\displaystyle=\sum_{k=0}^{N}\sum_{allj}\sum_{alli}\tau\sum_{h:l\in
h}f\left(N_{h}\left(R_{j}^{k}\right)\right)Y_{j}^{k}$
$\displaystyle=\sum_{k=0}^{N}\sum_{allj}\tau N_{SI}Y_{j}^{k}$
$\displaystyle=\tau[SI]$

The last two equality’s follow because there is a bijection between states
$R_{j}^{k}$ and $S_{j}^{k}$. This gives

$\frac{d[S]}{dt}=\gamma[I]-\tau[SI]$

Finally since $\mathbb{E}(S+I+D)=N$, then using linearity and differentiating,
$\frac{d[D]}{dt}=-\frac{d[I]}{dt}-\frac{d[S]}{dt}=\beta[I]$

Simulation Methods

So far our discussion has mostly general, however for the aim of simulating
epidemics, we choose to focus on Bi-Uniform Hypergraphs: Each vertex is
contained in exactly two edges of two distinct types which we call their
respective edge sets $H$ and $W$. Edges of type $H$ (respectively $W$ ) will
be disjoint and form a partition of $V$ and $(V,H\cup W)$ is a hypergraph. All
elements $h\in H$ have the same size $|h|=d_{h}$ and elements $w\in W$ also
have the same size $|w|=d_{w}.$ Random generation of such a hypergraph is done
as follows:

i) First generate a list $V$ of vertices. To select $k$ random vertices from
$V$, we randomly order $V$ and pick the first $k$ elements. i) Generating $H$
: Given $d_{h}$, we thus randomly pick $d_{h}$ elements from $V$, then $d_{h}$
elements from $V/h_{1}$ and iterate so that at the $i^{th}$ step we randomly
select $d_{h}$ elements from $V\backslash\left(\cup_{j<i}h_{j}\right)$.

ii) Generating $W$ : The same procedure as in generating $H$ is used.
Note that this requires $|V|$ to be divisible by $d_{h}$ and $d_{w}$ and we
arbitrarily choose $d_{w}\geq d_{h}$ $f$ is chosen to be non linear, as
otherwise it would be possible to get a very good approximation by using a
complete graph model of the hypergraph:

$f=d_{w}\arctan\left(\frac{\pi x}{2\left(d_{w}-1\right)}\right)$

We make various approximations throughout to estimate $[SI]$. If the number
total number of infected people is $I$, the number of infected nodes in $h\in
H$ is approximated as $\frac{d_{h}-1}{N}I$. Similarly the number of infected
nodes in $w\in W$ is approximated as $\frac{d_{w}-1}{N}I$. In the case of the
$SIR$ model this leads to:

$\displaystyle(t)$
$\displaystyle=\sum_{k=0}^{N}\sum_{j=1}^{b_{k}}\sum_{l:R_{j}^{k}(l)=S}\sum_{h:l\in
h}f\left(N_{h}\left(R_{j}^{k}\right)\right)Y_{j}^{k}$
$\displaystyle=\sum_{k=0}^{N}\sum_{j=1}^{b_{k}}\sum_{l:R_{j}^{k}(l)=S}f\left(\frac{d_{w}-1}{N}I\right)+f\left(\frac{d_{h}-1}{N}I\right)Y_{j}^{k}$
$\displaystyle=\sum_{k=0}^{N}\sum_{j=1}^{b_{k}}kf\left(\frac{d_{w}-1}{N}I\right)+f\left(\frac{d_{h}-1}{N}I\right)Y_{j}^{k}$
$\displaystyle=\mathbb{E}\left(S(t)\left(f\left(\frac{d_{w}-1}{N}I(t)\right)+f\left(\frac{d_{h}-1}{N}I(t)\right)\right)\right)$
$\displaystyle=[S]\left(f\left(\frac{d_{w}-1}{N}[I]\right)+f\left(\frac{d_{h}-1}{N}[I]\right)\right)$

In the last line we further approximated $S(t)$ to be independent of
$f\left(\frac{d_{w}-1}{N}I\right)+f\left(\frac{d_{h}-1}{N}I\right)$. These
approximations are only valid if they are sensible initially in the system
$(t=0)$ and if $N\gg d_{w}.$ Different approximation is similar for different
graphs and is done in (1).

In order to test the validity of our simulations, we perform stochastic
simulations of the models and compare them with solutions to our system of
differential equations. Label first the edges of $(V,H\cup W=E)$ as $e_{i}$
and nodes as $v_{i}$, and we define its incidence matrix $J$ as follows:

$J_{ij}=\begin{cases}1&\text{ if }v_{i}\in e_{j}\\\ 0&\text{ otherwise
}\end{cases}$

Let $x(t)=\left(x_{1}(t),\ldots,x_{N}(t)\right)$ be a vector such that
$x_{i}(t)=1$ if $v_{i}$ is infected at time $t$ and 0 otherwise. Define the
vector $y(t)$ similarly except that $y_{i}(t)=1$ if node $v_{i}$ is
susceptible at time $t$. Then $(xJ)_{j}\left((yJ)_{j}\right)$ gives the number
of infected (susceptible) nodes in hyperedge $e_{j}$. The probability of
infection of node $v_{i}$ in a time $\delta t$ is then
$1-\exp\left(-\tau\sum_{e_{i}\in E}J_{ij}f\left((xJ)_{j}\right)\delta
t\right.$. We then simulate transitions at each step by generating a random
$N$ dimensional vector $r$ such that each component $r_{i}\in(0,1)$. In a time
$\delta t$ an infection occurs at $v_{i}$ if
$r_{i}<1-\exp\left(-\tau\sum_{e_{i}\in E}J_{ij}f\left((xJ)_{j}\right)\delta
t\right)$ and $y_{i}=1$; a recovery occurs if $r_{i}<1-\exp(-\gamma\delta t)$,
and $x_{i}=1$; a deactivation occurs if $r_{i}<2-\exp(-\beta\delta
t)+\exp(-\gamma\delta t)$, $r_{i}>1-\exp(-\gamma\delta t)$ and $x_{i}=1$. The
number of deactivated nodes is then $N-\sum_{i}x_{i}(t)+y_{i}(t)$. This
generates a statistically accurate distribution []. Multiple simulations are
carried out with identical initial conditions ( with constant $\delta t$ and
length) and averaged to find $[I],[S]$ and $[D]$. We can make this simulation
arbitrarily accurate by taking $\delta t\rightarrow 0$.

Results

One important thing to note is that in (1) the model was not too accurate at
early times. This could be because the initial number of infected nodes was
not uniformly distributed. Indeed when we compare 1 and 2 , we see that the
predictions are significantly more accurate when the initial conditions are
uniformly distributed. This is because in 23 we made the approximation that
each hyperedge of a certain type contained the same number of infected nodes.
Comparing 3 with 1 and 2 we note that accuracy is significantly reduced when
$d_{h},d_{w}\ll N$. This is to be expected, as in this regime, it is unlikely
that each hyper edge of each type has the same number of infected nodes.
However when $d_{h},d_{w}$ are larger, this approximation is more likely to be
accurate [4]. Decreasing the value of $\tau$ shifts the location of the peak
to a lower time and flattens it as expected 7 . Taking $\gamma=0$ we recover
the $SIR$ model whose profile in 6 looks very similar to standard SIR models.
Similarly, taking $\beta=0$ we recover the SIS model 5. However for very small
values of $\beta$ do not seem to approach the SIS model. If we look at the
steady state solutions of $15,16,17$, we can see that $[I]=0,[S]=0,[D]=N$
always no matter what positive values of $\beta,\gamma$ are chosen. This
indicates a discontinuity in the phase space of this system.

[h!]

[h!]

Discussion

In this report, we generalised the SIS model on hypergraphs to a broader class
of spreading processes. The exact equations for the time evolution of expected
values were derived in a similar way to (1). A method to simulate the
corresponding stochastic process as described and numerical simulations were
compared to these. We found strong agreement was found as long as the
approximations and initial conditions were compatible and the size of each
hyperedge wasn’t small when compared to the number of nodes.

Our model predicts correct profiles but are only accurate in very specific
regimes due to our approximation being valid only in those regimes. All the
dependencies of our model were only hypothesized based on numerical
simulations. A proper phase space analysis is required to further investigate
the possible regimes of SISR model. Also, we have only looked at
approximations for one type of hypergraphs being used and these can perhaps be
generalised for hypergraphs with an arbitrary distribution of nodes (i.e. by
characterising the nodes hypergraph in terms of the probability that they
belong to some hyperedge). In order to apply this model to realistic
situations, the function $f$ would need to be fitted with known data. Further
depending on the needs, the parameters $\gamma,\beta$ should potentially
depend on other factors such that how many times a given node has been
infected. The advantage of the SISR model is that it takes into account that
not all elements transition into a 3rd state instantly. The model also
accounts for community structure as well as non-linear dependence of infection
pressure on the number of neighbours in each edge.

References.

1\. GYK Ágnes Bodó, Besson, PL Simon, Sis epidemic propagation on hypergraphs.
Bull. Math. Biol. 78, 713-735 (2016).

2\. PV Mieghem, Exact markovian sir and sis epidemics on networks and an upper
bound for the epidemic threshold (2014).

3\. GD Prato, Kolmogorov Equations for Stochastic PDEs. (Birkhãuser, Basel),
(2004).