\section{Introduction}

Inflation~\cite{inflation,Yokoyama:2014nga} plays a crucial role in cosmology of the very early Universe.
In particular, single-field slow-roll models of inflation generically
produce nearly scale-invariant, adiabatic, and Gaussian curvature perturbations
as the seeds for cosmic structure~\cite{fluctuation}.
The theoretical prediction
matches observational results e.g. of the Planck experiments fairly well~\cite{Ade:2013ktc,Ade:2013uln,Adam:2015rua,PlanckCollaboration2015b}.
During inflation primordial tensor modes (gravitational waves) are generated as well.
The tensor amplitude is conventionally parametrized by the tensor-to-scalar ratio, $r$,
and the Planck constraint on $r$ is given by $r<0.10$ (95\%
C.L.)~\cite{PlanckCollaboration2015b}.
Some of the single-field slow-roll models predict larger tensor modes,
and hence have been excluded by current observations.
One would then ask whether one can reduce the tensor amplitude somehow to save such models.
This is the question which we discuss in this paper.

General relativity is an underlying assumption of standard inflation models,
and nonstandard dynamics of the tensor modes can be obtained
by modifying this gravitational sector.
In doing so, one generically expects that
the behavior of the scalar perturbations is also modified.
This is however what we want to avoid, because
the standard inflationary predictions on the scalar perturbations are so successful.
In this paper,
we therefore explore the possibility of
modifying {\em only} the tensor modes and try to retain the same structure
of the scalar sector as in general relativity,
in order not to spoil the remarkable agreement between
the standard theoretical predictions of the scalar perturbations and observations.

It is natural to consider quadratic curvature terms in the action beyond general relativity
since such corrections are expected to arise as signatures of new physics at high energies.
Below we look for quadratic curvature terms that modify only the tensor sector
of cosmological perturbations without introducing any pathologies such as ghost instabilities.
It turns out that there are two independent combinations of the curvature tensors
fulfilling the above requirements.
Both combinations do not change the quadratic action for the scalar perturbations,
and one of them has no impact on the cubic action as well.
The resultant quadratic curvature terms are not of the form
${\cal R}^2$, ${\cal R}_{\mu\nu}{\cal R}^{\mu\nu}$, etc. which are familiar in the literature~\cite{Stelle:1976gc,Stelle1978},
but they have nontrivial coupling to the derivative of the inflaton field.
We study in detail how the tensor amplitude and tilt are modified,
and discuss the implications for observations.

The organization of this paper is as follows.
We determine the two possible quadratic curvature terms
satisfying our requirement in the next section.
In Sec.~III we evaluate the modified amplitude and tilt of the
primordial tensor modes, and then present the implications for observations.
We draw our conclusions in Sec.~IV.

\section{Construction of the Lagrangian}

The theory we consider is described by the action
\begin{eqnarray}
S=S_{\rm EH}+S_\phi +S_{\rm higher},
\end{eqnarray}
where $S_{\rm EH}$ is the Einstein-Hilbert term,
\begin{eqnarray}
S_{\rm EH}=\frac{1}{2\kappa}\int{\rm d}^4x\sqrt{-g}{\cal R},
\end{eqnarray}
$S_\phi$ is the action of the inflaton field,
\begin{eqnarray}
S_\phi =\int{\rm d}^4x\sqrt{-g}P(\phi, \partial^\mu \phi\partial_\mu\phi),
\end{eqnarray}
and $S_{\rm higher}$ represents higher curvature corrections,
\begin{eqnarray}
S_{\rm higher}=\frac{1}{\kappa} \int{\rm d}^4x\sqrt{-g}\left(
\frac{1}{M^2}{\cal R}_{\mu\nu\rho\sigma}{\cal R}^{\mu\nu\rho\sigma}+\cdots \right).\label{high-ex}
\end{eqnarray}
The simplest Lagrangian for the inflaton field would be of the form
$P=-\partial^\mu\phi\partial_\mu\phi/2-V(\phi)$,
but here we do not need to specify the concrete form of $P$.

It is known that
typical higher curvature terms like the one presented in Eq.~(\ref{high-ex})
give rise to new propagating degrees of freedom
which are plagued by (ghost) instabilities~\cite{Ostrogradski}.
In this paper, we carefully construct higher curvature terms
so that the resultant theory is free from such dangerous degrees of freedom.
Among such healthy theories we are interested in those in which
the dynamics of tensor perturbations is modified
while the scalar sector of cosmological perturbations is left unchanged.
To find the higher curvature terms having those properties,
we have to go beyond the familiar curvature invariants such as
${\cal R}_{\mu\nu\rho\sigma}{\cal R}^{\mu\nu\rho\sigma}$ and ${\cal R}_{\mu\nu}{\cal R}^{\mu\nu}$,
and consider the terms obtained by contracting with
the unit normal to constant $\phi$ hypersurfaces,
\begin{eqnarray}
u_\mu:=-\frac{\partial_\mu\phi}{\sqrt{-\partial^\nu\phi\partial_\nu\phi}},
\end{eqnarray}
and the induced metric,
\begin{eqnarray}
\gamma_{\mu\nu} = g_{\mu\nu}+u_\mu u_\nu,
\end{eqnarray}
e.g.,
${\cal R}_{\mu\nu\rho\sigma}{\cal R}_{\mu'\nu'\rho'\sigma'}\gamma^{\mu\mu'}\gamma^{\nu\nu'}\gamma^{\rho\rho'}u^\sigma u^{\sigma'}$.
This possibility was demonstrated
in the context of Weyl gravity in Ref.~\cite{Deruelle:2012xv}.

Focusing on quadratic curvature corrections,
we are going to identify the terms in the Lagrangian fulfilling the above requirements
in the following way. The basic idea here is along the same line as taken in Refs.~\cite{GLPV,Gao}.
We start by performing the Arnowitt-Deser-Misner (ADM) decomposition,
taking constant $\phi$ hypersurfaces as constant time hypersurfaces,
as the dynamics of cosmological perturbations is more transparent
in the ADM language.
The metric is written in terms of the ADM variables as
\begin{eqnarray}
{\rm d} s^2=-N^2{\rm d} t^2+\gamma_{ij}\left({\rm d} x^i+N^i{\rm d} t \right)\left({\rm d} x^j+N^j{\rm d} t \right).
\end{eqnarray}
The possible quadratic curvature terms in the Lagrangian are exhaustively
written in terms of the three-dimensional geometric quantities as
\begin{align}
\sqrt{\gamma}N &\;\times\; \left\{
K^4,\; K_{ij}K^{ij}K^2,\;\cdots,
\;
R^2,\;R_{ij}R^{ij},\right.
\notag \\
& \left.
K^2R,\; KK^{ij}R_{ij},\;\cdots,
\;
D_iK_{jk}D^iK^{jk},\;\cdots \right\} ,\label{termlist}
\end{align}
where $K_{ij}$ and $R_{ij}$ are
the extrinsic and intrinsic curvature tensors of
the constant $\phi$ hypersurfaces, respectively,
$D_i$ stands for the covariant derivative with respect to $\gamma_{ij}$,
and ellipses are used to indicate analogous terms whose indices are contracted in different ways.
We discard from the above candidates
the terms containing time derivatives of the extrinsic curvature,
because higher time derivatives of the metric imply the appearance of
additional propagating degrees of freedom other than $\phi$
and two tensor modes, signaling instabilities.\footnote{If
one has only $\dot K$ the theory is not necessarily unstable,
as is illustrated by the example of the ${\cal R}^2$ model.
This however adds an extra scalar degree of freedom
modifying the scalar sector of cosmological perturbations.
For this reason we avoid any time derivatives of $K_{ij}$.}

Let us consider cosmological perturbations,
\begin{align}
N=1+\delta N,\quad N_i=\partial_i\chi +\chi_i,\quad
\gamma_{ij}=a^2e^{2\zeta}\left(e^h \right)_{ij},
\end{align}
where $\zeta$ is the curvature perturbation on the uniform $\phi$ hypersurfaces,
$h_{ij}$ is the transverse and traceless tensor perturbation, and $\chi_i$
is the transverse vector perturbation.
Let us concentrate on the scalar sector for the moment.
To first order in perturbations,
the extrinsic curvature is given by
\begin{eqnarray}
K_i^{\;j} =H\delta_i^{\;j}+\frac{1}{3}\delta K\delta_i^{\;j} +
\delta \widetilde K_i^{\;j},
\end{eqnarray}
with
\begin{align}
\delta K =&\; -3H\delta N+3\dot\zeta-\frac{1}{a^{2}} \partial^2\chi,
\\
\delta\widetilde K_i^{\;j}=&\;-\frac{1}{a^{2}}
\left(\partial_i\partial^j-\frac{1}{3}\delta_i^{\;j}\partial^2\right)\chi,
\end{align}
and the intrinsic curvature is
\begin{eqnarray}
\delta R_i^{\;j}=-\frac{1}{a^2}\left(\partial_i\partial^j+\delta_i^{\;j}\partial^2 \right)\zeta.
\end{eqnarray}
The perturbation of the extrinsic curvature tensor has been decomposed into
its trace and traceless parts.
Using those quantities as building blocks,
one can construct
the following two combinations
of the form listed in Eq.~(\ref{termlist})
for which the scalar-type variables are canceled out
after integration by parts,
\begin{eqnarray}
2\partial_i\delta \widetilde K_{jk}\partial^i\delta\widetilde K^{jk}
-3\partial_i\delta\widetilde K^{ik}\partial^j\delta \widetilde K_{jk},
\label{quad-KK}
\end{eqnarray}
and
\begin{eqnarray}
\delta R_{ij}\delta R^{ij}-\frac{3}{8}\delta R^2,
\end{eqnarray}
at quadratic order in perturbations.
No other combinations can be found with vanishing scalar-type variables.
Now including vector and tensor perturbations we have
\begin{align}
&
2\partial_i\delta \widetilde K_{jk}\partial^i\delta\widetilde K^{jk}
-3\partial_i\delta\widetilde K^{ik}\partial^j\delta \widetilde K_{jk}
\notag \\
&\qquad\qquad\qquad
=\frac{1}{2a^2}\left(\partial_i\dot h_{jk}\right)^2+\frac{1}{4a^6}\left(\partial^2\chi_i \right)^2,
\\
&
\delta R_{ij}\delta R^{ij}-\frac{3}{8}\delta R^2 = \frac{1}{4a^4}\left(\partial^2h_{ij}\right)^2.
\label{RR-tensor}
\end{align}
Both of the two possible quadratic terms for $h_{ij}$ with four derivatives are obtained,
while we successfully exclude $\ddot h_{ij}^2$ which would cause Ostrogradski ghosts.
Since there is no kinetic term for $\chi_i$ here and in $S_{\rm EH}$,
the vector perturbation is not dynamical.
We therefore ignore the vector sector in this paper.

Having thus written the quadratic Lagrangian for perturbations
in terms of the geometric quantities, it is almost straightforward
to determine the full nonlinear Lagrangian in the ADM form as
\begin{align}
{\cal L}_1' =&\;
\frac{\sqrt{\gamma}N}{M^2}
\left(
2D_i \widetilde K_{jk}D^i\widetilde K^{jk}
-3D_i\widetilde K^{ik}D^j \widetilde K_{jk}
\right),
\\
{\cal L}_2 =&\;
\frac{\sqrt{\gamma}N }{M^2}\left(
R_{ij} R^{ij}-\frac{3}{8} R^2
\right),
\end{align}
where $\widetilde K_{ij}:=K_{ij}-(1/3)K\gamma_{ij}$
is the traceless part of the extrinsic curvature.
Note that ${\cal L}_1'$ is one of the candidates;
in fact, we have different choices
that reduce to Eq.~(\ref{quad-KK}) after integration by parts
at the level of the quadratic Lagrangian for perturbations.
Among them we adopt
\begin{align}
{\cal L}_1 =\frac{\sqrt{\gamma}N}{M^2}
\bigl(&
2D_i \widetilde K_{jk}D^i\widetilde K^{jk}
-D_i\widetilde K^{ik}D^j \widetilde K_{jk}
\notag\\ &
-2D_i\widetilde K_{jk}D^j \widetilde K^{ik}
\bigr)\label{L1true}
\end{align}
rather than ${\cal L}_1'$. What is particular to ${\cal L}_1$
is that it can be written as a square of some tensor as
\begin{eqnarray}
{\cal L}_1=\frac{\sqrt{\gamma}N}{M^2}W_{ijk}W^{ijk},
\end{eqnarray}
where
\begin{eqnarray}
W_{ijk}= 2D_{[i}\widetilde K_{j]k}+D_l\widetilde K^l_{\;[i}\gamma_{j]k}.
\end{eqnarray}
It is clear that the scalar perturbations do not participate in $W_{ijk}$
at first order.
This means that ${\cal L}_1$ does not modify
the scalar sector both
at quadratic and {\em cubic} order.
In other words, the prediction for non-Gaussianity of the curvature perturbation,
as well as that for the power spectrum, remains the same
in the presence of ${\cal L}_1$.
This is however not the case for ${\cal L}_1'$ and ${\cal L}_2$.

The covariant form of the Lagrangian can be recovered
by writing the extrinsic curvature as $K_{\mu\nu}=\gamma_\mu^{\;\rho}\gamma_\nu^{\;\sigma}\nabla_\rho u_\sigma$
and making use of the Gauss-Codazzi relations:
\begin{eqnarray}
\gamma^{\;\mu}_{\alpha}\gamma^{\;\nu}_{\beta}\gamma^{\;\gamma}_{\rho}\gamma^{\;\sigma}_{\delta}{\cal R}^{\rho}_{\;\sigma \mu \nu} &=& R^{\gamma}_{\;\delta \alpha \beta} + K^{\;\gamma}_{\alpha} K_{\delta \beta} - K^{\;\gamma}_{\beta} K_{\alpha \delta} ,\;\;\;\\
u^\mu\gamma_\alpha^{\;\nu}\gamma_\beta^{\;\rho}\gamma_\gamma^{\;\sigma}
{\cal R}_{\mu\nu\rho\sigma} &=& D_{\gamma} K_{\alpha\beta} - D_{\beta}K_{\alpha\gamma}.
\end{eqnarray}
We find
\begin{eqnarray}
{\cal L}_1 = \frac{\sqrt{-g}}{M^2}
C_{\mu\nu\rho\sigma}C_{\mu'\nu'\rho'\sigma'}\gamma^{\mu\mu'}\gamma^{\nu\nu'}\gamma^{\rho\rho'}u^\sigma u^{\sigma'},
\end{eqnarray}
where $C_{\mu\nu\rho\sigma}$ is the Weyl tensor.
Thus, it turns out that
${\cal L}_1$ reproduces the theory studied in Ref.~\cite{Deruelle:2012xv}.
One can repeat the same procedure also for ${\cal L}_2$ to derive its covariant form.
However, the covariant expression for ${\cal L}_2$ is messy and not so illuminating,
so that in this case it is better to work in the simpler ADM form.
It is worth noting that the covariant expression for ${\cal L}_2$ is also constructed
by contracting the Riemann curvature tensor with $u^\mu$ and hence it is Lorentz violating
in the same sense that ${\cal L}_1$ is.

In ${\cal L}_1$ and ${\cal L}_2$ in the ADM form one may consider time-dependent $M$.
This translates to the $\phi$-dependent coupling in the covariant language.
However, since $\phi$ is slowly rolling, it is reasonable to assume that
$M$ is only weakly time dependent. For simplicity we treat
$M$ as constant in the following.

\section{The tensor amplitude}

In the previous section we have identified the
two possible quadratic curvature terms
that make no contribution to the scalar sector of cosmological perturbations
at least at linear order.
Let us now investigate how
the amplitude of primordial tensor modes is modified due to those terms.
For clarity we study each term separately below.
Actually, we will find that a sizable modification
from ${\cal L}_2$ is prohibited because ${\cal L}_2$ also produces large non-Gaussianity.

\subsection{${\cal L}_1$}

First we consider
\begin{eqnarray}
S_{\rm higher}=\frac{1}{\kappa}\int{\rm d}^4x{\cal L}_1.
\end{eqnarray}
The quadratic action for the tensor perturbations is given by~\cite{Deruelle:2012xv}
\begin{align}
S=\frac{1}{8\kappa}\int{\rm d} t{\rm d}^3x\,a^3\left[
\dot h_{ij}^2-\frac{1}{a^2}(\partial_kh_{ij})^2
+\frac{4}{M^2a^2}(\partial_k\dot h_{ij})^2
\right].
\end{align}
Each Fourier mode of two polarization states, $h_k^\lambda(t)$ ($\lambda=+,\times$),
obeys a second-order evolution equation.
We use the canonically normalized variable
\begin{eqnarray}
f_k^\lambda(t) = \left(\frac{1}{4\kappa}\right)^{1/2}a^{3/2}\left(1+\frac{4k^2}{M^2a^2}\right)^{1/2} h_k^\lambda,
\end{eqnarray}
and omit $\lambda$ when unnecessary.
We then have
\begin{eqnarray}
\ddot f_k + \omega^2_k(t) f_k = 0,\label{evolution-f}
\end{eqnarray}
where
\begin{align}
\omega_k^2:= &\; -\frac{1}{4}\left(H^2+2\dot H \right)
+\frac{k^2/a^2-2H^2-\dot H}{1+4k^2/M^2a^2}
\notag \\
&\; -\frac{4H^2k^2/M^2a^2}{(1+4k^2/M^2a^2)^2}.\label{omega-eff}
\end{align}

We use the WKB solution,
\begin{eqnarray}
f_k\simeq \frac{1}{\sqrt{2\omega_k}}\exp\left[-\mathrm{i}\int^t\omega_k(t'){\rm d} t'\right],
\label{WKB-sol}
\end{eqnarray}
for the short wavelength modes with $k^2/a^2\gg H^2, M^2$.
In the short wavelength limit, Eq.~(\ref{omega-eff}) reduces to
\begin{eqnarray}
\omega_k^2 \simeq \frac{M^2}{4}-\frac{1}{4}\left(H^2+2\dot H \right).\label{omega-approx}
\end{eqnarray}
From this we see that during inflation the tensor perturbations are stable
if
\begin{eqnarray}
H< M.\label{stability}
\end{eqnarray}
We assume that the background evolution,
which is controlled by the inflaton sector $S_\phi$,
satisfies the condition~(\ref{stability}).
Since Eq.~(\ref{omega-approx}) gives the estimate
\begin{eqnarray}
\frac{\dot\omega_k}{\omega_k^2} \sim \frac{\epsilon H^3}{(M^2-H^2)^{3/2}},
\end{eqnarray}
where $\epsilon:=-\dot H/H^2\ll1$ is the slow-roll parameter,
the WKB approximation is justified as long as
$H$ is not too close to $M$.

In the long wavelength limit, $k^2/a^2\ll H^2, M^2$, we have
\begin{eqnarray}
\omega_k^2\simeq -\frac{9}{4}H^2-\frac{3}{2}\dot H=-\frac{(a^{3/2})\,\ddot{}}{a^{3/2}},
\end{eqnarray}
so that the standard result is recovered on superhorizon scales, $h_k\simeq$ const.

Let us compute the power spectrum of
the tensor modes,
\begin{eqnarray}
{\cal P}_T(k)=\frac{k^2}{\pi^2}\left|h_k \right|^2.
\end{eqnarray}
In general, Eq.~(\ref{evolution-f}) cannot be solved analytically,
and hence one needs numerical calculations to evaluate the power spectrum.
However, in the special case of the exact de Sitter background
one can solve Eq.~(\ref{evolution-f}) analytically
using the hypergeometric functions.
This was done in Ref.~\cite{Deruelle:2012xv}, and here
we quote their final result:
\begin{eqnarray}
{\cal P}_T=\frac{2\kappa H^2}{\pi^2}\Xi_1(H/M),\label{dS-formula}
\end{eqnarray}
where
\begin{eqnarray}
\Xi_1(x):=\frac{\cosh(\pi\nu/2)\coth(\pi\nu/2)|\Gamma(-1/4+\mathrm{i}\nu/4)|^4}{128\pi^2x^3},
\end{eqnarray}
with $\nu:=\sqrt{x^{-2}-1}$.
Based on this, one may expect that in the case where $H$ is varying
the power spectrum is given by
evaluating the de Sitter result~(\ref{dS-formula}) at horizon crossing,
\begin{eqnarray}
{\cal P}_T(k)=\left. \frac{2\kappa H^2}{\pi^2}\Xi_1(H/M)\right|_{k=aH},
\label{PT-at-hc}
\end{eqnarray}
as is commonly done in general relativity.

We numerically solved Eq.~(\ref{evolution-f}) in the case of power-law inflation, $a\propto t^p$,
using the initial condition~(\ref{WKB-sol}),
and verified that Eq.~(\ref{PT-at-hc}) reproduces
the numerical results very accurately, as shown in Fig.~\ref{fig:analyticsp.eps}.
We are thus allowed to use the formula~(\ref{PT-at-hc})
for slow-roll inflation.

\begin{figure}[tb]
  \begin{center}
    \includegraphics[keepaspectratio=true,height=52mm]{Weyl_inf_3_04-pic1-2.eps}
  \end{center}
  \caption{
The power spectrum ${\cal P}_T$ of tensor modes from power-law inflation with the correction from ${\cal L}_1$.
The upper line and points are for $a\propto t^{50}$ and the lower for $a\propto t^{400}$, respectively.
  Red points represent the numerical results, while the dashed line indicates the analytic estimate~(\ref{PT-at-hc}).
  The parameters for the upper line and points are given by $\sqrt{\kappa}H|_{t_{\rm end}}=10^{-4}$ and $H/M|_{t_{\rm end}}=0.16$,
  where $t_{\rm end}$ is the time at the end of inflation.
  The parameters for the lower ones are given by $\sqrt{\kappa}H|_{t_{\rm end}}=10^{-4}$ and $H/M|_{t_{\rm end}}=0.74$.
  }%
  \label{fig:analyticsp.eps}
\end{figure}

\begin{figure}[tb]
  \begin{center}
    \includegraphics[keepaspectratio=true,height=52mm]{consis_weyl2.eps}
  \end{center}
  \caption{${\cal D}$ as a function of $H/M$.
  }%
  \label{fig:consis_weyl.eps}
\end{figure}

The behavior of the function $\Xi_1(x)$ is as follows:
$\Xi_1\to 1$ as $x\to 0$, and $\Xi_1<1$ for $0<x\lesssim 0.95$.
Thus, the tensor amplitude is suppressed for $H\lesssim 0.95 M$.
The minimum of $\Xi_1$ is given by $\Xi_1\simeq 0.65$, which occurs at $x\simeq 0.74$, and
$\Xi_1$ diverges as $x\to 1$.

Since the power spectrum of the curvature perturbation remains unchanged
in our theory, the tensor-to-scalar ratio is given by
\begin{eqnarray}
r=16\epsilon \Xi_1. \label{eq-r}
\end{eqnarray}
The tensor tilt, $n_T:={\rm d}\ln{\cal P}_T/{\rm d}\ln k$, is evaluated as
\begin{eqnarray}
n_T=-\frac{2\epsilon}{1-\epsilon}
\left.\left[1+\frac{1}{2}\frac{{\rm d}\ln\Xi_1}{{\rm d}\ln(H/M)} \right]\right|_{k=aH}.\label{eq-nt}
\end{eqnarray}
We see that ${\rm d}\ln \Xi_1/{\rm d}\ln x < 0$ for $0<x\lesssim 0.74$
and its minimum value is given by ${\rm d}\ln \Xi_1/{\rm d}\ln x\simeq -0.46$ at $x\simeq 0.53$.
This shows that the tensor spectrum is always red.

From Eqs.~(\ref{eq-r})--(\ref{eq-nt}) it is clear that
the consistency relation~\cite{Lidsey} is violated.
The deviation from the standard consistency relation is characterized by
the following function,
\begin{eqnarray}
{\cal D}:=\left.\frac{1+(1/2){\rm d} \ln\Xi_1/{\rm d} \ln x}{\Xi_1} \right|_{x=H/M},
\end{eqnarray}
as $-8n_T/r\simeq {\cal D}|_{k=aH}$.
In Fig.~\ref{fig:consis_weyl.eps}, we plot ${\cal D}$ as a function of $H/M$.
We see that the violation depends on the scale $k$, and Fig.~\ref{fig:consis_weyl.eps}
tells us its scale dependence.

Figure~\ref{fig:planck_const6_1.eps} illustrates the observational implications of the ${\cal L}_1$ correction
by comparing the suppressed tensor amplitude with the Planck results.
The red stars in the figure indicate the case of power-law inflation,
assuming the maximal suppression ($\Xi_1=0.65$) at the observed scale.
Although original power-law inflation (represented by the dashed line)
is ruled out by observations, it can be within the 2$\sigma$ contour with the help of ${\cal L}_1$.
The same applies to other inflation models such as $V\propto \phi^2$.
Those models originally predict large tensor modes,
but the ${\cal L}_1$ correction can bring such models to
the observationally preferred region in the $n_s$-$r$ plane.

\begin{figure*}[htb]
  \begin{center}
    \includegraphics[keepaspectratio=true,height=90mm]{planck_const6_2.eps}
  \end{center}
  \caption{Reduced tensor amplitude and the Planck results~\cite{PlanckCollaboration2015b} in the $n_s$-$r$ plane.
  Red stars correspond to power-law inflation with the ${\cal L}_1$ correction, assuming the
  maximal suppression, $\Xi_1=0.65$.
  }%
  \label{fig:planck_const6_1.eps}
\end{figure*}

\subsection{${\cal L}_2$}

Next let us consider
\begin{eqnarray}
S_{\rm higher}=- \frac{1}{2\kappa}\int{\rm d}^4x{\cal L}_2.
\end{eqnarray}
Here we added a minus sign so that the
tensor perturbations are stable at high momenta.
The quadratic action for the tensor perturbations is
\begin{align}
S=\frac{1}{8\kappa}\int{\rm d} t{\rm d}^3x\,a^3\left[
\dot h_{ij}^2-\frac{1}{a^2}(\partial_kh_{ij})^2
-\frac{1}{M^2a^4}(\partial^2 h_{ij})^2
\right].
\end{align}
Using the conformal time defined by ${\rm d} \eta = {\rm d} t/a$ and
the canonically normalized variable $v_k^\lambda:=(4\kappa)^{-1/2}a h_k^\lambda$ in the Fourier space,
we obtain
\begin{eqnarray}
\frac{{\rm d}^2v_k}{{\rm d}\eta^2}+\omega^2_k(\eta)v_k = 0,\label{veom}
\end{eqnarray}
where
\begin{eqnarray}
\omega_k^2:=k^2+\frac{k^4}{M^2a^2}-\frac{1}{a}\frac{{\rm d}^2a}{{\rm d}\eta^2}.
\end{eqnarray}
This modified dispersion relation has been studied in detail
in the literature~\cite{Martin:2002kt,Ashoorioon:2011eg}.

The WKB solution
\begin{eqnarray}
v_k\simeq \frac{1}{\sqrt{2\omega_k}}\exp\left[-\mathrm{i}\int^\eta\omega_k(\eta'){\rm d}\eta'\right],
\label{vWKB}
\end{eqnarray}
may be used for the short wavelength modes, because $\omega_k^{-2}{\rm d}\omega_k/{\rm d}\eta \ll 1$
is always satisfied at large $k$.

Only in the case of exact de Sitter inflation for which
the scale factor is given by $a=1/H(-\eta)$,
Eq.~(\ref{veom}) can be solved analytically.
The solution that matches Eq.~(\ref{vWKB}) at large $k$ is obtained
in terms of the Whittaker function as~\cite{Ashoorioon:2011eg,Kobayashi:2015gga}
\begin{eqnarray}
v_k=\frac{e^{-\pi/8x}W_{\mathrm{i}/4x,3/4}(-\mathrm{i}xk^2\eta^2)}{\left(-2xk^2\eta\right)^{1/2}},
\end{eqnarray}
where $x=H/M$. This yields the power spectrum
\begin{eqnarray}
{\cal P}_T=\frac{2\kappa H^2}{\pi^2}\Xi_2(H/M),\label{P2dS}
\end{eqnarray}
where
\begin{eqnarray}
\Xi_2(x):=\frac{\pi}{4}\left[e^{\pi/(4x)}x^{3/2}\left|\Gamma(5/4+\mathrm{i}/(4x))\right|^2 \right]^{-1}.
\end{eqnarray}
In the case of slow-roll inflation, one may evaluate the de Sitter result~(\ref{P2dS})
at horizon crossing, $k=aH$. This can also be justified by a numerical calculation.

One sees that $\Xi_2$ is a monotonically decreasing function
and $\Xi_2\to 1$ as $x\to 0$. Therefore, also in this case
the tensor amplitude is suppressed relative to the standard result.
Since $\Xi_2\propto x^{-3/2}$ for $x\gg 1$, the tensor amplitude could potentially be
suppressed to a very small level.
However, as we see below, this possibility is hindered by
the generation of large non-Gaussianity.

In contrast to the case of ${\cal L}_1$, the cubic action for the curvature perturbation
is affected by ${\cal L}_2$. This implies that $M$ must be sufficiently large
in order to avoid large non-Gaussianities in $\zeta$.
Typically, ${\cal L}_2$ contains terms such as
\begin{eqnarray}
{\cal L}_2 \sim \frac{1}{ M^{2}} \zeta (\partial ^{2} \zeta)^{2}.
\end{eqnarray}
The non-Gaussianity generated by this term is estimated to be
\begin{eqnarray}
f_{\rm NL} \sim \frac{H^{2}}{\epsilon M^{2}} \,.
\end{eqnarray}
Requiring that $f_{\rm NL}\lesssim 1$, we have
\begin{equation}
\frac{H}{M} \lesssim \epsilon^{1/2} \ll 1.
\end{equation}
Therefore, in fact the suppression factor $\Xi_2$ cannot be much smaller than 1.
We conclude that the second Lagrangian ${\cal L}_2$ does not provide
an efficient way of suppressing the tensor amplitude.

One may consider a combination of the two Lagrangians, $a{\cal L}_1+b{\cal L}_2$.
Obviously, this does not change the quadratic Lagrangian for the scalar perturbations, and
to avoid large non-Gaussianities we must require $b\ll 1$.
Therefore, to suppress the tensor amplitude most effectively, essentially one can only use ${\cal L}_{1}$.

\section{Conclusions}

In this paper, we have studied inflationary predictions of theories with quadratic curvature corrections.
We began by looking for ghost-free quadratic curvature terms
that retain the same quadratic action for the curvature perturbation as in general relativity
while modifying the dynamics of tensor perturbations.
We have shown that such curvature terms can indeed be constructed, and
determined the two possible combinations (denoted by ${\cal L}_1$ and ${\cal L}_2$).
This was done by using the ADM formalism, and
recast in a covariant form those corrections contain
the curvature tensors contracted with the unit normal $u^\mu$ to
hypersurfaces on which the inflaton is homogeneous. It has turned out that
one of the two terms, ${\cal L}_1$, is in fact identical to the one introduced in
so-called Lorentz-violating Weyl gravity~\cite{Deruelle:2012xv}.
This term does not change the action of the curvature perturbation even at cubic order.
The other term, ${\cal L}_2$, in contrast, modifies the scalar sector at cubic order.

We have investigated the tensor amplitude in the presence of the
quadratic curvature corrections ${\cal L}_1$ and ${\cal L}_2$.
The analytic results were known only for exact de Sitter inflation,
and we have used the de Sitter formulas evaluated at horizon crossing
in the case of slow-roll inflation for which the Hubble parameter is varying.
The validity of the method has been checked by performing numerical calculations.
Both ${\cal L}_1$ and ${\cal L}_2$ reduce the amplitude of
primordial tensor perturbations. However, we have found that
${\cal L}_2$ could generate large non-Gaussianity of the curvature perturbation,
which places a stringent constraint on the amount of the suppression due to ${\cal L}_2$.
Since ${\cal L}_1$ does not change the cubic interaction of the curvature perturbation,
this evades the non-Gaussianity constraint.
The tensor power spectrum can be as small as 65\%
of the standard result due to ${\cal L}_1$, which
brings many inflation models with large tensor modes to the observationally preferred region
in the $n_s$-$r$ plane.
We have seen that the tensor tilt is also modified, though the spectrum can never be blue.

\acknowledgments
The authors thank Yuuiti Sendouda for helpful comments and discussion.
K.Y. was supported in part by Rikkyo University Special Fund for Research.
This work was supported in part by JSPS Grant-in-Aid for Young
Scientists (B) Grant No.~24740161 (T.K.).


\section{Introduction}
Theoretical investigation of Ruderman\cite{Ruder} leads to the conclusion that matter may be anisotropic in densities of order $10^{15}\;gm\;cm^{-3}$. The impact of anisotropy on steller configuration may be found in pioneering works of Bower and Liang\cite{BL} and Herrera and Santos\cite{HS}. Anisotropy may occur due to existance of type 3A superfluids\cite {Ruder,BL,KW}, or phase transition\cite{Soko}. Consenza \textit{et. al.}\cite{CHEW} developed the procedure to obtain anisotropic solutions from isotropic solutions of Einstein's field equations. Tikekar and Thomas\cite{TT} found exact solutions of Einstein's field equations for anisotropic fluid sphere on pseudo spheroidal spacetime. The key feature of their model is the high variation of density from centre to boundary of stellar configuration. The class of exact anisotropic solutions on spherically symmetric spacetime has been obtained by Mak and Harko\cite{MH}. Karmakar \textit{et. al.}\cite{KMSM} analysed the role of pressure anisotropy for Vaidya-Tikekar\cite{VT} model. Paul \textit{et. al}\cite{PCKT} developed anisotropic stellar model for strage star. A core-envelope model describing superdense stars with anisotropic fluid distribution has been obtained by Thomas and Ratanpal\cite{TR}, Thomas \textit{et. al.}\cite{TRV}  \& Tikekar and Thomas\cite{TT1}. Hence the study of anisotropic fluid distribution is important in general theory of relativity.\\\\
\noindent The study of Einstein-Maxwell system is carried out by several authors. Patel and Kopper\cite{PK} obtained charged analog of Vaidya-Tikekar\cite{VT} solution. The study of analytic models of quark stars have been carrried out by Komathiraj and Maharaj\cite{KM}, they found a class of solutions of Einstein-Maxwell system. Charged anisotropic matter with linear equation of state have been extensively studied by Thirukkanesh and Maharaj\cite{TM}.\\\\
\noindent Hence, both anisotropy and electromagnatic field are important in relativistic astrophysics. In this paper charged anisotropic model of a stellar configuration have been studied on the background of paraboloidal spacetime. Sec-2 describes, the field equations for charged static stallar configuration on paraboloidal spacetime and their solution is obtained. Sec-3 described the physical plausibility condition and Sec-4 contains discussion.

\section{Field Equations and Solution}
The interior of stellar configuration is described by static spherically symmetric paraboloidal spacetime metric,
\begin{equation}\label{ISpaceTime1}
	ds^{2}=e^{\nu(r)}dt^{2}-\left(1+\frac{r^{2}}{R^{2}} \right)dr^{2}-r^{2}\left(d\theta^{2}+\sin^{2}\theta d\phi^{2}\right),
\end{equation}
with the energy-momentum tensor for anisotropic charged fluid,
\begin{equation}\label{EMTensor}
	T_{ij}=diag\left(\rho+E^{2},\;p_{r}-E^{2},\;p_{t}+E^{2},\;p_{t}+E^{2} \right),
\end{equation}
where $\rho$ is the energy density, $p_{r}$ is the radial pressure, $p_{t}$ is the tangential pressure and $E$ is the electric field intensity. These quantities are measured relative to the comoving fluid velocity $u^{i}=e^{-\nu}\delta_{0}^{i}$. For the spacetime metric (\ref{ISpaceTime1}) and energy-momentum tensor (\ref{EMTensor}), the Einstein-Maxwell system takes the form,
\begin{equation}\label{rho1}
	\rho+E^{2}=\frac{3+\frac{r^{2}}{R^{2}}}{R^{2}\left(1+\frac{r^{2}}{R^{2}} \right)^{2}},
\end{equation}
\begin{equation}\label{pr1}
	p_{r}-E^{2}=\frac{\nu'}{r\left(1+\frac{r^{2}}{R^{2}}\right)}-\frac{1}{R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)},
\end{equation}
\begin{equation}\label{pt1}
	p_{t}+E^{2}=\frac{1}{1+\frac{r^{2}}{R^{2}}}\left[\frac{\nu''}{2}+\frac{\nu'^{2}}{4}+\frac{\nu'}{2r}\right]-\frac{\nu' r}{2R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)^{2}}-\frac{1}{R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)^{2}},
\end{equation}
\begin{equation}\label{Sigma1}
	\sigma=\frac{\left(r^{2}E\right)'}{r^{2}\sqrt{1+\frac{r^{2}}{R^{2}}}},
\end{equation}
where $\sigma$ is the proper charge density and prime denotes differentiation with respect to $r$. In field equations (\ref{rho1}) - (\ref{Sigma1}), velocity of light $c$ is taken as $1$ also $\frac{8\pi G}{c^{4}}=1$.\\\\
The anisotropic parameter $\Delta$ is defined as,
\begin{equation}\label{Delta1}
	\Delta=p_{t}-p_{r}.
\end{equation}
To solve the system (\ref{rho1}) - (\ref{Sigma1}) radial pressure is assumed to be of the form,
\begin{equation}\label{pr2}
	p_{r}=\frac{p_{0}\left(1-\frac{r^{2}}{R^{2}}\right)}{R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)^{2}},
\end{equation}
where $p_{0}>0$ is the model parameter and $\frac{p_{0}}{R^{2}}$ is central pressure. At the boundary of the star $r=R$, $p_{r}$ must vanish, which gives $r=R$ as the radius of the star.\\\\
This form of radial pressure is prescribed by Sharma and Ratanpal\cite{SR} to describe anisotropic stellar model admitting a quadratic equation of state on paraboloidal spacetime. Equations (\ref{pr2}) and (\ref{pr1}) gives,
\begin{equation}\label{NuDash1}
	\nu'=\frac{p_{0}r\left(1-\frac{r^{2}}{R^{2}}\right)}{R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)}+\frac{r}{R^{2}}-E^{2}r\left(1+\frac{r^{2}}{R^{2}}\right).
\end{equation}
We assume electric field intensity of the form,
\begin{equation}\label{E1}
	E^{2}=\frac{k\frac{r^{2}}{R^{2}}}{R^{2}\left(1+\frac{r^{2}}{R^{2}} \right)^{2}},
\end{equation}
where $k\geq 0$ is a model parameter, from equation (\ref{E1}) it is clear that $E$ is decreasing in radially outward direction. Equations (\ref{NuDash1}) and (\ref{E1}) leads to,
\begin{equation}\label{NuDash2}
	\nu'=\frac{\left(2p_{0}+k\right)r}{R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)}+\left(1-p_{0}-k\right)\frac{r}{R^{2}},
\end{equation}
and hence,
\begin{equation}\label{Nu}
	\nu=\log\left[C\left(1+\frac{r^{2}}{R^{2}}\right)^{\left(\frac{2p_{0}+k}{2}\right)}\right]+\left(\frac{1-p_{0}-k}{2}\right)\frac{r^{2}}{R^{2}},
\end{equation}
where $C$ is constant of integration. Therefore spacetime metric (\ref{ISpaceTime1}) is written as,
\begin{equation}\label{ISpaceTime2}
	ds^{2}=C\left(1+\frac{r^{2}}{R^{2}}\right)^{\left(\frac{2p_{0}+k}{2} \right)}e^{\left(\frac{1-p_{0}-k}{2}\right)\frac{r^{2}}{R^{2}}}dt^{2}-\left(1+\frac{r^{2}}{R^{2}} \right)dr^{2}-r^{2}\left(d\theta^{2}+\sin^{2}\theta d\phi^{2} \right).
\end{equation}
The spacetime metric (\ref{ISpaceTime2}) should continously match with Reissner-Nordstrom spacetime metric
\begin{equation}\label{ESpaceTime}
	ds^{2}=\left(1-\frac{2m}{r}+\frac{Q^{2}}{r^{2}}\right)dt^{2}-\left(1-\frac{2m}{r}+\frac{Q^{2}}{r^{2}}\right)^{-1}dr^{2}-r^{2}\left(d\theta^{2}+\sin^{2}\theta\right)d\phi^{2},
\end{equation}
at the boundary of the star $r=R$, where $\left(p_{r}\right)\left(r=R\right)=0$. This matching conditions gives
\begin{equation}\label{M}
	M=\frac{k+2R^{2}}{8R},
\end{equation}
and
\begin{equation}\label{C}
	C=\frac{e^{\left(\frac{p_{0}+k-1}{2}\right)}}{\left(2p_{0}+2+k\right)/2},
\end{equation}
where $M$ is the mass enclosed the spherical body of radius $R$, hence the electric field intensity parameter $k$ directly effects the mass of the star.
Equations (\ref{pr1}), (\ref{pt1}), (\ref{Delta1}) and (\ref{NuDash2}) gives anisotropic parameter $\Delta$ as,
\begin{equation}\label{Delta2}
	\Delta=\frac{\frac{r^{2}}{R^{2}}\left[X_{1}+Y_{1}\frac{r^{2}}{R^{2}}+Z_{1}\frac{r^{4}}{R^{4}}\right]}{4R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)^{3}},
\end{equation}
which vanishes at $r=0$, where $X_{1}=p_{0}^{2}-8p_{0}-12k+3$, $Y_{1}=-2p_{0}^{2}-2p_{0}k+2p_{0}-8k+4$, $Z_{1}=1+p_{0}^{2}+k^{2}-2p_{0}-2k+2p_{0}k$. Equations (\ref{rho1}) and (\ref{E1}) gives,
\begin{equation}\label{rho2}
	\rho=\frac{3+\left(1-k\right)\frac{r^{2}}{R^{2}}}{R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)^{2}},
\end{equation}
and from equations (\ref{Delta1}), (\ref{pr2}) and (\ref{Delta2}) we get expression of $p_{t}$ as,
\begin{equation}\label{pt2}
	p_{t}=\frac{4p_{0}+X_{1}\frac{r^{2}}{R^{2}}+Y_{1}\frac{r^{4}}{R^{4}}+Z_{1}\frac{r^{6}}{R^{6}}}{4R^{2}\left(1+\frac{r^{2}}{R^{2}}\right)^{3}}.
\end{equation}
Hence equations (\ref{rho2}), (\ref{pr2}), (\ref{pt2}), (\ref{E1}) and (\ref{Delta2}) describes matter density, radial pressure, tangential pressure, electric field intensity and measure of anisotropy respectively.
\section{Physical Plausibility Conditions}
Following Delgaty and Lake\cite{Delgaty} we impose the following condition on the sytem to make model physically acceptable.
\begin{itemize}
\item [(i)] $\rho(r),~p_{r}(r),~p_{t}(r) \geq 0 $ for $ 0 \leq r \leq R$.
\item [(ii)] $\rho-p_{r}-2p_{t} \geq 0$ for $ 0 \leq r \leq R$.
\item [(iii)] $\frac{d\rho}{dr},~\frac{dp_{r}}{dr},~\frac{dp_{t}}{dr} < 0$ for $0 \leq r \leq R$.
\item [(iv)] $0 \leq \frac{dp_{r}}{d\rho} \leq 1$; $0 \leq \frac{dp_{t}}{d\rho} \leq 1$, for $0 \leq r \leq R$.
\end{itemize}
From equation (\ref{rho2}), $\rho(r=0)=\frac{3}{R^{2}}>0$ and $\rho(r=R)=\frac{4-k}{4R^{2}}$,
therefore, $\rho>0$ for $0\leq r\leq R$ if $k\leq4$, i.e.
\begin{equation}\label{k1}
	0\leq k\leq 4.
\end{equation}
From equation (\ref{pr2}) $p_{r}(r=0)=\frac{p_{0}}{R^{2}}>0$ as $p_{0}>0$ and $p_{r}(r=R)=0$. Hence $p_{r}\geq 0$ for $0\leq r\leq R$. It is required that $p_{t}\geq 0$ for $0\leq r\leq R$ and further to get the simple bounds on $p_{0}$ and $k$, we assume $p_{r}=p_{t}$ at $r=R$.\\\\
from (\ref{pt2}), $p_{t}(r=0)=\frac{p_{0}}{R^{2}}>0$ as $p_{0}>0$, and $p_{t}(r=R)=\frac{k^{2}-22k+8-8p_{0}}{32R^{2}}$, at $r=R$, $p_{t}=p_{r}=0$ if,
\begin{equation}\label{p0}
	p_{0}=\frac{k^{2}-22k+8}{8},
\end{equation}
but $k$ should be chosen such that $p_{0}$ is positive, which restrict the value of $k$ as,
\begin{equation}\label{k2}
	k<0.3699.
\end{equation}
Hence,
\begin{equation}\label{k3}
	0\leq k<0.3699,\;\;\;\;\;p_{0}=\frac{k^{2}-22k+8}{8}
\end{equation}
which is condition for positivity of $p_{t}$.\\\\
Hence condition (i) is satisfied throughout the star. For the values of $k$ and $p_{0}$ specified in (\ref{k3}), programatically it has been varified that condition (ii) i.e. energy condition is satisfied throughout the star. From equation (\ref{rho2}),
\begin{equation}\label{drhodr}
	\frac{d\rho}{dr}=-\frac{2r}{R^{4}}\frac{\left[(5+k)+(1-k)\frac{r^{2}}{R^{2}}\right]}{\left(1+\frac{r^{2}}{R^{2}} \right)^{3}},
\end{equation}
from equation (\ref{drhodr}), $\left(\frac{d\rho}{dr}\right)(r=0)=0$ and $\left(\frac{d\rho}{dr}\right)(r=R)=-\frac{3}{2R^{3}}<0$. Hence $\rho$ is decreasing throught the star. From equation (\ref{pr2}),
\begin{equation}\label{dprdr}
	\frac{dp_{r}}{dr}=\frac{-2p_{0}r\left(3-\frac{r^{2}}{R^{2}}\right)}{R^{4}\left(1+\frac{r^{2}}{R^{2}}\right)^{3}}.
\end{equation}
Now, $\left(\frac{dp_{r}}{dr}\right)(r=0)=0$ and $\left(\frac{dp_{r}}{dr}\right)(r=R)=\frac{-p_{0}}{2R^{3}}<0$ as $p_{0}>0$. Hence $p_{r}$ is decreasing throughout the star. From equation (\ref{pt2}),
\begin{equation}\label{dptdr}
	\frac{dp_{t}}{dr}=\frac{r\left[X_{2}+Y_{2}\frac{r^{2}}{R^{2}}+Z_{2}\frac{r^{4}}{R^{4}}\right]}{2R^{4}\left(1+\frac{r^{2}}{R^{2}}\right)^{4}},
\end{equation}
where $X_{2}=p_{0}^{2}-20p_{0}-12k+3$, $Y_{2}=-6p_{0}^{2}+12p_{0}-4p_{0}k+8k+2$ and $Z_{2}=5p_{0}^{2}-4p_{0}+8p_{0}k+3k^{2}+2k-1$. Now $\left(\frac{dp_{t}}{dr}\right)(r=0)=0$ and $\left(\frac{dp_{t}}{dr}\right)(r=R)=\frac{-12p_{0}+3k^{2}-2k+4p_{0}k+4}{32R^{3}}$. Substituting $p_{0}$ from equation (\ref{p0}) in $\left(\frac{dp_{t}}{dr}\right)(r=R)$, $\left(\frac{dp_{t}}{dr}\right)(r=R)<0$ if $4k^{3}-76k^{2}+280k-64<0$. This further restrict value of $k$ as $k<0.2446$.\\\\
Hence, if
\begin{equation}\label{k4}
	0\leq k<0.2446,\;\;\;\;\;p_{0}=\frac{k^{2}-22k+8}{8},
\end{equation}
then $\frac{d\rho}{dr}$, $\frac{dp_{r}}{dr}$ and $\frac{dp_{t}}{dr}$ are decreasing in radially outward direction between $0\leq r\leq R$. From equations (\ref{drhodr}) and (\ref{dprdr}) we have,
\begin{equation}\label{dprdrho}
	\frac{dp_{r}}{d\rho}=\frac{p_{0}\left(3-\frac{r^{2}}{R^{2}}\right)}{\left[(5+k)+(1-k)\frac{r^{2}}{R^{2}}\right]}.
\end{equation}
At the centre of the star $\left(\frac{dp_{r}}{d\rho}\right)(r=0)<1$ if $k<24.8810$, which is consistent with condition (\ref{k4}) and at the boundary of the star $\left(\frac{dp_{r}}{d\rho}\right)(r=R)<1$ if $k<22.7047$, which is also consistent with conditon (\ref{k4}). From equations (\ref{drhodr}) and (\ref{dptdr}) we have,
\begin{equation}\label{dptdrho}
	\frac{dp_{t}}{d\rho}=\frac{-\left[X_{2}+Y_{2}\frac{r^{2}}{R^{2}}+Z_{2}\frac{r^{4}}{R^{4}}\right]}{4\left(1+\frac{r^{2}}{R^{2}}\right)\left[(5+k)+(1-k)\frac{r^{2}}{R^{2}}\right]}.
\end{equation}
Now, $\left(\frac{dp_{t}}{d\rho}\right)(r=0)<1$ if $k<19.4283$, which is consistent with condition (\ref{k4}) and $\left(\frac{dp_{t}}{d\rho}\right)(r=R)<1$ if $k<6.6371$, which is also consistent with condition (\ref{k4}). Hence for $0\leq k<0.2446$ and $p_{0}=\frac{k^{2}-22k+8}{8}$, all the physical plausibility conditions are satisfied.
\section{Dicussion}
Certain ascepts of charged relativistic star on paraboloidal spacetime have been discussed. It is observed that all the physical plausibility conditions are satisfied for $0\leq k\leq 0.2446$ and $p_{0}=\frac{k^{2}-22k+8}{8}$. The plots of $\rho$ (charged, uncharged), $p_{r}$ \& $p_{t}$ (charged), $p_{r}$ \& $p_{t}$ (uncharged), anisotropy $\Delta$ (charged, uncharged), $\rho-p_{r}-2p_{t}$ (charged, uncharged), $\frac{dp_{r}}{d\rho}$ \& $\frac{dp_{t}}{d\rho}$ (Charged), $\frac{dp_{r}}{d\rho}$ \& $\frac{dp_{t}}{d\rho}$ (uncharged) over $\frac{r^{2}}{R^{2}}$ for $R=10$, $k=0.2$ and taking $G=c^{2}=1$ are shown in figures 1 to 7. It is observed that energy condition is satisfied throughout the star. When $k=0$ the value of $p_{0}=1$ and the model reduces to the Sharma \& Ratanpal\cite{SR} model. Hence the model described here is charged generalization of particular case $p_{0}=1$ of uncharged Sharma \& Ratanpal\cite{SR} model.
\begin{figure}[h]
\begin{center}
\includegraphics[width=12cm]{rho.eps}\\
\caption{Variation of density ($\rho$) (charged \& uncharged) against $\frac{r^{2}}{R^{2}}$.}\label{rho}
\end{center}
\end{figure}

\pagebreak
\begin{figure}[h]
\begin{center}
\includegraphics[width=12cm]{pressureC.eps}\\
\caption{Variation of radial pressure ($p_{r}$) \& tangential pressure ($p_{t}$) (charged) against $\frac{r^{2}}{R^{2}}$.}\label{pressureC}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=12cm]{pressureU.eps}\\
\caption{Variation of radial pressure ($p_{r}$) \& tangential pressure ($p_{t}$) (uncharged) against $\frac{r^{2}}{R^{2}}$.}\label{pressureU}
\end{center}
\end{figure}

\pagebreak
\begin{figure}[h]
\begin{center}
\includegraphics[width=12cm]{anisotropy.eps}\\
\caption{Variation of anisotropy ($\Delta$) (charged \& uncharged) against $\frac{r^{2}}{R^{2}}$.}\label{anisotropy}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=12cm]{energycondition.eps}\\
\caption{Variation of energy condition ($\rho-p_{r}-2p_{t}$) (charged \& uncharged) against $\frac{r^{2}}{R^{2}}$.}\label{anisotropy}
\end{center}
\end{figure}

\pagebreak
\begin{figure}[h]
\begin{center}
\includegraphics[width=12cm]{soundspeedC.eps}\\
\caption{Variation of $\frac{dp_{r}}{d\rho}$ \& $\frac{dp_{t}}{d\rho}$ (charged) against $\frac{r^{2}}{R^{2}}$.}\label{soundspeedC}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=12cm]{soundspeedU.eps}\\
\caption{Variation of $\frac{dp_{r}}{d\rho}$ \& $\frac{dp_{t}}{d\rho}$ (uncharged) against $\frac{r^{2}}{R^{2}}$.}\label{soundspeedC}
\end{center}
\end{figure}

\pagebreak


\section{Introduction}

\label{sect:intro}

In recent years, attempts to understand the origin of the accelerated cosmic expansion \citep[][]{riess1998,perlmutter1999} have led to a large number of theoretical models \citep[][]{cst2006}. Apart from the current standard model, which assumes that the acceleration is caused by a cosmological constant $\Lambda$ (hence the name $\Lambda$-cold-dark-matter, or $\Lambda$CDM), these model can be roughly put in two categories: dark energy -- which replaces $\Lambda$ by some dynamical field, and modified gravity -- which assumes that there is no exotic matter species beyond the standard CDM model but gravity is not described by General Relativity (GR) on cosmological scales \citep[see, e.g.,][for some recent reviews]{clifton2012,joyce2015}. Some of the models in the latter class, such as the chameleon theory \citep{chameleon}, of which the well-known $f(R)$ gravity \citep{cddatt2005} is an example, have been active research topics in recent years.

Ultimately, any new cosmological model or theory of gravity should be put to rigorous tests against observational data. For this reason, a number of tests have been proposed or applied in the past to examine the viability of the models \citep[see, e.g.,][for some recent reviews in $f(R)$ gravity]{lombriser2014,demartino2015}. The present paper shall follow the same line to propose a test using observations of galaxy clusters.

Galaxy clusters are the largest gravitationally bound and virialised objects in our Universe. The most massive clusters observed today typically have masses in the range of $\sim10^{14}-10^{15}h^{-1}M_\odot$, of which the dominant component is dark matter. These are homes to galaxies, stars and eventually lives, which together hold the vast majority of the information that can be extracted from cosmological and astrophysical observations. In dark energy or modified gravity theories, the different cosmic expansion histories and gravitational laws between particles can have sizeable effects on how the clusters form and evolve. \citet{schmidt2009, cataneo2014},  based on this observation, have placed constraints on $f(R)$ gravity using cluster abundance data. In theories such as $f(R)$ gravity, massive and massless particles feel different strengths of gravity, thus allowing these theories to be constrained by comparing the so-called dynamical and lensing masses of clusters \citep{schmidt2010,zlk2011}. Combining with lensing observations, \citet{terukina2014,wilcox2015} obtained even stronger constraints on $f(R)$ gravity.

Given the abundant information associated with these rich objects, one expects that they will provide a wealth of other potential tests of new cosmological models. The arrival of the era of precision cosmology lends this perspective both more interest and more support. In this paper, we propose to utilise the observationally inferred mass fraction of hot X-ray emitting gas in galaxy clusters as a new test of gravity.

The gas fraction in clusters, $f_{\rm gas}$, is a well established and understood observable in the standard cosmological model, which can be used to place strong constraints on background cosmology \citep[see, e.g.,][and following-up works, for some examples]{allen2004,allen2008}. The basic assumptions (or approximations) are (i) clusters are in hydrostatic equilibrium between thermal pressure and gravity and (ii) as the largest objects in the Universe, the cluster baryon fraction, dominantly contributed by gas, is a faithful representation of the cosmological average baryon fraction $\Omega_{\rm b}/\Omega_{\rm m}$, in which $\Omega_{\rm b}$ and $\Omega_{\rm m}$ are respectively the fractional mass density of baryons and all matter \citep{white1993,eke1998}. Using (i), one can find, from the observed X-ray temperature and surface brightness profiles, the mass profiles of baryons and all matter inside a cluster, and consequently $f_{\rm gas}(r)$ -- the profile of gas fraction. Combined with (ii), this can have a say about $\Omega_{\rm m}$ provided that $\Omega_{\rm b}$ is measured elsewhere (e.g., from big bang nucleosynthesis or the cosmic microwave background (CMB)).

In modified gravity theories, the hydrostatic equilibrium inside clusters is changed as a result of the different law of gravity. Hence, a cluster can have a higher dynamical mass, with the same baryonic mass inside, leading to a $f_{\rm gas}^{\rm obs}$ that is lower than the true $f_{\rm gas}$ which is related to $\Omega_{\rm b}/\Omega_{\rm m}$. If a cosmologist wishes to infer $\Omega_{\rm b}/\Omega_{\rm m}$ from $f^{\rm obs}_{\rm gas}$ in a modified gravity universe, a correction has to be done to the end result, which can lead to inconsistency with other observational determinations of $\Omega_{\rm b}/\Omega_{\rm m}$, and therefore a constraint on the gravity theory. We shall demonstrate the potential constraints from this test using $f(R)$ gravity as an example.

This paper is organised as following: in Sect.~\ref{sect:fr_gravity} we briefly describe the $f(R)$ gravity theory and its equations which will be used in the discussion below. In Sect.~\ref{sect:gas_physics} we give a more detailed account of the physics related to the gas fraction test described above. Then we present a numerical example in Sect.~\ref{sect:numerical_examples}, which shows how current data of cluster gas fraction alone can give powerful constraints on gravity. We discuss the results and their implications in Sect.~\ref{sect:con}.

\section{The $f(R)$ gravity theory}

\label{sect:fr_gravity}

This section is devoted to a quick overview of $f(R)$ gravity. It will be kept brief and only include essential equations, given that there are already many papers in the literature covering this topic.

$f(R)$ gravity \citep{cddatt2005} is a simple generalisation of GR, by replacing the Ricci scalar, $R$, in the Einstein-Hilbert action with an algebraic function $f(R)$
\begin{eqnarray}\label{eq:fr_action}
S = \int{\rm d}^4x\sqrt{-g}\left\{\frac{1}{2}M_{\rm Pl}^2\left[R+f(R)\right]+\mathcal{L}_{\rm m}\right\},
\end{eqnarray}
where $M_{\rm Pl}$ is the reduced Planck mass, $M_{\rm Pl}^{-2}=8\pi G$ with $G$ being Newton's constant, $g$ is the determinant of the metric $g_{\mu\nu}$ and $\mathcal{L}_{\rm m}$ is the Lagrangian density for (normal plus dark) matter fields. The model is defined by specifying the functional form of $f(R)$.

The action in Eq.~(\ref{eq:fr_action}) leads to a modified Einstein equation
\begin{eqnarray}\label{eq:fr_einstein}
G_{\mu\nu} + f_RR_{\mu\nu} -\left[\frac{1}{2}f-\Box f_R\right]g_{\mu\nu}-\nabla_\mu\nabla_\nu f_R = 8\pi GT^{\rm m}_{\mu\nu},
\end{eqnarray}
in which $G_{\mu\nu}$, $R_{\mu\nu}$ are respectively the Einstein and Ricci tensors, $f_R\equiv {\rm d} f/{\rm d} R$, $\nabla_{\mu}$ the covariant derivative compatible to the metric $g_{\mu\nu}$, $\Box\equiv\nabla^\alpha\nabla_\alpha$ and $T^{\rm m}_{\mu\nu}$ is the matter energy momentum tensor. Eq.~(\ref{eq:fr_einstein}) can be considered as the standard Einstein equation of GR with an extra scalar field, $f_R$, whose dynamics is governed by
\begin{eqnarray}
\Box f_R = \frac{1}{3}\left(R-f_RR+2f+8\pi G\rho_{\rm m}\right),
\end{eqnarray}
where $\rho_{\rm m}$ is the mass density of baryons and dark matter. As we are interested in late times, photons and neutrinos will be neglected.

On scales well inside the Hubble radius, and for the models to be considered, it is safe to work with the quasi-static approximation \citep{bose2015}, in which the scalar equation becomes
\begin{eqnarray}\label{eq:fr_eqn_static}
\vec{\nabla}^2f_R = -\frac{1}{3}a^2\left[R(f_R)-\bar{R} + 8\pi G\left(\rho_{\rm m}-\bar{\rho}_{\rm m}\right)\right],
\end{eqnarray}
where $\vec{\nabla}$ denotes the three dimensional gradient, $a$ is the scale factor, and an overbar takes the background value of a quantity. Notice that $R$ can be expressed as a function of $f_R$ by inverting $f_R(R)$.

Similarly, the modified Poisson equation in this limit reads as
\begin{eqnarray}\label{eq:poisson_static}
\vec{\nabla}^2\Phi = \frac{16\pi G}{3}a^2\left(\rho_{\rm m}-\bar{\rho}_{\rm m}\right) + \frac{1}{6}a^2\left[R\left(f_R\right)-\bar{R}\right],
\end{eqnarray}
where $\Phi$ is the Newtonian potential.

Eq.~(\ref{eq:fr_eqn_static}) implies two limits of the behaviour of $f(R)$ gravity:

(i) When $f_R$ is small, or more accurately, when $|f_R|\ll|\Phi|$, it recovers the well-known GR solution $R=8\pi G\rho_{\rm m}$, and so Eq.~(\ref{eq:poisson_static}) reduces to GR as well. This is the {\it chameleon} \citep{chameleon} regime which any viable $f(R)$ model must be in to pass the stringent solar system and terrestrial tests of gravity.

(ii) When $|f_R|\sim\mathcal{O}(|\Phi|)$, the second term on the rhs of Eq.~(\ref{eq:poisson_static}) is negligible compared with the first term, so that we have a gravity that is $1/3$ stronger than in GR. This is ususally known as the {\it non-chameleon}, or {\it unscreened}, regime.

It is evident that the unscreened regime mostly happens where $\Phi$ is shallow, or in extensive regions of low density. On large scales, matter density is close to the cosmological average, and so the total gravity in $f(R)$ gravity is enhanced within scales comparable to the Compton wavelength of the scalar field $f_R$ (which in most models of interest is in the range $\mathcal{O}(1\sim10)h^{-1}$Mpc). This naturally leads to an enhanced large-scale structure formation, and features such as over-abundant and more massive galaxy clusters - a topic which has been extensively studied previously. This will also be the topic that we focus on in this paper.

\section{Cluster gas fraction}

\label{sect:gas_physics}

Galaxy clusters are the largest bound objects in the Universe, whose masses are dominated by the dark matter component, with the baryonic masses dominated by X-ray emitting intracluster gas, which is heated to temperatures of the order of keV during virialisation. It is the mass fraction of this gas component that we will employ to test the theory of gravity here.

In this section, we shall first give a brief overview of how the baryon fraction can be estimated observationally, and how it can be used to constrain cosmological models and their parameters. Then, we will discuss how this process might be affected if the underlying theory of gravity is modified. For simplicity, we shall neglect other baryonic components than the intracluster gas in our analysis
unless otherwise stated.

\subsection{The standard $\Lambda$CDM model}

\label{sect:standard_model}

In the standard cosmological scenario, halo density profiles can be universally described by the \citet[][NFW]{nfw} fitting formula, which is often expressed as
\begin{eqnarray}
\rho(r) = \frac{\rho_{\rm s}}{(r/r_{\rm s})\left(1+r/r_{\rm s}\right)^2},
\end{eqnarray}
in which $\rho(r)$ is the halo mass density as a function of the distance, $r$, from the halo's centre, $\rho_{\rm s}$ is a characteristic density and $r_{\rm s}$ is the scale radius. We shall assume that the halo is spherically symmetric and well relaxed throughout the analysis, unless otherwise stated.

The mass of the halo can be obtained by integrating the NFW profile from $r=0$ to $r=R_{\Delta}$, in which $R_{\Delta}$ is the edge of the halo and is defined as the radius within which the average mass density is $\Delta\times\rho_{\rm crit}(z)$, with $\rho_{\rm crit}(z)\equiv3H(z)^2/8\pi G$ the critical density at the redshift $z$ when the halo is identified. This leads to
\begin{eqnarray}
M_{\rm halo} = 4\pi\rho_{\rm s}r^3_{\rm s}\left[\ln\left(1+R_\Delta/r_{\rm s}\right)-\frac{R_\Delta/r_{\rm s}}{1+R_\Delta/r_{\rm s}}\right].
\end{eqnarray}

Observationally, the total and baryonic masses of a cluster
can be obtained by measuring its X-ray surface brightness profile, and the temperature profile of its X-ray gas. For a dynamically relaxed system that consists of dark matter and baryonic gas, a hydrostatic equilibrium can be achieved, which satisfies the following equation
\begin{eqnarray}\label{eq:hydrostatic}
\frac{1}{\rho_{\rm gas}(r)}\frac{{\rm d}}{{\rm d}r}P_{\rm gas}(r) = -\frac{GM_{\rm tot}(<r)}{r^2},
\end{eqnarray}
in which $M_{\rm tot}(<r)$ is the total mass of dark matter and gas within radius $r$, and $\rho_{\rm gas}(r)$, $P_{\rm gas}(r)$ are respectively the density and pressure of the gas at $r$. For simplicity, we neglect  non-thermal pressure in our discussion (the effects of non-thermal pressure, however, are taken into account in the error budget when modelling the relation between $f_{\rm gas}$ and $\Omega_{\rm b}/\Omega_{\rm m}$, cf.~Eq.~(\ref{eq:allen_model}) below).

For an ideal thermal gas, its pressure and density are related to its temperature, $T_{\rm gas}$, as
\begin{eqnarray}\label{eq:eos}
P_{\rm gas} = kn_{\rm gas}T_{\rm gas} = \frac{k}{\mu m_{\rm p}}\rho_{\rm gas}T_{\rm gas},
\end{eqnarray}
where $k$ is the Boltzmann constant, and the mass and number densities of the gas particles are connected by $\rho_{\rm gas}=\mu m_{\rm p}n_{\rm gas}$, where $m_{\rm p}$ is the proton mass and $\mu$ the mean molecular weight. Applying Eqs.~(\ref{eq:hydrostatic}, \ref{eq:eos}), we obtain
\begin{eqnarray}\label{eq:hydrostatic_spherical}
\frac{GM(<r)}{r} = -\frac{kT_{\rm gas}(r)}{\mu m_{\rm p}}\left[\frac{{\rm d}\ln\rho_{\rm gas}(r)}{{\rm d}\ln r}+\frac{{\rm d}\ln T_{\rm gas}(r)}{{\rm d}\ln r}\right],
\end{eqnarray}
which, if evaluated at the halo edge ($r=R_\Delta$), gives
\begin{eqnarray}\label{eq:halo_mass}
\frac{GM_{\rm halo}}{R_\Delta} = -\frac{kT_{\rm gas}(R_\Delta)}{\mu m_{\rm p}}\left[\frac{{\rm d}\ln\rho_{\rm gas}}{{\rm d}\ln r}+\frac{{\rm d}\ln T_{\rm gas}}{{\rm d}\ln r}\right]_{r=R_\Delta}.
\end{eqnarray}

In the mean time, the X-ray emission of galaxy clusters is produced mainly by thermal bremsstrahlung radiation, leading to a X-ray surface brightness profile that also depends on the gas density and temperature profiles, $\rho_{\rm gas}(r)$ and $T_{\rm gas}(r)$. In terms of gas density, the gas mass can be expressed as
\begin{eqnarray}\label{eq:gas_mass}
M_{\rm gas}(<r) = \int^r_04\pi r'^2{\rm d}r'\rho_{\rm gas}(r') \propto \rho_{\rm gas}(0)r_0^3,
\end{eqnarray}
where $\rho_{\rm gas}(0)$ is the gas density at the cluster centre, and we have assumed that $\rho_{\rm gas}(r)=\rho_{\rm gas}(0)g(r/r_0)$ with $g(x)$ some function describing the profile and $r_0$ a characteristic scale. In the isothermal $\beta$ model, for example, we have
\begin{eqnarray}
\rho_{\rm gas}(r) = \rho_{\rm gas}(0)\left[1+\left(\frac{r}{r_0}\right)^2\right]^{-3\beta/2},
\end{eqnarray}
where $\beta$ is a dimensionless constant.

The bolometric luminosity is given by \citep[e.g.,][]{sasaki1996}
\begin{eqnarray}\label{eq:gas_luminosity}
L_X(<r) \propto \int^r_04\pi r'^2{\rm d}r'T^{1/2}_{\rm gas}(r')\rho_{\rm gas}^2(r') \propto \rho_{\rm gas}^2(0)r_0^3,
\end{eqnarray}
where we have neglected the proportionality coefficients which are irrelevant here.

Combining Eqs.~(\ref{eq:gas_mass}, \ref{eq:gas_luminosity}), we have
\begin{eqnarray}\label{eq:temp1}
M_{\rm gas} \propto L_X^{1/2}r_0^{3/2} \propto d_{\rm L}f^{1/2}_X\Theta_0^{3/2}d_{\rm A}^{3/2} \propto (1+z)^2d^{5/2}_{\rm A},
\end{eqnarray}
where $d_{\rm L}$ and $d_{\rm A}$ are respectively the luminosity and angular diameter distances and are related by $d_{\rm L}=(1+z)^2d_{\rm A}$. Here we have used $L_X=4\pi d_{\rm L}^2f_X$ where $f_X$ is the X-ray flux, and $r_0=d_{\rm A}\Theta_0$, in which $\Theta_0$ denotes the angle spanned by $r_0$ at redshift $z$. Note that $f_X$ (or equivalently the surface brightness) and $\Theta_0$ are the observed quantities in this description. In real observations, one has the surface brightness and temperature profiles, or equivalently $L_X(<r)$ and $T_{\rm gas}(r)$, using which Eqs.~(\ref{eq:hydrostatic_spherical}, \ref{eq:gas_luminosity}) can be solved simultaneously to find $M_{\rm halo}(<r)$ and $M_{\rm gas}(<r)$ ({note that the innermost regions of clusters are often excluded due to complicated processes such as cooling flow. For example, the study of \citet{allen2008} uses only the data within $0.7\sim1.2 R_{2500}$ to measure $f_{\rm gas}$}). These provide the necessary information to find
\begin{eqnarray}\label{eq:xxx}
f_{\rm gas}(r) = \frac{M_{\rm gas}(<r)}{M_{\rm halo}(<r)}.
\end{eqnarray}
Often in observations, people quote the value of $f_{\rm gas}$ at $r=R_{2500}$ as the cluster gas fraction.  We note in passing that here $M_{\rm halo}$ enters the picture only through its gravitational effect on gas particles, as this fact is important for the discussion of modified gravity below.

To see how this can be used to constrain background cosmology, we note that, once the right-hand side of Eq.~(\ref{eq:halo_mass}) is known by observations,
we also have a fixed numerical value of the left-hand side and hence have the relation
\begin{eqnarray}\label{eq:temp2}
M_{\rm halo} \propto R_\Delta \propto d_{\rm A}.
\end{eqnarray}
Eqs.~(\ref{eq:temp1}, \ref{eq:temp2}) imply that the {measured} cluster gas fraction depends on redshift and angular diameter distance in the following specific way:
\begin{eqnarray}\label{eq:apparent_gas_fraction}
f_{\rm gas}(z) = \frac{M_{\rm gas}}{M_{\rm halo}} \propto (1+z)^2d_{\rm A}^{3/2}.
\end{eqnarray}
As mentioned earlier, we expect the cluster gas fraction to be a reasonably fair sample of the mean cosmological baryon fraction and is therefore roughly independent of $z$ for massive clusters at low $z$. This is a reasonable assumption which should hold regardless of the cosmological model/parameters. The measurement of the apparent gas fraction, as described above, involves the angular diameter distance $d_{\rm A}$ and is indeed dependent on both the cosmological model and its parameters. Therefore, with incorrect cosmological models or parameters, the constancy of the true gas fraction (hereafter $f_{\rm gas}^{\ast}$) is not guaranteed to be reflected in the observed value $f^{\rm obs}_{\rm gas}(z)$. This provides a powerful test \citep{allen2004} of {\it background} cosmology and can be used to constrain cosmological parameters.

Were galaxy clusters perfectly fair samples of the average matter components in the Universe, their baryon fraction would just be $\Omega_{\rm b}/\Omega_{\rm m}$. The true situation, however, is more complicated. To take the complexities into account, \citet{allen2008}, improving on the earlier work of \citet{allen2004}, propose the following model of the relation between cluster gas fractions and $\Omega_{\rm b}/\Omega_{\rm m}$,
\begin{eqnarray}\label{eq:allen_model}
f^\ast_{\rm gas} = \frac{K\gamma b(z)}{1+s(z)}\frac{\Omega_{\rm b}}{\Omega_{\rm m}},
\end{eqnarray}
in which:

$\bullet$ $K$ is a constant accounting for systematic effects such as the calibration of instrument and X-ray modelling -- it is assumed to be $K=1.0\pm0.1$ in \citet{allen2008};

$\bullet$ $\gamma$ models the non-thermal pressure support in galaxy clusters which can cause a bias in the estimate of $f^\ast_{\rm gas}$ of about $9$\%;

$\bullet$ $b(z)\equiv b_0(1+\alpha_{\rm b}z)$ is the so-called 'depletion' factor which is inspired by the observation that the baryon fraction at $R_{2500}$ in non-radiative simulations \citep{eke1998} is actually smaller than $\Omega_{\rm b}/\Omega_{\rm m}$, with $b_0=0.83\pm0.04$ and $\alpha_{\rm b}$ small indicating a weak redshift evolution below $z=1$;

$\bullet$ $s(z)\equiv s_0(1+\alpha_{\rm s}z)$ accounts for the fact that a small fraction of baryons can be in the form of stars, with $s_0=0.16\pm0.048$ and $-0.2<\alpha_{\rm s}<0.2$ describing its redshift evolution.

The model in Eq.~(\ref{eq:allen_model}) indeed has a weak redshift dependency. However, any additional dependency from the observed $f_{\rm gas}$ would imply that one is using the wrong background cosmology to extract data, cf.~Eq.~(\ref{eq:apparent_gas_fraction}). This is, to be clear, in the framework of standard GR.

\begin{figure*}
 \includegraphics[width=180mm]{simulation.eps}
\caption{(Colour online) {\it Left panel}: the halo gas temperature profiles from the non-radiative hydrodynamical simulations of \citet{he2015b} for the standard $\Lambda$CDM model (black solid line) and $f(R)$ gravity with $|\bar{f}_{R0}|=10^{-5}$ (blue dotted line for standard haloes and red dashed line for effective haloes. The same line styles are used for the other two panels). {\it Middle panel}: The same as the left panel, but for the gas density profiles. The red crosses are a rescaling of the $\Lambda$CDM curve by $3/4$ to take into account the fact that for a $\Lambda$CDM halo and a $f(R)$ effective halo of the same mass, their true masses differ by $1/3$. {\it Right panel}: The halo gas fraction profiles from the same simulations; the red crosses are the $\Lambda$CDM result scaled by $3/4$, which is almost identical to the red dashed line, cf.~Eq.~(\ref{eq:observer_wrong}). The profiles shown here are stacked results of the haloes in the mass range $10^{13}\sim10^{13.4}h^{-1}M_\odot$, for illustration purpose; we've checked that results from other halo mass bins follow the same trends.}
\label{fig:simulation}
\end{figure*}

\subsection{Modified gravity scenarios}

In many modified gravity theories, including $f(R)$ gravity, the way in which the trajectories of massive test bodies -- e.g., galaxies, stars and gas particles -- respond to the underlying matter distribution is different. This change of the dynamics of test bodies is sometimes described as the change of the {\it dynamical mass} of matter. Massless particles, such as photons, could behave differently: in some theories, such as the Galileon model \citep{galileon1, galileon2} and the K-mouflage model \citep{kmouflage1, kmouflage2}, photons can also feel a different mass of matter, but in other models, for example $f(R)$ gravity and the \citet[][DGP]{dgp} model, photon trajectories depend on the matter distribution in essentially the same way as in GR, since the conformal coupling does not affect geodesics of massless particles. For distinction, the mass felt by photons is usually called the {\it lensing mass}. The differences in the dynamical and lensing masses of galaxy clusters have been used to constrain $f(R)$ gravity in, e.g., \citet[][]{terukina2014,wilcox2015}

As mentioned above, gas particles, like dark matter particles and galaxies, do feel the dynamical mass of a cluster. In $f(R)$ gravity, the same cluster can have a dynamical mass $4/3$ times its value in GR. This maximum enhance factor of 1/3, however, is not necessarily realised in all clusters, because of the chameleon screening \citep{chameleon}. The screening helps to reduce the difference between the dynamical and lensing masses, especially for more massive clusters. Consequently, constraints relying on the dynamical masses of clusters are in general weaker than those coming from astrophysical considerations. Nevertheless, they have cleaner physics than that of astrophysical observables -- which can often depend on whether the considered astrophysical system lives inside a screened cluster -- and are amongst the tightest constraints obtained using cosmological data \citep{terukina2014}.

Evidently, without a reliable measurement of the lensing mass, it is difficult to tell whether one observes a cluster of true mass $M$ in GR, or one with a smaller mass in $f(R)$ gravity, since both have the same dynamical masses. This has motivated \citet{he2015a} to propose the concept of {\it effective haloes}. Briefly speaking, the idea is to redefine the right-hand side of the modified Poisson equation in $f(R)$ gravity, Eq.~(\ref{eq:poisson_static}),
so that it can be rewritten as
\begin{eqnarray}\label{eq:poisson_eff}
\vec{\nabla}^2\Phi = 4\pi Ga^2\delta\rho_{\rm m,eff}.
\end{eqnarray}
In the above, $\Phi$ is the Newtonian potential that determines the dynamics of massive bodies and $\delta\rho_{\rm m}\equiv\rho_{\rm m}-\bar{\rho}_{\rm m}$ is the density perturbation of (dark plus baryonic) matter. $\delta\rho_{\rm m,eff}$ is the {\it effective density field}, with which the Poisson equation takes exactly the same form as in GR [cf.~Eq.~(\ref{eq:poisson_eff})]. In this way, the complicated new physics in $f(R)$ gravity is absorbed into $\rho_{\rm m,eff}$, and with that solved (e.g., in numerical simulations) one can in principle proceed assuming GR as the true theory of gravity. \citet{he2015b}, for example, show with hydrodynamical simulations that cluster gas temperatures depend only on the masses of the corresponding effective haloes, and that with certain rescaling depending on the effective halo mass the scaling relations -- such as the $L_X$-$M$ relation with $L_X$ the X-ray luminosity -- in $f(R)$ gravity can be derived reliably using existing knowledge of GR.

In the left panel of Fig.~\ref{fig:simulation}, we present the gas temperature profiles for standard haloes in GR and effective haloes in $f(R)$ gravity, both in the mass bin $10^{13}\sim10^{13.4}h^{-1}M_\odot$. Though there are differences in the inner regions -- which could be due to different halo density profiles or screening -- we notice that beyond $\sim100h^{-1}$kpc the two agree very well. \citet{he2015b} find that the average gas temperatures in the two also show very good agreement, and indeed the temperature-mass relation is barely distinguishable in the two models, provided that effective haloes are used in $f(R)$ gravity.

The fact that the cluster gas temperatures depend on the mass of the effective haloes is as expected, since for relaxed systems the virial temperature depends on the Newtonian potential, which does not distinguish between standard (GR) and effective ($f(R)$ gravity) density fields. For a polytropic gas with an equation of state $P_{\rm gas}\propto\rho_{\rm gas}^\Gamma$ in which the constant $\Gamma\geq1$, the hydrostatic equation implies that \citep[see, e.g.,][]{mvw-book} the temperature can be analytically expressed as a function of the potential $\Phi$.

Let us now consider two haloes, one identified in the standard dark matter field in a model with GR as the gravity theory, another from the effective density field in $f(R)$ gravity. The profiles of the two haloes are the same so that Eq.~(\ref{eq:poisson_eff}) sees no difference in them. Since gravity only enters the picture through Eq.~(\ref{eq:hydrostatic_spherical}), we make the following two observations/predictions:

$\bullet$ the gas temperature profiles are the same in these two haloes;

$\bullet$ the two sides of the spherical hydrostatic equation, Eq.~(\ref{eq:hydrostatic_spherical}), are the same for the two haloes.

\begin{figure*}
 \includegraphics[width=182mm]{obom.eps}
\caption{The inferred $\Omega_{\rm b}/\Omega_{\rm m}$ (triangles with 1$\sigma$ error bars) using the $f^{\rm obs}_{\rm gas}$ for the 42 clusters from \citet{allen2008}, as a function of the cluster redshift. 4 cases are shown, with different underlying  models of gravity: GR (upper left), and Hu-Sawicki $n=1$ $f(R)$ model with $|\bar{f}_{R0}|=3\times10^{-5}$ (upper right), $5\times10^{-5}$ (lower left), $10^{-4}$ (lower right). The horizontal solid and dotted lines are respectively the mean and 1$\sigma$ range of $\Omega_{\rm b}/\Omega_{\rm m}$ from Planck CMB data. A good match between $f_{\rm gas}$ and CMB for the GR case, and progressively worse matches for the $f(R)$ models, can be seen by a quick inspection by eye. For simplicity, all haloes are assumed to have a mass of $7.5\times10^{14}h^{-1}M_\odot$ and a concentration of $3.3$ when determining the effect of the chameleon screening.}
\label{fig:obom}
\end{figure*}

Since $T_{\rm gas}(r)$ and ${\rm d}\ln T_{\rm gas}(r)/{\rm d}\ln r$ in Eq.~(\ref{eq:hydrostatic_spherical}) are the same, we conclude that ${\rm d}\ln\rho_{\rm gas}(r)/{\rm d}\ln r$ is also the same in the haloes. This, however, does not necessarily mean that the haloes have identical gas density profiles, because we can rescale $\rho_{\rm gas}$ by a constant factor without changing ${\rm d}\ln\rho_{\rm gas}(r)/{\rm d}\ln r$. To confirm this, in the middle panel of Fig.~\ref{fig:simulation} we compare the gas density profile of $f(R)$ effective haloes with that of $\Lambda$CDM haloes of the same mass, and find that the two show a constant shift by 1/4 beyond $r\sim100h^{-1}$ kpc. Note that the simulated haloes in the plot do not have perfectly identical total -- standard or effective -- mass profiles, which is why in the middle panel of Fig.~\ref{fig:simulation} the red crosses and red dashed line do not agree on scales below $\sim100h^{-1}$kpc. However, as mentioned above Eq.~(\ref{eq:xxx}), in real observations, such innermost regions are not used in the determination of $f_{\rm gas}$ anyway.

As a result, to obtain the gas density profile, we need further, independent, information to fix its normalisation -- as opposed to its shape -- which brings us back to the measurements of cluster X-ray surface brightness. Inspecting the equation for the cluster X-ray luminosity, Eq.~(\ref{eq:gas_luminosity}), we notice that the luminosity density (i.e., the integrand) depends on (i) the {\it physical} gas density $\rho_{\rm gas}$, and (ii) the gas temperature $T_{\rm gas}$ which, as we have seen above, depends on the total mass of the {\it effective} halo. Consequently, should the physical gas densities be the same for $f(R)$ effective and $\Lambda$CDM haloes of the same mass, there would be no difference in their X-ray surface brightness profiles.

However, despite the standard (GR) and effective ($f(R)$ gravity) haloes above having identical gas temperature and halo mass, their {\it actual} (physical, or lensing) masses are different, and it is important to remember that $\Omega_{\rm m}$ characterises the amount of the {\it actual} mass in the Universe. If, as we have assumed so far, the gas fraction in clusters is a fair sample of the cosmological value, it would be the ratio of the gas mass and {\it actual} halo's mass that satisfies Eq.~(\ref{eq:allen_model}). Gas fractions inferred observationally, in the way described in the previous subsection, are in fact the ratio of the gas mass and that of the {\it effective} halo. If we denote the ratio of the effective and actual masses of a halo by $\eta$, then $1\leq\eta\leq4/3$, depending on the actual mass and environment of the halo, its redshift, as well as the $f(R)$ model parameters\footnote{More accurately speaking, effective haloes are identified from the effective density field, $\rho_{\rm m,eff}$ in Eq.~(\ref{eq:poisson_eff}), and standard haloes are identified from the physical density field $\rho_{\rm m}$. They do not necessarily share the same physical particles. Here, for simplicity, when talking about the effective and actual masses of some halo, we mean the masses of the effective and standard haloes that would be considered as matched haloes in the two catalogues.}. Here, as we are interested in the most massive clusters, with $M_{\rm halo}\sim10^{14}-10^{15}h^{-1}M_{\odot}$, we can for simplicity neglect the impacts of the halo's environment, so that $\eta$ mainly depends on $z$, i.e., $\eta=\eta(z)$, for a given halo mass.

Consider the extreme case in which $\eta=4/3$ as example. The apparent $f_{\rm gas}$ inferred from X-ray cluster observations would be
\begin{eqnarray}\label{eq:observer_wrong}
f^{\rm obs}_{\rm gas} = \frac{M_{\rm gas}}{M_{\rm halo, eff}} = \frac{3}{4}\frac{M_{\rm gas}}{M_{\rm halo,actual}} = \frac{3}{4}f^\ast_{\rm gas}.
\end{eqnarray}
This is also confirmed by hydro simulations, as shown in the right panel of Fig.~\ref{fig:simulation}. In that plot, we see that if haloes are defined using their actual masses in $f(R)$ gravity, then they share the same $f_{\rm gas}$ as $\Lambda$CDM haloes of the same mass (black solid vs. blue dotted lines). For effective $f(R)$ haloes, on the other hand, their $f_{\rm gas}$ profiles are a constant downward shift by $1/4$ from the results of $\Lambda$CDM haloes of the same masses (red dashed line vs.~red crosses), which is just what Eq.~(\ref{eq:observer_wrong}) predicts.

It is worthwhile to pause for a moment and try to understand the physics behind the behaviour of Fig.~\ref{fig:simulation}. It may seem surprising that, although the gas density profiles in $\Lambda$CDM and $f(R)$ standard haloes are significantly different within $r\sim100h^{-1}$kpc, their gas fraction profiles are very close to each other. This suggests that the two also have different dark matter (or total) mass profiles to cancel the differences in $\rho_{\rm gas}(r)$. In other words, $\rho_{\rm gas}(r)$ follows $\rho_{\rm DM}(r)$, which is the {\it physical} dark matter mass density, in the same way in $\Lambda$CDM and $f(R)$ standard haloes, even though they have different mass and even more different potential profiles. To understand this, we note that in Eq.~(\ref{eq:hydrostatic_spherical}) the $G$ on the left-hand side and $T_{\rm gas}$ on the right-hand side are both modified in $f(R)$ standard haloes. Indeed, if one assumes hydrostatic equilibrium and $P_{\rm gas}\propto\rho^\Gamma_{\rm gas}$, then the gas density profile can be written as \citep{komatsu2001}:
\begin{eqnarray}\label{eq:xxx}
\frac{\rho^{\Gamma-1}_{\rm gas}(r)}{\rho^{\Gamma-1}_{\rm gas}(0)} = 1- \frac{\Gamma-1}{\Gamma}\frac{G\mu m_{\rm p}M_{\rm halo}}{R_{\rm vir}kT_{\rm gas}(0)}\frac{c}{m(c)}\int^{\frac{r}{r_{\rm s}}}_0\frac{m(x)}{x^2}{\rm d}x,
\end{eqnarray}
in which $m(r/r_{\rm s})\equiv M(<r)/4\pi\rho_{\rm s}r^3_{\rm s}$. Eq.~(\ref{eq:xxx}) is derived under the assumption of self-similarity of gas density profiles, but the key point therein, that the modified gravity effects on $G$ and $T_{\rm gas}$ can be cancelled out, is not affected by this assumption. If this cancellation happens, then the gas density profile is determined by the total mass profile under the assumption of hydrostatic equilibrium, regardless of the theory of gravity\footnote{Note in Eq.~(\ref{eq:xxx}) it is $\rho_{\rm gas}(r)/\rho_{\rm gas}(0)$ that is determined by $m(r/r_{\rm s})$. The normalisation of $\rho_{\rm gas}$ will then be fixed by the total gas fraction inside the halo.}.

In clusters, gas is heated by accretion shocks during the assembly of the halo, a process which involves the conversions of energy from gravitational to kinetic (that of the cold accreted gas) and then to thermal (via shocks). Assuming a complete thermalisation, the post-shock gas temperature is proportional to $v_{\rm infall}^2$, with $v_{\rm infall}$ the infall speed of the accreted gas \citep[e.g.,][]{mvw-book}. Consequently, energy conservation implies that the final gas temperatures in the central regions will be affected in the same way as $v_{\rm infall}^2$ of the cold gas and hence $G$ in modified gravity.
Of course, this is only an approximation, and the cancellation of the effects of modified gravity on $G$ and $T_{\rm gas}$ depend on various factors including the screening and formation history of a cluster, which is not expected to be complete. However, Fig.~\ref{fig:simulation} suggests that it works pretty well for the haloes we use here. We checked explicitly that it works slightly less well for more massive haloes, for which the agreement between the $f_{\rm gas}(r)$ in $\Lambda$CDM and $f(R)$ standard haloes is slightly less perfect -- this may be because those haloes became unscreened only very recently.

The argument above in theory also applies to effective haloes, for which $G$ is the same as in GR, but the effects of modified gravity are incorporated in $M_{\rm halo}$. However, in the effective halo case the normalisation is different because of the different total gas fraction (see footnote 2) -- although the shape is the same -- hence the nearly constant rescaling of the red dashed curve compared with the black solid one in Fig.~\ref{fig:simulation}.

Coming back to the discussion prior to the previous three paragraphs, our result suggests two possible tests of $f(R)$ gravity:

(i) If an observer actually lives in a universe shaped by $f(R)$ gravity, then the true cluster gas fraction is given by $f^{\ast}_{\rm gas}=\eta f^{\rm obs}_{\rm gas}$. Assuming that Eq.~(\ref{eq:allen_model}) still holds for $f^\ast_{\rm gas}$, the observer will need to do the following transformation to get the true $\Omega_{\rm b}/\Omega_{\rm m}$:
\begin{eqnarray}
f^{\rm obs}_{\rm gas} = \frac{1}{\eta}\frac{K\gamma b(z)}{1+s(z)}\left[\frac{\Omega_{\rm b}}{\Omega_{\rm m}}\right]_{\rm true} \Rightarrow \left[\frac{\Omega_{\rm b}}{\Omega_{\rm m}}\right]_{\rm true} = \eta\left[\frac{\Omega_{\rm b}}{\Omega_{\rm m}}\right]_{\rm obs}.
\end{eqnarray}
As $(\Omega_{\rm b}/\Omega_{\rm m})_{\rm obs}$ depends only on the actual observational data, the observer will obtain the same value as an observer in a standard GR universe would do. The resulting $(\Omega_{\rm b}/\Omega_{\rm m})_{\rm true}$ might then be too large to be compatible with other constraints, such as the one from the CMB.

(ii) Alternatively, if one takes the $\Omega_{\rm b}/\Omega_{\rm m}$ measured by other probes as the true value and starts from there, then Eq.~(\ref{eq:observer_wrong}) implies that the observed cluster gas fraction $f^{\rm obs}_{\rm gas}$ will be smaller than what the $\Lambda$CDM model and simulations predict. Because of the time dependence of $\eta(z)$ (see above), if the $f(R)$ model parameters happen to take the values for $\eta$ to evolve from 1 to $4/3$ between $z=1$ and the present for the clusters of interest, there may also be an apparent decrease of $f^{\rm obs}_{\rm gas}(z)$ as $z$ decreases, by a maximum of $25\%$.

\section{Numerical Examples}

\label{sect:numerical_examples}

In this section, we use a simplified example to illustrate the power of the cluster gas fraction test proposed above. For this, we will use the gas fraction data of the 42 clusters studied by \citet[][Table 3]{allen2008}. As described above, these $f_{\rm gas}$ data are obtained by fitting the gas temperature and X-ray surface brightness profiles of these clusters simultaneously, assuming NFW profiles for the total mass in clusters. We have also found that, in the context of $f(R)$ gravity, as long as we use effective haloes, the dynamics of gas particles can be calculated using standard gravity theory. Therefore, in this work {\it we can directly take the data of} \citet{allen2008} {\it as $f^{\rm obs}_{\rm gas}$}, bearing in mind that the cluster mass inferred therein would be the effective mass and therefore $f^{\rm obs}_{\rm gas}$ can be different from $f^{\ast}_{\rm gas}$ for unscreened clusters, cf.~Eq.~(\ref{eq:observer_wrong}).

To obtain an estimation of the mean and standard deviation of $\left(\Omega_{\rm b}/\Omega_{\rm m}\right)_{\rm obs} = \frac{1+s(z)}{K\gamma b(z)}f^{\rm obs}_{\rm gas}$ from each cluster, random samples of size $10^5$ are drawn for each parameter or data: $K, \gamma, b_0, \alpha_{\rm b}, s_0, \alpha_{\rm s}$ and $f^{\rm obs}_{\rm gas}$. Of these, $f^{\rm obs}_{\rm gas}$ is taken, for a given cluster, from Table 3 of \citet{allen2008}, and is assumed to satisfy a Gaussian distribution with mean and standard deviation given by \citet{allen2008}. The other parameters and their distributions are shown in Table~\ref{tab:parameters}. We therefore obtain $10^5$ realisations of $\left(\Omega_{\rm b}/\Omega_{\rm m}\right)_{\rm obs}$, from which its mean and standard deviation can be calculated. This procedure is repeated for all 42 clusters.

\begin{table}
\caption{The assumed ranges and distributions of the model parameters in Eq.~(\ref{eq:allen_model}). For more details the readers are referred to Sect.~\ref{sect:standard_model} or \citet{allen2008}.}
\begin{tabular}{@{}llll}
\hline\hline
param & physical effect described & mean $\pm$ stddev & prior \\
\hline
$K$ & overall calibration & $1.000\pm0.100$ & Gaussian \\
$\gamma$ & non-thermal pressure & $1.050\pm0.050$ & Uniform \\
$b_0$ & gas bias: normalisation & $0.825\pm0.175$ & Uniform \\
$\alpha_{\rm b}$ & gas bias: evolution & $0.000\pm0.100$ & Uniform \\
$s_0$ & stellar fraction: normalisation & $0.160\pm0.048$ & Gaussian \\
$\alpha_{\rm s}$ & stellar fraction: evolution & $0.000\pm0.200$ & Uniform \\
\hline
\end{tabular}
\label{tab:parameters}
\end{table}

The estimation of the effects of modified gravity, i.e., the factor $\eta(z)$, is more complicated, since it depends on the cluster mass, density profile, environment, redshift, as well as the $f(R)$ parameters. Because the main purpose of this paper is to illustrate the basic idea, we shall leave a full analysis using real cluster data for future work, and instead adopt a simplified modelling. The cluster masses are assumed to be the same, with $M_{\rm halo}=7.5\times10^{14}h^{-1}M_\odot$, for all 42 clusters, since this is a typical value for massive X-ray clusters. The cluster's radius ($R_{200}$ or $R_{\rm vir}$) is taken to be $1.5h^{-1}$Mpc and its concentration parameter, $c\equiv R_{200}/r_{\rm s}$ (or $R_{\rm vir}/r_{\rm s}$) is $3.3$. The cluster is assumed to live on the cosmological background, so that the ratio of its effective and actual masses can be approximated as \citep[see, e.g.,][]{lzk2012}
\begin{eqnarray}\label{eq:etaz}
\eta(z)=\frac{M_{\rm halo,eff}}{M_{\rm halo,actual}} = \min\left\{1+\frac{\bar{f}_R(z)}{2\Phi_{\rm N}},\frac{4}{3}\right\}
\end{eqnarray}
where $\bar{f}_R(z)$ is the background value of $f_R$ at redshift $z$, and
\begin{eqnarray}
\Phi_{\rm N} = -\frac{GM_{\rm halo}}{R_{\rm halo}}\frac{\ln(1+c)}{\ln(1+c)-c/(1+c)} \approx -5\times10^{-5},
\end{eqnarray}
is the Newtonian potential at the edge of an NFW halo. We adopt the $f(R)$ model by \citet{hs2007} with $n=1$, for which
\begin{eqnarray}
\bar{f}_R(z) = \left[\frac{\Omega_m+4\Omega_\Lambda}{\Omega_m(1+z)^3+4\Omega_\Lambda}\right]^2\bar{f}_{R0}.
\end{eqnarray}
The values $\Omega_m=0.316$ and $\Omega_\Lambda=0.684$ are taken from the latest results of \citet{planck2015}. As a result, the physics of modified gravity is completely governed by $\bar{f}_{R0}$, which is the present-day value of $\bar{f}_R$. Once this is specified, we can obtain $\eta(z)$, and therefore infer $\left(\Omega_{\rm b}/\Omega_{\rm m}\right)_{\rm true}$ given $\left(\Omega_{\rm b}/\Omega_{\rm m}\right)_{\rm obs}$ and $z$ of a cluster.

In Fig.~\ref{fig:obom} we show the $\left(\Omega_{\rm b}/\Omega_{\rm m}\right)_{\rm true}$ result obtained from $f^{\rm obs}_{\rm gas}$ for 4 different cases: standard $\Lambda$CDM (upper left), and $f(R)$ gravity models with $|\bar{f}_{R0}|=3\times10^{-5}$ (upper right), $5\times10^{-5}$ (lower left) and $10^{-4}$ (lower right). For comparison, we have also, in each panel, plotted the mean value (solid line) and 1$\sigma$ confidence level (dotted) of $\Omega_{\rm b}/\Omega_{\rm m}$ from \citet{planck2015}\footnote{Note that the $f(R)$ models studied here have practically identical CMB power spectra as the $\Lambda$CDM model with the same $\Omega_{\rm m}$ and $\Omega_\Lambda$. As a result, using CMB data only, the constraints on cosmological parameters such as $\Omega_{\rm m}$ would be the same in all these models. Because of this, the CMB constraints are less model-dependent.}. The results from the 42 clusters, with 1$\sigma$ errors, are shown as symbols.

A quick naked-eye inspection shows that the $f_{\rm gas}$ method and the CMB observation give compatible $\Omega_{\rm b}/\Omega_{\rm m}$ if one assumes the $\Lambda$CDM paradigm (upper left). The $f(R)$ model with $|\bar{f}_{R0}|=10^{-4}$ (lower right), on the other hand, leads to a significantly higher value of $\Omega_{\rm b}/\Omega_{\rm m}$ than what CMB says, and is therefore inconsistent. The other two cases are more interesting: for $|\bar{f}_{R0}|=5\times10^{-5}$ (lower left), $\eta(z)$ increases to $4/3$ at $z\sim0.45$, while for $|\bar{f}_{R0}|=3\times10^{-5}$ (upper right) $\eta(z)$ only increases to $\sim1.28$ at $z=0$. In both cases, however, the inferred values of $\Omega_{\rm b}/\Omega_{\rm m}$ are still substantially larger than the Planck result, especially for the low-$z$ clusters. This shows that cluster gas fraction can be a potentially powerful test of gravity, using X-ray observations only. In such tests, lensing data can be a useful addition, but is not necessary.

\begin{figure}
 \includegraphics[width=88mm]{fgas_evol.eps}
\caption{(Colour online) The evolution of $f_{\rm gas}$ in 4 models with the same cosmic value of $\Omega_{\rm b}/\Omega_{\rm m}$ from \citet{planck2015}. The models are, respectively, $\Lambda$CDM (black solid line), and Hu-Sawicki $n=1$ $f(R)$ model with $|\bar{f}_{R0}|=3\times10^{-5}$ (blue dotted line), $5\times10^{-5}$ (green dashed line) and $10^{-4}$ (red dot-dashed line). Black triangles with error bars are the $f_{\rm gas}$ values of the 42 clusters used in \citet[][Tab 3]{allen2008}. The blue shaded region denotes the standard deviation around the mean $f_{\rm gas}$ for the model with $|\bar{f}_{R0}|=3\times10^{-5}$, for illustration, which shows that the theoretical uncertainty is roughly of the same order as current observational errors in $f_{\rm gas}$; therefore, the constraining power can be further improved if either of these uncertainties is reduced in the future.}
\label{fig:fgas_evol}
\end{figure}

The test can be done in an alternative way. For this, we assume the value of $\Omega_{\rm b}/\Omega_{\rm m}$ obtained from CMB observations, and check what value of $f^{\rm obs}_{\rm gas}$ an observer would have found  if living in a $f(R)$ universe. The idea is that, if this value differs too much from what {\it our} observers have told us \citep[e.g., in][]{allen2008}, then it would place a constraint on the extent to which the assumed $f(R)$ model can deviate from standard $\Lambda$CDM. As in the previous case, we draw random samples of size $10^5$ for the parameters $K, \gamma, b_0, \alpha_{\rm b}, s_0, \alpha_{\rm s}$ from which we find $10^5$ realisations of $K\gamma b(z)/(1+s(z))$. Then, by modelling modified gravity effects using Eq.~(\ref{eq:etaz}), we compute the mean $f_{\rm gas}$ for the 4 models shown in Fig.~\ref{fig:obom}, and these are shown as curves in Fig.~\ref{fig:fgas_evol} together with the observed values of $f_{\rm gas}$ from \citet{allen2008}. Again, we note that current data favour $\Lambda$CDM over all three variants of $f(R)$ gravity. As clusters are less screened at late times, we find that low-$z$ data is more useful in constraining the model than high-$z$ data.

\section{Discussion and Conclusions}

\label{sect:con}

In this paper, we proposed a new cosmological test of gravity, by inferring the cosmic baryon fraction from the apparent gas fractions of massive clusters, and comparing with the results from other, less model-dependent, measurements such as the CMB. In theories with a stronger gravity, the apparent gas fraction is smaller than that in $\Lambda$CDM for a fixed $\Omega_{\rm b}/\Omega_{\rm m}$. Reversely, if the observed value $f^{\rm obs}_{\rm gas}$ is fixed, we would find a higher value of $\Omega_{\rm b}/\Omega_{\rm m}$ than in GR, that can be inconsistent with the model-independent measurements. Taking the Hu-Sawicki $f(R)$ model as an example: our quick calculation shows that model parameters $|\bar{f}_{R0}|\sim3-5\times10^{-5}$ are in tension with the gas fraction data of the 42 clusters from \citet{allen2008}, though a more rigorous constraint will be left for future work.

$f_{\rm gas}$ has been a rather widly-used observable \citep[e.g.,][]{white1993}, and its power in constraining cosmology -- in particular dark energy models -- is convincingly demonstrated in various previous works \citep[e.g.,][]{sasaki1996,allen2004,allen2008}. The inclusion of baryons opens a new dimension for tests of gravity, since ultimately most cosmological observables can be tracked back to lights emitted by interactions involving baryons. In the mean time, the physics of the X-ray-emitting hot gas in massive clusters is relatively clean, making it easier both for the modelling and to use the observational data. As an example, the assumption that gas temperature depends on the gravitational potential and our main conclusions that (i) $f^\ast_{\rm gas}$ -- the true gas fraction -- is {\it unchanged} with modified gravity while (ii) $f^{\rm obs}_{\rm gas}$ is {\it changed} are supported by hydrodynamical simulations in $f(R)$ gravity \citep[e.g.,][]{he2015b}. Some uncertainties remain in relating $f_{\rm gas}$ to $\Omega_{\rm b}/\Omega_{\rm m}$, but these have been included in the error budget estimate above. Furthermore, within our current state of understanding, slightly changing its modelling (e.g., from \citet{allen2004} to \citet{allen2008}) does not change results drastically.

Here, we would like to emphasise the use of effective haloes \citep{he2015a} in our analysis. Though the idea has a similar origin as that of the dynamical mass of halo \citep[e.g.,][]{schmidt2010}, there are fundamental differences. Dynamical mass is a certain attribute of a given halo which is defined in the standard way, while effective halo is a completely new way to define and identify haloes. Given an effective halo, all gravitational effect can be calculated from GR, and in particular this means that the way in which $f_{\rm gas}$ is currently extracted from observational data -- and the resulting $f_{\rm gas}$ results -- can be directly used for our purpose. Thus, with a little extra effort from the people who generate a halo catalogue, the analyses of end users can be made much more straightforward, and this provides an efficient bridge between simulators, theorists and observers.

One may naturally wonder about the generality of this method. As a cosmological test, it relies on galaxy clusters being totally or partially unscreened. Because we are talking about massive clusters which tend to be better screened, this test, {like most other cosmological ones}, will probably not be able to constrain $\bar{f}_{R0}$ to substantially smaller than the quoted values here. However, it does provide a fairly clean test -- with good observational data available -- that has the potential to place one of the strongest constraints from cosmology on $f(R)$ gravity. Also, one can always combine $f_{\rm gas}$ and other observables, such as lensing \citep{terukina2014}, cluster scaling relations \citep[e.g.][]{arnold2014}, and cluster gas pressure profiles \citep[][]{demartino2014}, to place joint, and likely stronger, constraints. In principle, the test would be more powerful if observational data for smaller galaxy clusters (e.g., those in the mass range $10^{13}\sim10^{14}h^{-1}M_\odot$) and galaxy groups are included, because these objects are less screened and so gravity deviates more from GR in general. This, however, requires a better understanding of the feedbacks in different models, which are not well studied so far.

The test can be applied not only to $f(R)$ gravity and the more general chameleon theory, but also to similar models such as dilatons \citep{dilaton} and symmetrons \citep{symmetron}. These models are all featured by a universal coupling of all matter species to a scalar field, that effectively enhances the gravity for all particles (at least in unscreened regimes). There are models in which only certain matter species, e.g.. dark matter, experiences the scalar coupling: therein, baryons can still feel a different gravity depending on how the dark matter particle mass evolves with time, in which case the proposed test does apply. In addition to these theories, we have mentioned above the DGP, Galileon and k-mouflage models. In the first two classes, the deviations from GR are strongly suppressed inside dark matter haloes \citep[][]{barreira2013,barreira2014b}, and so we do not expect the new test to work. For K-mouflage, as is for the so-called non-local gravity \citep{nonlocal1,nonlocal2,barreira2014a}, there can be a time evolution of Newton's constant inside clusters, making it possible to use this test. However, one needs to bear in mind that in many of these theories the background evolution history is also modified, and that can affect the $f_{\rm gas}$ test (whether it leads to degeneracies or stronger constraints can only be told by a case-by-case study in the future).

As mentioned earlier, the aim of this paper is to illustrate the main idea of using $f_{\rm gas}$ as a test  of gravity theories, and therefore we have made a simplified estimate and have not quoted any numerical results on the confidence levels of the constrained $f_{R0}$. A more complete and rigorous analysis will require one to relax the simplification that all observed clusters share the same mass, radius and concentration, and use the real observational results of these for all clusters. If the cluster mass is obtained from its dynamical effects, we also need to account for the fact that different clusters may have experienced different degrees of screening, and so a more accurate modelling of the screening is needed to compute the cluster mass profile. These will be left for future work. We note that hydrodynamical simulations for modified gravity theories started to appear recently \citep[e.g.,][]{arnold2014,hammami2015,he2015b}, and such works will be useful for improving the constraining power of this test in the future.

\section*{Acknowledgments}

We thank Sownak Bose, Vince Eke, Claudio Llinares and Gongbo Zhao for discussions and comments. The work has used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.  BL acknowledges support by the UK STFC Consolidated Grant No.~ST/L00075X/1 and No.~RF040335. JHH is supported by the Italian Space Agency (ASI) via contract agreement I/023/12/0. LG acknowledges support from {NSFC} grants Nos. 11133003 and 11425312, the Strategic Priority Research Program {\it The Emergence of Cosmological Structure}  of the Chinese Academy of Sciences (No.~XDB09000000), {MPG} partner Group family, and an {STFC} and Newton Advanced Fellowship.