# Fully Connected Reconfigurable Intelligent Surface Aided Rate-Splitting
Multiple Access for Multi-User Multi-Antenna Transmission

Tianyu Fang∗†‡ , Yijie Mao∗, Shanpu Shen§, Zhencai Zhu‡, Bruno Clerckx¶
This work was sponsored by Shanghai Sailing Program under Grant 22YF1428400.
∗School of Information Science and Technology, ShanghaiTech University,
Shanghai, China
†University of Chinese Academy of Sciences, Beijing, China
‡Innovation Academy for Microsatellites, Chinese Academy of Sciences,
Shanghai, China
§Department of Electronic and Computer Engineering, The Hong Kong University
of Science and Technology, Hong Kong
¶Department of Electrical and Electronic Engineering, Imperial College London,
United Kingdom
Email: {fangty<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Rate-splitting multiple access (RSMA) has been recognized as a promising and
powerful multiple access (MA) scheme, non-orthogonal transmission framework
and interference management strategy for 6G. Inspired by the appealing
spectral efficiency gain achieved by RSMA over conventional MA schemes in
multi-user multi-antenna transmission, in this paper we introduce RSMA to
reconfigurable intelligent surface (RIS)-aided multiple-input single-out
(MISO) broadcast channel (BC). To further enhance the spectral efficiency, a
more generalized RIS architecture called fully connected RIS is considered. By
jointly optimizing the scattering matrix of the fully connected RIS and the
transmit beamformers to maximize the sum-rate, we show that the proposed fully
connected RIS aided RSMA transmission scheme significantly improves the
spectral efficiency compared with the conventional single connected RIS
schemes and the schemes without RIS. It acts as a new benchmark for linearly
precoded multi-user multi-antenna networks.

###### Index Terms:

fully connected, multiple-user multi-antenna network, rate-splitting multiple
access, reconfigurable intelligent surface, spectral efficiency.

##  I Introduction

The past few years have witnessed the development of rate-splitting multiple
access (RSMA) for multi-antenna networks. It has been recognized as a
promising physical-layer non-orthogonal transmission strategy, a powerful
interference management approach, and a candidate of the multiple access
technique for the sixth generation wireless network (6G) [1, 2, 3]. By
splitting the user messages into common and private parts, encoding the common
parts into the common streams to be decoded by multiple users, and encoding
the private parts respectively into the private streams to be decoded by the
corresponding users only, RSMA enables a flexible interference management
capability of partially decoding the interference and partially treating the
interference as noise [4, 5]. Existing works have shown that RSMA outperforms
other multiple access techniques (including linearly precoded space division
multiple access–SDMA, power-domain non-orthogonal multiple access–NOMA,
orthogonal multiple access-OMA, and multicasting) in terms of spectral
efficiency [4, 6, 5, 7, 8], energy efficiency [9, 7], max-min fairness [8,
10], and robustness towards inaccuracies of channel state information at the
transmitter (CSIT) [6, 11].

Reconfigurable intelligent surfaces (RISs), consisting of a large number of
reconfigurable scattering elements, is another promising technology for 6G
[12]. RISs can smartly reconfigurable the wireless propagation environment so
as to effectively enhance the spectral and energy efficiency. Most of existing
RIS research relies on using a simple architecture referred to as single
connected RIS, which is characterized by a diagonal matrix with constant
modulus entries [13, 14, 15, 12]. To enhance the RIS performance, a more
general architecture referred to as fully connected RIS has recently been
proposed in [16]. The fully connected RIS is characterized by a complex
symmetric unitary matrix, which achieves a better performance than the single
connected RIS.

The appealing performance benefits of RSMA and RIS have motivated the study on
the integration of them [17, 18, 19, 20, 21]. The interplay of RSMA and RIS is
first investigated in [17] for multi-user multi-antenna networks, where RIS
aided RSMA transmission model achieves higher energy efficiency than NOMA and
orthogonal frequency division multiple access (OFDMA). RIS aided RSMA
transmission is further shown to enhance the fairness among users [18], reduce
the transmit power [19], and achieve superior outage performance [20, 21] over
conventional RIS-aided networks. However, all the above works only use the
single connected RIS architecture. Therefore, we would like to further enhance
the performance by leveraging the more advanced fully connected RIS
architecture. To the best of our knowledge, there is no existing work that
investigates the sum-rate achieved by the fully connected RIS aided RSMA
transmission networks.

In this work, we propose a fully connected RIS aided RSMA downlink multi-
antenna multi-user transmission model. The scattering matrix of RIS and the
transmit beamformers of RSMA are jointly optimized in order to maximize the
sum-rate. To solve the problem, we propose an alternative optimization
framework where the scattering matrix of RIS and the beamformers of RSMA are
iteratively optimized. Numerical results show that by synergizing RSMA and
fully connected RIS, the proposed scheme significantly improves the spectral
efficiency in multi-user multiple-input single-output (MU-MISO) MU-MISO
networks. The proposed scheme explores a larger achievable sum-rate than the
conventional single connected RIS aided schemes and the schemes without RIS.
It acts as a new benchmark for linearly precoded multi-user multi-antenna
networks.

Notations: Vectors and matrices are denoted by bold lower and upper case
letters. $\mathbb{R}^{m\times n}$ and $\mathbb{C}^{m\times n}$ represent the
real-valued and complex-valued spaces with dimension $m\times n$. $|x|$
indicates the magnitude of a complex number $x$ and $\Re(x)$ denotes its real
part. For a vector $\mathbf{x}$, $\|\mathbf{x}\|$ denotes its Euclidean norm.
$\mathbb{E}\\{\cdot\\}$ denotes the statistical expectation operator for a
random variable. $(\cdot)^{H}$, $(\cdot)^{T}$ and $\mathrm{tr}(\cdot)$
respectively denote the conjugate transpose, transpose and trace operators.
$\mathrm{diag}\left(x_{1},x_{2},\cdots,x_{n}\right)$ is a diagonal matrix with
$(x_{1},x_{2},\cdots,x_{n})$ being its diagonal elements. $\mathbf{I}$ is the
identity matrix. $\mathcal{CN}(\mu,\sigma^{2})$ denotes the circularly
symmetric complex Gaussian (CSCG) distribution with mean $\mu$ and variance
$\sigma^{2}$.

##  II System Model and Problem Formulation

We consider a MU-MISO communication network consisting of one base station
(BS) equipped with $M$ antennas, one RIS with a set of $N$ passive reflecting
elements indexed by $\mathcal{N}=\\{1,2,\ldots,N\\}$, and a set of $K$ single-
antenna users indexed by $\mathcal{K}=\\{1,2,\ldots,K\\}$. As illustrated in
Fig. 1, the BS simultaneously serves the $K$ users with the assistance of a
fully-connected RIS (as in Fig. 1 (a)). The reconfigurable impedance network
of RIS is adjusted and determined by a smart controller attached to the RIS,
which also acts as a gateway to exchange the information between the BS and
the RIS. The channels from the BS to the users, from the RIS to the users, and
from the BS to the RIS are denoted as $\mathbf{g}_{k}\in\mathbb{C}^{M\times
1}$, $\mathbf{h}_{k}\in\mathbb{C}^{N\times 1}$, $k\in\mathcal{K}$ and
$\mathbf{G}\in\mathbb{C}^{N\times M}$, respectively. All channels are assumed
to be invariant during one transmission block and perfect CSI is available at
the BS. Although the assumption of perfect CSI is ideal, the proposed scheme
explores a larger achievable sum-rate than the conventional schemes, which
therefore acts as a new benchmark for multi-user multi-antenna networks as
well as future study for the corresponding imperfect CSIT settings. At the BS,
message $W_{k}$ intends to user $k$ is split into a common part $W_{c,k}$ and
a private part $W_{p,k}$. The common parts of all users are combined and
encoded into a common stream $s_{0}$ while the private parts are independently
encoded into the private streams $s_{1},\cdots,s_{K}$. Denote
$\mathbf{s}=[s_{0},s_{1},\cdots,s_{K}]^{T}$ and
$\mathbf{W}=[\mathbf{w}_{0},\mathbf{w}_{1},\cdots,\mathbf{w}_{K}]\in\mathbb{C}^{M\times(K+1)}$
as the data stream vector and beamforming matrix for all streams,
respectively. We assume that each stream $s_{k},$ $k\in\mathcal{K}\cup\\{0\\}$
has zero mean and unit variance, i.e.,
$\mathbb{E}\\{\mathbf{s}\mathbf{s}^{H}\\}=\mathbf{I}$. The transmitted signal
at the BS is

$\mathbf{x}=\sum\limits_{k=0}^{K}\mathbf{w}_{k}s_{k},$ (1)

and the transmit power constraint is

$\mathrm{tr}(\mathbf{WW}^{H})\leq P_{t},$ (2)

where $P_{t}$ refers to the maximum transmit power of the BS. The signal is
transmitted through the direct signal path from the BS to the users as well as
the RIS-aided path. At user $k$, the total received signal is

$y_{k}=(\mathbf{g}_{k}^{H}+\mathbf{h}_{k}^{H}\mathbf{\Theta}\mathbf{G})\sum\limits_{i=0}^{K}\mathbf{w}_{i}s_{i}+z_{k},$
(3)

where $\mathbf{\Theta}\in\mathbb{C}^{N\times N}$ refers to the scattering
matrix of the $N$-port reconfigurable impedance network in the $N$-element
RIS, and $z_{k}\sim\mathcal{CN}(0,\sigma_{k}^{2})$ is the additive white
Gaussian noise (AWGN).

Figure 1: A multi-antenna multi-user transmission network with the assistance
of a (a) fully connected RIS, (b) single connected RIS.

###  II-A Fully Connected Reconfigurable Intelligent Surface

In this work, we focus on using a fully connected RIS [16] to enhance the
spectral efficiency. An example of a 4-element fully connected RIS is
illustrated in Fig. 1 (a). In the reconfigurable impedance network of fully
connected RIS, each port is connected with other ports through a
reconfigurable reactance. Accordingly, the scattering matrix of a fully
connected RIS $\mathbf{\Theta}$ satisfies the constraints

$\displaystyle\mathbf{\Theta}^{H}\bm{\Theta}=\mathbf{I},$ (4a)
$\displaystyle\mathbf{\Theta}=\mathbf{\Theta}^{T}.$ (4b)

As per [16, 22], constraint (4) is equivalent to

$\displaystyle\mathbf{\Theta}$
$\displaystyle=(j\mathbf{X}+Z_{0}\mathbf{I})^{-1}(j\mathbf{X}-Z_{0}\mathbf{I}),$
(5a) $\displaystyle\mathbf{X}$ $\displaystyle=\mathbf{X}^{T},$ (5b)

where $Z_{0}$ refers to the reference impedance and $\mathbf{X}$ is a
symmetric real matrix referring to the reactance matrix of the reconfigurable
impedance network in RIS. With constraint (5), a closed-form expression for
scattering matrix $\mathbf{\Theta}$ satisfying constraint (4) is obtained by
introducing an unconstrained symmetrical real matrix $\mathbf{X}$.

Remark 1. When each port is disconnected with other ports in the
reconfigurable impedance network, the fully connected RIS reduces to the
single connected RIS [23] as illustrated in Fig. 1 (b). The single connected
RIS has been widely used in existing works [17, 18, 20, 19, 21], where the
scattering matrix $\mathbf{\Theta}$ satisfies the constraint

$\bm{\Theta}=\mathrm{diag}\left(e^{j\theta_{1}},e^{j\theta_{2}},\cdots,e^{j\theta_{N}}\right),$
(6)

where $\theta_{n}\in[0,2\pi)$ denotes the phase of the scattering parameter of
the $n$-th port in reconfigurable impedance network. Accordingly, it can be
also equivalently transformed to

$\displaystyle\bm{\Theta}$
$\displaystyle=(j\mathbf{X}+Z_{0}\mathbf{I})^{-1}(j\mathbf{X}-Z_{0}\mathbf{I}),$
(7a) $\displaystyle\mathbf{X}$
$\displaystyle=\mathrm{diag}\left(x_{1},x_{2},\cdots,x_{N}\right),$ (7b)

where $x_{n}\in\mathbb{R}$ is the reconfigurable reactance component connected
to the $n$-th port.

It should be noted that a $N$-port fully connected RIS given in (5) requires
to tune $N(N+1)/2$ scattering parameters while a $N$-port single connected RIS
given in (7) only requires to tune $N$ scattering parameters. The fully
connected RIS therefore brings a larger searching space for the optimal RIS
design.

###  II-B Problem Formulation

At user sides, each user first decodes the common stream by treating all
private streams as interference. Thus, the signal-to-interference-plus-noise
ratio (SINR) of $s_{0}$ at user $k$ is

$\gamma_{0,k}=\frac{|(\mathbf{g}_{k}^{H}+\mathbf{h}_{k}^{H}\mathbf{\Theta}\mathbf{G})\mathbf{w}_{0}|^{2}}{\sum\limits_{i=1}^{K}|(\mathbf{g}_{k}^{H}+\mathbf{h}_{k}^{H}\mathbf{\Theta}\mathbf{G})\mathbf{w}_{i}|^{2}+\sigma_{k}^{2}},$
(8)

and the corresponding transmission rate is
$r_{0,k}=\log_{2}\left(1+\gamma_{0,k}\right)$. To ensure common message
$s_{0}$ is successfully decoded by all users, the achievable rate of $s_{0}$
should satisfy $r_{0}=\min_{k\in\mathcal{K}}r_{0,k}$. After decoding the
common stream $s_{0}$, each user employs SIC to remove the common stream from
the received signal, and then decodes the intended private stream with the
SINR

$\gamma_{k}=\frac{|(\mathbf{g}_{k}^{H}+\mathbf{h}_{k}^{H}\mathbf{\Theta}\mathbf{G})\mathbf{w}_{k}|^{2}}{\sum\limits_{i=1,i\neq
k}^{K}|(\mathbf{g}_{k}^{H}+\mathbf{h}_{k}^{H}\mathbf{\Theta}\mathbf{G})\mathbf{w}_{i}|^{2}+\sigma_{k}^{2}}.$
(9)

The rate of decoding private message is
$r_{k}=\log_{2}\left(1+\gamma_{k}\right)$. User $k$ then reconstructs its
message by combining the submessages $W_{c,k}$ and $W_{p,k}$ respectively
decoded from the common and private streams [4].

In this work, we aim at jointly optimizing the scattering matrix of RIS
$\mathbf{\Theta}$ and the beamforming matrix $\mathbf{W}$ of RSMA to maximize
the sum-rate of the system. The sum-rate problem for the downlink fully
connected RIS aided RSMA network can be formulated as:

$\displaystyle(\mathcal{P}_{1})\,\,\max\limits_{\bm{\Theta},\mathbf{W}}\,\,$
$\displaystyle\sum\limits_{k=0}^{K}r_{k}$ (10a) s.t.
$\displaystyle\mathrm{tr}(\mathbf{WW}^{H})\leq P_{t},$ (10b)
$\displaystyle\bm{\Theta}^{H}\bm{\Theta}=\bm{I},$ (10c)
$\displaystyle\bm{\Theta}=\bm{\Theta}^{T}.$ (10d)

Constraint (10b) is the transmit power constraint at the BS. (10c) and (10d)
show that the reconfigurable impedance network in fully connected RIS is a
lossless and reciprocal circuit network.

When constraints (10c) and (10d) for scattering matrix $\mathbf{\Theta}$ are
replaced by constraint (6), $\mathcal{P}_{1}$ reduces to the sum-rate problem
of the single connected RIS aided RSMA [24]. When the power allocated to
$\mathbf{w}_{0}$ is fixed to zero, the problem reduces to the sum-rate problem
for the conventional single connected RIS aided SDMA [14].

##  III Alternative Optimization Framework

Problem $\mathcal{P}_{1}$ is a joint beamforming matrix and RIS scattering
matrix optimization problem. It is non-convex and the beamformers are coupled
with the scattering matrix in multiple fractional SINR expressions. Following
existing works [13, 14, 17, 18, 20, 19, 21], we propose an alternative
optimization (AO) framework to solve $\mathcal{P}_{1}$. Specifically, the
problem is first decomposed into the subproblems of beamforming design and
scattering matrix design. The former is solved by a weighted minimum mean
square error (WMMSE)-based approach while the latter is solved by the quasi-
Newton algorithm. The two subproblems are solved iteratively until
convergence. In the following subsections, the proposed optimization algorithm
is delineated.

###  III-A Beamforming Optimization

With a given scattering matrix $\bm{\Theta}$, the channel responses from the
RIS to the users are fixed. To ease notations, we denote the effective channel
from the BS and the RIS to user $k$ as

$\widetilde{\mathbf{g}}_{k}^{H}=\mathbf{g}_{k}^{H}+\mathbf{h}_{k}^{H}\mathbf{\Theta}\mathbf{G}.$
(11)

And $\mathcal{P}_{1}$ reduces to

$\displaystyle(\mathcal{P}_{2})\,\,\max\limits_{\mathbf{W}}\,\,$
$\displaystyle\sum\limits_{k=0}^{K}r_{k}$ (12a) s.t.
$\displaystyle\mathrm{tr}(\mathbf{WW}^{H})\leq P_{t},$ (12b)

which can be solved by the WMMSE algorithm [6] as briefly introduced below.

Denote the equalizers to estimate $s_{0}$ and $s_{k}$ as $e_{0,k}$ and
$e_{k}$, respectively. $\hat{s}_{0,k}=e_{0,k}y_{k}$ is the estimate of
$s_{0}$, and
$\hat{s}_{k}=e_{k}(y_{k}-\widetilde{\mathbf{g}}_{k}^{H}\mathbf{w}_{0}\hat{s}_{0,k})$
is the estimate of $s_{k}$ at user $k$. The mean square errors (MSEs) of
decoding $s_{0}$ and $s_{k}$ are calculated as

$\begin{split}\varepsilon_{0,k}&\triangleq\mathbb{E}\\{|\hat{s}_{0,k}-s_{0,k}|^{2}\\}=|e_{0,k}|^{2}T_{0,k}-2\Re\\{e_{0,k}\widetilde{\mathbf{g}}_{k}^{H}\mathbf{w}_{0}\\}+1,\\\
\varepsilon_{k}&\triangleq\mathbb{E}\\{|\hat{s}_{k}-s_{k}|^{2}\\}=|e_{k}|^{2}T_{k}-2\Re\\{e_{k}\widetilde{\mathbf{g}}_{k}^{H}\mathbf{w}_{k}\\}+1,\end{split}$
(13)

where
$T_{0,k}=\textstyle\sum_{i=0}^{K}|\widetilde{\mathbf{g}}_{k}^{H}\mathbf{w}_{i}|^{2}+\sigma_{k}^{2}$,
and
$T_{k}=\textstyle\sum_{i=1}^{K}|\widetilde{\mathbf{g}}_{k}^{H}\mathbf{w}_{i}|^{2}+\sigma_{k}^{2}$
are the average power of the received signal and the signal after removing the
common stream, respectively. By setting $\partial\varepsilon_{0,k}/\partial
e_{0,k}$ and $\partial\varepsilon_{k}/\partial e_{k}$ to zero respectively,
the minimum MSE (MMSE) equalizers are given by

$e_{0,k}^{\text{MMSE}}=\mathbf{w}_{0}^{H}\widetilde{\mathbf{g}}_{k}T_{0,k}^{-1},\,\,e_{k}^{\text{MMSE}}=\mathbf{w}_{k}^{H}\widetilde{\mathbf{g}}_{k}T_{k}^{-1}.$
(14)

Substituting (14) into (13), the MMSEs are given by
$\varepsilon_{0,k}^{\text{MMSE}}=(T_{0,k}-|\widetilde{\mathbf{g}}_{k}^{H}\mathbf{w}_{0}|^{2})T_{0,k}^{-1},\,\,\varepsilon_{k}^{\text{MMSE}}=(T_{k}-|\widetilde{\mathbf{g}}_{k}^{H}\mathbf{w}_{k}|^{2})T_{k}^{-1}.$
Then, the SINRs coresponding to $s_{0}$ and $s_{k}$ can be transformed to
$\gamma_{0,k}=1/\varepsilon_{0,k}^{\text{MMSE}}-1,\,\gamma_{k}=1/\varepsilon_{k}^{\text{MMSE}}-1.$
The rates of common and private streams become
$r_{0,k}=-\log_{2}(\varepsilon_{0,k}^{\text{MMSE}})$ and
$\,\,r_{k}=-\log_{2}(\varepsilon_{k}^{\text{MMSE}}).$ However, the logarithmic
rate-MMSE relationships above cannot be used directly for the sum-rate
problem. To tackle the issue, the augmented MMSEs are introduced as follows

$\xi_{0,k}\triangleq\lambda_{0,k}\varepsilon_{0,k}-\log_{2}(\lambda_{0,k}),\,\,\xi_{k}\triangleq\lambda_{k}\varepsilon_{k}-\log(\lambda_{k}),$
(15)

where $\lambda_{0,k}$ and $\lambda_{k}$ are auxiliary variables (also known as
weights) for the rate-WMMSE relationships of $r_{0,k}$ and $r_{k}$,
respectively. By
calculating$\frac{\partial\xi_{0,k}}{\partial\lambda_{0,k}}=0,\,\frac{\partial\xi_{k}}{\partial\lambda_{k}}=0$,
we obtain the optimum weights given by

$\lambda_{0,k}^{\text{MMSE}}=(\varepsilon_{0,k}^{\text{MMSE}})^{-1},\,\,\lambda_{k}^{\text{MMSE}}=(\varepsilon_{k}^{\text{MMSE}})^{-1}.$
(16)

Substituting (13) and (16) into (15), the rate-WMMSE relationships are
established as
$\xi_{0,k}^{\text{MMSE}}=1-r_{0,k},\,\,\xi_{k}^{\text{MMSE}}=1-r_{k}.$ With
the rate-WMMSE relationships above, $\mathcal{P}_{2}$ is equivalently
transformed into the WMMSE problem

$\displaystyle(\mathcal{P}_{3})\,\min\limits_{\mathbf{W},\bm{\lambda},\mathbf{e}}\,$
$\displaystyle\sum\limits_{k=0}^{K}\xi_{k}$ (17a) s.t.
$\displaystyle\mathrm{tr}(\mathbf{WW}^{H})\leq P_{t},$ (17b)

where
$\bm{\lambda}=[\lambda_{0,1},\cdots,\lambda_{0,K},\lambda_{1},\cdots,\lambda_{K}]^{T}$
is the weight vector and
$\mathbf{e}=[e_{0,1},\cdots,e_{0,K},e_{1},\cdots,e_{K}]^{T}$ is the equalizer
vector. $\xi_{0}=\max_{k\in\mathcal{K}}\xi_{0,k}$. $\mathcal{P}_{3}$ is still
non-convex, an AO framework is applied to decompose it into three convex
subproblems. For each block, one of $\mathbf{W},\bm{\lambda},\bm{e}$ is
optimized by fixing the other two blocks. Algorithm 1 specifies the procedure
of the WMMSE method to optimize the beamforming vectors. Readers are referred
to [25] for the details of the convergence proof for Algorithm 1.

1 Initialize: $t\leftarrow 0$, $\epsilon$, $\mathbf{W}^{[t]}$,
$\mathrm{SR}^{[t]}$;

2 repeat

3 t $\leftarrow$ t+1;

4 update $\mathbf{e}^{[t]},\mathbf{\lambda}^{[t]}$ by (14), (16);

5 update $\mathbf{W}^{[t]}$ by solving problem $\mathcal{P}_{3}$ using
$\mathbf{e}^{[t]},\mathbf{\lambda}^{[t]}$ ;

6

7until  _$|\mathrm{SR}^{[t]}-\mathrm{SR}^{[t-1]}| <\epsilon$_;

Algorithm 1 WMMSE algorithm for beamforming design

###  III-B Scattering Matrix Optimization

Similarly, with a given beamforming design $\mathbf{W}$, problem
$\mathcal{P}_{1}$ is simplified as

$\displaystyle(\mathcal{P}_{4})\,\,\max\limits_{\bm{\Theta}}\,\,$
$\displaystyle\sum\limits_{k=0}^{K}r_{k}$ (18a) s.t.
$\displaystyle\bm{\Theta}^{H}\bm{\Theta}=\bm{I},$ (18b)
$\displaystyle\bm{\Theta}=\bm{\Theta}^{T}.$ (18c)

However, it is challenging to transform problem $\mathcal{P}_{4}$ into a
convex problem due to the non-convex matrix equality constraints. Hence, we
apply equality (5) to equivalently reformulate $\mathcal{P}_{4}$ as

$\displaystyle(\mathcal{P}_{5})\,\max\limits_{\bm{X}}\,$
$\displaystyle\sum\limits_{k=0}^{K}r_{k}$ (19a) s.t.
$\displaystyle\mathbf{\Theta}=(j\mathbf{X}+Z_{0}\mathbf{I})^{-1}(j\mathbf{X}-Z_{0}\mathbf{I}),$
(19b) $\displaystyle\mathbf{X}=\mathbf{X}^{T}.$ (19c)

Substituting (19b) into (19a) , and removing (19c) (by defining a symmetric
matrix variable), problem $\mathcal{P}_{5}$ becomes an unconstrained
optimization problem. Moreover, the matrix variable $\mathbf{X}$ is a real
symmetric matrix in which $N(N+1)/2$ variables are adjustable. Such
unconstrained optimization problem can be directly solved by the quasi-Newton
method [16].

###  III-C Alternative Optimization Algorithm

1 Initilize: $n\leftarrow 0$, $\epsilon$, $\mathbf{W}^{[n]}$,
$\mathbf{X}^{[n]}$, $\mathrm{SR}^{[n]}$;

2 repeat

3 $n\leftarrow n+1$ ;

4 calculate $\mathbf{\Theta}^{[n-1]}$ by (19b);

5 Given $\mathbf{\Theta}^{[n-1]}$, update $\mathbf{W}^{[n]}$ by solving
$\mathcal{P}_{3}$ with Algorithm 1;

6 Given $\mathbf{W}^{[n]}$, update $\mathbf{X}^{[n]}$ by solving
$\mathcal{P}_{5}$ with the quasi-Newton method using $\mathbf{X}^{[n-1]}$ as
the initial point;

7

8until  _$|\mathrm{SR}^{[n]}-\mathrm{SR}^{[n-1]}| <\epsilon$_;

Algorithm 2 Alternative Optimization to solve ($\mathcal{P}_{1}$)

The proposed AO algorithm to jointly maximize the scattering matrix and the
beamforming matrix is specified in Algorithm 2. Starting with a feasible
beamforming matrix $\mathbf{W}^{[0]}$ and a symmetric reactance matrix
$\mathbf{X}^{[0]}$, in $n$-th iteration, we first calculate the scattering
matrix $\mathbf{\Theta}^{[n-1]}$ with reactance matrix $\mathbf{X}^{[n-1]}$
from the last iteration. Next, with a fixed scattering matrix
$\mathbf{\Theta}^{[n-1]}$, the beamforming matrix $\mathbf{W}^{[n]}$ is
updated by Algorithm 1. For a given $\mathbf{W}^{[n]}$, the reactance matrix
$\mathbf{X}^{[n]}$ is then updated based on the quasi-Newton method using
$\mathbf{X}^{[n-1]}$ as the initial point. The sum-rate $\mathrm{SR}^{[n]}$ is
then calculated based on the updated $\mathbf{W}^{[n]}$ and
$\mathbf{X}^{[n]}$. The process is repeated until convergence.

Convergence Analysis: In each iteration $[n]$, the solution of
$\mathcal{P}_{1}$ is also a feasible solution of $\mathcal{P}_{1}$ for the
next iteration. Hence, the sum-rate $\mathrm{SR}^{[n+1]}$ is larger than or
equal to $\mathrm{SR}^{[n]}$. Moreover, the non-decreasing sequences generated
by Algorithm 2 is bounded above by the transmit power constraint. Therefore,
the proposed AO algorithm is guaranteed to converge with a given tolerance
$\epsilon$.

##  IV Numerical Results

In this section, we evaluate the performance of the proposed system model and
the proposed algorithm. The following six schemes are compared:

  * •

Fully RIS RSMA: This is the scheme proposed in Section II.

  * •

Fully RIS SDMA: This is a special case of the proposed scheme when the power
allocated to the common stream is fixed to zero.

  * •

Single RIS RSMA: This is the single connected RIS aided RSMA scheme, as
studied in [24].

  * •

Single RIS SDMA: This is the conventional single connected RIS aided SDMA
scheme, as studied in [14].

  * •

no RIS RSMA: This is the conventional RSMA scheme without using RIS, as
studied in [4, 6, 5, 7, 8].

  * •

no RIS SDMA: This is the conventional multi-user linearly precoded SDMA scheme
without using RIS, as studied in [4, 26].

The sum-rate maximization problems of fully/single RIS RSMA and fully/single
RIS SDMA are solved by Algorithm 2 while the corresponding problems of no RIS
RSMA and no RIS SDMA are solved by WMMSE directly. Problem $\mathcal{P}_{3}$
is solved by the CVX toolbox [27] with the interior-point method, and problem
$\mathcal{P}_{5}$ is solved by the optimization toolbox in Matlab with the
quasi-Newton method.

The setting of the simulation follows [13], which is a planar RIS-aided
network as shown in Fig. 2. The BS and RIS are located at $(0,0)$ and
$(50,50)$, respectively. In addition, there are $K=4$ users randomly generated
in a circle centered at $(150,0)$ meters with a diameter of 20 meters. The
path loss of the channels are modeled as $P(d)=L_{0}d^{-\alpha}$, where
$L_{0}=-30$ dB is the reference path loss at $d=1$ m, $d$ refers to the link
distance, and $\alpha$ denotes the path loss exponent. Assuming that the
location of RIS is chosen carefully, we set the path loss exponent of BS to
users, BS to RIS, and RIS to users are 3.5, 2, 2.2, respectively [13]. For
simplicity, the small-scale fading of all channels are modeled as Rayleigh
fading. Hence, the channels are given as
${\mathbf{g}}_{k}\sim\mathcal{CN}(0,P(d_{k}^{g})\mathbf{I})),\mathbf{h}_{k}\sim\mathcal{CN}(0,P(d_{k}^{h})\mathbf{I}),$
and $\mathbf{G}\sim\mathcal{CN}(0,P(d^{G})\mathbf{I})$, where
$d_{k}^{g},d_{k}^{h}$ and $d^{G}$ respectively denote the distance between the
BS and user $k$, the distance between the RIS and user $k$, and the distance
between the BS and RIS. Besides, the reference impedance of RIS is $Z_{0}=50$
$\Omega$, the convergence tolerance is $\epsilon=10^{-3}$, and the noise at
user $k$ is $\sigma_{k}^{2}=1$. The transmit singal-to-noise ratio (SNR)
SNR$\triangleq P_{t}/\sigma_{k}^{2}$ is therefore equal to the transmit power
numerically. All simulation results are averaged over 100 random channel
realizations.

Figure 2: The Simulated RIS-aided $K$-user MISO transmission scenario. Figure
3: Sum rate versus the transmit power, when $M=4$, $K=4$, and $N=32$.

Fig. 3 shows the sum-rate of different strategies versus the transmit power
when $M=4$, $K=4$, and $N=32$. It shows that the proposed fully connected RIS
aided RSMA scheme outperforms all other baseline schemes. The relative sum-
rate gain of fully RIS RSMA over single RIS RSMA and no RIS RSMA are at least
$4.6\%$ and $16.5\%$ respectively when SNR is $30$ dB. By using the fully
connected RIS aided RSMA model, the sum-rate of the multi-user multi-antenna
network increases significantly. Moreover, the single connected RIS aided RSMA
achieves approximately the same sum-rate as the fully connected RIS aided SDMA
in the high SNR regime.

Figure 4: Sum rate versus the number of RIS Elements, when $M=4$, $K=4$, and
SNR is $30$ dB.

Fig. 4 shows the impact of the number of passive reflecting elements at the
RIS (i.e., $N$) to the sum-rate of different strategies when $M=4$, $K=4$, and
SNR is $30$ dB. For both fully RIS RSMA and SDMA schemes, the sum-rate
increases faster than the corresponding single RIS RSMA and SDMA as $N$
increases. Particularly, the single RIS RSMA achieves a higher sum-rate than
the fully RIS SDMA when the number of elements is less than $16$. In this
regime, the gain obtained by RSMA scheme is more significant than the gain
obtained by fully connected RIS.

Figure 5: Convergence of the algorithms in one channel realization.

Fig. 5 illustrates the convergence of Algorithm 2 for the fully RIS RSMA
scheme and other baseline schemes (no RIS schemes are not included since they
only adopt Algorithm 1) when $M=4$, $K=4$, $N=32$ and SNR is $25$ dB. It can
be observed that single connected RIS aided schemes converge faster than the
fully connected RIS aided schemes due to the smaller number of variables in
the RIS scattering matrix. In general, the algorithm can converge with 100
iterations.

##  V Conclusion

In this work, we propose a fully connected RIS aided RSMA downlink
transmission network. The beamforming vectors at the BS and the scattering
matrix of the fully connected RIS are jointly designed to maximize the sum-
rate of the network. To solve this problem, we propose an effective algorithm
that alternatively optimizes the beamforming and scattering matrices.
Simulation results show the outstanding spectral efficiency of the proposed
fully connected RIS aided RSMA scheme over the existing transmission schemes.
It acts as a new benchmark for linearly precoded multi-user multi-antenna
networks. Moreover, we show that the single connected RIS aided RSMA can
achieve approximately the same spectral efficiency as the fully connected RIS
aided SDMA in the high SNR regime. Therefore, we conclude that by marrying
RSMA and RIS, the spectral efficiency can be enhanced significantly.

## References

  * [1] O. Dizdar, Y. Mao, Y. Xu, P. Zhu, and B. Clerckx, “Rate-splitting multiple access for enhanced URLLC and eMBB in 6G,” _Proc. Int. Symp. Wirel. Commun. Syst. (ISWCS)_ , 2021.
  * [2] H. Tataria, M. Shafi, A. F. Molisch, M. Dohler, H. Sjoland, and F. Tufvesson, “6G wireless systems: Vision, requirements, challenges, insights, and opportunities,” _Proc. IEEE_ , vol. 109, no. 7, pp. 1166–1199, 2021.
  * [3] Y. Mao, O. Dizdar, B. Clerckx, R. Schober, P. Popovski, and H. V. Poor, “Rate-splitting multiple access: Fundamentals, survey, and future research trends,” _arXiv:2201.03192_ , 2021.
  * [4] Y. Mao, B. Clerckx, and V. O. Li, “Rate-splitting multiple access for downlink communication systems: Bridging, generalizing, and outperforming SDMA and NOMA,” _Eurasip J. Wireless Commun. Networking_ , 2018.
  * [5] B. Clerckx, H. Joudeh, C. Hao, M. Dai, and B. Rassouli, “Rate splitting for MIMO wireless networks: A promising PHY-layer strategy for LTE evolution,” _IEEE Commun. Mag._ , vol. 54, no. 5, pp. 98–105, 2016.
  * [6] H. Joudeh and B. Clerckx, “Sum-rate maximization for linearly precoded downlink multiuser MISO systems with partial CSIT: A rate-splitting approach,” _IEEE Trans. Commun._ , vol. 64, no. 11, pp. 4847–4861, 2016\.
  * [7] Y. Mao, B. Clerckx, and V. O. K. Li, “Rate-splitting for multi-antenna non-orthogonal unicast and multicast transmission: Spectral and energy efficiency analysis,” _IEEE Trans. Commun._ , vol. 67, no. 12, pp. 8754–8770, 2019.
  * [8] B. Clerckx, Y. Mao, R. Schober, E. Jorswieck, D. J. Love, J. Yuan, L. Hanzo, G. Y. Li, E. G. Larsson, and G. Caire, “Is NOMA efficient in multi-antenna networks? A critical look at next generation multiple access techniques,” _IEEE open J. Commun. Soc._ , 2021.
  * [9] Y. Mao, B. Clerckx, and V. O. Li, “Energy efficiency of rate-splitting multiple access, and performance benefits over SDMA and NOMA,” in _Proc. Int. Symp. Wirel. Commun. Syst. (ISWCS)_ , 2018, pp. 1–5.
  * [10] Y. Mao, B. Clerckx, J. Zhang, V. O. Li, and M. A. Arafah, “Max-min fairness of k-user cooperative rate-splitting in MISO broadcast channel with user relaying,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 10, pp. 6362–6376, 2020.
  * [11] Y. Mao and B. Clerckx, “Beyond dirty paper coding for multi-antenna broadcast channel with partial CSIT: A rate-splitting approach,” _IEEE Trans. Commun._ , vol. 68, no. 11, pp. 6775–6791, 2020.
  * [12] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface aided wireless communications: A tutorial,” _IEEE Trans. Commun._ , 2021.
  * [13] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, 2019.
  * [14] H. Guo, Y.-C. Liang, J. Chen, and E. G. Larsson, “Weighted sum-rate maximization for reconfigurable intelligent surface aided wireless networks,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 5, pp. 3064–3076, 2020.
  * [15] X. Yu, V. Jamali, D. Xu, D. W. K. Ng, and R. Schober, “Smart and reconfigurable wireless communications: From IRS modeling to algorithm design,” _arXiv preprint arXiv:2103.07046_ , 2021.
  * [16] S. Shen, B. Clerckx, and R. Murch, “Modeling and architecture design of reconfigurable intelligent surfaces using scattering parameter network analysis,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2021.
  * [17] Z. Yang, J. Shi, Z. Li, M. Chen, W. Xu, and M. Shikh-Bahaei, “Energy efficient rate splitting multiple access (RSMA) with reconfigurable intelligent surface,” in _IEEE Int. Conf. Commun. Workshops (ICC Workshops)_ , 2020, pp. 1–6.
  * [18] H. Fu, S. Feng, and D. W. Kwan Ng, “Resource allocation design for IRS-aided downlink MU-MISO RSMA systems,” in _IEEE Int. Conf. Commun. Workshops (ICC Workshops)_ , 2021, pp. 1–6.
  * [19] K. Weinberger, A. A. Ahmad, and A. Sezgin, “On dynergistic benefits of tate dplitting in IRS-assisted cloud radio access networks,” in _IEEE Int. Conf. Commun. (ICC)_ , 2021, pp. 1–6.
  * [20] A. Bansal, K. Singh, and C.-P. Li, “Analysis of hierarchical rate splitting for intelligent reflecting surfaces-aided downlink multiuser MISO communications,” _IEEE open J. Commun. Soc._ , vol. 2, pp. 785–798, 2021\.
  * [21] A. Bansal, K. Singh, B. Clerckx, C.-P. Li, and M.-S. Alouini, “Rate-splitting multiple access for intelligent reflecting surface aided multi-user communications,” _IEEE Trans. Veh. Technol_ , vol. 70, no. 9, pp. 9217–9229, 2021.
  * [22] D. M. Pozar, “Microwave engineering USA: John Wiley & Sons,” 2009.
  * [23] N. K. Kundu, Z. Li, J. Rao, S. Shen, M. R. McKay, and R. Murch, “Optimal grouping strategy for reconfigurable intelligent surface assisted wireless Ccommunications,” _arXiv preprint arXiv:2111.10550_ , 2021.
  * [24] A. Jolly, S. Biswas, and K. Singh, “An analysis on rate-splitting multiple access for IRS aided 6G communication,” _arXiv preprint arXiv:2106.04418_ , 2021.
  * [25] Q. Shi, M. Razaviyayn, Z.-Q. Luo, and C. He, “An iteratively weighted MMSE approach to distributed sum-utility maximization for a MIMO interfering broadcast channel,” _IEEE Trans. Signal Processing_ , vol. 59, no. 9, pp. 4331–4340, 2011.
  * [26] S. S. Christensen, R. Agarwal, E. De Carvalho, and J. M. Cioffi, “Weighted sum-rate maximization using weighted MMSE for MIMO-BC beamforming design,” _IEEE Trans. Wireless Commun._ , vol. 7, no. 12, pp. 4792–4799, 2008.
  * [27] M. Grant, S. Boyd, and Y. Ye, “CVX: Matlab software for disciplined convex programming,” 2008.


# Non-minimal coupling inspires the Dirac cosmological model

H<EMAIL_ADDRESS>H.
<EMAIL_ADDRESS>A. H<EMAIL_ADDRESS>U. K<EMAIL_ADDRESS>1 Research Institute for Astronomy and
Astrophysics of Maragha (RIAAM), University of Maragheh, P.O. Box 55136-553,
Maragheh, Iran
2 Physics Department, Faculty of Sciences, University of Sistan and
Baluchestan, Zahedan, Iran
3 Department of Mathematics, Institute of Applied Sciences and Humanities, GLA
University, Mathura-281406, Uttar Pradesh, India

###### Abstract

In the framework of the generalized Rastall theory (GRT), we study the ability
of a non-minimal coupling between geometry and matter fields in order to
provide a setting which allows for a variable $G$ during the cosmic evolution.
In this regard, the compatibility of this theory with Dirac hypothesis on the
variations of $G$ is investigated, and additionally, the possibility of
obtaining the current accelerated universe is also addressed. In summary, our
study indicates that, in GRT, having in hand the $G$ profile, one may find the
corresponding non-minimal coupling between the energy source and geometry and
vise versa, in a compatible way with the current accelerated universe.

##  I Introduction

The idea that $G$ (the Newtonian gravitational coupling) has probably
experienced diverse values during the cosmic evolution has many motivations.
It began with Dirac’s proposal dir1 ; dir2 ; dir3  which states that, the
ubiquitousness of certain large dimensionless numbers (LDN’s), arising in
combinations of physical constants and cosmological quantities WZE  was not a
coincidence but an outcome of an underlying relationship between them BTB . In
his proposal, Dirac pointed out that the electrical force between proton and
electron within a hydrogen atom i.e., $F_{e}=e^{2}/4\pi\epsilon_{0}r^{2}$ is a
large number being 40 orders of magnitude greater than their gravitational
force $F_{G}=Gm_{p}m_{e}/r^{2}$, i.e.,

$\displaystyle{\rm
LN}_{1}=\frac{F_{e}}{F_{G}}=\frac{e^{2}}{4\pi\epsilon_{0}Gm_{p}m_{e}}\approx
10^{40},$ (1)

where $m_{e},e,m_{p},\epsilon_{0}$ and $G$ are the mass and charge of
electron, the proton mass, the vacuum permittivity and gravitational constant,
respectively. On the other side, the ratio of the age of the universe and the
time for light to traverse an electron is also nearly of the same size, i.e.,

$\displaystyle{\rm LN}_{2}=\frac{t}{e^{2}/4\pi\epsilon_{0}m_{e}c^{3}}\approx
10^{40}.$ (2)

Dirac then suggested that the above two quantities are equal. As a result of
such a relationship, some of the fundamental constants cannot remain constant
for ever since ${\rm LN}_{2}$ varies with the age of the universe. According
to Dirac’s hypothesis, atomic parameters cannot change with time and thus $G$
should change inversely with time, i.e., $G\propto t^{-1}$ CHK , see also
DIRACREV  for recent reviews. Since the advent of this idea, it has led to
interesting implications within theoretical physics, and has attracted a great
deal of attention during the past decades ras2 ; vin ; sab ; bap ; bee ; wu ;
deg ; bar1 ; bar2 ; bar3 ; mans ; gaz ; clif ; bro ; sol ; uza1 ; uza2 ; smo ;
fri ; les . Moreover, it has even interesting power to justify baryogenesis
les , the current and primary accelerated universes ell  and can support the
de Sitter spacetime uza1 ; uza2 .

In Newtonian gravity one is allowed to write an explicit time variation of $G$
without the need of satisfying any further constraint. However, the situation
is different in GR as there are further constraints to be satisfied. Consider
the Einstein field equation $G^{\mu}_{\,\nu}=8\pi GT^{\mu}_{\,\,\nu}$ with the
assumption of $G=G(t)$ and $c\equiv 1$. If one takes the covariant divergence
of this equation the left hand side vanishes as a result of Bianchi identity.
Then, if the ordinary energy-momentum conservation law (OCL) is assumed to
hold, i.e., $T^{\mu}_{\,\,\nu;\mu}=0$, one finds that $G$ must be a constant
with respect to spacetime coordinates, i.e., $\partial G/\partial x^{\mu}=0$
always. In this respect, GR does not allow for any variation in the
gravitational coupling $G$ owing to the fact that the Einstein tensor is
divergence free and the divergence of energy-momentum tensor is also zero.
Hence, in the light of Dirac’s proposal, some modifications of GR field
equation are essential. This is because, if we simply let $G$ to be a variable
then the OCL is violated CanutoAdams . In this respect, investigating the
effects of a varying $G$ can be performed only through modified field
equations along with modified conservation laws. From these arguments, one may
intuitively imagine that a varying $G$ could contribute as a new degree of
freedom within the OCL. As in GR, $G$ denotes mutual relation between geometry
and matter fields, hence, variations of $G$ together with the violation of OCL
may be considered as a signal for the idea that another relation between
geometry and matter fields may exist that connects their changes to each
other. However, there are modifications of GR with a varying $G$ that respect
the OCL such as Brans-Dicke theory, in which, the dynamical scalar field
$\phi$ can be considered as the origin of gravitational coupling and thus it
varies as $G\propto\frac{1}{\phi}$ 11 ; 12 ; 13 ; 14 ; bar2 .

OCL, as one of the cornerstones of GR pois , is not respected in all modified
gravity theories, for example, it is broken in the non-minimal curvature
matter coupling theories od1 ; all ; koi ; bert ; hark ; car ; boh ; ras1 ;
mor1 ; mora . Rastall gravity is a pioneering theory in this area ras1  in
accordance with various observations li ; raw1 ; raw2 ; raw3 ; maj ; arb ;
rah1 ; rah2 ; rah3 ; mor2 ; man ; ortiz2020 ; shabooni2020  and its
corresponding cosmology avoids the age and entropy problems arisen in the
framework of the standard cosmology fab . In fact, this theory can even
provide a better platform for describing the matter dominated era compared to
the Einstein theory raw2 . A generalized form of this theory allows us to
relate the current and primary accelerated universe to the ability of the
spacetime to couple with the energy-momentum sources, filling the background,
and in fact, introduces this coupling as a candidate for dark energy and
inflaton field mor1 .

In addition to inflationary models powered by employing varying $G$ theories
ell , there are also other models to describe inflation without considering an
inflaton field jos ; mor1 ; wat ; gam . In Ref. mor1 , it has been shown that
while the existence of an interaction between the geometry and matter fields
may model the primary and current inflationary eras, it does not necessarily
lead to the break-down of OCL. In fact, if geometry has the ability of being
non-minimally coupled with the matter fields, then this ability may support
the primary inflationary era and the current accelerated phase mor1 . To
obtain these results, authors focus on the special case of $T^{\mu\nu}_{\ \
;\mu}=0$, and find out the form of non-minimal coupling in each cosmic era.

The study of various non-minimal couplings can at least make us familiar with
their consequences and properties which may finally lead to a better
understanding of spacetime that helps us provide better predictions about its
behavior and nature. In GRT, cosmological scenarios jos ; mor1 ; das ; lin
imply the power of non-minimal coupling in $i$) providing both singular and
non-singular universes, $ii$) describing the particle production process,
$iii$) avoiding the coincidence problem, and $iv$) playing the role of dark
energy (unlike the original Rastall theory mor1 ; batis ). In this regard,
thermodynamically it has also been shown that the confusion in defining energy
and some of its outcomes which may lead to the OCL generalization (or
equivalently, the breakdown of OCL) could make the cosmos dark mor3 ; mor4 .

Since in Rastall gravity, the gravitational coupling is a constant, but
differs from those of GR and Newtonian gravity (NG) mor5 ; ras1 , Rastall
theory (and indeed a mutual non-minimal coupling between the geometry and
matter fields in the Rastall way) cannot provide a theoretical basis for the
probable variations of $G$ during the cosmic evolution. These points will be
reopened in more details in the next section.

Motivated by the above arguments, it is reasonable to $i$) examine the ability
of non-minimal coupling between geometry and matter fields in producing a non-
constant $G$, and also $ii$) study the results of a non-constant $G$ in the
framework of Rastall theory. The latter is tried to be answered by some
authors in Ref. ref , by combining Rastall and Brans-Dicke theories with each
other. In the present study, the changes in $G$ is not originated by the
Rastall theory meaning that the first part is still unsolved and debateable.
We therefore focus on GRT to describe the compatibility of a non-minimal
coupling with Dirac’s idea on evolution of $G$. Indeed, we are eager to show
that, at least phenomenologically, a non-minimal coupling may itself change
$G$ and play the role of dark energy.

The present work is then arranged as follows. In Sects. II and III, a brief
review on the Rastall theory and its generalization mor1  has been provided,
and some of their predictions about the variations of $G$ are addressed. Sect.
IV includes our survey on the possibility of explaining a well-known Dirac
cosmological model, previously introduced by other authors, within the
framework of GRT. To show the ability of non-minimal coupling in satisfying
Dirac hypothesis and describing the cosmic evolution, simultaneously, a new
model is also introduced in Sect. V. Sect. VI is devoted to concluding
remarks. Here, we use $c=\hbar=1$ units.

##  II Rastall theory and a model for varying $G$

Originally, P. Rastall argued that the OCL may not be valid in a curved
spacetime leading to ras1

$\displaystyle T^{\mu\nu}_{\ \ ;\mu}\neq 0,$ (3)

in the non-flat spacetimes. From the mathematical point of view,
$T^{\mu\nu}_{\ \ ;\mu}$ is a ranked one tensor field written as $T^{\mu\nu}_{\
\ ;\mu}=Q^{,\nu}$ where $Q$ is an unknown scalar function found out from other
parts of physics, mathematics and observations ras1 . Since $Q$ is a scalar
and Rastall hypothesis admits the violation of OCL in a curved spacetime
(where Ricci scalar is not always zero), therefore Ricci scalar, R, can be
considered as a suitable suggestion for $Q$, and thus ras1

$\displaystyle T^{\mu\nu}_{\ \ ;\mu}=\lambda^{\prime}R^{;\nu},$ (4)

where $\lambda^{\prime}$ is called the Rastall constant parameter. Using the
Bianchi identity, it is easy to get

$\displaystyle
G_{\mu\nu}+\kappa^{\prime}\lambda^{\prime}g_{\mu\nu}R=\kappa^{\prime}T_{\mu\nu},$
(5)

which $\kappa^{\prime}$ is a constant ras1  called the Rastall gravitational
coupling constant. Applying the Newtonian limit on this result, we obtain ras1

$\displaystyle\frac{\kappa^{\prime}}{4\kappa^{\prime}\lambda^{\prime}-1}\left(3\kappa^{\prime}\lambda^{\prime}-\frac{1}{2}\right)=\kappa_{G},$
(6)

where $\kappa_{G}\equiv 4\pi G$. Hence, since $\kappa^{\prime}$ and
$\lambda^{\prime}$ are constants, $G$ should also be a constant as well (the
current value of $G$, namely $G_{0}$, is proper option leading to
$\kappa_{G}\equiv\kappa_{G_{0}}=4\pi G_{0}$). We therefore conclude that,
since the left hand side of (6) is a constant then a mutual non-minimal
interaction between the geometry and matter fields within the framework of
original version of Rastall gravity does not support the varying $G$ theories.
Eq. (6) also reveals that the Rastall gravitational coupling constant
($\kappa^{\prime}$) differs from that of GR ($2\kappa_{G}=8\pi G$) and only if
$\lambda^{\prime}=0$ then they will be equal.

It is also useful to note that one may use Eq. (5) in order to introduce the
generalized energy-momentum tensor
$\Theta_{\mu\nu}=T_{\mu\nu}-(\kappa^{\prime}\lambda^{\prime})/(4\kappa^{\prime}\lambda^{\prime}-1)Tg_{\mu\nu}$
which finally leads to the GR counterpart form of the Rastall field equations,
given as $G_{\mu\nu}=\kappa^{\prime}\Theta_{\mu\nu}$. In this manner, although
the obtained field equations are similar to those of GR, their solutions for
$T_{\mu\nu}$ differ in general from those of GR mor4 ; dar , a result
confirmed by various observational data, see e.g., li ; mor2 ; dar  and
references therein).

One can also generalize the Rastall theory by considering
$\lambda^{\prime}\rightarrow\lambda$, where $\lambda$ is a varying parameter.
Therefore Eq. (4) is extended as follows mor1

$\displaystyle T^{\mu\nu}_{\ \ \ ;\mu}=\left(\lambda R\right)^{;\nu},$ (7)

which finally leads to

$\displaystyle G_{\mu\nu}+\kappa\lambda g_{\mu\nu}R=\kappa T_{\mu\nu},$ (8)

where $\kappa$ is again a constant but $\lambda$ can change over time. Using
the trace of Eq. (8), one can also rewrite this equation as

$G_{\mu\nu}+\tau Tg_{\mu\nu}=\kappa T_{\mu\nu},$ (9)

in which

$\tau=\frac{\kappa^{2}\lambda}{4\kappa\lambda-1}.$ (10)

Now, since $\kappa$ is constant, the covariant derivative of Eq. (9) leads to

$\tau^{,\nu}T+\tau T^{,\nu}=\kappa T^{\nu\,\,\,;\mu}_{\,\,\,\mu},$ (11)

meaning that even if OCL is respected and until $\tau\neq constant$ (or
equally, $\lambda\neq constant$), the non-minimal coupling affects the
evolution of the energy-momentum source and vice versa mor1 . Therefore,
unlike the Rastall theory, OCL can be met in this framework even in the
presence of non-minimal coupling. In this regard, it is shown that, in the
framework of Eq. (8), even if OCL is met, the accelerated universe can be
explained under the shadow of $\lambda$ without resorting to a dark energy
source mor1 .

Now, considering the Newtonian limit (ignoring the pressure of $T_{\mu\nu}$
and utilizing relation $R_{00}=\nabla^{2}\phi$, in which $\phi$ denotes the
Newtonian potential mor6 ), one can easily find

$\displaystyle\frac{\kappa}{4\kappa\lambda-1}\left(3\kappa\lambda-\frac{1}{2}\right)=\kappa_{G}.$
(12)

Due to the similarity of Eqs. (8) and (5), one could expect that the Newtonian
limit of field equations (8) is obtainable by replacing $\kappa^{\prime}$ and
$\lambda^{\prime}$ with $\kappa$ and $\lambda$, respectively, in Eq. (6). Eq.
(12) also indicates that $G$ (or equally $\kappa_{G}$) does not necessarily
remain constant in this theory. Therefore, this generalization of Rastall
theory provides a basis for theories including a varying $G$ dir1 ; dir2 ; vin
; sab ; bap ; bee ; wu ; mans ; gaz ; bar1 ; bar2 ; bar3 ; clif ; uza1 ; uza2
; smo ; fri ; les . In fact, this equation tells that a non-minimal coupling
between the geometry and matter fields can make $G$ variable mot  meaning that
such coupling can be considered as a theoretical basis for varying $G$
theories.

##  III Newtonian limit, a model for running $G$, and the value of $\kappa$

Now, using Eq. (12), and following Ref. mor1 , in which
$\kappa\lambda\equiv\beta=[4+\theta(1+z)^{3}]^{-1}$, where $\theta$ is an
unknown constant and $z$ denotes the redshift, one can obtain

$\displaystyle\kappa_{G}=\frac{\kappa}{2}\left[1-\frac{2}{\theta(1+z)^{3}}\right],$
(13)

finally leading to

$\displaystyle\kappa_{G}=\frac{\kappa}{2}\left[1-\frac{2}{\theta}\right]\equiv\kappa_{G_{0}},$
(14)

and

$\displaystyle\kappa_{G}=\frac{\kappa}{2},$ (15)

for $z\rightarrow 0$ and $z\rightarrow\infty$, respectively. Based on Ref.
mor1 , whenever $0<\theta\leq 1/2$ (leading to $\beta>0$), the current
accelerated universe is explainable in the presence of OCL, and without
considering a dark energy-like source. Moreover, expression
$\beta=[4+\theta(1+z)^{3}]^{-1}$ is present in both of the matter dominated
era (MDE) and the current accelerated universe mor1 . Hence, Eq. (15) can be
considered as the value of $G$ at the beginning of MDE whereas the value of
$\kappa$ is obtainable by using Eq. (14)

$\displaystyle\kappa=\frac{8\pi G_{0}}{1-\frac{2}{\theta}},$ (16)

combined with Eq. (15) to see that $\kappa$, and thus $\kappa_{G}$ are
negative at the beginning of MDE. Therefore, in the model proposed in Ref.
mor1  which still respects OCL in the framework of (8), $G$ is not always
positive during the cosmic evolution. Negative values of $\kappa$ provide a
setting for baryonic matters to support traversable wormholes in the Rastall
framework mor5 . Moreover, in the framework of GRT, it has been shown that
negative values of $\kappa$ could have their own effects on matter
perturbations and formation of structures in large scale universe AHH2020 . In
this regard, overdense and underdense regions in the universe could form
periodically so that both large scale structures and voids could form as the
universe evolves from MDE to present time. Also, emergence of structures in a
class of alternative theories of gravity has been reported in Lohiya1996 ,
where the authors considered a non-minimally coupled scalar field in addition
to an induced negative gravitational constant and studied structure formation
with repulsive gravitation on the large scale. In the framework of general
scalar tensor-theories, a cosmological mechanism has been proposed in which it
is possible for $G$ to change sign from a positive branch (attracting) to a
negative branch (repulsive gravity) and vice versa Nunez2019 . It is also
worth mentioning that negative values of $G$ have previously been reported in
some other approaches studying the variations of $G$ bar2 ; uza1 ; uza2 .

Beside the effects of repulsive gravity (represented by a universal negative
coupling) on the evolution of perturbations and formation of structures, the
study of possible consequences of $\kappa<0$ on the stability of the model is
of particular importance. In this regard, from the viewpoint of perturbative
analysis, the existence of a repulsive gravity phase in the evolution of the
universe could lead to growing models with respect to scalar perturbations
producing then, large inhomogeneities. Hence a repulsive phase may destroy
homogeneity and in this sense it may be unstable Batista2001 . In Star1981 ,
it has been discussed that a transition from positive gravitational coupling
$G$ to negative one results in an instability, in such a way that, small
deviations from isotropy and homogeneity within the gravitational field will
grow unboundedly, leading to a true cosmological singularity at the boundary
between gravity and anti gravity. Also, investigating classical stability of
the model through dynamical system approach is of long-standing interest and
significance. Work along this line has been carried out for a class of GRT
models Lin2020 , where the authors have shown that the eventual fate of the
universe ends in late time attractors which are classically stable. However,
investigating these issues for the present model needs a deeper analysis with
more scrutiny and future studies will be reported elsewhere. Finally, we note
that, since $\dot{G}$ does not decrease with time for $0<\theta\leq 1/2$
($\dot{G}>0$ in this manner), this model does not respect the Dirac’s
hypothesis claiming that $G$ should decrease as a function of time dir1 ; dir2
; vin ; bap . Hence, more comprehensive non-minimal couplings are needed to
provide settings for Dirac hypothesis and also to model the cosmic evolution
without considering a mysterious fluid (dark energy), simultaneously.

###  III.1 Another possibility

In Ref. das , choosing $\lambda=(1+d_{0}H)/[3\kappa(w+1)]$, in which $w\equiv
p/\rho$ (where $p$ and $\rho$ denote the pressure and energy density of the
cosmic fluid, respectively), it has been shown that non-singular cosmic
evolution is obtainable in GRT. In this case $d_{0}$ is a free parameter, and
some outcomes of this proposal in various cosmic eras have also been studied
in Ref. das . Accepting this proposal along with considering the unit
$\kappa=8\pi G_{0}$ and also assuming $G(H_{0})=G_{0}$ (which helps us in
finding $d_{0}$), one easily reaches

$\displaystyle G(H)=G_{0}\frac{3(1-w)H_{0}-6H}{(1-3w)H_{0}-4H},$ (17)

where $H_{0}$ is the current value of $H$ and use has been made of Eq. (12).

##  IV Dirac cosmological model

As in the present model there is no evolution equation for the variation of
$G$ which is promoted as a dynamical field, one then has to impose a suitable
ansatz on the behavior of this parameter. Based on Dirac hypothesis, $G$
should decrease with time, i.e, $G\propto t^{-1}$ CHK . In general, one may
consider $G=G_{0}f$, in which $f$ is a decreasing function of time dir1 ; dir2
; vin ; bap ; clif ), in order to preserve Dirac hypothesis. Now, combining
Eq. (12) with $\kappa=8\pi G_{0}\alpha$ raw2 , along with Eqs. (7) and (8) for
a flat FLRW universe, one finds

$\displaystyle\gamma\equiv\lambda\kappa=\frac{f-\alpha}{4f-6\alpha},$ (18)
$\displaystyle
3\int(\rho+p)\frac{da}{a}=\frac{1}{2\alpha}\Big{[}(f-3\alpha)\rho-3(f-\alpha)p\Big{]},$
$\displaystyle H^{2}=\frac{1}{6}\Big{[}(3\alpha-f)\rho+3(f-\alpha)p\Big{]},$
$\displaystyle
q=-1-\frac{\dot{H}}{H^{2}}=-1+\frac{3\alpha(\rho+p)}{\rho(3\alpha-f)+3(f-\alpha)p},$

whenever a fluid with energy density $\rho$ and pressure $p$ fills the
background. We note that $\gamma$ is a varying parameter and, $q$ and $a$
denote deceleration parameter and scale factor, respectively, and we also have
assumed $8\pi G_{0}=1$.

Figure 1: The evolution of $q$ and state parameter $w$ versus $z$ for
$H(z=0)=67$ dom . Upper panels are provided for the case (i) and the lower
ones are depicted for case (ii) discussed in Sect. IV. The model parameters
used to draw the curves of $w$ are the same as those of $q$ diagrams.

The case with $f=a^{-n}$ leads to a decreasing function of time whenever $n>0$
gaz ; smo . In this manner, assuming $w\equiv p/\rho=0$, together with using
Eqs. (18), one easily finds $q=(3\alpha-1)^{-1}$, and
$\rho=\rho_{0}a^{n}(1-3\alpha a)^{-(n+2)/n}$, where $\rho_{0}$ is the
integration constant. These results indicate that, at limit $a\rightarrow 1$,
the obtained pressureless fluid can accelerate the universe expansion with
$q\leq-1/2$ for $-1/3\leq\alpha<1/3$. Consequently, the non-minimal coupling
$\gamma=[(1+z)^{n}-\alpha]/[4(1+z)^{n}-6\alpha]$ allows $G$ to vary as
$G=G_{0}(1+z)^{n}$ gaz , where we used the $1+z=1/a$ relation. It is also easy
to see that the universe described by this model has begun from a primary
inflationary phase ($q=-1$) corresponding to the $a\rightarrow 0$ point. In
fact, in this limit, we also have $\gamma=1/4$, a value that supports an
inflationary phase for even an empty universe mor1 .

Now, let us consider two more comprehensive cases i.e., $i$)
$p=k\rho^{1+1/m}$, where $m$ and $k$ are unknown constants to be evaluated
later, and $ii$) $p=\sigma\rho/(a-b\rho)-c\rho^{2}$ in which $\sigma$, $a$,
$b$ and $c$ are unknown coefficients. In this manner, as it is obvious from
Fig. 1, a proper behavior is obtainable for the cosmos. Here, $w\equiv p/\rho$
denotes the equation of state of cosmic fluids. Depending on the values of
unknown parameters, the universe can also experience a transition at $z_{t}$
which can even take values smaller than $1$. Clearly, both fluids behave as
dark energy sources, and the corresponding non-minimal coupling can not be
considered as a dark energy source.

##  V A new proposal for $\lambda$ parameter

Now, let us consider a flat FRW universe filled by a pressureless fluid with
energy density $\rho$ when $\lambda R=\zeta H^{n}$ in which $\zeta$ and $n$
are unknown constants. In this manner, the $\lambda$ parameter takes the form

$\displaystyle\lambda=\zeta\frac{H^{n}}{R}=\frac{\zeta}{6}\frac{H^{n}}{\dot{H}+2H^{2}},$
(19)

whence, the corresponding Friedmann equations read

$\displaystyle H^{2}-\frac{\kappa\zeta}{3}H^{n}=\frac{\kappa}{3}\rho,$
$\displaystyle H^{2}+\frac{2}{3}\dot{H}-\frac{\kappa\zeta}{3}H^{n}=0.$ (20)

Defining $\Omega=8\pi G\rho/3H^{2}$, while $\Omega_{0}$ denotes its current
value macq , the evolution of $q$ and $G/G_{0}$ have been plotted in Fig. (2).
For the employed parameters, transition redshift ($z_{t}$) lies within the
range of $0.4\leq z_{t}\leq 0.88$. The sensitivity of diagrams to the values
of $\Omega_{0}$ and $H_{0}$ is so weak compared with those of $\zeta$ and $n$
and $\kappa$. Indeed, although we only consider a baryonic source for current
density parameter $\Omega_{0}=0.049$ macq , and $H_{0}=67.66$ agh , the
obtained behaviors are also achievable for other candidates of $\Omega$ (such
as dark matter) and also the other values of $H_{0}$, reported in the
literature. Hence, suitable behavior of $q$ is obtainable by only considering
the baryonic content of the universe, meaning that the $\zeta H^{n}$ term may
play the role of the unknown parts (dark components) of cosmos. Dirac
hypothesis is also respected during the cosmic evolution. Remarkably, $G$ will
take negative values in future meaning that gravity will become repulsive
which speeds the universe expansion rate up more i.e., $q$ decreases. All
these happen under the shadow of the existence of non-minimal coupling
$\lambda$ which varies during the evolution of the universe. In Fig. (3),
$H(z)$ far  and the distance modulus ama  are plotted for the $\Lambda$CDM
model and also our model.

Figure 2: The evolution of $q$ and $G/G_{0}$ assuming $w=0$, for the case
discussed in Sec.V. The diagrams for $G/G_{0}$ are plotted using the same
model parameters as of $q$ diagrams.

Figure 3: The evolution of $H(z)$ and $\mu(z)$ whenever $w=0$, for the case
discussed in Sec.V. The same values of parameters as of Fig. 2 have been used.
The black dashed lines show $H(z)$ and $\mu(z)$ for $\Lambda$CDM model.

The negative value of $G$ is the direct result of the assumed $\lambda$, and
changes in the values of model parameters do not affect this result. There are
also other works that predict negative values for $G$ bar2 ; uza1 ; uza2 .
Theoretically, our model shows that a non-minimal coupling between geometry
and matter fields can accelerate the universe expansion and has an ability to
satisfy Dirac hypothesis.

##  VI concluding remarks

After addressing some properties of previously cosmological models introduced
in the framework of GRT mor1 ; das , the implications of GRT on obtaining
varying $G$ has been studied through considering the Newtonian limit of the
field equations. Thereinafter, following a proposal of Dirac hypothesis
introduced in gaz ; smo , the required non-minimal coupling needed to support
Dirac model was also obtained. Our results show that the dark sectors of
cosmos can be unified into one cosmic fluid which behaves as a pressureless
fluid in high redshift limit, and also accelerates the universe in line with
the current observations (Fig. 1). We also proposed a non-minimal coupling
(Sec. V) which can play the role of dark side of cosmos satisfying Dirac
hypothesis. Indeed, the present study addresses a deep connection between non-
minimal coupling (between the matter fields and geometry) and the idea of
variable $G$. This translates into saying that one may find the footprints of
non-minimal coupling between the matter fields and geometry by having the
observationally confirmed profile of $G$ and conversely.

Although relying on the Rastall hypothesis on relation between the changes in
spacetime curvature and violation of OCL we only focused on the implications
of the violation of OCL in cosmology and its connection with Dirac hypothesis,
the OCL violation can also be allowed due to the quantum considerations such
as uncertainty principle, and in the framework of unimodular gravity producing
significant cosmological outcomes jos . Indeed, even in the framework of GR
and thanks to the Bianchi identity, OCL is violated as the result of the
existence of a non-constant $G$. In summary, it was our goal to address $i$)
probable connection between Dirac hypothesis and non-minimal couplings, and
simultaneously, $ii$) the ability of such couplings in being responsible for
the unknown parts (dark sides) of cosmos. Therefore, such couplings need to be
further studied from both of the theoretical and observational viewpoints.
Finally, we would like to mention that, though Rastall gravity and its
generalizations provide interesting results, cosmological models based on this
theory need to be accurately tested by observations. In the present model, we
tried to explore theoretical consequences of a varying G cosmology based on
GRT and also briefly examined observational aspects of the theory. However, a
full observational treatment of the present model, e.g., in light of
Akarsu2020 , needs to be done and work along this line can be considered as an
interesting subject for future studies and developments.

## References

  * (1) Dirac P. A. M., Nature 139 (1937), 323
  * (2) Dirac P. A. M.,  Proc. Roy. Soc. London, Ser. A 165 (1938), 199.
  * (3) Dirac P. A. M., Nature, 139, (1937), 1001.
  * (4) Weyl H. Ann. Phys., 59, 129 (1919);
Zwicky F., Phys. Rev., 55, 726 (1939);
Eddington A. S., The Mathematical Theory of Relativity, Cambridge University
Press, London (1923).

  * (5) Barrow J. D., Tipler F. J., The Anthropic Cosmological Principle, Oxford University Press, Oxford, (1986);
J. D. Barrow, Varying G and Other Constants, In: S$\acute{a}$nchez N.,
Zichichi A. (eds) Current Topics in Astrofundamental Physics: Primordial
Cosmology. NATO ASI Series (Series C: Mathematical and Physical Sciences), vol
511. Springer, Dordrecht (1998).

  * (6) Chandrasekhar S., Nature 139 (1937), 757;
Kothari D. S., Nature, 142 (1938), 354.

  * (7) S. Ray, U. Mukhopadhyay, S. Ray, A. Bhattacharjee, Int. Journal Mod. Phys. D 28 (2019), 1930014.
  * (8) Rastall P., Can. J. Phys. 54 (1976), 66
  * (9) Vinti J. P., Celestial Mechanics 16 (1977), 391
  * (10) De Sabbata V., Acta Cosmologica Zesz. 9 (1980), 63
  * (11) Baptista J. P., Batista A. B., Fabris J. C., Revista Brasileira de Fisica. 14 (1984), 208
  * (12) Beesham A., Int. J. Theo. Phys. 25 (1986), 1295
  * (13) Wu Y. S., Wang Z., Phys. Rev. Lett. 57 (1986), 16
  * (14) Degl’Innocenti S. et al., A&A 312 (1996), 345
  * (15) Barrow J. D., Mon. Not. R. Astron. Soc. 282 (1996), 1397
  * (16) Barrow J. D., 1997, arXiv:gr-qc/9711084.
  * (17) Barrow J. D., The Constants of Nature, (Vintage Books, London, 2002)
  * (18) Mansouri R., Nasseri F., Khorrami M., Phys. Lett. A 259 (1999), 194
  * (19) Gaztañaga E. et al., Phys. Rev. D 65 (2001), 023506
  * (20) Clifton T., Mota D., Barrow J. D, Mon. Not. R. Astron. Soc. 358 (2005), 601
  * (21) Bronnikov K. A., Kononogov S. A., Metrologia 43 (2006), 1
  * (22) Solà J.,  J. Phys. A: Math. Theor. 41 (2008), 164066
  * (23) Uzan J. P., Rev. Mod. Phys. 75 (2003), 403
  * (24) Uzan J. P., Liv. Rev. Relativ. 14 (2011), 2
  * (25) Smolin L., Class. Quantum Grav. 33 (2016), 025011
  * (26) Fritzsch H., Solà J., Nunes R. C., Eur. Phys. J. C 77 (2017), 193
  * (27) Leszczyńska K., Da̧browski M. P., Denkiewicz T., Eur. Phys. J. C 79 (2019), 222
  * (28) Ellis G. F. R., Maartens R., Maccallum M. A. H., The Constants of Nature, (Cambridge University Press, UK, 2012).
  * (29) Canuto, V., Adams, P. J., Hsieh, S. H., Tsiang, E., Phy. Rev. D 16, 6 (1977);
Wesson, P., Goodson, R. E., Observ. 101, 105 (1981).

  * (30) C. Brans, R. H. Dicke, Phys. Rev. 124, 925 (1961).
  * (31) R. H. Dicke, Phys. Rev. 125, 2163 (1962).
  * (32) R. H. Dicke, Rev. Mod. Phys. 29, 355 (1957).
  * (33) R. H. Dicke, Nature 192, 440 (1961).
  * (34) Poisson E., A Relativist’s Toolkit, (Cambridge University Press, UK, 2004).
  * (35) Nojiri S., Odintsov S. D., Phys. Lett. B 599 (2004), 137
  * (36) Allemandi G. et al., Phys. Rev. D 72 (2005), 063505
  * (37) Koivisto T.,Class. Quant. Grav. 23 (2006), 4289
  * (38) Bertolami O. et al., Phys. Rev. D 75 (2007), 104016.
  * (39) Harko T., Lobo F. S. N., Galaxies 2 (2014), 410.
  * (40) Carloni S., Phys. Lett. B 766 (2017), 55.
  * (41) Boehmer C. G., Carloni S.,Phys. Rev. D 98 (2018), 024054.
  * (42) Rastall P., Phys. Rev. D 6 (1972), 3357.
  * (43) Moradpour H. et al., The European Physical Journal C, 77 (2017), 259.
  * (44) De Moraes W. A. G., Santos A. F., Gen. Relativ. Gravit. 51 (2019), 167.
  * (45) Li R. et al., Mon. Not. R. Astron. Soc. 486 (2019), 2407.
  * (46) Al-Rawaf A. S., Taha O. M., Phys. Lett. B 366 (1996), 69.
  * (47) Al-Rawaf A. S., Taha O. M., Gen. Relat. Gravit. 28 (1996), 935.
  * (48) Al-Rawaf A. S., Int. J. Mod. Phys. D 14 (2005), 1941.
  * (49) Majernik V., Gen. Relat. Gravit. 35 (2003), 1007.
  * (50) Arbab A. I., J. Cosmol. Astropart. Phys. 05 (2003), 008.
  * (51) Abdel-Rahman A. M. M., Astrophys. Space. Sci. 278 (2001), 383.
  * (52) Abdel-Rahman A. M. M., Hashim M. H. A., Astrophys. Space. Sci. 298 (2005), 519.
  * (53) Abdel-Rahman A. M. M., Riad I. F., Astron. J. 134 (2007), 1931.
  * (54) Moradpour H. et al., Phys. Rev. D. 96 (2017), 123504.
  * (55) Manna T., Rahaman F., Mondal M., Mod. Phys. Lett. A 35 (2020), 2050034.
  * (56) S. K. Maurya and F. T.-Ortiz, Phys. Dark Univ. 29 (2020), 100577.
  * (57) H. Shabani and A. H. Ziaie, Europhysics Letters 129, (2020) 20004.
  * (58) Fabris J. C., Kerner R., Tossa J.,  Int. J. Mod. Phys. D 9 (2000), 111.
  * (59) Josset T., Perez A.,  Phys. Rev. Lett. 118 118 (2017), 021102.
  * (60) Watson S. et al., J. Cosmol. Astropart. Phys. 07 (2017), 11.
  * (61) Gamboa J. et al., Phys. Rev. D 96 (2017), 083534.
  * (62) Das D., Dutta S., Chakraborty S., Eur. Phys. J. C 78 (2018), 810.
  * (63) Lin K., Qian W. L., Eur. Phys. J. C 80 (2020), 561.
  * (64) C. E. M. Batista, M. H. Daouda, J. C. Fabris, O. F. Piattella, D. C. Rodrigues, Phys. Rev. D 85, (2012), 084008.
  * (65) Moradpour H. et al., Mod. Phys. Lett. A 32 (2017), 1750078
  * (66) Moradpour H. et al., Adv. High Energy Phys. 2018 (2018), 7124730
  * (67) Moradpour H., Sadeghnezhad N., Hendi S. H., Can. J. Phys. 95 (2017), 1257
  * (68) T. R. P. Carames. et al., Eur. Phys. J. C74 (2014) 3145.
  * (69) Darabi F. et al., Eur. Phys. J. C 78 (2018), 25
  * (70) Moradpour H. et al., Mod. Phys. Lett. A 33 (2019), 1950096
  * (71) Mota C.E., et al., arXiv:2007.01968.
  * (72) A. H. Ziaie, H. Moradpour, H. Shabani, Eur. Phys. J. Plus 135 (2020), 916.
  * (73) D. Lohiya, A. Batra, S. Mehra, S. Mahajan and A. Mukherjee, Astron. Astrophys. Trans. 14 (1997), 199.
  * (74) I. Ayuso, J. P. Mimoso and N. J. Nunes, Galaxies, 7 (2019), 38.
  * (75) A. B. Batista, J. C. Fabris and S. V. B. Goncalves, Class. Quant. Grav. 18 (2001), 1389.
  * (76) A. A. Starobinskij, Pisma v Astronomicheskii Zhurnal, 7, (1981), 67; Soviet Astronomy Letters, 7, (1981), 36, Translation.
  * (77) K. Lin and W.-L. Qian, Eur. Phys. J. C 80 (2020), 561.
  * (78) Domínguez A. et al., Astrophys. J. 885 (2019), 137.
  * (79) Macquart J. et al., Nature 581 (2020), 391.
  * (80) Farooq O. et al., Astrophys. J. 835 (2017), 26.
  * (81) Aghanim N. et al., A&A 641 (2020), A6.
  * (82) Amanullah et al., Astrophys. J. 716 (2010), 712.
  * (83) O. Akarsu, N. Katirci, S. Kumar, R. C. Nunes, B. Ozturk and S. Sharma, Eur. Phys. J. C 80 (2020), 1050.