# Horizontal and Vertical Collaboration for VR Delivery in MEC-Enabled Small-
Cell Networks

Zhuojia Gu  Hancheng Lu  and Chenkai Zou

###### Abstract

Due to the large bandwidth, low latency and computationally intensive features
of virtual reality (VR) video applications, the current resource-constrained
wireless and edge networks cannot meet the requirements of on-demand VR
delivery. In this letter, we propose a joint horizontal and vertical
collaboration architecture in mobile edge computing (MEC)-enabled small-cell
networks for VR delivery. In the proposed architecture, multiple MEC servers
can jointly provide VR head-mounted devices (HMDs) with edge caching and
viewpoint computation services, while the computation tasks can also be
performed at HMDs or on the cloud. Power allocation at base stations (BSs) is
considered in coordination with horizontal collaboration (HC) and vertical
collaboration (VC) of MEC servers to obtain lower end-to-end latency of VR
delivery. A joint caching, power allocation and task offloading problem is
then formulated, and a discrete branch-reduce-and-bound (DBRB) algorithm
inspired by monotone optimization is proposed to effectively solve the
problem. Simulation results demonstrate the advantage of the proposed
architecture and algorithm in terms of existing ones.

###### Index Terms:

VR delivery, horizontal and vertical collaboration, MEC, end-to-end latency.

##  I introduction

On-demand virtual reality (VR) applications require ultra-large bandwidth,
ultra-low latency and abundant computation resources in delivery networks,
which hinders its widespread popularity in existing 5G networks [1]. By
providing caching and computation capabilities at the network edge, mobile
edge computing (MEC) can significantly reduce the transmission and computation
delays, which are the main components of the end-to-end latency of VR delivery
[2, 3].

Existing works on MEC can be categorized into two types: 1) vertical
collaboration (VC)-based MEC, and 2) horizontal collaboration (HC)-based MEC.
The VC-based MEC typically adopts a three-tier hierarchical architecture
including the terminal tier, the edge tier, and the cloud tier to conduct
collaborative caching and computation offloading strategies [4]. In the
scenario of VR delivery, the authors in [5, 6] proposed the joint
communications, caching and computation (3C) modeling of a hierarchical VR
delivery network architecture, and revealed the 3C resources trade-off for
saving communication bandwidth while meeting the low latency requirement. On
the other hand, considering that mobile edge servers (MESs) are usually
deployed sinking to the densely deployed base stations (BSs) in 5G wireless
networks [7, 8], HC-based MEC has been proposed for higher utilization
efficiency of edge resources, facilitated by the higher speed wireless
connections in dense small-cell networks than that of the backhaul link [9].
The authors in [10, 11] investigated the HC strategies of caching or
computation resources among MESs, respectively, and the performance gains
using HC in terms of the service delay were validated.

From the discussion above, it is evident that both VC and HC are effective
methods to improve the performance of delay-sensitive applications in wireless
networks. Intuitively, combining HC with VC in MEC-enabled small-cell networks
is a potential solution to further meet the stringent latency requirement of
computation-intensive VR applications. However, there is a lack of studies
focusing on the joint HC and VC-assisted MEC for VR video delivery.
Considering the specific characteristic of VR video delivery, the design of
joint HC and VC in MEC-enabled small-cell networks is non-trivial. First,
since the VR field of views (FoVs) needs to be projected and rendered from
monocular videos (MVs) into stereoscopic videos (SVs) before it can be
presented at VR head-mounted displays (HMDs), caching MVs in MESs can save the
cache capacity, but at the cost of consuming more computation resources.
Considering limited caching and computation resources in MESs, it is a
critical issue to design HC-based caching and task offloading strategies among
MESs to reduce the transmission and computation delays while considering the
VC task offloading aspects. Second, since HC and VC are to be considered in
collaboration to further improve the end-to-end latency performance for VR
delivery, it will incur a trade-off issue between HC and VC for collaborative
caching and task offloading. In addition, power allocation in dense small-cell
networks has a great impact on the transmission rates of wireless links. In
the case of limited power resources at BSs, power allocation should be
reconsidered in coordination with joint HC and VC-based caching and task
offloading strategies to achieve lower total latency for VR delivery.

In this letter, we propose a joint HC and VC architecture in MEC-enabled
small-cell networks for VR delivery. A novel task offloading problem with
consideration of collaborative caching and power allocation is formulated with
the aim of minimizing the end-to-end latency for VR delivery. A joint caching,
power allocation and task offloading algorithm envisioned by discrete monotone
optimization is proposed to effectively solve the problem. Numerical
simulations have been conducted to justify the effectiveness of the joint HC
and VC architecture for VR delivery and show that the proposed JCPT algorithm
outperforms the benchmark algorithms.

##  II System Model

###  II-A Network Architecture

Figure 1: A joint horizontal and vertical collaboration architecture in MEC-
enabled small-cell networks for VR delivery.

The proposed network architecture is shown in Fig. 1. There are $M$ small-cell
base stations (SBSs) equipped with MEC servers (MESs) in the network denoted
by an index set $\mathcal{M}=\\{1,...,M\\}$, with each MES $m\in\mathcal{M}$.
The VR head-mounted devices (HMDs) are denoted by an index set
$\mathcal{U}=\\{1,...,U\\}$, with each HMD $u\in\mathcal{U}$. The set of all
viewpoints of VR videos is denoted by $\mathcal{N}=\\{1,...,N\\}$, and each
viewpoint $i\in\mathcal{N}$ corresponds to a MV and a SV. MVs should be
projected and rendered into SVs by using the computation resources at MESs or
HMDs, thus the size of the SV version of viewpoint $i$ should be
$\alpha(\alpha>2)$ times larger than that of the corresponding MV version.

Due to the densification of small-cell networks, a HMD can be within the
coverage of multiple SBSs, thus the viewpoints can be retrieved from multiple
SBSs caching the required MVs or SVs, and the computation task of projecting
MVs into SVs can also be processed on multiple MESs. The multiple connections
of a HMD to SBSs can be realized by adopting the multi-connectivity technology
proposed in 3GPP specification [12, 13]. Let $T_{imu}^{M}\in\\{0,1\\}$ denote
whether the computation task of the $i$-th MV on HMD $u$ is offloaded to MES
$m$. Similarly, let $T_{imu}^{C}$ denote whether the computation task of the
$i$-th MV on HMD $u$ is offloaded to the cloud service via MES $m$. Assume
that the offloading task for each viewpoint is performed at one server, thus
we have

$T_{imu}^{M}+T_{imu}^{C}\leq 1,\forall i\in\mathcal{N},\forall
m\in\mathcal{M},\forall u\in\mathcal{U},$ (1) $\sum_{m=1}^{M}T_{imu}^{M}\leq
1,\forall i\in\mathcal{N},\forall u\in\mathcal{U},$ (2)
$\sum_{m=1}^{M}T_{imu}^{C}\leq 1,\forall i\in\mathcal{N},\forall
u\in\mathcal{U}.$ (3)

###  II-B 3C Resource Model

SBSs and HMDs are equipped with caching capacity of $C_{u}^{M}$ and
$C_{u}^{V}$, respectively. Let $c_{im}^{M,M},c_{im}^{M,S}\in\\{0,1\\}$ denote
whether the MV or SV version of viewpoint $i$ is cached on MES $m$,
respectively. Similarly, let $c_{iu}^{V,M},c_{iu}^{V,S}\in\\{0,1\\}$ denote
whether the MV or SV version of viewpoint $i$ is cached on VR HMD $u$. The
size of the MV version of viewpoint $i$ is denoted as $d_{i}$, then the size
of the corresponding SV is $\alpha d_{i}$. Considering the caching capacity
constraints at MESs and HMDs, we have

$\sum_{i=1}^{N}(c_{im}^{M,M}+\alpha c_{im}^{M,S})d_{i}\leq C_{m}^{M},\forall
m\in\mathcal{M},$ (4) $\sum_{i=1}^{N}(c_{iu}^{V,M}+\alpha
c_{iu}^{V,S})d_{i}\leq C_{u}^{V},\forall u\in\mathcal{U}.$ (5)

MESs and VR HMDs are all equipped with computation resources. According to
[14], if the computation process is executed on an MES, the consumed
computation energy to project MV $i$ into SV $i$ can be calculated as
$e_{i}^{M}=k_{M}f_{M}^{2}d_{i}w_{i}$, where $k_{M}$ is the energy efficiency
coefficient related to the computing process unit (CPU) architecture of the
MES, $f_{m}$ is the CPU frequency of the MES, and $w_{i}$ is the required CPU
cycles of projection MV $i$ into SV $i$ for 1 bit. Then the computation delay
on the MES to process MV $i$ can be calculated as
$\tau_{i}^{M}=\frac{d_{i}w_{i}}{f_{M}}$. On the other hand, if the computation
process is executed on a VR HMD, the consumed computation energy is
$e_{i}^{V}=k_{V}f_{V}^{2}d_{i}w_{i}$, and the corresponding computation delay
is $\tau_{i}^{V}=\frac{d_{i}w_{i}}{f_{V}}$. Assume that the requested
probabilities of all the viewpoints conforms to the Zipf distribution, i.e.,
the requested probability of viewpoint $i$ is
$q_{i}=\frac{1/\mathrm{i}^{\lambda}}{\sum_{\mathrm{j}=1}^{N}1/\mathrm{j}^{\lambda}}$,
where $\lambda$ reflects the preference skewness of the Zipf distribution. To
ensure that the average energy consumption of each MES or HMD does not exceed
the threshold $E_{m}^{M}$ or $E_{u}^{V}$, respectively, we have

$\sum_{u=1}^{U}\sum_{i=1}^{N}q_{i}(1-c_{im}^{M,S})e_{i}^{M}\mathbb{I}_{\big{\\{}\sum_{m=1}^{M}T_{imu}^{M}=1\big{\\}}}\leq
E_{m}^{M},\forall m\in\mathcal{M},$ (6)
$\sum_{i=1}^{N}q_{i}(1-c_{iu}^{V,S})e_{i}^{V}\mathbb{I}_{\big{\\{}\sum_{m=1}^{M}T_{imu}^{M}+T_{imu}^{C}=0\big{\\}}}\leq
E_{u}^{V},\forall u\in\mathcal{U},$ (7)

where $\mathbb{I}_{\\{\cdot\\}}$ is an indicator function, i.e.,
$\mathbb{I}_{\\{z\\}}=1$ if condition $z$ is true, and zero otherwise.

Each SBS allocates transmission power to the associated HMDs, denoted by
$p_{mu}$. Assuming universal frequency reuse among SBSs, according to
Shannon’s formula, the downlink wireless link rate $R_{mu}$ from SBSs $m$ to
HMD $u$ can be expressed as

$R_{mu}=W\log_{2}(1+\gamma_{mu}(\mathbf{p})),$ (8)

where $W$ is the bandwidth allocated to HMD $u$,
$\gamma_{mu}(\mathbf{p})\triangleq\frac{p_{mu}h_{mu}}{I_{Inter-cell}+\ \zeta
I_{Intra-cell}+n_{0\ }}$, $h_{mu}$ is the channel gain between SBS $m$ and HMD
$u$, $n_{0}$ is the additive white Gaussian noise (AWGN) power. The inter-cell
and intra-cell interferences can be calculated as $I_{Inter-cell}=\sum_{b\neq
m}^{b\in M}\sum_{v\neq u}^{v\in U}p_{bv}h_{bv}$ and $I_{Intra-
cell}=\sum_{v\neq u}^{v\in U}p_{mv}h_{mv}$, respectively, where
$\zeta\in[0,1]$ is the orthogonality factor that represents the capability of
intra-cell interference cancelation at the HMD. Considering the total power
constraint $P_{T}$ at SBS $m$, we have

$\sum_{u=1}^{U}p_{mu}\leq P_{T},\forall m\in\mathcal{M}.$ (9)

###  II-C HC and VC Task Offloading Model

A HMD can be associated to multiple SBSs simultaneously to fetch the cached
viewpoints and utilize the computation resources on multiple MESs, which is
referred to as horizontal collaboration (HC) among MESs. Assume that the end-
to-end latency of VR delivery mainly includes the transmission delay and the
computation delay. Given the task offloading variable $T_{imu}^{M}$ and the
caching status of viewpoint $i$ at MES $m$, the end-to-end latency for HMD $u$
requesting viewpoint $i$ from MES $m$ can be calculated as

$\tau_{iu}^{M}=\sum_{m=1}^{M}\bigg{(}\frac{c_{im}^{M,M}d_{i}w_{i}}{f_{M}}+\frac{\alpha
d_{i}}{R_{mu}}\bigg{)}T_{imu}^{M}.$ (10)

Note that the task can be offloaded to MES $m$ can only when the MV or SV is
cached at MES $m$. Thus we have the coupled caching and task offloading
constraint,

$T_{imu}^{M}(1-c_{im}^{M,M})(1-c_{im}^{M,S})<1,\forall i\in\mathcal{N},\forall
m\in\mathcal{M},\forall u\in\mathcal{U}.$ (11)

On the other hand, computation task can also be implemented locally at the
HMD, or be offloaded to cloud server, which is referred to as vertical
collaboration (VC). Accordingly, if the computation task is implemented at the
HMD, the end-to-end latency just includes the computation delay, which can be
calculated as

$\tau_{iu}^{V}=\mathbb{I}_{\big{\\{}\sum_{m=1}^{M}T_{imu}^{M}+T_{imu}^{C}=0\big{\\}}}\frac{c_{iu}^{V,M}d_{i}w_{i}}{f_{V}}.$
(12)

Note that $\sum_{m=1}^{M}T_{imu}^{M}+T_{imu}^{C}=0$ can hold only when
$c_{iu}^{V,M}=1$ or $c_{iu}^{V,S}=1$, i.e., we have

$T_{imu}^{M}+T_{imu}^{C}+c_{iu}^{V,M}+c_{iu}^{V,S}>0,\forall
i\in\mathcal{N},\forall m\in\mathcal{M},\forall u\in\mathcal{U}.$ (13)

Otherwise, the requested viewpoint should be fetched via the cloud server, and
it would incur an additional backhaul retrieving delay $\tau_{b}$, thus the
end-to-end latency can be calculated as

$\tau_{iu}^{C}=\sum_{m=1}^{M}\bigg{(}\frac{\alpha\
d_{i}}{R_{mu}}+\tau_{b}\bigg{)}T_{imu}^{C}.$ (14)

Taking consideration of collaborated HC and VC of the MESs and the HMDs in the
network, the end-to-end latency for HMD $u$ requesting viewpoint $i$ can be
written as

$\tau_{iu}=\tau_{iu}^{M}+\tau_{iu}^{V}+\tau_{iu}^{C}.$ (15)

##  III Problem Formulation and Solution

In this letter, we aim to minimize the average end-to-end latency of VR
delivery based on the collaborated HC and VC architecture by jointly
considering the caching, power allocation and task offloading strategies.
Thus, the proposed optimization problem can be formulated as follows,

$\displaystyle\min_{\begin{subarray}{c}\\{c_{iu}^{V,M},c_{iu}^{V,S},c_{im}^{M,M},c_{im}^{M,S},\\\
p_{mu},T_{imu}^{M},T_{imu}^{C}\\}\end{subarray}}\quad$
$\displaystyle\sum_{u=1}^{U}\sum_{i=1}^{N}q_{i}\tau_{iu}$ (16a)
$\displaystyle\text{s.t.}\qquad c_{iu}^{V,M},$ $\displaystyle
c_{iu}^{V,S},c_{im}^{M,M},c_{im}^{M,S}\in\\{0,1\\},$ (16b) $\displaystyle
T_{imu}^{M},$ $\displaystyle T_{imu}^{C}\in\\{0,1\\},$ (16c)
$\displaystyle\qquad\qquad(\ref{eq1})-$
$\displaystyle(\ref{eq9}),(\ref{coupling-constraint}),(\ref{local-coupling-
constraint}).$

Problem (16) is a mixed integer nonlinear programming problem (MINLP), which
is hard to solve in general. In the following, we exploit the hidden
monotonicity of the problem and utilize the monotonic optimization method to
solve the problem effectively.

###  III-A Monotonicity Analysis and Problem Transformation

It can be observed that the objective function (16a) monotonically increases
with the increase of the caching variables $c_{iu}^{V,M},c_{im}^{M,M}$ and the
task offloading variables $T_{imu}^{M},T_{imu}^{C}$. However, (16a) is not a
monotonic function of the power allocation variables $p_{mu}$. Nevertheless,
inspired by the concept of exploiting the hidden monotonicity in the objective
function using the general linear fractional programming (GLFP) [15], (16a) is
actually a monotonically decreasing function of the positive-valued function
$\gamma_{mu}(\mathbf{p})$. More specifically, for fixed caching and task
offloading variables, Problem (16) is equivalent to

$\displaystyle\min_{\mathbf{y}}\quad$
$\displaystyle\sum_{u=1}^{U}\sum_{i=1}^{N}q_{i}\tau_{iu}(\mathbf{y})$ (17a)
s.t. $\displaystyle\quad\mathbf{y}\in\mathcal{G},$ (17b)

where $\mathcal{G}=\\{\mathbf{y}|0\leq
y_{mu}\leq\gamma_{mu}(\mathbf{p}),\forall m\in\mathcal{M},\forall
u\in\mathcal{U},(\ref{eq9})\\}$. Then, we should transform Problem (16) into
the form of a monotonic optimization problem.

Definition (Monotonic Optimization) [16]: An optimization problem can be
classified as a monotonic optimization (MO) problem if it can be written in
the following canonical form,

$\max_{\mathbf{x}}\\{f(\mathbf{x})\mid\mathbf{x}\in\mathcal{G}\cap\mathcal{H}\\},\mathbf{x}\in[\mathbf{a},\mathbf{b}],$
(18)

where $\mathbf{x}\in\mathbb{R}_{+}^{n}$,
$[\mathbf{a},\mathbf{b}]\subset\mathbb{R}_{+}^{n}$,
$f:\mathbb{R}_{+}^{n}\rightarrow\mathbb{R}$ is a monotonic increasing function
of $\mathbf{x}$, $\mathcal{G}\subset\mathbb{R}_{+}^{n}$ is a normal set with
nonempty interior, and $\mathcal{H}$ is a conormal set in $[0,\mathbf{b}]$.

Based on the above monotonicity analysis of Problem (16), it is observed that
the monotonicity of the objective function (16a) with respect to the caching
and task offloading variables is the opposite of that with respect to the
power allocation variables. According to the definition of monotonic
optimization, the objective function to be minimized should be a monotonically
decreasing function with respect to the optimization variables. Thus, variable
substitutions are made for the caching and task offloading variables to meet
the monotonic requirements. Specifically, let $\bar{c}_{iu}^{V,M}\triangleq
1-c_{iu}^{V,M}$, and $\bar{c}_{im}^{M,M}\triangleq 1-c_{im}^{M,M}$ for the
caching variables, and similar substitutions can be done to the task
offloading variables $T_{imu}^{M},T_{imu}^{C}$. Besides, it can be inferred
that the constraints of Problem (16) forms the intersection of a normal set
and a conormal set. Then Problem (16) can be equivalently transformed into a
monotonic optimization problem.

###  III-B Proposed Algorithm Based on DBRB

For the monotonic optimization problem, the most widely adopted algorithm is
the polyblock outer approximation (POA) approach [16], which attempts to
construct a sequence of polyblocks on the boundary of the constraints to
approximate the optimal solution of continuous variables. However, Problem
(16) cannot be solved using the POA approach since it involves the discrete
caching and task offloading optimization variables. A discrete branch-reduce-
and-bound (DBRB) algorithm was proposed in [17] to deal with the discrete
variables in a monotonic optimization problem. The DBRB algorithm draws on the
idea of the standard branch-and-bound algorithm that using box111As defined in
[16], a box $[\mathbf{a},\mathbf{b}]$ is the set of all
$\mathbf{x}\in\mathbb{R}^{n}$ that satisfies
$\mathbf{a}\leq\mathbf{x}\leq\mathbf{b}$. A box is also referred to as a
hyper-rectangle. division based on three basic operations at each iteration:
branching, reduction, and bounding, to converge to a near-optimal solution
rapidly. In this letter, we propose an improved DBRB algorithm that adapts to
the specific properties of Problem (16).

The proposed joint caching, power control and task offloading (JCPT) algorithm
based on the DBRB framework is outlined in Algorithm 1, where $B_{n}$,
$f_{max}^{(n)}$ and $f_{min}^{(n)}$ represents the set of boxes, the upper
bound and the lower bound, respectively. $r(\cdot)$ represents the reduction
operation for a box. in the $nth$ iteration. On the basis of three basis
operations at each iteration, we have made two improvements considering the
properties of the problem to accelerate the solution searching process. First,
the branching and reduction operations in steps 6-8 can be improved
considering that the caching and task offloading variables are all binary
variables. Specifically, given a Box
$[\mathbf{a},\mathbf{b}]\subset\\{0,1\\}^{n}$, if $b_{k}-a_{k}=1$, we can
branch on its $k$th dimension and obtain two new Box
$[\mathbf{a},\mathbf{a^{\prime}}]$ and $[\mathbf{b^{\prime}},\mathbf{b}]$,
where the value of each dimension of $\mathbf{a^{\prime}}$ is the same as that
of $\mathbf{b}^{\prime}$ except that the $k$th dimension of $a_{k}^{\prime}$
is 0, and the $kth$ dimension $b_{k}^{\prime}$ is 1. Thus, the dimension of
each branched box is reduced by one, which reduces the searching space. Also,
the reduction operation can also be improved using the same technique. Second,
the bounding operation in step 11 can be improved considering the binary
feature of the task offloading variables. In the standard DBRB framework, the
bounding operation is to find a projection point of the upper right vertex of
a box on the constraint curve, then the projection point serves as the lower
bound of the solution, and the L-shaped area around the projection point
serves as the upper bound of the solution. Nevertheless, the binary task
offloading variables make the projection point not be a feasible solution.
Alternatively, a heuristic algorithm can be adopted to find a tighter upper
and lower bound of the solution, e.g., by randomly assigning the MEC task
offloading solution and obtain the joint caching and power allocation solution
using a greedy algorithm. Thus, a feasible solution utilized to compute the
upper or lower bound can be obtained.

Algorithm 1 The JCPT algorithm

1: Input: $\mathbf{a}$, $\mathbf{b}$;

2: Initialization: $\mathcal{B}_{1}=\\{r[\mathbf{a},\mathbf{b}]\\}$,
$n\leftarrow 1$, $f_{min}^{(1)}\leftarrow 0$,
$f_{max}^{(1)}\leftarrow+\infty$, $\gamma\leftarrow+\infty$

3: while not converged do

4: $n\leftarrow n+1;$

5: Select Box $\mathbf{B}_{s}$ that has the maximum upper bound of objective
function in $\mathcal{B}_{n-1}$;

6: Branch on Box $\mathbf{B}_{s}$ to obtain Box $\mathbf{B}_{s}^{(1)}$ and
$\mathbf{B}_{s}^{(2)}$;

7: if $\mathbf{B}_{s}^{(k)}$ is feasible and contains a better solution then

8: Compute $r(\mathbf{B}_{s}^{(k)})$, $k\in\\{1,2\\}$;

9: end if

10: if $r(\mathbf{B}_{s}^{(k)})$ is feasbile then

11: Adopt a heuristic algorithm to compute the upper and lower bounds in this
iteration;

12: Update $f_{min}^{(n)}$ and $f_{max}^{(n)}$, $k\in\\{1,2\\}$;

13: end if

14: if $f_{max}^{(n)}<\gamma$ then

15: Update $\gamma\leftarrow f_{max}^{(n)}$ and record the solution $x^{*}$.

16: end if

17: Add the feasible box in $r(\mathbf{B}_{s}^{(1)})$ and
$r(\mathbf{B}_{s}^{(2)})$ into $\mathcal{B}_{n-1}$ and remove Box
$\mathbf{B}_{s}$;

18: $\mathbf{B}_{n}\leftarrow\mathbf{B}_{n-1}$;

19: end while

20: Output $(x^{*},\gamma)$;

##  IV Performance Evaluation

Simulations are carried out to validate the proposed joint HC and VC network
architecture for VR delivery. The simulations are conducted in a square area
of 100 m $\times$ 100 m, where 40 SBSs and 100 HMDs are distributed in this
area unless otherwise stated. There are a total of 100 VR viewpoints ranging
in size from 10 Mb to 30 Mb. The maximum transmit power of SBSs is 30 dBm, and
the downlink bandwidth is set to 1 GHz. For the sake of comparison, four
benchmark algorithms are listed as follows. 1) Nearest Offloading (NO): a task
is offloaded to the nearest MES if the local caching or computation resources
cannot meet the requirements. In other words, HC is not adopted in this
algorithm; 2) Power Equally Allocation (PEA): transmit power is equally
allocated among HMDs within coverage, then the caching and task offloading
strategies are jointly optimized; 3) Popularity-first (PF): MESs and HMDs
cache the most popular viewpoints, then the task offloading and power
allocation strategies are jointly optimized; 4) Least Recently Used (LRU): a
classical caching placement algorithm that substitutes the least recently
requested MVs or SVs.

Figure 2: Sum of delays versus (a) caching capacity on a MES, and (b) number
of SBSs.

Fig. 2 shows the sum of delays with respect to the MES caching capacity and
the number of SBSs using different algorithms. It is observed that the
proposed JCPT algorithm achieves the lowest sum of delays in various network
parameters. The NO algorithm achieves the second best delay performance, while
the PF algorithm achieves the worst performance. This is because the proposed
algorithm comprehensively considers the HC and VC in wireless and edge
networks, especially making full use of the overlapping coverage
characteristics of dense small cell networks, so that the utilization of
communication, caching and computing resources for VR delivery can be
maximized. In contrast, the NO algorithm is based on the classical nearby
association and task offloading principles, which does not make full use of
the HC among multiple MESs, and only considers VC to offload local tasks to
the MESs or the cloud server, thus resulting in performance degradation. The
PEA algorithm does not take into consideration the coupling relationship
between wireless transmit power and heterogeneous caching and computation
resources on the MESs, so the sum of delays is even higher. The PF algorithm
is a greedy heuristic algorithm which does not consider the characteristics of
VR videos, so it fails to obtain a better strategy in the decision of caching
MVs or SVs of different viewpoints, which results in the highest the delay
performance.

Figure 3: Cache Hit ratio versus (a) caching capacity on a MES, and (b) the
zipf parameter.

Fig. 3 shows the cache hit ratio with respect to the MES cache capacity and
the zipf viewpoint popularity parameter with different algorithms. It can be
seen that the proposed JCPT algorithm achieves the highest cache hit ratio.
Combining with the observation in Fig. 2, it can be inferred that a higher
cache hit ratio will lead to a lower sum of delays, which indicates the
importance of an edge caching hit event to reduce the end-to-end latency of VR
delivery over MEC-enabled small cell networks. Since the proposed JCPT
algorithm exploits the multi-connectivity advantage of a HMD associated to
multiple SBSs, the HMD can obtain the requested viewpoints directly from
different MESs. Therefore, the JCPT algorithm takes advantage of this
characteristic in making caching strategies, which can be viewed as a
collaborative caching strategy among MESs. Meanwhile, the computation
requirements of VR videos have been jointly considered with the computation
resources on MESs and wireless resources of SBSs using the JCPT algorithm,
thus providing a higher cache hit ratio and reducing the end-to-end latency of
VR delivery.

##  V Conclusion

In this letter, we design a collaborative HC and VC architecture in MEC-
enabled small cell networks for low-latency on-demand VR delivery. A joint
horizontal and vertical task offloading model is introduced with consideration
of caching and wireless communication resources to fully exploit the advantage
of small cell networks for VR delivery. The hidden monotonicity of the
formulated end-to-end latency minimization problem is analyzed and an improved
DBRB-based algorithm is designed to efficiently solve the problem. Simulation
results validate the effectiveness of the proposed JCPT algorithm in
comparison with the existing benchmark algorithms and show great promise of
the proposed joint HC and VC architecture for improving the quality of VR
delivery over MEC-enabled small cell networks.

## References

  * [1] M. Hosseini, “View-aware tile-based adaptations in 360 virtual reality video streaming,” in 2017 IEEE Virtual Reality (VR), Los Angeles, CA, 2017, pp. 423–424.
  * [2] X. Liu and Y. Deng, “Learning-based prediction, rendering and association optimization for MEC-enabled wireless virtual reality (VR) network,” IEEE Transactions on Wireless Communications, doi: 10.1109/TWC.2021.3073623.
  * [3] Y. Qian, R. Wang, J. Wu, B. Tan, and H. Ren, “Reinforcement learning-based optimal computing and caching in mobile edge network,” IEEE Journal on Selected Areas in Communications, vol. 38, no. 10, pp. 2343–2355, 2020.
  * [4] N. Abbas, Y. Zhang, A. Taherkordi, and T. Skeie, “Mobile edge computing: A survey,” IEEE Internet of Things Journal, vol. 5, no. 1, pp. 450–465, 2018.
  * [5] Y. Sun, Z. Chen, M. Tao, and H. Liu, “Communications, caching, and computing for mobile virtual reality: Modeling and tradeoff,” _IEEE Transactions on Communications_ , vol. 67, no. 11, pp. 7573–7586, 2019.
  * [6] T. Dang and M. Peng, “Joint radio communication, caching, and computing design for mobile virtual reality delivery in fog radio access networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 37, no. 7, pp. 1594–1607, 2019.
  * [7] H. Wu and H. Lu, “Delay and power tradeoff with consideration of caching capabilities in dense wireless networks,” IEEE Transactions on Wireless Communications, vol. 18, no. 10, pp. 5011–5025, 2019.
  * [8] Z. Gu, H. Lu, M. Zhang, H. Sun, and C. W. Chen, “Association and caching in relay-assisted mmWave networks: From a stochastic geometry perspective,” IEEE Transactions on Wireless Communications, doi: 10.1109/TWC.2021.3091815.
  * [9] K. Wang, H. Yin, W. Quan, and G. Min, “Enabling collaborative edge computing for software defined vehicular networks,” IEEE Network, vol. 32, no. 5, pp. 112–117, 2018.
  * [10] Y. M. Saputra, D. T. Hoang, D. N. Nguyen, and E. Dutkiewicz, “A novel mobile edge network architecture with joint caching-delivering and horizontal cooperation,” IEEE Transactions on Mobile Computing, vol. 20, no. 1, pp. 19–31, 2021.
  * [11] P. Dai, K. Hu, X. Wu, H. Xing, F. Teng, and Z. Yu, “A probabilistic approach for cooperative computation offloading in MEC-assisted vehicular networks,” IEEE Transactions on Intelligent Transportation Systems, pp. 1–13, 2020.
  * [12] 3GPP, “Evolved universal terrestrial radio access (E-UTRA) and NR; Multi-connectivity,” (Release 16), 3GPP TS 37.340, V16.6.0, 2021.
  * [13] A. Wolf, P. Schulz, M. Dorpinghaus, J. C. S. S. Filho, and G. Fettweis, “How reliable and capable is multi-connectivity?,” IEEE Transactions on Communications, vol. 67, no. 2, pp. 1506–1520, 2019.
  * [14] Y. Mao, C. You, J. Zhang, K. Huang, and K. B. Letaief, “A survey on mobile edge computing: The communication perspective,” IEEE Communications Surveys $\&$ Tutorials, vol. 19, no. 4, pp. 2322–2358, 2017.
  * [15] Y. J. A. Zhang, L. Qian, and J. Huang, “Monotonic optimization in communication and networking systems,” Foundations and Trends in Networking, vol. 7, no. 1, pp. 1–75, 2013.
  * [16] H. Tuy, “Monotonic optimization: Problems and solution approaches,” SIAM Journal on Optimization, vol. 11, no. 2, pp. 464–494, 2000.
  * [17] H. Tuy, M. Minoux, and N. T. Hoai-Phuong, “Discrete monotonic optimization with application to a discrete location problem,” SIAM Journal on Optimization, vol. 17, no. 1, pp. 78–97, 2007.


# Minimal theory of massive gravity in the light of CMB data and the $S_{8}$
tension

José C. N. de Araujo<EMAIL_ADDRESS>Divisão de Astrofísica,
Instituto Nacional de Pesquisas Espaciais, Avenida dos Astronautas 1758, São
José dos Campos, 12227-010, SP, Brazil  Antonio De Felice
<EMAIL_ADDRESS>Center for Gravitational Physics, Yukawa
Institute for Theoretical Physics, Kyoto University, 606-8502, Kyoto, Japan
Suresh Kumar<EMAIL_ADDRESS>Department of Mathematics, Indira Gandhi
University, Meerpur, Haryana-122502, India  Department of Mathematics, BITS
Pilani, Pilani Campus, Rajasthan-333031, India  Rafael C. Nunes
<EMAIL_ADDRESS>Divisão de Astrofísica, Instituto Nacional de Pesquisas
Espaciais, Avenida dos Astronautas 1758, São José dos Campos, 12227-010, SP,
Brazil

###### Abstract

We investigate the Minimal Theory of Massive Gravity (MTMG) in the light of
different observational data sets which are in tension within the $\Lambda$CDM
cosmology. In particular, we analyze MTMG model, for the first time, with the
Planck-CMB data, and how these precise measurements affect the free parameters
of the theory. The MTMG model can affect the CMB power spectrum at large
angular scales and cause a suppression on the amplitude of the matter power
spectrum. We find that on adding Planck-CMB data, the graviton has a small,
positive, but non-zero mass at 68% confidence level, and from this
perspective, we show that the tension between redshift space distortions
measurements and Planck-CMB data in the parametric space $S_{8}-\Omega_{m}$
can be resolved within the MTMG scenario. Through a robust and accurate
analysis, we find that the $H_{0}$ tension between the CMB and the local
distance ladder measurements still remains but can be reduced to $\sim
3.5\sigma$ within the MTMG theory. The MTMG is very well consistent with the
CMB observations, and undoubtedly, it can serve as a viable candidate amongst
other modified gravity theories.

††preprint: YITP-21-55

##  I Introduction

Over the past two decades, a large volume of cosmological information has been
obtained by several data surveys, making the estimates of cosmological
parameters increasingly accurate. The understanding of the cosmological probes
has been well modeled through the standard $\Lambda$-Cold Dark Matter scenario
(the $\Lambda$CDM scenario). The $\Lambda$CDM cosmological model provides a
wonderful fit to the current cosmological data, but recently a few tensions
and anomalies became statistically significant while analyzing different data
sets. The most discussed tension in the literature is in the estimation of the
Hubble constant, $H_{0}$, between the CMB and the direct local distance ladder
measurements. Assuming the $\Lambda$CDM scenario, Planck-CMB data analysis
provides $H_{0}=67.4\pm 0.5$ km s-1Mpc-1 [1], which is in $4.4\sigma$ tension
with a cosmological model-independent local measurement $H_{0}=74.03\pm 1.42$
km s-1Mpc-1 [2] from the Hubble Space Telescope (HST) observations of 70 long-
period Cepheids in the Large Magellanic Cloud. Additionally, a combination of
time-delay cosmography from H0LiCOW lenses and the distance ladder
measurements are in $5.2\sigma$ tension with the Planck-CMB constraints [3].
Recently, based on a joint analysis from several geometrical probes, a tension
of $\sim 6\sigma$ was observed in [4] under the consideration of minimal
theoretical assumptions. Motivated by these observational discrepancies, it
has been widely discussed in the literature whether a new physics beyond the
standard cosmological model can solve the $H_{0}$ tension (see [5, 6, 7] and
reference therein for a review). In addition to the $H_{0}$ disagreement, a
significant tension between Planck data with redshift surveys data and weak
lensing measurements has been reported, about the value of the matter density
$\Omega_{m}$, and the amplitude of growth of the structures, $\sigma_{8}$,
also quantified in terms of the parameter:
$S_{8}=\sigma_{8}\sqrt{\Omega_{m}/0.30}$. A higher $S_{8}$ value is estimated
from CMB data assuming $\Lambda$CDM, namely, $S_{8}=0.834\pm 0.016$ [1] from
Planck data and $S_{8}=0.840\pm 0.030$ from ACT+WMAP joint analysis [8]. The
cosmic shear and redshift space distortions (RSD) measurements have predicted
a lower value of $S_{8}$. This tension is above the 2$\sigma$ level with
KiDS-450 ($S_{8}=0.745\pm 0.039$) [9], KiDS+VIKING-450
($S_{8}=0.737^{+0.040}_{-0.036}$) [10] and DES ($S_{8}=0.783\pm 0.021)$ [11].
The KiDS-1000 team reported a 3$\sigma$ tension with Planck-CMB [12]. The
tension becomes 3.2$\sigma$ if we consider the combination of VIKING-450 and
DESY1 [13] and 3.4$\sigma$ for BOSS+VIKING-450 ($S_{8}=0.728\pm 0.026)$ [14].
Also, in agreement with a lower value, there is an estimate from the BOSS
Galaxy Power Spectrum $S_{8}=0.703\pm 0.045$ [15]. See [16] and references
therein for a mini-review and additional information on $S_{8}$ estimations.
Although this tension could be due to systematic errors, it is worthwhile to
investigate the possibility of new physics beyond the standard model to
explain the $S_{8}$ tension [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 23, 30, 31]. Additionally, the authors in [24] identified a large tension
between RSD and CMB measurements. Disagreement between CMB and combinations of
RSD measurements with other datasets, including the $E_{G}$ statistic, is
discussed in [32], pointing out a tension up to 5$\sigma$, depending on the
datasets.

On the other hand, there are theoretical and observational reasons to believe
that general relativity (GR) should be modified when gravitational fields are
strong and/or on large scales. From an observational point of view, the
physical mechanism responsible for accelerating the expansion of Universe at
late times is still an open question, and it has been intensively investigated
whether modified gravity scenarios can explain such an accelerated stage, as
well as to fit very well with the observational data from different sources
(see [33] for review). Also, theories beyond GR can serve as alternatives to
explain the current $H_{0}$ tension [34, 35, 36, 37, 38, 39]. One of the most
interesting possibilities for modification of gravity is to give a mass to the
graviton (see [40] for review about massive gravity theories). In this work,
we investigate in detail the observational feasibility of the Minimal Theory
of Massive Gravity (MTMG) [41, 42], due to its infrared Lorentz violations
measurable at cosmological scales. We robustly constrain the MTMG framework
using observational probes of the $f\sigma_{8}$ growth rate of cosmological
perturbations, which is a useful bias free statistical cosmic test, as well as
we combine these growth rate data with several geometric observations such as
the most recent data from Type Ia supernovae and baryon acoustic oscillations.
We derive and discuss observational perspectives from how the MTMG scenario
predict the statistical plane for $S_{8}-\Omega_{m}$ in direct comparison with
the $\Lambda$CDM model, as well as the potential of the model in alleviating
the $H_{0}$ tension. Furthermore, we investigate for the first time how the
Planck-CMB data constrains the MTMG baseline model. Our results also provide
new constraint on the graviton mass from Planck-CMB data, obtained in the MTMG
context. We find that the use of Planck-CMB data selects parameter space for
MTMG which is compatible with a non-zero but positive value of
$\mu^{2}/H_{0}^{2}=0.25^{{+}0.16}_{{-}0.10}$, $\mu$ being the mass of the
graviton. This is probably the newest and most interesting result from this
study.

This paper is structured as follows. In Section II, we review and introduce
the MTMG scenario. In Section III, we present the data sets and methodology
used in this work. In section IV, we discuss the main results of our analysis.
In section V, we summarize main findings of this study.

##  II Minimal Theory of Massive Gravity

The model we study here is a model which has been introduced in order to fix
the cosmology of massive gravity, i.e., to have a stable FLRW background
without ghosts or strong-coupling issues. The model has been introduced in
order to study a standard cosmological phenomenology to the theories with a
massive graviton. In order for this goal to be achieved, the theory in its
simplest form (i.e., in vacuum), has only two gravitational degrees of freedom
(corresponding to the two polarization gravitational waves). Furthermore, on
FLRW, MTMG shares the same background of dRGT [43, 44]. This fact was
implemented by construction in MTMG. Therefore also MTMG is endowed, as dRGT,
with two branches which are named as 1) the self-accelerated branch; and 2)
the normal branch.

These two branches have, in general, different FLRW background and
perturbation dynamics. In particular, the self-accelerated branch has a
phenomenology which is identical in everything to the standard $\Lambda$CDM
except for the propagation of tensor modes which are now massive (compared to
GR). On the other hand, in the normal branch we have the possibility of giving
a general dynamics for the graviton mass. This allows for many
phenomenological possibilities. However, in order to study this branch, it is
somehow convenient to study its simplest realization, namely the one for which
the effective energy density of the MTMG component becomes a constant during
the whole evolution.

More in detail, MTMG relies on the presence and choice of a three-dimensional
fiducial metric which, in unitary gauge, will be written here as
${\tilde{\gamma}}_{ij}=\tilde{a}^{2}\,\delta_{ij}$, where
$\tilde{a}=\tilde{a}(t)$ is the fiducial scale factor. On top of that, there
is need of the fiducial lapse, $M=M(t)$ and of another 3D tensor
$\tilde{\zeta}^{i}{}_{j}$, which represents the time variation of the fiducial
vielbein (for details, see [41, 42]). Out of the fiducial and physical
metrics, one builds up the following quantities:

$\displaystyle\mathcal{K}^{k}{}_{n}$ $\displaystyle=$
$\displaystyle(\sqrt{{\tilde{\gamma}}^{-1}\gamma})^{k}{}_{n}\,,$ (1)
$\displaystyle\mathfrak{K}^{k}{}_{n}$ $\displaystyle=$
$\displaystyle(\sqrt{{\gamma}^{-1}{\tilde{\gamma}}})^{k}{}_{n}\,,$ (2)
$\displaystyle\mathfrak{K}^{k}{}_{n}\mathcal{K}^{n}{}_{m}$ $\displaystyle=$
$\displaystyle\delta^{k}{}_{m}=\mathcal{K}^{k}{}_{n}\mathfrak{K}^{n}{}_{m}\,.$
(3)

Out of these quantities, one defines the following symmetric tensor

$\displaystyle\Theta^{ij}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\tilde{\gamma}}}{\sqrt{\gamma}}\\{c_{1}(\gamma^{il}\mathcal{K}^{j}{}_{l}+\gamma^{jl}\mathcal{K}^{i}{}_{l})$
(4) $\displaystyle+$ $\displaystyle
c_{2}[\mathcal{K}(\gamma^{il}\mathcal{K}^{j}{}_{l}+\gamma^{jl}\mathcal{K}^{i}{}_{l})-2{\tilde{\gamma}}^{ij}]\\}+2c_{3}\gamma^{ij}\,,$

where $c_{1,2,3}$ are constants. We then need the following building blocks

$\displaystyle\mathcal{C}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,m^{2}\,M\,K_{ij}\Theta^{ij}$ (5) $\displaystyle-$
$\displaystyle
m^{2}M\left\\{\frac{\sqrt{\tilde{\gamma}}}{\sqrt{\gamma}}[c_{1}\tilde{\zeta}+c_{2}(\mathcal{K}\zeta-\mathcal{K}^{m}{}_{n}{\tilde{\zeta}}^{n}{}_{m})]\right.$
$\displaystyle+$
$\displaystyle\left.c_{3}\mathfrak{K}^{m}{}_{n}{\tilde{\zeta}}^{n}{}_{m}\right\\},$
$\displaystyle\mathcal{C}^{n}{}_{i}$ $\displaystyle=$
$\displaystyle-m^{2}M\left\\{\frac{\sqrt{\tilde{\gamma}}}{\sqrt{\gamma}}\bigl{[}\tfrac{1}{2}(c_{1}+c_{2}\mathcal{K})(\mathcal{K}^{n}{}_{i}+\gamma^{nm}\mathcal{K}^{l}{}_{m}\gamma_{li})\right.$
(6) $\displaystyle-$
$\displaystyle\left.c_{2}{\tilde{\gamma}}^{nl}\gamma_{li}\bigr{]}+c_{3}\delta^{n}{}_{i}\right\\},$

where $\gamma$ and $\tilde{\gamma}$ represent the determinants of
$\gamma_{ij}$ and ${\tilde{\gamma}}_{ij}$, respectively. Furthermore $K_{ij}$
is the extrinsic curvature of $\gamma_{ij}$, namely
$K_{ij}=\frac{1}{2N}({\dot{\gamma}}_{ij}-\mathcal{D}_{i}N_{j}-\mathcal{D}_{j}N_{i})$,
where $N$ is the lapse, $N^{i}$ the shift vector and $\mathcal{D}$ is the
covariant derivative compatible with the 3D-metric $\gamma_{ij}$. Finally,
$\mathcal{K}\equiv\mathcal{K}^{n}{}_{n}$ and
${\tilde{\zeta}}\equiv{\tilde{\zeta}}^{n}{}_{n}$. We are now ready to define
the action of MTMG as follows

$\displaystyle S$ $\displaystyle=$ $\displaystyle S_{\rm pre}+\frac{M_{\rm
P}^{2}}{2}\int
d^{4}xN\sqrt{\gamma}\left(\frac{m^{2}}{4}\frac{M}{N}\,\lambda\right)^{\\!2}(\Theta_{ij}\Theta^{ij}-\tfrac{1}{2}\Theta^{2})$
(7) $\displaystyle-$ $\displaystyle\frac{M_{\rm P}^{2}}{2}\int
d^{4}x\sqrt{\gamma}[\lambda\mathcal{C}_{0}-(\mathcal{D}_{n}\lambda^{i})\mathcal{C}^{n}{}_{i}]+S_{\rm
mat}\,,$

where $S_{\rm mat}$ represents the standard matter action. We still need to
define the last bits of the theory, which are given by

$\displaystyle S_{\rm pre}$ $\displaystyle=$ $\displaystyle S_{\rm
GR}+\frac{M_{\rm P}^{2}}{2}\sum_{i=1}^{4}\int d^{4}x\mathcal{S}_{i}\,,$ (8)
$\displaystyle S_{\rm GR}$ $\displaystyle=$ $\displaystyle\frac{M_{\rm
P}^{2}}{2}\int d^{4}xN\sqrt{\gamma}[{}^{(3)}R+K^{ij}K_{ij}-K^{2}]\,,$ (9)
$\displaystyle\mathcal{S}_{1}$ $\displaystyle=$
$\displaystyle-m^{2}c_{1}\sqrt{\tilde{\gamma}}(N+M\mathcal{K})\,,$ (10)
$\displaystyle\mathcal{S}_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\,m^{2}c_{2}\sqrt{\tilde{\gamma}}(2N\mathcal{K}+M\mathcal{K}^{2}-M{\tilde{\gamma}}^{ij}\gamma_{ji})\,,$
(11) $\displaystyle\mathcal{S}_{3}$ $\displaystyle=$
$\displaystyle-m^{2}c_{3}\sqrt{\gamma}(M+N\mathfrak{K})\,,$ (12)
$\displaystyle\mathcal{S}_{4}$ $\displaystyle=$
$\displaystyle-m^{2}c_{4}\sqrt{\gamma}N\,.$ (13)

The theory has been defined on the unitary gauge (for the Stuckelberg fields),
so we cannot make any further gauge choice. In what follows, we describe the
background and perturbation analysis by using the following variables. For the
scalar modes, the flat 3D metric is then written according to

$ds_{3}^{2}=\gamma_{ij}dx^{i}dx^{j}=[a(t)^{2}\,(1+2\zeta)\,\delta_{ij}+2\partial_{i}\partial_{j}E]\,dx^{i}dx^{j}\,,$
(14)

$a$ being the scale factor of the physical metric. Then we can define the
lapse and shift vector as

$\displaystyle N$ $\displaystyle=$ $\displaystyle N(t)\,(1+\alpha)\,,$ (15)
$\displaystyle N_{i}$ $\displaystyle=$ $\displaystyle
N(t)\,\partial_{i}\chi\,.$ (16)

We also need to give

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\lambda(t)+\delta\lambda\,,$ (17) $\displaystyle\lambda^{i}$
$\displaystyle=$
$\displaystyle\frac{1}{a^{2}}\delta^{ij}\partial_{j}\delta\lambda_{V}\,.$ (18)

Finally, we need to give all matter actions together with their own variables.

It is possible to show that on the background, $\lambda(t)=0$. In this case,
the dynamics of MTMG imposes the following constraint

$\left(HX-H\,\frac{M}{N}+\frac{\dot{X}}{N}\right)({X}^{2}c_{{1}}+2\,Xc_{{2}}+c_{{3}})=0\,,$
(19)

where we have defined the variable $X\equiv\tilde{a}/a$, being the ratio
between the fiducial and the physical scale factors. The function
$H=\dot{a}/(aN)$ is the Hubble parameter, and the coefficients $c_{1,2,3}$ are
instead constant coefficients. On solving Eq. (19), we can see that in this
setup only two solutions are possible: 1) setting
${X}^{2}c_{{1}}+2\,Xc_{{2}}+c_{{3}}=0$, which defines the self-accelerating
branch and leads to $X={\rm const.}$; 2) the normal branch for which
$M=NX+\dot{X}/H$ where the time dependence of $X$ is still not imposed. Here
onwards, we will only focus on the normal branch and, in particular, we fix
the dynamics of $X$ so that $X=X_{0}={\rm const.}$

As a consequence of this choice, the background equations of MTMG read as
follows:

$\displaystyle 3M_{\rm P}^{2}H^{2}$ $\displaystyle=$
$\displaystyle\sum_{I}\rho_{I}+\rho_{\rm MTMG}\,,$ (20) $\displaystyle 2M_{\rm
P}^{2}\,\frac{\dot{H}}{N}$ $\displaystyle=$
$\displaystyle-\sum_{I}(\rho_{I}+P_{I})\,,$ (21)
$\displaystyle\frac{\dot{\rho}_{I}}{N}$ $\displaystyle=$
$\displaystyle-3\,H\,(\rho_{I}+P_{I})\,,$ (22)

where $\rho_{\rm MTMG}\equiv\frac{{m}^{2}M_{\rm
P}^{2}}{2}(c_{1}X_{0}^{3}+3c_{2}X_{0}^{2}+3c_{3}X_{0}+c_{4})$ is a constant.
Here the sum over the index $I$ is over all the standard matter components. In
this case, we can see that the background equations of motion exactly reduce
to the ones of $\Lambda$CDM (endowed with an effective cosmological constant,
i.e., the model is self-accelerating).

Let us now describe the dynamics of the perturbation equations of motion. We
will follow here a procedure similar to the one introduced in [45] and [46,
34]. We will have the perturbation equations of motion coming from the gravity
sector (the ones for $\alpha$, $\chi$, $\zeta$, $E$, $\delta\lambda$ and
$\delta\lambda_{V}$) together with the ones coming from each matter component.
We then find the equations of motion for each of the fields we have. In the
following, we will (and can) set $N(t)=a(t)$ as to have dynamics in terms of
the conformal time. Although we cannot fix any gauge for the perturbations, we
can still perform field redefinitions. In particular we will introduce the
following ones:

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\psi-\frac{\dot{\chi}}{a}+\frac{1}{a}\,\frac{d}{dt}\\!\left[a\,\frac{d}{dt}\\!\left(\frac{E}{a^{2}}\right)\right],$
(23) $\displaystyle\zeta$ $\displaystyle=$
$\displaystyle-\phi-H\,\chi+aH\frac{d}{dt}\\!\left(\frac{E}{a^{2}}\right).$
(24)

These definitions correspond to the usual gauge-invariant definitions for the
Bardeen potentials $\phi$ and $\psi$. Then we will use the Newtonian gauge
invariant fields, $\phi$ and $\psi$, as to rewrite the modified Einstein
equations. For each matter component, labeled with $I$, we also introduce
gauge invariant combinations $\delta_{I}$’s and $\theta_{I}$’s (in place of
the density contrast $\delta\rho_{I}/\rho_{I}$ and the scalar speed
$u_{Ii}=\partial_{i}v_{I}$, respectively) as follows:

$\displaystyle\frac{\delta\rho_{I}}{\rho_{I}}$ $\displaystyle=$
$\displaystyle\delta_{I}-\frac{\dot{\rho}_{I}}{a\rho_{I}}\,\chi+\frac{\dot{\rho}_{I}}{\rho_{I}}\,\frac{d}{dt}\\!\left(\frac{E}{a^{2}}\right),$
(25) $\displaystyle v_{I}$ $\displaystyle=$
$\displaystyle-\frac{a}{k^{2}}\,\theta_{I}+\chi-a\,\frac{d}{dt}\\!\left(\frac{E}{a^{2}}\right).$
(26)

In terms of these variables the equations of motion for each matter component,
for the lowest multipoles, read as follows:

$\displaystyle\dot{\delta}_{I}$ $\displaystyle=$ $\displaystyle
3a(w_{I}-c_{sI}^{2})H\delta_{I}-(1+w_{I})\,(\theta_{I}-3\dot{\phi})\,,$ (27)
$\displaystyle\dot{\theta}_{I}$ $\displaystyle=$ $\displaystyle
aH(3c_{sI}^{2}-1)\,\theta_{I}+{k}^{2}\psi+{\frac{c_{sI}^{2}k^{2}}{1+w_{I}}}\,\delta_{I}-{k}^{2}\sigma_{I}\,,$
(28)

where $w_{I}\equiv P_{I}/\rho_{I}$, and $c_{sI}^{2}=\left(\frac{\partial
p_{I}}{\partial\rho_{I}}\right)_{s}$ is the speed of propagation for each
matter species (i.e., $\dot{p}_{I}/\dot{\rho}_{I}$, which vanishes for dust
and equals 1/3 for photons). This result merely shows that the dynamical
equations of motion for the matter components exactly coincide with the same
ones in $\Lambda$CDM. This is not a surprise as modifications of gravity only
enter in the gravity sector.

The equations of motion for $\delta\lambda_{V}$ impose

$\chi=-\frac{\phi}{H}+\frac{\dot{E}}{a}-2HE\,.$ (29)

We can solve the equations of motion for $\chi$ in terms of $\delta\lambda$
(and other fields), and the equations of motion for $\alpha$ in terms of $E$.
We can consider a linear combination of the equations of motion for $\zeta$
and $E$ in order to find a solution for $\delta\lambda_{V}$.

At this point, the equation of motion for $\delta\lambda$ gives

$\displaystyle\mathcal{E}_{0}$ $\displaystyle\equiv$
$\displaystyle\left[{\frac{(Y\theta-2){k}^{2}}{{a}^{2}}}+{\frac{9\,\theta\,Y\Gamma}{2}}\right]\phi-3\sum_{I}\varrho_{I}\delta_{I}$
(30) $\displaystyle+$
$\displaystyle{\frac{9aH(Y\theta-2)}{2{k}^{2}}}\,\sum_{I}\Gamma_{I}\theta_{I}=0\,,$

where we have introduced $\rho_{I}=3M_{\rm P}^{2}\varrho_{I}$, $P_{I}=3M_{\rm
P}^{2}p_{I}$, $\Gamma_{I}\equiv\varrho_{I}+p_{I}$,
$\Gamma=\sum_{I}\Gamma_{I}$, $Y=H_{0}^{2}/H^{2}$ and
$\theta\equiv\frac{m^{2}X_{0}}{2H_{0}^{2}}(c_{1}X_{0}^{2}+2c_{2}X_{0}+c_{3})$,
which represents the only extra parameter which enters into the equations of
motion in addition to the background effective cosmological constant. In the
$\Lambda$CDM limit, whenever $|\theta Y|\ll 1$, we recover the GR standard
evolution.

The equation of motion for $E$, the last one available in the gravity sector,
gives

$\displaystyle\mathcal{E}_{1}$ $\displaystyle\equiv$
$\displaystyle\dot{\phi}+{\frac{3\,Y\theta\,a\Gamma\,\phi}{2\,H\left(Y\theta-2\right)}}+aH\psi+{\frac{3\,{a}^{2}}{{k}^{2}(Y\theta-2)}}\sum_{I}\Gamma_{I}\theta_{I}$
(31) $\displaystyle=$ $\displaystyle 0\,.$

On taking the time-derivative of $\mathcal{E}_{0}$ and removing $\dot{\phi}$
from $\mathcal{E}_{1}$, we find the following derived equation of motion:

$\displaystyle\mathcal{E}_{2}$ $\displaystyle\equiv$
$\displaystyle\psi+{\frac{9\,{a}^{2}}{2\,{k}^{2}}}\sum_{I}\Gamma_{I}\sigma_{I}-{\frac{9{a}^{2}\,\theta\,Y}{{k}^{2}\left(2\,Y\theta-4\right)}}\sum_{I}c_{sI}^{2}\varrho_{I}\delta_{I}$
(32) $\displaystyle-$
$\displaystyle\left[1+{\frac{3\theta\,Y\Gamma}{\left(2\,Y\theta-4\right){H}^{2}}}\right.$
$\displaystyle+$ $\displaystyle\left.{\frac{27Y\theta
a^{2}}{2\,{k}^{2}(2-Y\theta)}\left(\sum_{J}c_{sJ}^{2}\Gamma_{J}-{\frac{{\Gamma}^{2}}{2{H}^{2}}}\right)}\right]\phi$
$\displaystyle-$
$\displaystyle{\frac{27\,\theta\,Y\Gamma\,{a}^{3}}{4\,{k}^{4}(Y\theta-2)H}}\sum_{I}\Gamma_{I}\theta_{I}=0\,,$

where we have also used the matter equations of motion in order to replace
$\dot{\delta}_{I}$, etc. In the limit $|\theta Y|\to 0$, the above equations
reduce to the standard GR form. We will use, as in GR, both $\mathcal{E}_{1}$
and $\mathcal{E}_{2}$ as the dynamical equations in the gravity sector. The
fact that we do not have any additional dynamical equation of motion is a
consequence of the fact that MTMG does not add any new dynamical degree of
freedom in the theory.

Although the previous equations, never appeared before in the literature,
completely define the behavior of gravity in MTMG, at any time and scale, we
finish the description of the theory by considering the behavior of the theory
under the influence of a single dust fluid at late times (i.e., neglecting
radiation), in order to understand the effect of MTMG for the effective
gravitational constant. Since this last step has already been performed in the
literature (see e.g. [42, 47, 48]), this last calculation can then be
considered to be a check of the calculations presented so far in this paper.
Then on considering a dust fluid (mimicking CDM or late-time baryons), Eqs.
(27) and (28) reduce to

$\displaystyle\mathcal{E}_{D}$ $\displaystyle\equiv$
$\displaystyle{\dot{\delta}}_{c}+\theta_{c}-3\dot{\phi}=0\,,$ (33)
$\displaystyle\mathcal{E}_{V}$ $\displaystyle\equiv$
$\displaystyle{\dot{\theta}}_{c}+aH\theta_{c}-k^{2}\psi=0\,.$ (34)

We can then solve Eq. (32) for $\psi$ in terms of $\phi$ and $\theta_{c}$
(dust having no shear, and vanishing $c_{s}^{2}$). On replacing, by using Eq.
(31), $\dot{\phi}$ into Eq. (33), we can then solve $\mathcal{E}_{D}=0$ for
$\theta_{c}$ in terms of $\phi$ and ${\dot{\delta}}_{c}$. We can now replace,
$\theta_{c}$, $\psi$ and $\dot{\phi}$ (coming from ${\dot{\theta}}_{c}$)
inside Eq. (34), and solve $\mathcal{E}_{V}=0$ for $\phi$ in terms of
${\dot{\delta}}_{c}$ and ${\ddot{\delta}}_{c}$. Finally, on replacing $\phi$
and $\theta_{c}$ in Eq. (30), we find a closed second order differential
equation for $\delta_{c}$. On evaluating this last equation for subhorizon
scales, that is for $k/(aH)\gg 1$, we find

${\ddot{\delta}}_{c}+aH{\dot{\delta}}_{c}-\frac{3}{2}\,\frac{G_{\rm
eff}}{G_{N}}\,\varrho_{c}\,a^{2}\,\delta_{c}=0\,,$ (35)

where $M_{\rm P}^{2}=(8\pi G_{N})^{-1}$, $\Omega_{c}=\varrho_{c}(t)/H^{2}$,
and

$\frac{G_{\rm eff}}{G_{N}}=\frac{2}{2-\theta Y}-\frac{3\theta
Y\Omega_{c}}{(2-\theta Y)^{2}}\,,$ (36)

matching the result found, e.g. in [47]. The same dynamical equation holds
also for the late-time baryonic fluid. Since we are going to use bounds on
ISW-galaxy cross correlations, it is worthy to find the expression for the
ISW-field, defined as $\psi_{\rm ISW}\equiv\phi+\psi$, in the subhorizon
scales approximation. We find that

$\displaystyle\psi_{\rm ISW}$ $\displaystyle=$
$\displaystyle-\frac{3H_{0}^{2}\Omega_{m0}}{k^{2}}\,\frac{\Sigma\,\delta_{m}}{a}\,,$
(37) $\displaystyle\Sigma$ $\displaystyle\equiv$
$\displaystyle\frac{8-\left(4+3\Omega_{m}\right)\theta Y}{2\left(2-\theta
Y\right)^{2}},$ (38)

where here $\Omega_{m}=\Omega_{m}(t)=\rho_{m}/(3M_{\rm P}^{2}H^{2})$, and
$\Omega_{m0}=\Omega_{m}(z=0)$ (in the next sections, unless specified
otherwise, $\Omega_{m}$ will be used to rewrite $\Omega_{m0}$). This result
for $\psi_{\rm ISW}$ agrees with the one found in [48].

Let us finally discuss the equations of motion for the gravitational waves
propagating on this flat FLRW background (the vector modes do not show any
deviation from GR). In fact, the contribution of MTMG to the tensor modes
consists, by construction (see e.g. [41, 42]), of introducing a non-zero
squared mass $\mu^{2}=H_{0}^{2}\theta$ to them. The equations of motion for
the two tensor modes polarizations ($h_{f}$ labeled with
$f\in\\{{+},{\times}\\}$) become

${\ddot{h}}_{f}=-2\,\frac{\dot{a}}{a}\,{\dot{h}}_{f}-(k^{2}+\theta\,H_{0}^{2}\,a^{2})\,h_{f}\,,$
(39)

where in the right hand side a source term can further be added.

Figure 1: Deviations of the CMB TT power spectrum for $\Lambda$CDM from MTMG
for some values of $\theta$.

Figure 2: Left panel: Deviations of the matter power spectrum, $\Delta
P(k,z)=|P(k,z)^{MTMG}-P(k,z)^{GR}|/(P(k,z)^{MTMG}+P(k,z)^{GR})$, for
$\Lambda$CDM from MTMG with $\theta=-0.2$ for some values of $z$. Right panel:
Same as in left panel, but assuming $\theta=0.2$.

We implemented the model described above in CLASS code [49]. In Figure 1, we
quantify the MTMG model affects on the CMB temperature power spectrum, i.e.,
$C_{l}^{TT}$, with respect to GR-$\Lambda$CDM prediction. We can see that
increasing $\theta$ enhances the CMB temperature power spectrum at low $l$.
This effect is well known for dark energy and MG models, and it is due to the
integrated Sachs-Wolfe effect, which affects the CMB spectrum at low $l$, but
has no significant effect at large scales. For scales larger than $l>100$, we
note $C_{l}^{GR,TT}\approx C_{l}^{MTMG,TT}$. In Figure 2, we show the
difference on the $P(k,z)$ at four different redshift values on the scales
$k\in[10^{-4},1]$ h/Mpc. These $z$ values are chosen because our RSD sample
covers the range $z\in[0.02,1.94]$. For a quantitative example, taking
$\theta=-0.2$, we find that the difference with respect to $\Lambda$CDM is
$\sim 25$% at $k=10^{-3}$ h/Mpc, while for $k>10^{-2}$ h/Mpc the difference is
$<0.2$%, for all $z$ values. For the case with $\theta=0.2$, we note that for
$k>10^{-3}$ h/Mpc, we have a similar difference as in the previous case, but
with different $P(k,z)$ behavior at different $z$ values, as expected, since
the dynamics depends on the signal in $\theta$. For scales $k<10^{-3}$,
significant deviations are noted. It is important to mention that on very
large scale, the cosmic variance effects can be significant. In these
simulations, the predictions of the non-linear effects are performed using the
HMcode code [50]. Non-linear dynamics of the MTMG framework using N-body
simulations were investigated in detail in [51].

Table 1: Constraints at 95% CL on $\Omega_{m}$, $\sigma_{8}$, $S_{8}$ and $\theta=\mu^{2}/H_{0}^{2}$, inferred from different data set combinations, in the MTMG model. data set | $\theta$ | $\Omega_{m}$ | $\sigma_{8}$ | $S_{8}$
---|---|---|---|---
RSD | $-5.2^{+7.1}_{-4.8}$ | $0.37^{+0.17}_{-0.16}$ | $0.808^{+0.086}_{-0.084}$ | $0.89^{+0.19}_{-0.20}$
RSD + BAO + Pantheon | $-1.8^{+3.3}_{-4.4}$ | $0.305^{+0.041}_{-0.039}$ | $0.796^{+0.091}_{-0.073}$ | $0.802^{+0.10}_{-0.096}$
RSD + BAO + Pantheon + ISW | $-0.12^{+0.28}_{-0.26}$ | $0.293^{+0.018}_{-0.018}$ | $0.775^{+0.055}_{-0.055}$ | $0.766^{+0.057}_{-0.055}$

Figure 3: Left panel: 68% and 95% confidence levels on the $S_{8}-\Omega_{m}$
plane from various data sets in $\Lambda$CDM and MTMG models. Right panel: 68%
and 95% confidence levels on the $\theta-\Omega_{m}$ plane from various data
sets in MTMG model. GR is recovered for $\theta=0$.

##  III Data and Methodology

In this section, we present the data sets and methodology used to obtain the
observational constraints on the model parameters by performing Bayesian MCMC
analysis. In order to constrain the parameters, we use the following data
sets.

CMB: We use the full Planck-2018 [1] CMB temperature and polarization data
which comprise of the low-l temperature and polarization likelihoods at $l\leq
29$, temperature (TT) at $l\geq 30$, polarization (EE) power spectra, and
cross correlation of temperature and polarization (TE) while also including
the Planck-2018 CMB lensing power spectrum likelihood [52] in our analysis.

RSD: The growth rate data comprise of $f(z)\sigma_{8}(z)$ measurements from a
variety of redshift space distortion surveys. The current measurements of
$f\sigma_{8}(z)$ come from a plethora of different surveys with different
assumptions and systematics, thus an approach to study the statistical
properties and robustness of the data is imperative. In [53], an internal
robustness analysis was used to validate a subset of $f\sigma_{8}(z)$
measurements. In our analysis, we use the compilation of the $f\sigma_{8}(z)$
measurements presented in Table I in [53], which are minimally affected by
systematic-contaminated data-points. We refer this data set to as the data set
of RSD measurements. To build the likelihood function, we follow the same
methodology as presented in [54]. It is important to highlight that the RSD
sample is still fiducial to some cosmology, for example fixed on a
$\Lambda$CDM baseline, different assumptions on the reference value of
$\Omega_{m}$, and non-linearities modeling of which starts to play an
important role on smaller scales and at later epochs. These points should be
taken into account to accurately estimate the cosmological parameters. Our
sample is minimally affected by systematic-contaminated $f\sigma 8$ data
points, checked through Bayesian model comparison framework described in [53],
identifying potential outliers as well as subsets of data affected by
systematics or new physics.

BAO: The Baryon Acoustic Oscillations (BAO) provide an important cosmological
probe, which can trace expanding spherical wave of baryonic perturbations from
acoustic oscillations at recombination time through the large-scale structure
correlation function. In this work, we consider the most recent BAO data
compilation comprised of the $D_{V}(z)/r_{d}$, $D_{M}(z)/r_{d}$, and
$D_{H}(z)/r_{d}$ measurements compiled in Table 3 in [55], regarding BAO-only
data.

Pantheon: The Supernovae Type Ia have traditionally been one of the most
important astrophysical tools in establishing the so-called standard
cosmological model. For the present analysis, we use the Pantheon compilation,
which consists of 1048 SN Ia distributed in the redshift range
$z\in[0.01,2.3]$ [56].

ISW: The late-time Integrated Sachs-Wolfe effect (ISW) on the CMB is an effect
imprinted in the angular pattern of the CMB in the presence of a time varying
cosmological gravitational potential, which can be due to a non-flat universe
[57], as well as for a flat one in the presence of dark energy or modified
gravity theories [58, 59, 60]. A non-zero ISW necessarily implies the presence
of a generating physical source for the accelerated expansion of the universe
at late times. For the present analysis, we use the cross correlation of the
CMB with galaxy surveys that derive constraints on the ISW as obtained in
[61]. We use the five catalogs of extragalactic sources as presented in [61],
namely, the 2MASS Photometric Redshift catalog, the WISE × SuperCOSMOS photo-z
catalog, the NRAO VLA Sky Survey radio sources catalog, the SDSS DR12 and SDSS
DR6 QSO photometric catalogs.

We use the Metropolis-Hastings mode in CLASS+MontePython code [49, 62, 63] to
derive the constraints on cosmological parameters using various data
combinations from the data sets described above, ensuring a Gelman-Rubin
convergence criterion of $R-1<10^{-3}$. In what follows, we describe our main
results.

##  IV Results and discussions

We divide the analysis into two parts. First, we consider RSD and RSD + BAO +
Pantheon data, that is, the growth data and its combination with the geometric
data. Thus, we can quantify the dynamics of the model and the constraints on
$\theta$ without CMB data influence. It is important to check the prediction
of the model about $\sigma_{8}$ or $S_{8}$ in the absence of the CMB data.
Then, we add ISW information to these data. In the second round of analysis,
we analyze the model with CMB data, and discuss the potential of the model to
solving and/or alleviating the $H_{0}$ and $S_{8}$ tensions.

Table 2: Summary of the 1$\sigma$ constraints on the baseline parameters of the MTMG scenario from CMB and its combination with several other data sets. In the last column, the Full Joint Analysis means CMB + BAO + Pantheon + RSD + ISW. The parameter $H_{\rm 0}$ is measured in the units of km s-1 Mpc-1. Parameter | CMB | CMB + BAO + Pantheon | CMB + BAO + Pantheon + RSD | Full Joint Analysis
---|---|---|---|---
$10^{2}\omega_{b}$ | $2.242\pm 0.015$ | $2.245\pm 0.013$ | $2.252\pm 0.013$ | $2.247\pm 0.013$
$\omega_{\rm cdm}$ | $0.1195\pm 0.0012$ | $0.11903\pm 0.00090$ | $0.11819\pm 0.00090$ | $0.11886\pm 0.00095$
$100\theta_{\rm s}$ | $1.04195\pm 0.00029$ | $1.04197\pm 0.00028$ | $1.04202\pm 0.00028$ | $1.04197\pm 0.00028$
$\ln 10^{10}A_{\rm s}$ | $3.042\pm 0.014$ | $3.042\pm 0.014$ | $3.035\pm 0.014$ | $3.036\pm 0.015$
$n_{\rm s}$ | $0.9663\pm 0.0044$ | $0.9672\pm 0.0036$ | $0.9692\pm 0.0039$ | $0.9681\pm 0.0037$
$\tau_{\rm reio}$ | $0.0538\pm 0.0071$ | $0.0542\pm 0.0071$ | $0.0520\pm 0.0070$ | $0.0517\pm 0.0075$
$\theta$ | $0.19^{+0.14}_{-0.10}$ | $0.21^{+0.13}_{-0.10}$ | $0.26^{+0.14}_{-0.10}$ | $0.129^{+0.13}_{-0.088}$
$\Omega_{\rm m}$ | $0.3115\pm 0.0077$ | $0.3089\pm 0.0054$ | $0.3038\pm 0.0054$ | $0.3018\pm 0.0056$
$H_{0}$ | $67.65\pm 0.57$ | $67.84\pm 0.40$ | $68.22\pm 0.42$ | $68.44\pm 0.43$
$\sigma_{8}$ | $0.8071\pm 0.0064$ | $0.8057\pm 0.0065$ | $0.8001\pm 0.0060$ | $0.8169\pm 0.0068$
$S_{8}$ | $0.822\pm 0.014$ | $0.818\pm 0.011$ | $0.805\pm 0.010$ | $0.819\pm 0.011$

Figure 4: Left panel: 68% and 95% CLs in the $S_{8}-H_{0}$ plane from various
data sets in the MTMG model. Right panel: Same as the left panel, but in the
$\theta-\Omega_{m}$ plane. GR is recovered for $\theta=0$.

###  IV.1 Analysis with the growth, geometric and ISW data

Model-wise, we consider $\Lambda$CDM+$\theta$ model baseline, spanned by the 5
parameters: the Hubble constant $H_{0}$ or equivalently the reduced Hubble
constant $h\equiv H_{0}/(100\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1})$, the
physical baryon density $\omega_{b}\equiv\Omega_{b}h^{2}$, the physical cold
dark matter density $\omega_{c}\equiv\Omega_{c}h^{2}$, $\sigma_{8}$ and
$\theta$, quantifying deviations from GR. The matter density parameter today
$\Omega_{m}$ is derived by $\Omega_{m}=\Omega_{b}+\Omega_{c}$. Another
important derived parameter is $S_{8}\equiv\sigma_{8}\sqrt{\Omega_{m}/0.3}$.
To constrain the physical baryon density, we adopt a Gaussian prior on
$\omega_{b}$ from Big Bang Nucleosynthesis (BBN): $100\omega_{b}=2.233\pm
0.036$ [64]. In order that the tensor modes not to develop instability whose
time scale is shorter than the age of the universe, we impose the prior
$\theta\geqslant-10$. In the limit $\theta=0$, we recover the evolution
equation for the perturbations as in GR. Note that the $\Lambda$CDM and MTGM
models are indistinguishable at the background evolution. Thus, the data
combination BAO + Pantheon is indistinguishable for both scenarios. We use BAO
+ Pantheon to break the statistical degeneracy on $\Omega_{m}$, that may be
present when considering RSD data only. Since both $\Lambda$CDM and MTGM
models can predict different constraints on $\Omega_{m}$ in light of only the
RSD data, breaking this background degeneracy with BAO + Pantheon can
significantly improve the results. Finally, we also consider ISW sample data.
Table 1 summarizes the main results of our statistical analyses.

Figure 3, on the left panel, shows the parameter space in the
$S_{8}-\Omega_{m}$ plane at 68% and 95% CL with various data sets in
$\Lambda$CDM and MTMG models. First, we note that, in direct comparison with
the constraints provided with $\Lambda$CDM, the MTMG model can provide an
enlargement in both $S_{8}$ and $\Omega_{m}$ estimations. The right panel
shows the regions at 68% and 95% CL in the $\theta-\Omega_{m}$ plane with RSD
and RSD + BAO + Pantheon data in the MTMG model. We find
$\theta=-5.2^{+7.1}_{-4.8}$, $\theta=-1.8^{+3.3}_{-4.4}$ and
$\theta=-0.12^{+0.28}_{-0.26}$ at 95% CL from RSD, RSD + BAO + Pantheon and
RSD + BAO + Pantheon + ISW data, respectively. The addition of the ISW data
improves the constraints on $\theta$ because $\theta$ can significantly affect
the late-time ISW effect [48]. So, ISW data improve $\theta$, and consequently
by correlation, improvements on the other parameters in the baseline are also
reached, which are already well constrained from other data-sets. The joint
analysis with RSD + BAO + Pantheon + ISW data provides robust and accurate
constraints on the theory using only late time probes. As we will show in the
next section, the constraints are further improved with the addition of CMB
data. All these constraints are consistent with the GR prediction, i.e.,
$\theta=0$ even at 68% CL. For comparison, assuming $\Lambda$CDM cosmology,
from RSD data, we find $\Omega_{m}=0.293^{+0.10}_{-0.083}$ and
$S_{8}=0.758^{+0.082}_{-0.073}$. From the joint analysis with RSD + BAO +
Pantheon data, we find $\Omega_{m}=0.286^{+0.015}_{-0.015}$ and
$S_{8}=0.765^{+0.055}_{-0.054}$. We can clearly see that the addition of BAO +
Pantheon improves the constraint by breaking down the degeneracy on
$\Omega_{m}$, in both models, i.e., $\Lambda$CDM and MTMG. As we will see in
the next section, in MTMG model, there is no tension in the $S_{8}-\Omega_{m}$
plane, when analyzed from CMB and RSD data. Therefore, MTMG brings an
agreement between these data.

###  IV.2 MTMG in the light of CMB data

For the first time, here we find the constraints on the MTMG model from the
full Planck-CMB data set alone and its combination with several other data.
The baseline seven free parameters set of the MTMG model is given by:

$\mathcal{P}=\Big{\\{}\omega_{\rm b},\,\omega_{\rm
cdm},\,\theta_{s},\,\ln(10^{10}A_{s}),\,n_{s},\,\tau_{\rm
reio},\,\theta\Big{\\}},$

where the first six parameters are the baseline parameters of the standard
$\Lambda$CDM model, namely: $\omega_{\rm b}$ and $\omega_{\rm cdm}$ are
respectively the dimensionless densities of baryons and cold dark matter;
$\theta_{s}$ is the ratio of the sound horizon to the angular diameter
distance at decoupling; $A_{s}$ and $n_{s}$ are respectively the amplitude and
spectral index of the primordial curvature perturbations, and $\tau_{\rm
reio}$ is the optical depth to reionization. As commented before, the
parameter $\theta$ quantifies deviations from GR induced by the MTMG
framework.

Table 2 displays summary of the our statistical analyses using the full
Planck-CMB data set and its combination with several other data sets. In the
last column, the Full Joint Analysis stands for CMB + BAO + Pantheon + RSD +
ISW. Figure 4, on the left panel (right panel), shows the parameter space in
the $S_{8}-H_{0}$ ($\theta-\Omega_{m}$) plane, at 68% and 95% CL from Planck-
CMB data set and its combination with other data.

In all the analyses carried out here, we note that $\theta$ can be non-null at
68% CL, but consistent with $\theta=0$ at 95% CL. Thus, we do not find any
significant evidence for deviations from GR. In the Full Joint Analysis, we
obtain $\theta=0.25^{+0.16}_{-0.10}$ at 68% CL. The constraints on MTMG using
the CMB data alone are well consistent with the $\Lambda$CDM baseline. Despite
$\theta$ being non-null at 68% CL, all other parameters do not show
significant deviation from $\Lambda$CDM baseline. When BAO + Pantheon and BAO
+ Pantheon + RSD are added, we notice a minor shift of $\Omega_{m}$ to low
values and $H_{0}$ to high values. This behavior is clear in Fig. 4. We find
that RSD and CMB data are not in tension with each other in MTMG, see Fig. 6.
On the other hand, these data are known to be in tension ($\sim 3\sigma$)
within $\Lambda$CDM [32]. Thus, we can combine CMB and RSD in MTMG scenario.
In other words, there is no tension in the $S_{8}-\Omega_{m}$ plane for the
MTMG model.

Figure 5: Left panel: Scatterplot quantifying the correlations for some
parameters of interest in MTMG model from the RSD + BAO + Pantheon + ISW data.
Right panel: Same as the left panel, but for CMB + RSD + BAO + Pantheon + ISW
data.

We know that MTMG affects the CMB spectrum predominantly at low $l$, and
practically no effects are observed at high $l$. Thus, we also consider ISW
data in combination with the data-sets described above. We find
$H_{0}=68.44\pm 0.43$ km s-1Mpc-1 from CMB + BAO + Pantheon + RSD + ISW data.
This constraint minimizes the $H_{0}$ tension up to 3.5$\sigma$ compared to
cosmological model-independent local measurement of $H_{0}$ from the HST
observations [2]. From our joint analysis, we have $\Omega_{\rm m}=0.301\pm
0.0056$ and $\sigma_{8}=0.8169\pm 0.0068$, which represents a deviation to low
(high) values on $\Omega_{\rm m}$ ($\sigma_{8}$), respectively, compared to
the CMB only constraints.

Figure 5 shows a scatterplot for the correlations between the parameters of
interest in our model from the RSD + BAO + Pantheon + ISW and CMB + RSD + BAO
+ Pantheon + ISW data. We note that in both cases, $\theta$ exhibits negative,
positive and negative correlations with $\Omega_{m}$, $H_{0}$ and $S_{8}$,
respectively. We found similar correlations in other analyses as well. Despite
all the constraints on $\theta$, from the different analyses performed here be
statistically consistent with each other, it is interesting to note that
$\theta$ shows a best fit preference to $\theta<0$ from RSD + BAO + Pantheon +
ISW, and $\theta>0$ when we add the CMB data. We know that $\theta<0$ causes a
suppression on the amplitude of the matter perturbations. This suppression
behavior is necessary to fit well, for instance, with the RSD data and its
combination with BAO + Pantheon sample. On the other hand, in light of these
data, the parameter $\theta$ is degenerate with the other baseline parameters.
Thus, it is expected that in the presence of RSD data with BAO + Pantheon +
ISW data, preference exists for $\theta<0$ values while not excluding the
possibility of $\theta>0$ values. Next, the addition of CMB data is expected
to yield a tight constraint on $\theta$. We find in this case that the
correlation of $\theta$ with $\Omega_{m}$, $H_{0}$ and $S_{8}$, increases. It
is quantified in Fig. 5 comparing both panels. That increase especially occurs
on $\Omega_{m}$ and $H_{0}$, where the correlation with $H_{0}$ increases
400%. Note that $H_{0}$ is strongly constrained using CMB data, in MTMG model,
$H_{0}$ is constrained to 0.6% accuracy. The global constraint on $\theta$ is
also very tight, and the best fit on $\theta$ changes the signal due to the
correlation change in the analysis RSD + BAO + Pantheon + ISW to CMB + RSD +
BAO + Pantheon + ISW. That is caused by the addition of the CMB data. We also
note from the joint analysis CMB + RSD + BAO + Pantheon + ISW that $\theta>0$
at 68% CL. Changes in the other baseline parameters are understandable looking
in Figure 5.

Figure 6: Left panel: 68% and 95% CLs on the parametric space
$S_{8}-\Omega_{m}$ from CMB and RSD data sets under the $\Lambda$CDM and MTMG
model assumptions. Right panel: 68% and 95% CLs on the parametric space
$\mu^{2}-\Omega_{m}$ from several different data sets.

It is interesting to note that CMB data can break the degeneracy present in
RSD and BAO + Pantheon + RSD data, providing a tight constraint on $\theta$
when all these data are combined. It is clear from a comparison between the
tables 1 and 2. Using the relation $\mu^{2}=H_{0}^{2}\theta$, and constraints
on $H_{0}$, we can infer direct constraints on the graviton mass. Figure 6 on
the right panel shows the constraint on the graviton mass squared from CMB
alone, CMB + BAO + Pantheon, CMB + BAO + Pantheon + RSD and CMB + BAO +
Pantheon + RSD + ISW. If we take an inference with prior $\mu^{2}\geq 0$, we
find the bound: $\mu\lessapprox 8.8\times 10^{-67}$ ${\rm eV}^{2}$,
$\lessapprox 9.2\times 10^{-67}$ ${\rm eV}^{2}$ eV, $\lessapprox 9.9\times
10^{-67}$ ${\rm eV}^{2}$ eV, $\lessapprox 7.2\times 10^{-67}$ ${\rm eV}^{2}$
eV at 95% CL from CMB alone, CMB + BAO + Pantheon, CMB + BAO + Pantheon + RSD
and CMB + BAO + Pantheon + RSD + ISW, respectively. Without loss of
generality, this inference on $\mu$ can be performed for the other data sets.
We chose to perform only using the most accurate measurements.

Figure 7: Left panel: 68% and 95% CLs on the parametric space $S_{8}-\mu^{2}$
from our full joint analysis combined with a Gaussian prior on $S_{8}$ from
DES results [65]. Right panel: Same as in left panel, but the $S_{8}$ prior
from KiDS-1000 [66].

Our constraints are consistent with the bound $\mu<10^{-23}$ eV on the
graviton mass set by the LIGO collaboration [67]. Bounds from the solar system
tests provide $\mu<10^{-24}$ eV [68]. Analyzing galactic dynamics under n-body
simulations, it was shown that $\mu<10^{-26}$ eV [69]. Constraints from
observation-derived energy condition bounds show that $\mu<10^{-31}$ eV [70].
Recently, using 14 well-timed binary pulsars, from their intrinsic orbital
decay rates, the authors in [71] found $\mu<10^{-28}$ eV. Primordial
Gravitational Waves modeled with a massive graviton will induce extra effects
on the B-mode polarization of the CMB at low $l$ scales, which can place a
bound on the massive graviton to $\mu<10^{-30}$ eV [72]. A more strong bound
can be imposed from weak gravitational lensing observations [73] and lunar
laser ranging experiments [74], viz., $\mu<10^{-32}$ eV. In Table I in [75],
the authors provide a compilation of several bounds on the graviton mass.
Although found for a given (but realistic and stable) model describing a
massive graviton, our constraint on the graviton mass is complementary and
some orders of magnitude stronger than the previous ones. We refer to [47] for
other bounds on $\mu$ in MTMG.

##  V Final Remarks

Despite of the $\Lambda$CDM successes, some challenges at both the theoretical
and the observational levels have placed the $\Lambda$CDM cosmology in
crossroads. At the observational level, it faces in particular the $H_{0}$ and
growth tensions. The growth tension is based on the fact that the observed
growth of cosmological perturbations at low $z$ is weaker than the growth
predicted by the standard Planck-$\Lambda$CDM parameter values. In this work,
we have investigated in detail the MTMG model in the light of different
observational data sets. In particular, we have derived constraints on the
MTMG using Planck-CMB data for the first time. From our full joint analysis,
we have found $\mu^{2}/H_{0}^{2}=0.25^{{+}0.16}_{{-}0.10}$, $\mu$ being the
mass of the graviton. It represents a non-null measurement on $\mu$ at 68% CL.

From a theoretical side, in obtaining our results, we have considered the
scalar and tensor perturbative effects under a FLRW metric, while the
background evolution is equivalent to $\Lambda$CDM model as described in
section II. We know that the $H_{0}$ value from CMB data depends on the
angular scale $\theta_{*}=d_{s}^{*}/D_{A}^{*}$, where $d_{s}^{*}$ is the sound
horizon at decoupling (the distance traveled by a sound wave from the big bang
to the epoch of the CMB-baryons decoupling) and $D_{A}^{*}$ is the angular
diameter distance at decoupling, which in turn depends on the expansion
history, $H(z)$, after decoupling, controlled also by the ratio
$\Omega_{m}/\Omega_{de}$ and $H_{0}$ mainly. In its simplest form, the MTMG
scenario does not lead to changes in $H(z)$ evolution with regard to
$\Lambda$CDM model. Thus, significant changes on $H_{0}$ that are obtained in
other models in the literature, are not observed here. All changes on $H_{0}$
as well as on the full baseline are due to the changes on the evolution of the
scalar potentials $\phi$ and $\psi$ in the Einstein’s field equations. Within
this framework, we have found that the current $H_{0}$-tension at 4.4$\sigma$
can be reduced to 3.5$\sigma$. It does not completely solve this tension, but
leads us to conclude that the possibility for a non-null $\mu$ can affect the
cosmological parameters estimation to CMB level too. Going beyond the simplest
implementation of the MTMG model, namely changing also $H(z)$, may address
this point more significantly. It may be interesting, in a future
communication, to verify massive theories with change on the background
evolution, which certainly should reduce the $H_{0}$ tension significantly. We
also conclude that the well-known tension on the $S_{8}$ parameter can be
solved in MTMG model, as shown in the previous section. During the final stage
of preparation of this work, the weak lensing and galaxy clustering
measurements from the Dark Energy Survey (DES) have been updated, providing
new and robust estimates on $S_{8}$, viz., $S_{8}=0.776^{+0.017}_{-0.017}$ and
$S_{8}=0.812^{+0.012}_{-0.012}$ (DES + CMB) [76, 65], within a $\Lambda$CDM
baseline. It is interesting to interpret that the MTMG scenario, without using
the CMB data, predicts $S_{8}=0.802\pm 0.06$, while enlarging the error bar a
bit, in general, on the $S_{8}-\Omega_{m}$ plane compared to $\Lambda$CDM.
Thus, comparing with DES constraint within $\Lambda$CDM baseline results, we
can postulate that MTMG can also solve the $S_{8}$ tension between DES and
CMB, which is around 2.3$\sigma$.

We simulate an addition of the DES and Kilo-Degree Survey (KiDS-1000) data on
our Full Joint Analysis by adding Gaussian priors on $S_{8}$ to our analysis.
Figure 7, on the left panel, shows the parameter space in the $S_{8}-\mu^{2}$
plane under the prior $S_{8}=0.776\pm 0.017$ from the DES results [65]. On the
right panel, we show the results assuming the prior
$S_{8}=0.759^{+0.024}_{-0.021}$ from the KiDS-1000 [66]. In both cases, we
note that the amplitude in $S_{8}$ decreases, and slightly extends the
expectation of the graviton mass, but consistent with all other results
developed here. Of course, an analysis within MTMG baseline using the full DES
and KiDS likelihoods must confirm this prediction. We hope to provide such
investigation in a future communication.

Therefore, in this work, we conclude that the MTMG is very well consistent
with the CMB observations as well as to other observational ones that are
sensitive to the cosmological perturbation theory. Undoubtedly, the MTMG model
is a viable candidate among the modified gravity theories.

###### Acknowledgements.

J.C.N.A. thanks FAPESP (2013/26258-4) and CNPq (308367/2019-7) for partial
financial support. The work of A.D.F. was supported by Japan Society for the
Promotion of Science Grants-in-Aid for Scientific Research No. 20K03969. S.K.
gratefully acknowledges the support from SERB-DST project No. EMR/2016/000258.
R.C.N. would like to thank the agency FAPESP for financial support under the
project No. 2018/18036-5. Part of the numerical computation in this work was
carried out at the Yukawa Institute Computer Facility.

## References

  * Aghanim  _et al._ [2020a] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, and et al., Planck 2018 results, Astronomy and Astrophysics 641, A6 (2020a).
  * Riess  _et al._ [2019] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, Large magellanic cloud cepheid standards provide a 1% foundation for the determination of the hubble constant and stronger evidence for physics beyond lcdm, The Astrophysical Journal 876, 85 (2019).
  * Wong  _et al._ [2019] K. C. Wong, S. H. Suyu, G. C.-F. Chen, Rusu, and et al., H0licow – xiii. a 2.4 per cent measurement of h0 from lensed quasars: 5.3$\sigma$ tension between early- and late-universe probes, Monthly Notices of the Royal Astronomical Society 498, 1420–1439 (2019).
  * Bonilla  _et al._ [2021] A. Bonilla, S. Kumar, and R. C. Nunes, Measurements of $h_{0}$ and reconstruction of the dark energy properties from a model-independent joint analysis, The European Physical Journal C 81, 10.1140/epjc/s10052-021-08925-z (2021).
  * Di Valentino and Anchordoqui [2021] E. Di Valentino and L. A. e. a. Anchordoqui, Snowmass2021 - letter of interest cosmology intertwined ii: The hubble constant tension, Astroparticle Physics 131, 102605 (2021).
  * Valentino  _et al._ [2021] E. D. 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 (2021), arXiv:2103.01183 [astro-ph.CO] .
  * Perivolaropoulos and Skara [2021] L. Perivolaropoulos and F. Skara, Challenges for $\lambda$cdm: An update (2021), arXiv:2105.05208 [astro-ph.CO] .
  * Aiola  _et al._ [2020] S. Aiola, E. Calabrese, L. Maurin, S. Naess, B. L. Schmitt, M. H. Abitbol, G. E. Addison, P. A. R. Ade, D. Alonso, M. Amiri, and et al., The atacama cosmology telescope: Dr4 maps and cosmological parameters, Journal of Cosmology and Astroparticle Physics 2020 (12), 047–047.
  * Joudaki  _et al._ [2017a] S. Joudaki, A. Mead, C. Blake, A. Choi, J. de Jong, T. Erben, I. Fenech Conti, R. Herbonnet, C. Heymans, H. Hildebrandt, and et al., Kids-450: testing extensions to the standard cosmological model, Monthly Notices of the Royal Astronomical Society 471, 1259–1279 (2017a).
  * Joudaki  _et al._ [2017b] S. Joudaki, C. Blake, A. Johnson, and et al., Kids-450 + 2dflens: Cosmological parameter constraints from weak gravitational lensing tomography and overlapping redshift-space galaxy clustering, Monthly Notices of the Royal Astronomical Society 474, 4894–4924 (2017b).
  * Troxel  _et al._ [2018] M. Troxel, N. MacCrann, J. Zuntz, and et al., Dark energy survey year 1 results: Cosmological constraints from cosmic shear, Physical Review D 98, 10.1103/physrevd.98.043528 (2018).
  * Heymans  _et al._ [2020] C. Heymans, T. Tröster, M. Asgari, and et al., Kids-1000 cosmology: Multi-probe weak gravitational lensing and spectroscopic galaxy clustering constraints (2020), arXiv:2007.15632 [astro-ph.CO] .
  * Asgari  _et al._ [2020] M. Asgari, T. Tröster, C. Heymans, and et al., Kids+viking-450 and des-y1 combined: Mitigating baryon feedback uncertainty with cosebis, Astronomy and Astrophysics 634, A127 (2020).
  * Tröster  _et al._ [2020] T. Tröster, A. G. Sánchez, and et al., Cosmology from large-scale structure, Astronomy and Astrophysics 633, L10 (2020).
  * Ivanov  _et al._ [2020] M. M. Ivanov, M. Simonović, and M. Zaldarriaga, Cosmological parameters from the boss galaxy power spectrum, Journal of Cosmology and Astroparticle Physics 2020 (05), 042–042.
  * Di Valentino  _et al._ [2021] E. Di Valentino, L. A. Anchordoqui, and et al., Cosmology intertwined iii: fs8 and s8, Astroparticle Physics 131, 102604 (2021).
  * Kumar  _et al._ [2019a] S. Kumar, R. C. Nunes, and S. K. Yadav, Dark sector interaction: a remedy of the tensions between cmb and lss data, The European Physical Journal C 79, 10.1140/epjc/s10052-019-7087-7 (2019a).
  * Kumar  _et al._ [2019b] S. Kumar, R. C. Nunes, and S. K. Yadav, Testing the warmness of dark matter, Monthly Notices of the Royal Astronomical Society 490, 1406–1414 (2019b).
  * Lambiase  _et al._ [2019] G. Lambiase, S. Mohanty, A. Narang, and P. Parashari, Testing dark energy models in the light of $\sigma_{8}$ tension, The European Physical Journal C 79, 10.1140/epjc/s10052-019-6634-6 (2019).
  * Barros  _et al._ [2020] B. J. Barros, T. Barreiro, T. Koivisto, and N. J. Nunes, Testing $f(q)$ gravity with redshift space distortions, Physics of the Dark Universe 30, 100616 (2020).
  * Abellan  _et al._ [2020] G. F. Abellan, R. Murgia, V. Poulin, and J. Lavalle, Hints for decaying dark matter from $s_{8}$ measurements (2020), arXiv:2008.09615 [astro-ph.CO] .
  * Heimersheim  _et al._ [2020] S. Heimersheim, N. Schöneberg, D. C. Hooper, and J. Lesgourgues, Cannibalism hinders growth: Cannibal dark matter and the s8 tension, Journal of Cosmology and Astroparticle Physics 2020 (12), 016–016.
  * Choi  _et al._ [2021] G. j. Choi, T. T. Yanagida, and N. Yokozaki, A model of interacting dark matter and dark radiation for h0 and $\sigma_{8}$ tensions, Journal of High Energy Physics 2021, 10.1007/jhep01(2021)127 (2021).
  * Skara and Perivolaropoulos [2020] F. Skara and L. Perivolaropoulos, Tension of the eg statistic and redshift space distortion data with the planck-lcdm model and implications for weakening gravity, Physical Review D 101, 10.1103/physrevd.101.063521 (2020).
  * Marra and Perivolaropoulos [2021] V. Marra and L. Perivolaropoulos, A rapid transition of $g_{\rm eff}$ at $z_{t}\simeq 0.01$ as a solution of the hubble and growth tensions (2021), arXiv:2102.06012 [astro-ph.CO] .
  * Kumar [2021] S. Kumar, Remedy of some cosmological tensions via effective phantom-like behavior of interacting vacuum energy, Physics of the Dark Universe 33, 100862 (2021), arXiv:2102.12902 [astro-ph.CO] .
  * Lucca [2021] M. Lucca, Dark energy-dark matter interactions as a solution to the $s_{8}$ tension (2021), arXiv:2105.09249 [astro-ph.CO] .
  * Yang  _et al._ [2021] W. Yang, E. Di Valentino, S. Pan, and O. Mena, Emergent dark energy, neutrinos and cosmological tensions, Physics of the Dark Universe 31, 100762 (2021).
  * Yang  _et al._ [2020] W. Yang, E. Di Valentino, O. Mena, and S. Pan, Dynamical dark sectors and neutrino masses and abundances, Physical Review D 102, 10.1103/physrevd.102.023535 (2020).
  * Choi and Yanagida [2021] G. Choi and T. T. Yanagida, Projecting cosmological problems on the story of gravitinos (2021), arXiv:2104.02958 [hep-ph] .
  * Kazantzidis and Perivolaropoulos [2018] L. Kazantzidis and L. Perivolaropoulos, Evolution of the f$\sigma_{8}$ tension with the planck15/lcdm determination and implications for modified gravity theories, Physical Review D 97, 10.1103/physrevd.97.103503 (2018).
  * Nunes and Vagnozzi [2021] R. C. Nunes and S. Vagnozzi, Arbitrating the $s_{8}$ discrepancy with growth rate measurements from redshift-space distortions (2021), arXiv:2106.01208 [astro-ph.CO] .
  * Ishak [2018] M. Ishak, Testing general relativity in cosmology, Living Reviews in Relativity 22, 10.1007/s41114-018-0017-4 (2018).
  * De Felice  _et al._ [2021] A. De Felice, S. Mukohyama, and M. C. Pookkillath, Addressing $H_{0}$ tension by means of VCDM, Phys. Lett. B 816, 136201 (2021), arXiv:2009.08718 [astro-ph.CO] .
  * Ballardini  _et al._ [2020] M. Ballardini, M. Braglia, F. Finelli, D. Paoletti, A. A. Starobinsky, and C. Umiltà, Scalar-tensor theories of gravity, neutrino physics, and the $h_{0}$ tension, Journal of Cosmology and Astroparticle Physics 2020 (10), 044–044.
  * Nunes [2018] R. C. Nunes, Structure formation in $f(t)$ gravity and a solution for $h_{0}$ tension, Journal of Cosmology and Astroparticle Physics 2018 (05), 052–052.
  * D’Agostino and Nunes [2020] R. D’Agostino and R. C. Nunes, Measurements of $h_{0}$ in modified gravity theories: The role of lensed quasars in the late-time universe, Physical Review D 101, 10.1103/physrevd.101.103505 (2020).
  * Felice  _et al._ [2020] A. D. Felice, C.-Q. Geng, M. C. Pookkillath, and L. Yin, Reducing the $h_{0}$ tension with generalized proca theory, Journal of Cosmology and Astroparticle Physics 2020 (08), 038–038.
  * Alestas  _et al._ [2020] G. Alestas, L. Kazantzidis, and L. Perivolaropoulos, $h0$ tension, phantom dark energy, and cosmological parameter degeneracies, Physical Review D 101, 10.1103/physrevd.101.123516 (2020).
  * de Rham [2014] C. de Rham, Massive gravity, Living Reviews in Relativity 17, 10.12942/lrr-2014-7 (2014).
  * De Felice and Mukohyama [2016] A. De Felice and S. Mukohyama, Minimal theory of massive gravity, Phys. Lett. B 752, 302 (2016), arXiv:1506.01594 [hep-th] .
  * Felice and Mukohyama [2016] A. D. Felice and S. Mukohyama, Phenomenology in minimal theory of massive gravity, Journal of Cosmology and Astroparticle Physics 2016 (04), 028–028.
  * de Rham and Gabadadze [2010] C. de Rham and G. Gabadadze, Generalization of the fierz-pauli action, Physical Review D 82, 10.1103/physrevd.82.044020 (2010).
  * de Rham  _et al._ [2011] C. de Rham, G. Gabadadze, and A. J. Tolley, Resummation of massive gravity, Physical Review Letters 106, 10.1103/physrevlett.106.231101 (2011).
  * Aoki  _et al._ [2020] K. Aoki, A. De Felice, S. Mukohyama, K. Noui, M. Oliosi, and M. C. Pookkillath, Minimally modified gravity fitting Planck data better than $\Lambda$CDM, Eur. Phys. J. C 80, 708 (2020), arXiv:2005.13972 [astro-ph.CO] .
  * Pookkillath  _et al._ [2019] M. C. Pookkillath, A. De Felice, and S. Mukohyama, Baryon Physics and Tight Coupling Approximation in Boltzmann Codes, Universe 6, 6 (2019), arXiv:1906.06831 [astro-ph.CO] .
  * De Felice and Mukohyama [2017] A. De Felice and S. Mukohyama, Graviton mass might reduce tension between early and late time cosmological data, Physical Review Letters 118, 10.1103/physrevlett.118.091104 (2017).
  * Bolis  _et al._ [2018a] N. Bolis, A. De Felice, and S. Mukohyama, Integrated Sachs-Wolfe-galaxy cross-correlation bounds on the two branches of the minimal theory of massive gravity, Phys. Rev. D 98, 024010 (2018a), arXiv:1804.01790 [astro-ph.CO] .
  * Blas  _et al._ [2011] D. Blas, J. Lesgourgues, and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes, JCAP 07, 034, arXiv:1104.2933 [astro-ph.CO] .
  * Mead  _et al._ [2015] A. J. Mead, J. A. Peacock, C. Heymans, S. Joudaki, and A. F. Heavens, An accurate halo model for fitting non-linear cosmological power spectra and baryonic feedback models, Monthly Notices of the Royal Astronomical Society 454, 1958–1975 (2015).
  * Hagala  _et al._ [2021] R. Hagala, A. De Felice, D. F. Mota, and S. Mukohyama, Non-linear dynamics of the minimal theory of massive gravity, Astronomy and Astrophysics 10.1051/0004-6361/202040018 (2021).
  * Aghanim  _et al._ [2020b] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, and et al., Planck2018 results, Astronomy and Astrophysics 641, A8 (2020b).
  * Sagredo  _et al._ [2018] B. Sagredo, S. Nesseris, and D. Sapone, Internal Robustness of Growth Rate data, Phys. Rev. D 98, 083543 (2018), arXiv:1806.10822 [astro-ph.CO] .
  * Arjona  _et al._ [2020] R. Arjona, J. García-Bellido, and S. Nesseris, Cosmological constraints on nonadiabatic dark energy perturbations, Physical Review D 102, 10.1103/physrevd.102.103526 (2020).
  * Alam  _et al._ [2021] S. Alam, M. Aubert, S. Avila, C. Balland, J. E. Bautista, M. A. Bershady, D. Bizyaev, M. R. Blanton, A. S. Bolton, J. Bovy, and et al., Completed sdss-iv extended baryon oscillation spectroscopic survey: Cosmological implications from two decades of spectroscopic surveys at the apache point observatory, Physical Review D 103, 10.1103/physrevd.103.083533 (2021).
  * Scolnic  _et al._ [2018] D. M. Scolnic _et al._ , The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample, Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Kinkhabwala and Kamionkowski [1999] A. Kinkhabwala and M. Kamionkowski, New constraint on open cold-dark-matter models, Physical Review Letters 82, 4172–4175 (1999).
  * Giacomello  _et al._ [2019] F. Giacomello, A. D. Felice, and S. Ansoldi, Bounds from isw-galaxy cross-correlations on generalized covariant galileon models, Journal of Cosmology and Astroparticle Physics 2019 (03), 038–038.
  * Bolis  _et al._ [2018b] N. Bolis, A. De Felice, and S. Mukohyama, Integrated sachs-wolfe-galaxy cross-correlation bounds on the two branches of the minimal theory of massive gravity, Physical Review D 98, 10.1103/physrevd.98.024010 (2018b).
  * Song  _et al._ [2007] Y.-S. Song, W. Hu, and I. Sawicki, Large scale structure of $f(r)$ gravity, Physical Review D 75, 10.1103/physrevd.75.044004 (2007).
  * Stölzner  _et al._ [2018] B. Stölzner, A. Cuoco, J. Lesgourgues, and M. Bilicki, Updated tomographic analysis of the integrated sachs-wolfe effect and implications for dark energy, Physical Review D 97, 10.1103/physrevd.97.063506 (2018).
  * Audren  _et al._ [2013] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, Conservative Constraints on Early Cosmology: an illustration of the Monte Python cosmological parameter inference code, JCAP 02, 001, arXiv:1210.7183 [astro-ph.CO] .
  * Brinckmann and Lesgourgues [2019] T. Brinckmann and J. Lesgourgues, MontePython 3: boosted MCMC sampler and other features, Phys. Dark Univ. 24, 100260 (2019), arXiv:1804.07261 [astro-ph.CO] .
  * Mossa  _et al._ [2020] V. Mossa _et al._ , The baryon density of the Universe from an improved rate of deuterium burning, Nature 587, 210 (2020).
  * Collaboration [2021a] D. Collaboration, Dark energy survey year 3 results: Cosmological constraints from galaxy clustering and weak lensing (2021a), arXiv:2105.13549 [astro-ph.CO] .
  * Asgari  _et al._ [2021] M. Asgari, C.-A. Lin, B. Joachimi, B. Giblin, C. Heymans, H. Hildebrandt, A. Kannawadi, B. Stölzner, T. Tröster, J. L. van den Busch, and et al., Kids-1000 cosmology: Cosmic shear constraints and comparison between two point statistics, Astronomy and Astrophysics 645, A104 (2021).
  * Collaboration and the Virgo Collaboration [2020] T. L. S. Collaboration and the Virgo Collaboration, Tests of general relativity with binary black holes from the second ligo-virgo gravitational-wave transient catalog (2020), arXiv:2010.14529 [gr-qc] .
  * Will [2018] C. M. Will, Solar system versus gravitational-wave bounds on the graviton mass, Classical and Quantum Gravity 35, 17LT01 (2018).
  * Brandao and de Araujo [2012] C. S. S. Brandao and J. C. N. de Araujo, A recipe to probe alternative theories of gravitation vian-body numerical simulations. i. spiral galaxies, The Astrophysical Journal 750, 29 (2012).
  * Alves  _et al._ [2018] M. E. S. Alves, F. C. Carvalho, J. C. N. de Araujo, M. Penna-Lima, and S. D. P. Vitenti, Cosmological constant constraints from observation-derived energy condition bounds and their application to bimetric massive gravity, The European Physical Journal C 78, 10.1140/epjc/s10052-018-6190-5 (2018).
  * Shao  _et al._ [2020] L. Shao, N. Wex, and S.-Y. Zhou, New graviton mass bound from binary pulsars, Physical Review D 102, 10.1103/physrevd.102.024069 (2020).
  * Dubovsky  _et al._ [2010] S. Dubovsky, R. Flauger, A. Starobinsky, and I. Tkachev, Signatures of a graviton mass in the cosmic microwave background, Physical Review D 81, 10.1103/physrevd.81.023523 (2010).
  * CHOUDHURY [2004] S. CHOUDHURY, Probing large distance higher dimensional gravity from lensing data, Astroparticle Physics 21, 559–563 (2004).
  * Dvali  _et al._ [2003] G. Dvali, A. Gruzinov, and M. Zaldarriaga, The accelerated universe and the moon, Physical Review D 68, 10.1103/physrevd.68.024012 (2003).
  * de Rham  _et al._ [2017] C. de Rham, J. T. Deskins, A. J. Tolley, and S.-Y. Zhou, Graviton mass bounds, Reviews of Modern Physics 89, 10.1103/revmodphys.89.025004 (2017).
  * Collaboration [2021b] D. Collaboration, Dark energy survey year 3 results: Cosmology from cosmic shear and robustness to data calibration (2021b), arXiv:2105.13543 [astro-ph.CO] .