# Pseudo-labelling and Meta Reweighting Learning for Image Aesthetic Quality
Assessment

Xin Jin, Hao Lou, Huang Heng, Xiaodong Li*, Shuai Cui, Xiaokun Zhang, Xiqiao
Li This work is partially supported by the National Natural Science Foundation
of China (62072014), the Beijing Natural Science Foundation (L192040), and the
CAAI-Huawei MindSpore Open Fund (CAAIXSJLJJ-2021-022A). Xin Jin is with the
Department of Cyber Security, Beijing Electronic Science and Technology
Institute, Beijing, China, and also with the Beijing Institute for General
Artificial Intelligence (BIGAI), Beijing, China. E-mail:
<EMAIL_ADDRESS>Lou, Heng Huang, Xiaokun Zhang and Xiqiao Li are
with the Department of Cyber Security, Beijing Electronic Science and
Technology Institute, Beijing, China.Xiaodong Li is with the Department of
Cyber Security, Beijing Electronic Science and Technology Institute, Beijing,
China. Xiaodong Li* is the corresponding author. E-mail:
<EMAIL_ADDRESS>Cui was with the Department of Philosophy and
Mathematics, University of California, Davis, CA, USA.

###### Abstract

In the tasks of image aesthetic quality evaluation, it is difficult to reach
both the high score area and low score area due to the normal distribution of
aesthetic datasets. To reduce the error in labeling and solve the problem of
normal data distribution, we propose a new aesthetic mixed dataset with
classification and regression called AMD-CR, and we train a meta reweighting
network to reweight the loss of training data differently. In addition, we
provide a training strategy acccording to different stages, based on pseudo
labels of the binary classification task, and then we use it for aesthetic
training acccording to different stages in classification and regression
tasks. In the construction of the network structure, we construct an aesthetic
adaptive block (AAB) structure that can adapt to any size of the input images.
Besides, we also use the efficient channel attention (ECA) to strengthen the
feature extracting ability of each task. The experimental result shows that
our method improves 0.1112 compared with the conventional methods in SROCC.
The method can also help to find best aesthetic path planning for unmanned
aerial vehicles (UAV) and vehicles.

###### Index Terms:

image aesthetic quality assessment, meta reweighting network, pseudo labels,
aesthetic adaptive block, path planning.

##  I Introduction

Aesthetics was raised by the German philosopher Bavmgarten in 1750. It focuses
on the aesthetic relationship between the public and the world. The aesthetic
research subjects are the aesthetic activities, an abstract concept, with
strong subjectivity. Only the experts can evaluate the abstractly aesthetic
consensus in real life. However, the public can not give the accurately
aesthetic evaluation without professional training, which leads to
computational aesthetics [1].

The goal of computational aesthetics is aesthetic intelligence that computes
and robots are able to recognize, produce, and create beauty. That can make
computers and robots focus on human cognition and the inner mechanism of
recognizing beauty. Computational visual aesthetics [2] is an intelligent
comprehension about visual information by visual computation techniques.
Nowadays, visual computation aesthetics depends on deep learning. In the
related researches, the visual computation aesthetics, mainly by training
massive data, get neural network modal. So, the model can give the evaluation
of aesthetic quality. Although the generally existing aesthetic benchmark
dataset like AVA [3], AADB [4], PCCD [5] and others have numeric score labels,
most of them are not over 20,000 images. AVA is the largest dataset with
25,5530 images, but the image of over 7 and under 3 only take 4.4% in the
overall score. In the training process, that can trigger seriously
underfitting problems. AADB includes 9,958 images, and each image have 11
aesthetic attribute scores. However, the aesthetic labels and amount are too
few to have good labeling quality. CHUKPQ [6] only have the simply two-value
label instead of the numeric label. The above datasets have obvious
disadvantages. Moreover, comparing with target detection, sematic
segmentation, semantic classification, and other tasks, it is more challenging
to extract aesthetic features in detail. That brings a lot of challenges for
evaluation of aesthetic quality. It is hard to extract lighting, color,
composition, and other features in aesthetics. That brings serious problems
for computational aesthetics.

In order to solve those problems, we propose a method of aesthetic quality
evaluation based on pseudo-labelling and meta-learning. Specifically, we make
the following four contributions:

  1. •

We filtering and constructing an Aesthetic Mixed Dataset with Classification
and Regression (AMD-CR) with an appropriate scale. AMD-CR have both more
reasonable data distribution and better aesthetic quality labels than the
classic aesthetic dataset.

  2. •

We introduce meta-learning in aesthetic sample re-weight. We use meta
reweighting network to provide appropriate weights of the training samples.

  3. •

We propose Aesthetic Adaptive Block (AAB) to compose aesthetic network model
that is used to compare images, with ratio between length and width, for
effectively evaluation.

  4. •

To the best of our knowledge, we are the first to use binary classification
training pseudo-label in the aesthetic tasks; and, our method shows good
results in image aesthetic quality evaluation and path planning for unmanned
aerial vehicles (UAV) and vehicles.

The rest of this paper is organized as follows: in Section II we will
introduce the related work of image aesthetic quality assessment and meta
learning. In Section III we will show the network model and AMD-CR. In Section
IV, we will describe the experiments in detail and analyze the results.
Finally, Section V will focus on conclusion and future work.

##  II Related work

Figure 1: The architecture of the network. The input image size is arbitrary
and the main network contains EfficienNet-B4 [8]. AAB is the aesthetic
adaptive block, GAP is the global average pooling, $\sigma$ is the sigmoid
function, meta reweighting network contains a hidden layer of 100 nodes, and
each layer will be activated by the relu function.

We divide the discussions of related work into the following three
subsections.

###  II-A Image aesthetic datasets

The most common benchmark aesthetic datasets are AVA [3] and AADB [4]. There
are over 250,000 images in AVA. Each image has one numeric aesthetic quality
evaluation label, 66 sematic labels, and 14 typed labels. AADB comes from
University of California, Irvine. Kong and others designed a new image
aesthetic attribute dataset with total 9,958 images. This dataset includes
images from both professional photographers and ordinary workers. In fact,
there are obvious disadvantages in both datasets. AVA dataset both high score
and low score images take less than 5%. Although AADB includes 11 aesthetic
attributes, including balancing elements, colorful harmony, and others, its
labels are few and subjective. The limitations of datasets stop the
development of image aesthetic quality evaluation.

###  II-B Image aesthetic quality assessment

Based on the definition of the classification task in basic aesthetic
evaluation tasks, there are two categories of images: high quality and low
quality [10]. However, the judgements are usually based on distinguished
aesthetic common sense. If two images are both ordinary images, it will be
hard to distinguish them by binary values. So, we need continuously numeric
value to evaluate them. In the early stages, researchers designed many
handcrafted aesthetics features [11, 12, 13, 14], Luo et al. [6] divided
1,7613 images to 7 categories based on semantics; they design a series image
detection salient regions method according to different image material. They
use machine learning algorithm to train classification to extract features and
predict classified results.

Nowadays, image aesthetic features mainly use deep learning network model [15]
to extract automatically. The following researches [16, 17] directly or
indirectly model the extra information, such as style; the classification of
this model is more accurate than the model based on handcraft. Then,
researchers use probability distribution [18] to describe aesthetics and
explain the subjectivity of aesthetic evaluation. However, the above single
image aesthetic vector is hard to quantize complex aesthetic evaluation.
Recently, aesthetic attributes evaluation field proposes a method based on
hierarchical multi-task network [19], incremental learning of multi-task
network [20]. This method quantize aesthetic evaluation by contrasting
composition, light, color, matching, exposure, depth of field, camera use.
However, all researches are limited by the amount of data and data
distribution. Our research [9] has great improvement in binary classification
of aesthetic evaluation and numeric evaluation of regression tasks. Beside the
application of basic aesthetic evaluation, we propose a multistream network,
composed by spatial, motion, and structural streams, to gain the multimodal
features of path planning. That inspires us to explore path planning based on
image aesthetic quality evaluation.

###  II-C Meta learning

In deep learning, the modal has to be retrained by changing scene. Meta
learning, learn to learn, can learn the ability of human learning new things
by limit data [21]. Meta learning mainly includes memory [22], prediction of
gradient [23], attention mechanisms [24], LSTM [25], reinforcement learning
[26] and others. By the inspiration of meta learning [25, 21], the recent
researches study reweight cases adaptatively from data sample. That shows
learning is more automatic and dependent. This idea can effectively solve the
longtail problem caused by few data taggers.

##  III Resarch methodology

In this section, we mainly introduce the structure and the algorithm of the
aesthetic assessment model in detail. The network structure is constructed in
Fig.1.

###  III-A Meta reweighting network

The loss functions of conventional deep learning are generally cross entropy
or mean square error (MSE) loss function, which calculate the loss of average
on batch size samples; this method is easily influenced by the sample
distribution. In the classificational task, if the proportion of the positive
samples is high in samples, then the final transmission of the loss value must
be influenced by most of the positive sample loss. Although a few negative
samples have higher loss, the network will be less influenced by these
negative samples. The idea of the meta reweighting in this paper is to learn
the loss weight of the data automatically, optimize the normal data
distribution, and enhance the model robustness.

By the meta reweighting method mentioned in [27], we apply this method to the
aesthetic quality assessment tasks. The idea of meta reweighting is to
minimize the loss $\frac{1}{N}\sum_{i=1}^{N}\ell\left(y_{i},{f}(x_{i})\right)$
in the training set and extract the best classifier parameters $w^{*}$. In the
biased training set, each sample weight is expressed as
$\mathcal{V}\left(L_{i}^{train}(w);\mathrm{\Theta}\right)$,
$L_{i}^{train}\left(w\right)=\ \ell\left(y_{i},{f}(x_{i})\right)$, and it
means the loss of _i_ th training sample; _$\Theta$_ represents the parameter
of the meta reweighting network(MRN). The best parameter $w^{*}$ is calculated
as Equation (1).

$w^{\ast}(\Theta)=\frac{1}{N}\sum_{i=1}^{N}\mathcal{V}\left(L_{i}^{\mathrm{train}}\left(w\right);\Theta\right)L_{i}^{\mathrm{train\
}}\left(w\right)$ (1)

In order to learn hyperparametric $\Theta$ automatically, we generate a
multilayer perceptron based on the hidden layer in the method of meta
learning, shown as the yellow network in Fig.1. The input of the meta
reweighting network is the loss of samples, and the output is the weight of
the loss of samples. We extract a few high quality images from the training
set _S train_ to construct the meta set _S meta_. _N_ is the amount of _S
train_, and _M_ is the amount of _S meta_, _N_ $\gg$ _M_. The advantages of _S
meta_ are correct labels and equal data distribution. According to the score
labels of aesthetic quality assessment datasets, we divide the data into 10
segments equally; we extract 200 images from each segment artificially. If the
amount of any segment is less than 200, we will choose appropriate amount of
images from the adjacent segments to replenish the lack images. Finally, we
filter 2000 images fo construct _S meta_. There are two stages in the training
process; we update the _$\Theta$_ and $w$ through a single optimization cycle.
We firstly update the parameter _$\Theta$_ of the MRN, and then update the
parameter $w$ of the main network in each iteration optimization. The update
equation for the first stage is defined as Equation (2), _$\beta$_ is the step
size of the MRN.

$\Theta^{(t+1)}=\Theta^{(t)}-\left.\beta\frac{1}{m}\sum_{i=1}^{m}\nabla_{\Theta}L_{i}^{meta}\left({\hat{w}}^{(t)}(\Theta)\right)\right|_{\theta^{(t)}}$
(2)

$\left.\nabla_{\Theta}L_{i}^{meta}\left({\hat{w}}^{(t)}(\Theta)\right)\right|_{\theta^{(t)}}$
represents the gradient of the MRN on meta set,
${\hat{w}}^{(t)}\left(\theta\right)$ is calculated as Equation (3), _$\alpha$_
is the step size of the main network.

$\displaystyle{\hat{w}}^{\left(t\right)}\left(\theta\right)=$ (3)
$\displaystyle
w^{\left(t\right)}-\left.\alpha\frac{1}{n}\sum_{i=1}^{n}v\left(L_{i}^{train}\left(w^{\left(t\right)}\right);\Theta\right)\nabla_{\mathbf{w}}L_{i}^{\mathrm{train}}\left(w\right)\right|_{w^{\left(t\right)}}$

In the first stage, we need to calculate the weighted loss of _S train_ on the
main network. In order to distinguish real training from meta-learning
training, we copy the parameters of the main network at each single cycle
training. The parameters of MRN are optimized by Equation (2).

In the second stage, we optimize and update the main network according to the
loss of _S train_ and the weight from MRN on the meta network, it is as
Equation (4).

$\displaystyle w^{(t+1)}=w^{(t)}-$ (4)
$\displaystyle\left.\alpha\frac{1}{n}\sum_{i=1}^{n}\nu\left(L_{i}^{train}\left(w^{\left(t\right)}\right);\Theta^{\left(t+1\right)}\right)\nabla_{\mathbf{w}}L_{i}^{train}\left(w\right)\right|_{w^{\left(t\right)}}$

###  III-B Aesthetic adaptive block

Figure 2: (a) is the image for cropping, and (b) is for resizing. Figure 3:
The structure of the aesthetic adaptive block, AAP is the adaptive average
pooling layer.

For the pretraining of images, there are three ways: cropping, resizing and
adaptive padding. The best input size of images in EfficientNet-B4 is 380×380.
For cropping, we adjust any size images to make the short edge reach 380
proportionally, and then crop the center of images into the best input size.
For resizing, we stretch any size images to 380×380 instead of cropping. Both
methods may lose a part of the aesthetic features. We propose a structure,
aesthetic adaptive block (AAB), based on adaptive padding; this method can
adapt any input size without cropping and resizing. In AAB, we adjust the long
edge of the image to 800 proportinally, and then we fill the short edge with 0
to make it reach 800. So, we construct an 800x800 input image. Then, for the
first convolution operation of the EfficientNet-B4, we use an adaptive average
pooling layer to reshape the feature map of 48×190×190; 48 is the amount of
the feature channel. Finally, the feature map enters the remaining network
layers in EfficientNet-B4. Three methods are shown as Fig.2 and Fig.3.

###  III-C Eifficent channel attention

Recently, the channel attention mechanism [28] shows a great potential in
improving the performance of deep learning network models. However, most of
the previous focused on building complex attention modules to achieve good
performance. In order to reduce the complexcity of models, this paper cites
the ECA module [7]. It adopts a local cross-channel attention mechanism and
the adaptive one-dimensional convolution kernel, which can improve performance
obviously.

The ECA captures local cross-channel attention information by interacting with
each channel and its _k_ neighbors. This can be effectively implemented by a
_1D_ convolution with the kernel size of _k_ , where _k_ represents the local
cross-channel attention and means the number of neighbor nodes involved in the
attention prediction of this channel. In order to avoid adjusting _k_ in the
validation phase, the ECA proposes a method that _k_ can be adaptively
proportional to the channel dimensions, _k_ is expressed as Equation (5).

$k=\psi(C)\
=\left|\frac{\log_{2}(C)}{\gamma}+\frac{b}{\gamma}\right|_{\mathrm{odd}}$ (5)

${|t|}_{odd}$ means the nearest odd number of _t_ , the value of _$\gamma$_ is
2, the value of _b_ is 1. We can get _k_ = 7 if _C_ = 1792.

###  III-D Multi-task network

Both the aesthetic multi-classification and regression tasks share the same
main network parameters to extract features. The multi-classification network
and the regression network structure are similar as shown in Fig.1. The main
network outputs the _N_ ×1792×11×11 feature map. _N_ is for the batch size.
The number of channels is first reduced to 448 using a convolution layer with
the 3×3 kernel, then the feature map is transformed into fully connected layer
in 448 dimensions by an adaptive average pooling layer. Finally, the last
layer has 10 nodes in the multi-classification network and 1 node in the
regression network. Loss function formula of multi-classification network is
expressed as Equation (6).

${Loss}_{class}=\sum_{x}p(x)logq(x)$ (6)

In this equation, _p(x)_ represents the predicted output probability, _q(x)_
represents the ground-truth probability, and the probability is obtained by
the following Equation (7).

$y_{i}=\frac{e^{i}}{\mathrm{\Sigma}_{j}e^{j}}$ (7)

Loss function formula of regression network is expressed as Equation (8).

$Loss_{regre}=\frac{1}{N}\sum_{i}{(s_{i}-\widehat{s_{i}})}^{2}$ (8)

_$s_{i}$_ is the output score, _$\hat{s_{i}}$_ is the ground-truth score.

###  III-E Piecewise strategy based on pseudo-labelling of the binary
classification task

The pseudo-labelling was firstly proposed in the [29], which are not
completely consistent with the ground-truth labels; and. the ground-truth
labels are often used as the hard labels. Pseudo labels can be considered as
soft labels after training. We propose a piecewise training method based on
pseudo labels of the binary classification task. For first stage, we train a
binary classifier, and the images are divided into _0_ and _1_ based on the
boundary of 5 points, while 0 represents low aesthetic quality and 1
represents high aesthetic quality. TABLE I shows the accuracy and error rates
of each segment in the testing set.

TABLE I: Correctness and error rates for classification Score Segment | Correctness rate | Error Rate
---|---|---
0.0-1.0 | 100% | 0%
1.0-2.0 | 100% | 0%
2.0-3.0 | 98.67% | 1.33%
3.0-4.0 | 92.85% | 7.15%
4.0-5.0 | 70.06% | 29.94%
5.0-6.0 | 75.73% | 24.27%
6.0-7.0 | 91.42% | 8.58%
7.0-8.0 | 96.01% | 3.99%
8.0-9.0 | 97.01% | 2.99%
9.0-10.0 | 100% | 0%

Although there are still a few false classification samples in 2.0-3.0 and
8.0-9.0 intervals, the error rate in the middle interval 4.0-6.0 is higher; it
is because, in the middle interval, there are no huge difference in the scores
of images in all perspectives of aesthetic quality assessment. However, this
does not influence the final assessment of image aesthetics by the network
model. Through continuous score regression, the assessment measures of the
model can be measured by quantitative errors. If the scores of data is
subjective labels, the result of the classification task by the binary
classifier can be regarded as pseudo labels of the model. In the second stage,
the pseudo-labelling redivides the datasets into two datasets. We performed
multi-classification training and fine-grained regression training separately
to extract better image aesthetic features. The experimental results in
Section IV show that the piecewise strategy helps to improve the accuracy of
the regression.

###  III-F Data

In order to construct an aesthetic dataset with reasonable distribution and
high quality of aesthetic annotation, we filter and reconstruct a dataset
called the aesthetic mixed dataset with classification and regression (AMD-
CR). There are two datasets in AMD-CR: aesthetic mixed dataset with
classification (AMD-C) applied to the binary classification training task, and
aesthetic mixed dataset with regression (AMD-R) applied to the regression
task.

####  III-F1 Construction of AMD-CR

The part of aesthetic mixed dataset with classification and regression (AMD-
CR) comes from image aesthetic benchmark datasets: DPChallenge, Photo.net, AVA
[3], CUHKPQ [6] and SPAQ [30]; another part comes from self-built datasets:
GLAMOUR, OUTDOOR, NG and PSA. AVA and DPChallenge come from
www.dpchallenge.com. So, we need to delete the duplicate images. Photo.net is
from www.photo.net, GLAMOUR is from www.glamour-photos.org, OUTDOOR is from
www.outdoor-photos.com, NG is from www.dili360.com, PSA is from
www.psachina.org. We obtain a portion of the public images from these websites
and integrate them into our aesthetic dataset. We distribute continuous labels
to [0, 10.0] through standardized fractional labeling and map binary labels to
_0_ and _1_. _0_ represents the images with low aesthetic quality and _1_
represents the images with high aesthetic quality. Discrete datasets only have
binary labels, while continuous datasets both have binary labels and contiuous
labels. We focused on selecting images of high and low quality, and we
performed a preliminary selection of the images with general aesthetic quality
to balance the data distribution. TABLE II shows the constitution of the
aesthetic mixed dataset with classification and regression (AMD-CR).

TABLE II: The constitution of AMD-CR Datasets | Volume | labels
---|---|---
DPChallenge + AVA | 61461 | Continuous
CUHKPQ | 29690 | Discrete
Photo.net | 8985 | Continuous
NG | 4140 | Discrete
SPAQ | 4102 | Continuous
PSA | 3460 | Discrete
GLAMOUR | 2701 | Discrete
OUTDOOR | 1916 | Discrete

####  III-F2 Construction of AMD-R

We fliter the data with the continuous labels in the aesthetic mixed dataset
with regression (AMD-CR). The images of low and high segments with the scores
less than 4.0 points or more than 6.0 points are reserved, and the middle
segment images (4.0-6.0) are randomly sampled. Then we get 59,371 images form
AMD-R. The segment distribution of AMD-R is shown in Fig.4.

Figure 4: The Segment Distribution Of AMD-R

####  III-F3 Construction of AMD-C

We removed the images ranging between 4.0 points and 6.0 points in the AMD-CR
to ensure that the ratio of positive and negative samples is 1:1. Finally, we
obtain the aesthetic mixed dataset with classification (AMD-C) containing
61,660 images.

##  IV Experiment

###  IV-A Training details

We set the classification batch size to be 32, the regression batch size to be
64 and the learning rate to be 0.0001. We use Adam as the optimizer; betas are
set as (0.98, 0.999); weight decay is set as 0.0001. If either the accuracy
rate of classification is not improved or regression loss is not decreased in
two consecutive rounds in validation datasets, the learning rate will multiply
by 0.5. Our running environment is in Pytorch 1.5.0 and Nvidia TITAN XPs.

The datasets used in this paper are AMD-C and AMD-R. We divide the dataset
into three sets; the ratio of training set and validation set and testing set
is 8:1:1. EfficientNet-B4 is the main network in the model. After we first
trained the binary classifier, we will get the training model as _C 2_. Then
the training set and validation set are divided into _S train0_, _S train1_,
_S valid0_, and _S valid1_ according to the classifier, _0_ represents the low
quality and _1_ represents the high quality in aesthetics. Then we use _S
train0_, _S valid0_ and _S train1_, _S valid1_ to train classification and
regression models for the next stage.

In the aesthetic assessment training. We start with a 10-class classification
training, with [0,10.0] divided into ten average segments, corresponding to
the category _A_ when the scoring label is in the (_A_ , _A_ +1.0] interval,
_A_ $\in$ (0,1,2,3,4,5,6,7,8,9). If there is an image with 0 score, it is
classified as the _0_ category. When we perform the ten-class trainning task
with the images of the same pseudo labels obtained by the binary classifier,
the parameters of the regression network should not be returned. We use cross-
entropy function to get classification loss. Then, we release the regression
network and freeze the main network, classification network, and MRN to train
the regression model, we can get models _R 0_ and _R 1_ corresponding to the
_S train0_ and _S train1_ respectively. In the same method, we can get
regression model _R all_ training on all the AMD-R. Finally we can get the
result score like this: if $C_{2}(i)$ = 0,
$score=\frac{R_{0}(i)+R_{all}(i)}{2}$,else
$score=\frac{R_{1}(i)+R_{all}(i)}{2},i\in S_{test}$.

Figure 5: Test samples. _G_ for the ground-truth score, _M_ for the score of
MRN+AAB+ECA(PCR), _A_ for the score of AAB+ECA(PCR), _R_ for the score of
Resizing+ECA(CR) and _C_ for the score of Cropping+ECA(CR). Figure 6: Path
planning samples. The better path is identified in bold and the best capture
in this path is identified by a yellow star.

###  IV-B Analysis of experimental results

We use these indicators for assessment: first, we calculate mean square error
(MSE): calculating the loss between the output score and the ground-truth
score, it is as Equation (9):

$MSE=\frac{1}{N}\sum_{i}{(r_{i}-\widehat{r_{i}})}^{2}$ (9)

Mean absolute error (MAE) measures the differences between ground-truth and
regressed scores directly, and then find the mean of differences, it is as
Equation (10):

$MAE=\frac{1}{N}\sum_{i}\left|r_{i}-\hat{r}_{i}\right|$ (10)

We caluculate Spearman’s rank order correlation coefficient (SROCC): the
correlation between the output score and the ground-truth score, it is as
Equation (11):

$SROCC=1-\frac{6\sum_{i}{(r_{i}-\widehat{r_{i}})}^{2}}{N^{3}-N}$ (11)

Accuracy indicates whether the prediction score is consistent with the ground-
truth score in the binary classification task when the dividing line is _5_ ,
it is as Equation (12):

$ACCURACY=\frac{TP+TN}{P+N}$ (12)

$ACCURACY_{\left|error\right|\leq 1}$ indicates whether the absolute value of
the error between the output score and the ground-truth score is within 1
point, it can be expressed as Equation (13):

$ACCURACY_{\left|error\right|\leq 1}=\frac{N_{\left|error\right|\leq 1}}{N}$
(13)

Several experiments were conducted as follows: for training strategies, we
record _R_ for the single regression training method, _CR_ for the
classification before regression training method, and _PCR_ for classification
before regression based on pseudo-labelling training in TABLE III. We verify
the goodness of guidance for classification to regression firstly. We use the
cropping method in Section III; we can obviously compare _Cropping(R)_ with
_Cropping(CR)_ to find that if classification before regression will get
better result on the testing set. Then we verify that the ECA plays a great
role in the aesthetic assessment task. The comparison experiments of
_Cropping(CR)_ and _Cropping+ECA(CR)_ can verify that the ECA can effectively
improve the accuracy of the assessment. Besides, we can verify the training
method of pseudo labels can improve the result by _AAB+ECA(PCR)_. Finally, we
can also verify the meta reweighting network can improve the result by
_MRN+AAB+ECA(PCR)_ , all the ablation experiments results are shown in Table
III.

TABLE III: Ablation experimental results Method |  MSE↓ | SROCC↑ |  Acc↑ |  $Acc_{\left|error\right|\leq 1}$↑
---|---|---|---|---
MRN+AAB+ECA(PCR) |  0.6691 | 0.8471 |  86.53% | 80.61%
AAB+ECA(PCR) |  0.6888 | 0.8409 |  86.48% | 79.80%
AAB+ECA(CR) |  0.7575 | 0.8238 |  84.58% | 78.73%
Resizing+ECA(CR) |  0.7684 | 0.8213 |  84.82% | 78.23%
Cropping+ECA(CR) |  0.8213 | 0.8117 |  84.43% | 76.54%
Cropping(CR) |  1.0325 | 0.7743 |  79.23% | 69.58%
Cropping(R) |  1.2520 | 0.7359 |  76.39% | 63.78%
TABLE IV: Comparisons on AADB dataset Method | MAE↓ | SROCC↑ | Acc↑
---|---|---|---
Reg-Net [4] | 0.1268 | 0.678 | -
RGNet [31] | - | 0.710 | -
Lee _et al._ [32] | 0.1141 | 0.879 | -
Ours | 0.1016 | 0.7185 | 80.5%
Figure 7: Comparison of scatter plots. 6(a) represents Cropping(CR), 6(b)
represents Cropping + ECA(CR), 6(c) represents Resizing + ECA(CR), 6(d)
represents AAB + ECA(CR), 6(e) represents AAB + ECA(PCR), and 6(f) represents
MRN + AAB + ECA(PCR). The horizontal axis represents the prediction scores,
and the vertical axis represents the ground-truth scores.

Fig.5 compares the test score results for four main methods. We select some
typical samples in each score segment: _G_ represents the ground truth; _M_
represents the score of _MAE(PCR)_ ; _A_ represents the score of _AE(PCR)_ ;
_R_ represents the score of _RE(CR)_ ; _C_ represents the score of _CE(CR)_.
It can be found that each segment of the dataset has obvious distinction, and
_MAE(PCR)_ method has the best ability for asssessment. In TABLE IV, by
comparison, our method has better performance in MAE and accuracy than others.

Fig.6 shows the path planning samples; we choose three different scenarios,
and each scenary has two different path planning strategies. It is easy for
our model to select one of the path with high aesthetic quality assessment.
Besides, we choose the image with the highest aesthetic score for each
selected path. These images are all better in color, light, composition and
other aesthetic photography features.

We can get the prediction of the distribution with different methods from the
scatter plots of Fig.7. Fig.7 shows that the best model can predict the scores
between 1 and 9 points, which is very close to the actual label score
interval; the prediction results of the high and low segments are ideal; (f)
shows a long and straight spinning cone, that means (f) has the best
distribution for aesthetic regression.

Figure 8: Images scaling in different sizes TABLE V: Prediction results of different image preprocessing methods Image size | Origin | AAB | Resizing | Cropping
---|---|---|---|---
1000×582 | 8.12 | 8.14 | 8.21 | 6.78
1600×582 | 8.12 | 7.67 | 8.21 | 6.97
1000×1782 | 8.12 | 7.80 | 8.29 | 7.33
250×146 | 8.12 | 7.83 | 8.02 | 7.23

Comparing with three different image preprocessing methods in Section III.B,
we can find that AAB method can save the information about the length and the
width in the original image. Resizing can save the whole picture information,
but, when the ratios between length and width are either extremly high or
extremely low, resizing may lead to the loss of aesthetic features. Cropping
only cuts the center into 380×380, which may lose some important image feature
information. As shown in the Fig.8.

Table V shows that the original image is stretched horizontally to 1600×592,
longitudinally to 1000×1792, and isometric scaling to 250×146. AAB model has
different output in different ratios of length and width. The scores of
resizing are closer, and the cropping model has low scores. It is because both
methods may neglect of local features.

##  V Conclusions

It is a challenge task to construct a new dataset in image aesthetic field. We
construct AMD-CR with reasonable classification regression by mixing and
filtering massive datasets; we proposed a new aesthetic evaluation model based
on meta reweight network and binary classification pseudo labels. This model
can effectively improve ability of predict image aesthetic. Moreover, we
design AAB structure to use in any input size images. We are also the first
team use ECA in aesthetics. This attention can improve ability of extracting
features without changing the number of channels. The best test set, SROCC, is
0.8471 that is higher than classic deep learning regression model.

Now, we only evaluate images generally/ In the future, we will improve
application aesthetic quality evaluation, Interpretability, by more image
attributes. We also consider how to use emotions, arts, themes, and other high
semantic information to help aesthetic evaluation tasks. Moreover, we hope to
explore more in the area of aesthetic guidance.

## References

  * [1] R. Datta, D. Joshi, J. Li, and J. Z. Wang, “Studying aesthetics in photographic images using a computational approach,” in _European conference on computer vision_. Springer, 2006, pp. 288–301.
  * [2] A. Brachmann and C. Redies, “Computational and experimental approaches to visual aesthetics,” _Frontiers in computational neuroscience_ , vol. 11, p. 102, 2017.
  * [3] N. Murray, L. Marchesotti, and F. Perronnin, “Ava: A large-scale database for aesthetic visual analysis,” in _2012 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 2012, pp. 2408–2415.
  * [4] S. Kong, X. Shen, Z. Lin, R. Mech, and C. Fowlkes, “Photo aesthetics ranking network with attributes and content adaptation,” in _European Conference on Computer Vision_. Springer, 2016, pp. 662–679.
  * [5] K.-Y. Chang, K.-H. Lu, and C.-S. Chen, “Aesthetic critiques generation for photos,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 3514–3523.
  * [6] W. Luo, X. Wang, and X. Tang, “Content-based photo quality assessment,” in _2011 International Conference on Computer Vision_. IEEE, 2011, pp. 2206–2213.
  * [7] Q. Wang, B. Wu, P. Zhu, P. Li, W. Zuo, and Q. Hu, “Eca-net: Efficient channel attention for deep convolutional neural networks, 2020 ieee,” in _CVF Conference on Computer Vision and Pattern Recognition (CVPR). IEEE_ , 2020.
  * [8] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 6105–6114.
  * [9] Q. Kuang, X. Jin, Q. Zhao, and B. Zhou, “Deep multimodality learning for uav video aesthetic quality assessment,” _IEEE Transactions on Multimedia_ , vol. 22, no. 10, pp. 2623–2634, 2019.
  * [10] R. W. Sperry, M. S. Gazzaniga, and J. E. Bogen, “Interhemispheric relationships: the neocortical commissures; syndromes of hemisphere disconnection,” 1969.
  * [11] H. Tong, M. Li, H.-J. Zhang, J. He, and C. Zhang, “Classification of digital photos taken by photographers or home users,” in _Pacific-Rim Conference on Multimedia_. Springer, 2004, pp. 198–205.
  * [12] Y. Luo and X. Tang, “Photo and video quality evaluation: Focusing on the subject,” in _European Conference on Computer Vision_. Springer, 2008, pp. 386–399.
  * [13] M. Nishiyama, T. Okabe, I. Sato, and Y. Sato, “Aesthetic quality classification of photographs based on color harmony,” in _CVPR 2011_. IEEE, 2011, pp. 33–40.
  * [14] L. Marchesotti, F. Perronnin, D. Larlus, and G. Csurka, “Assessing the aesthetic quality of photographs using generic image descriptors,” in _2011 international conference on computer vision_. IEEE, 2011, pp. 1784–1791.
  * [15] X. Lu, Z. Lin, H. Jin, J. Yang, and J. Z. Wang, “Rapid: Rating pictorial aesthetics using deep learning,” in _Proceedings of the 22nd ACM international conference on Multimedia_ , 2014, pp. 457–466.
  * [16] Y. Kao, C. Wang, and K. Huang, “Visual aesthetic quality assessment with a regression model,” in _2015 IEEE International Conference on Image Processing (ICIP)_. IEEE, 2015, pp. 1583–1587.
  * [17] X. Lu, Z. Lin, X. Shen, R. Mech, and J. Z. Wang, “Deep multi-patch aggregation network for image style, aesthetics, and quality estimation,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 990–998.
  * [18] O. Wu, W. Hu, and J. Gao, “Learning to predict the perceived visual quality of photos,” in _2011 International Conference on Computer Vision_. IEEE, 2011, pp. 225–232.
  * [19] X. Jin, L. Wu, X. Zhou, G. Zhao, X. Zhang, X. Li, and S. Ge, “Predicting aesthetic radar map using a hierarchical multi-task network,” in _Chinese Conference on Pattern Recognition and Computer Vision (PRCV)_. Springer, 2018, pp. 41–50.
  * [20] X. Jin, X. Zhou, X. Li, X. Zhang, H. Sun, X. Li, and R. Liu, “Incremental learning of multi-tasking networks for aesthetic radar map prediction,” _IEEE Access_ , vol. 7, pp. 183 647–183 655, 2019.
  * [21] C. Finn, P. Abbeel, and S. Levine, “Model-agnostic meta-learning for fast adaptation of deep networks,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 1126–1135.
  * [22] A. Santoro, S. Bartunov, M. Botvinick, D. Wierstra, and T. Lillicrap, “Meta-learning with memory-augmented neural networks,” in _International conference on machine learning_. PMLR, 2016, pp. 1842–1850.
  * [23] M. Andrychowicz, M. Denil, S. Gomez, M. W. Hoffman, D. Pfau, T. Schaul, B. Shillingford, and N. De Freitas, “Learning to learn by gradient descent by gradient descent,” _arXiv preprint arXiv:1606.04474_ , 2016.
  * [24] O. Vinyals, C. Blundell, T. Lillicrap, K. Kavukcuoglu, and D. Wierstra, “Matching networks for one shot learning,” _arXiv preprint arXiv:1606.04080_ , 2016.
  * [25] S. Ravi and H. Larochelle, “Optimization as a model for few-shot learning,” 2016\.
  * [26] N. Schweighofer and K. Doya, “Meta-learning in reinforcement learning.” _Neural Networks_ , vol. 16, no. 1, pp. 5–9, 2003.
  * [27] J. Shu, Q. Xie, L. Yi, Q. Zhao, S. Zhou, Z. Xu, and D. Meng, “Meta-weight-net: Learning an explicit mapping for sample weighting,” _arXiv preprint arXiv:1902.07379_ , 2019.
  * [28] J. Hu, L. Shen, and G. Sun, “Squeeze-and-excitation networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 7132–7141.
  * [29] G. Hinton, O. Vinyals, and J. Dean, “Distilling the knowledge in a neural network,” _arXiv preprint arXiv:1503.02531_ , 2015.
  * [30] Y. Fang, H. Zhu, Y. Zeng, K. Ma, and Z. Wang, “Perceptual quality assessment of smartphone photography,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 3677–3686.
  * [31] D. Liu, R. Puri, N. Kamath, and S. Bhattacharya, “Composition-aware image aesthetics assessment,” in _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision_ , 2020, pp. 3569–3578.
  * [32] J.-T. Lee and C.-S. Kim, “Image aesthetic assessment based on pairwise comparison a unified approach to score regression, binary classification, and personalization,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 1191–1200.

|  XIN JIN was born in Anhui, China. He received the Ph.D. degree from Beihang
University, China. He is currently an Associate Professor with the Department
of Cyber Security at Beijing Electronic and Science Technology Institute and a
Visiting Scholar at Beijing Institute for General Artificial Intelligence
(BigAI). He received his Ph.D. degrees in Computer Science from Beihang
University, China, in 2013. He was a visiting student at Lotus Hill Institute,
China and a visiting scholar at Tsinghua University, Beijing China. His
research interests include computational aesthetics, computer vision and
artificial intelligence security. He served as the program committee members
and session chairs of multiple conferences including AAAI 2020-2022, IJCAI
2020, ACM MM 2021, ISAIR 2017-2022, ICVRD2017-2020. He has received the Best
Paper Awards at ROSENET 2018/2019.
---|---
|  HAO LOU was born in Zhejiang, China. He received the bachelor’s degree from
Northwest Agriculture and Forestry University, China. He is currently pursuing
the master’s degree with the Beijing Electronic Science and Technology
Institute, China. His research interest includes visual media aesthetics.
---|---
|  HUANG HENG was born in Fujian, China. He received the bachelor’s degree
from Shangdong University, China. He is currently pursuing the master’s degree
with the Beijing Electronic Science and Technology Institute, China. His
research interest includes visual media aesthetics.
---|---
|  XIAODONG LI was born in Henan, China. He received the Ph.D. degree from
Northwestern Polytechnic University, China. He is currently an Associate
Professor with the Beijing Electronic Science and Technology Institute, China.
His researches are focused on information security and visual media security.
---|---
|  SHUAI CUI was born in Beijing, China. She received Bachelor of Arts and
Sciences degree in philosophy and mathematics from University of California,
Davis, CA, USA in 2021. Now, she is applying master programs. Her research
interests include logic, metaphysics, epistemology, and Artificial
Intelligence.
---|---
|  XIAOKUN ZHANG was born in Gansu, China. He is currently a Professor with
the Department of Cyber Security, Beijing Electronic Science and Technology
Institute, China. His research interest includes computer science and
technology.
---|---
|  XIQIAO LI was born in Shanxi, China. He received the master’s degree from
Beijing Electronic Science and Technology Institute, Beijing, China. His
research interest includes visual media aesthetics.
---|---


# Cosmological future singularities in massive gravity and massive bigravity

M. Mousavi<EMAIL_ADDRESS>K. Atazadeh<EMAIL_ADDRESS>Department of Physics, Azarbaijan Shahid Madani University, Tabriz, 53714-161
Iran

###### Abstract

We study the future cosmological singularities in the framework of massive
gravity and minimal massive bigravity theory. In this regards, we consider the
possible classes of finite-time future singularities such as sudden, big rip,
big freeze and big brake singularities in the massive universe. In dRGT model
with an open expanding universe we obtain the sudden singularity in the future
at a finite-time which generally without taking account of any particular
realistic equation of state, is not avoidable and except the fluid density,
all dynamical physical quantities such as pressure approach to infinity. To
complete our study, we search the future cosmological singularities in the
context of minimal massive bigravity theory and we find that the cosmology of
this theory suffers from the sudden and big brake singularities, in which we
can see that the parameters of the model approaches to zero the sudden
singularity can be removed.

##  I Introduction

The discovery of current cosmic acceleration via the type Ia supernovas data 1
; 2 , the cosmic microwave background (CMB) radiation 3  and the large scale
structure4 ; 5 , has constructed a renewed interest in theories that modify
standard gravity. Apart from current explanations as the cosmological constant
scenario and the dynamical dark energy option, there is also the belief that
the theory of general relativity (GR) is no longer working on the cosmological
scale and needs a modification. As a result, many modified theories have been
introduced to justify the present expanding acceleration of our universe
without the need of introducing the unknown dark energy. Giving a tiny mass to
graviton is an alternative modification of general relativity which can cause
a small cosmological constant term to cope with the rationalizing this picture
of accelerated expanding universe. About the old history of massive gravity as
a comparably successful modified theory we should refer to the year 1939, when
Fierz and Pauli proposed their linear model 6 . Actually, the Fierz-Pauli
linear action dose not eventuate the linearized Einstein gravity action in the
zero-mass limit besides cannot satisfy the solar system tests resulted of the
vDVZ discontinuity 7 . The familiar Vainshtein mechanism says that nonlinear
interactions can avoid this discontinuity 8  and successively, the nonlinear
terms create a Boulware-Deser ghost 9 ; 10 . In the end, a ghost-free
nonlinear massive gravity theory was prepared successfully by de Rham,
Gabadadze and Tolly (dRGT) 11 . In this manner, by studding 12  we can review
the levels leading to the modern approach.
In the scenario of massive gravity, in addition to the principal metric
$g_{\mu\nu}$ a second tensor field $f_{\mu\nu}$ plays a key role. In this
respect, the dRGT theory was presented as a ghost-free one 13 . It is notable
that to say the cosmology in this model dose not support the flat Friedmann-
Robertson-Walker (FRW) solution 14 , thus it may be ruled out by cosmological
observations. It is noticeable that open cosmological solutions have been
obtained with two classifications of solutions in which the massive terms in
the field equations lead to an effective cosmological constant 15 .
Additionally, about the birth of the massive bigravity theory as a modified
massive gravity theory one can refer to 16 ; 17 , in which the authors convert
the second tensor field to a dynamical one, similarly as principal metric
$g_{\mu\nu}$, although the only difference in this work is that the matter
source just couples to the auxiliary metric. Thus, the resulted ghost-free
massive bigravity theory allows cosmologically viable solutions. As discussed,
two dynamical metrics by which we construct this theory behave in a completely
symmetric way and then breaks the belief in aether-like concept of the
principal metric in massive gravity. These mentioned features of massive
gravity and massive bigravity models encourages people to go more through the
deep study of current cosmological issues in them. In the several literatures,
the cosmological solutions in massive gravity and massive bigravity theories
have considered in Refs. 18 ; 19 ; 20 ; 21 ; 22 ; 23 .
Actually, the discovery of the accelerated expanding universe creates much
enthusiasms to do research about the features of the cosmic fluid which causes
this accelerating expansion named dark energy picturing the overabundance of
new types of singularities which are completely different than the initial big
bang singularity, the so called “finite-time future singularities”.
Empathically, about the initial cosmological singularity in the early universe
we are supposed to refer to the emergent universe scenario in which the
universe emerges from a static state as the Einstein static universe (ESU)
indicating a non-zero spatial volume with a positive curvature by which the
initial big bang singularity will be avoided 24 . Around our discussion, the
ESU in massive gravity model has been studied in 25 , in which we see that all
curvature spaces $\kappa=\pm 1$ and $\kappa=0$ are allowed to describe the
ESU. Also, the initial singularity problem and the ESU in the massive
bigravity models have been considered in 26 , in which the authors have
obtained that the flat FRW universe is not capable to explain the ESU. As a
result in the following study we take $\kappa=-1$ and choose an open universe.

Many works have studied about the theoretical possibility that expanding
universes end up a violent final destiny at a future finite-time without
meeting an expansion maximum or sequent collapse to a “big crunch”
singularity. These new types of singularities are characterized by the
violation of some of the energy conditions. As a result, we see that this
violation causes divergencies in some of the physical quantities involving the
scale factor, energy density and pressure profile. As examples, one may
mention big rip singularity 27 , sudden singularity 28 , big freeze
singularity 29  and other classifications for the possible future
singularities 30 . Additionally, there are some works about cosmological
singularity issues which have been studied beyond the Einstein’s general
relativity 30r ; 300 . Thus, regarding the interesting aspects of the massive
gravity and massive bigravity theories, we have a great motivation to study
the features of the future cosmological singularities in these modified self
accelerating gravity models.

The paper is organized as follows. In section II, we review briefly the
massive gravity theory and its cosmological equations. In section III, we
study the possible finite-time future singularities in the massive universe.
In section IV, we consider the massive bigravity theory with its modified
Friedmann equations and then we address the possible future finite-time
cosmological singularities in massive bigravity theory, in section V. Finally,
we close the paper whit a conclusion in section VI.

##  II Friedman equations in the Massive Gravity Theory

The theory of massive gravity is introduced by the following action as 14

$S=M_{g}^{2}\int d^{4}x\sqrt{-{\rm
det}g}\left(-\frac{1}{2}R+m^{2}\mathcal{L}_{Int}\right)+S_{m},$ (1)

where $M_{g}$, $R$ and $m$ are the Planck mass of the metric $g_{\mu\nu}$,
Ricci scalar and the graviton mass, respectively. We should note that $S_{m}$
is just describing the ordinary matter and $\mathcal{L}_{Int}$, non-derivative
massive graviton terms defined as

$\mathcal{L}_{Int}=\frac{1}{2}\left(S^{2}-S^{A}_{B}S^{B}_{A}\right)+\frac{c_{3}}{3!}\epsilon_{MNPQ}\epsilon^{ABCQ}S^{M}_{A}S^{N}_{B}S^{P}_{C}+\frac{c_{4}}{4!}\epsilon_{MNPQ}\epsilon^{ABCQ}S^{M}_{A}S^{N}_{B}S^{Q}_{D},$
(2)

that we have $S=S^{A}_{A}$ and also $c_{3}$ and $c_{4}$ are two constants
beside that $\epsilon_{MNPQ}$ is the Levi-Civita tensor density. Moreover

$S^{A}_{B}=\delta^{A}_{B}-\gamma^{A}_{B},$ (3)

in which $\gamma^{A}_{B}$ is defined by

$\gamma^{A}_{B}\gamma^{C}_{A}=g^{AC}f_{AB}.$ (4)

Here $f_{AB}$ is a symmetric tensor field. The original dRGT theory respects
the Poincarè symmetry in the field space and thus the reference metric is
Minkowski, i.e. $f_{AB}=\eta_{AB}$$=$ diag$(-1,1,1,1)$.
The spatially homogeneous and isotropic FRW metric can be written as

$ds^{2}=-dt^{2}+a^{2}(t)\left(\frac{dr^{2}}{1-\kappa
r^{2}}+r^{2}d^{2}\Omega\right),$ (5)

here $a(t)$ is the cosmic scale factor, $t$ is the cosmic time and
$\kappa=0,\pm 1$ plays the role of the three dimensional space constant
curvature. According to 31  the first Friedmann equation with dimensionless
constants $c_{3}$ and $c_{4}$ by imposing the Bianchi constraints is as
follows

$3\frac{\dot{a}^{2}+\kappa}{a^{2}}=m^{2}\left(4c_{3}+c_{4}-6+3\mu\frac{3-3c_{3}-c_{4}}{a}+3\mu^{2}\frac{c_{4}+2c_{3}-1}{a^{2}}-\mu^{3}\frac{c_{3}+c_{4}}{a^{3}}\right)+\frac{\rho}{M_{g}^{2}}.$
(6)

Where $\mu$ is an integrating constant. As a result of minimal matter coupling
to gravity, the matter equation of motion reads

$\dot{\rho}+3H\left(\rho+P\right)=0,$ (7)

that $H=\frac{\dot{a}}{a}$. If we assume an equation of state parameter,
namely $\omega$ then we have $P=\omega\rho$. According to 32  the simplest
choice that presents a successful Vainshtein effect in the weak field limit is
reducing the parameter space to the subset $c_{3}=-c_{4}$. Using this
simplification, the first and second Friedmann equations are extracted as
follows

$H^{2}=\frac{\rho}{3M_{g}^{2}}-\frac{\kappa}{a^{2}}+\frac{m^{2}}{3}\left(h_{1}+\frac{h_{2}}{a}+\frac{h_{3}}{a^{2}}\right),$
(8)

$\dot{H}=\frac{-\rho}{2M_{g}^{2}}\left(1+\omega\right)+\frac{\kappa}{a^{2}}-\frac{m^{2}}{6}\left(\frac{h_{2}}{a}+\frac{2h_{3}}{a^{2}}\right).$
(9)

Where $\frac{a}{\mu}\rightarrow a$ and also

$\displaystyle h_{1}=$ $\displaystyle-3c_{4}-6,$ $\displaystyle h_{2}=$
$\displaystyle 3\left(3+2c_{4}\right),$ $\displaystyle h_{3}=$
$\displaystyle-3\left(1+c_{4}\right).$

As we have $\dot{H}=\frac{\ddot{a}}{a}-\frac{\dot{a}^{2}}{a^{2}}$, thus we can
write

$\frac{\ddot{a}}{a}-\frac{\dot{a}^{2}}{a^{2}}+\frac{\rho}{2M_{g}^{2}}\left(1+\omega\right)-\frac{\kappa}{a^{2}}+\frac{m^{2}}{6}\left(\frac{h_{2}}{a}+\frac{2h_{3}}{a^{2}}\right)=0.$
(11)

Rearranging the firs Friedmann equation (8) and also the above equation, we
have

$\frac{\rho}{3M_{g}^{2}}-\frac{\dot{a}^{2}}{a^{2}}-\frac{\kappa}{a^{2}}+\frac{m^{2}}{3}\left(-3c_{4}-6+\frac{1}{a}\left(3\left(3+2c_{4}\right)\right)+\frac{1}{a^{2}}\left(-3\left(1+c_{4}\right)\right)\right)=0,$
(12)

and

$\frac{\rho}{2M_{g}^{2}}\left(1+\omega\right)-\frac{\kappa}{a^{2}}+\frac{\ddot{a}}{a}-\frac{\dot{{}^{2}}}{a^{2}}+\frac{m^{2}}{6}\left(\frac{1}{a}\left(3\left(3+2c_{4}\right)\right)+\frac{1}{a^{2}}\left(-6\left(1+c_{4}\right)\right)\right)=0.$
(13)

One finds $\rho(t)$ and $P(t)$ terms

$\rho(t)=\frac{3M_{g}^{2}}{a^{2}}\left(\dot{a}^{2}+\kappa+m^{2}\left(\left(c_{4}+2\right)a^{2}-\left(3+2c_{4}\right)a+1+c_{4}\right)\right).$
(14)

$P(t)=-\frac{M_{g}^{2}}{a^{2}}\left(2\ddot{a}a+\dot{a}^{2}+\kappa+m^{2}\left(3\left(2+c_{4}\right)a^{2}-2\left(3+2c_{4}\right)a+1+c_{4}\right)\right).$
(15)

In the next section we try to study the possible types of finite-time
singularities may be occurred in the future of an open FRW universe in the
context of massive gravity theory.

##  III Future cosmological singularities in Massive Gravity Theory

As mentioned in the introduction the possible types of finite-time future
singularities in GR have been studied by Barrow 28 ; 34 . We are also
interested in analyzing massive gravity theory in this research area to obtain
all allowed types of finite-time future singularities to imagine different
kinds of late time cosmological evolutions may happen in this model.
To start considering the finite-time future singularities in context of the
massive gravity and massive bigravity, first we give short review of this
issue in the GR case. Thus, by considering the standard FRW universe with the
spatial curvature parameter ($\kappa$) and $8\pi G=c=1$, the Einstein field
equations reduce to

$3H^{2}=\rho-\frac{\kappa}{a^{2}},$ (16)

$\frac{\ddot{a}}{a}=-\left(\frac{\rho+3P}{6}\right),$ (17)

Note that to solve the above equations we need the continuity equation, i.e.
equation (7).

So, to search the finite-time future cosmological singularities in the time
interval $0<t<t_{s}$, we can write the following explicit expression for the
scale factor which it has been studied in 28

$a(t)=1+\alpha t^{q}+\beta\left(t_{s}-t\right)^{n}.$ (18)

Where $\alpha$, $\beta$, $q$ and $n$ are free positive constants which should
be found. Also, the most general form of this solution in the vicinity of the
singularity at $t=t_{s}$ for Friedmann equations has been considered in 34 .
The suggested expression for the scale factor (18) fulfils our purposes about
its behaviors in the cosmic time interval $0<t<t_{s}$ in which we expect the
scale factor starts from zero then grows to a $a_{s}$ value when $t\rightarrow
t_{s}$. By fixing the initial condition $a(0)=0$ at $t\rightarrow 0$, we get
$\beta=-\frac{1}{t^{n}_{s}}$. Thus with the assumption $q$, $n$$>0$ we obtain
the following expression for the scale factor

$a(t)=1+\left(a_{s}-1\right)\left(\frac{t}{t_{s}}\right)^{q}-\left(1-\frac{t}{t_{s}}\right)^{n},$
(19)

where $a_{s}\equiv a\left(t_{s}\right)$. Then we find

$\dot{a}(t)=\frac{\left(-1+a_{s}\right)q\left(\frac{t}{t_{s}}\right)^{q}}{t}+\frac{n\left(1-\frac{t}{t_{s}}\right)^{-1+n}}{t_{s}}.$
(20)

And so

$\ddot{a}(t)=\frac{\left(-1+a_{s}\right)\left(-1+q\right)q\left(\frac{t}{t_{s}}\right)^{q}}{t^{2}}-\frac{\left(-1+n\right)n\left(1-\frac{t}{t_{s}}\right)^{-2+n}}{t^{2}_{s}}.$
(21)

If $t\rightarrow t_{s}$ the above extracted equation reads

$\ddot{a}\rightarrow q\left(q-1\right)\alpha
t^{q-2}-\frac{n\left(n-1\right)}{t_{s}^{2}\left(1-\frac{t}{t_{s}}\right)^{2-n}}\rightarrow-\infty.$
(22)

According to (19), for $t\rightarrow t_{s}$ we should choose the intervals
$1<n<2$ and $0<q<1$ we get $a\rightarrow a_{s}$, also $H_{s}$ and $\rho_{s}$
are finite but $P_{s}\rightarrow\infty$ because there is not any upper bound
on the pressure profile. Moreover, for $2<n<3$ we have a finite $\ddot{a}$ but
$\dddot{a}$ and also $P_{s}$ remains finite while
$\dot{P_{s}}\rightarrow\infty$. In other words, there is an initial
singularity when $t\rightarrow 0$ with $\rho\rightarrow\infty$ and
$P\rightarrow\infty$. From (22) and (17) it can be seen that during the time
interval $0<t<t_{s}$ the strong energy conditions hold ($\rho>0$ and
$\rho+3P>0$). Therefor, this especial solution in an expanding universe
indicates that it is contingent to extend a “sudden” singularity at a finite-
time future (type II), regardless of the sign of the spatial curvature
parameter 34 .

Additionally, the “big rip” singularity (type I) can arise in GR and modified
gravity theories by the effect of matter fields which do not respect the
dominant energy condition such as phantom or ghost field as $\rho+P<0$ 27 .

it is worthwhile to note that, the extreme phantom equation of state
($\rho+P<0$ ) for generating a future singularity in a non-contracting
universe is sufficient but not necessary, according to the first reference in
28 . Therefore, this finite-time future singularity can appear in the
expanding Friedmann universe without taking the condition $\rho+P<0$. Thus,
the solution (18) indicates that for an expanding universe, there is a
possibility to develop a “big-rip” singularity (type I) in a finite-time in
the future even the matter fields satisfy the strong energy conditions.
Other types of cosmological singularities (type III and IV) are different from
the sudden finite-time future singularity in the sense that $\rho$ diverges.
Usually, in these singularities the pressure profile and density have a non-
realistic equation of state, e.g. $P=-\rho-f(\rho)$ which they are studied in
35 .

In Table $1$ we summarize the types of possible future singularities in an
expanding FRW universe considering the behavior of $a(t)$, $\rho(t)$ and
$P(t)$ in a finite cosmological time 35 .

Table 1:  The possible finite-time future singularities for $t\rightarrow
t_{s}$
Singularity Type $a(t)$ $\rho(t)~{}~{}~{}$ $P(t)~{}$ I $a\rightarrow\infty$
$\rho\rightarrow\infty$ $|P|\rightarrow\infty$ II $a\rightarrow a_{s}$
$\rho\rightarrow\rho_{s}$ $|P|\rightarrow\infty$ III $a\rightarrow a_{s}$
$\rho\rightarrow\infty$ $|P|\rightarrow\infty$ IV $a\rightarrow a_{s}$
$\rho\rightarrow 0$ $|P|\rightarrow 0$

According to table $1$, the type I singularity is known as “big rip or cosmic
doomsday” explained in 27  beside the type II, III and IV singularities are
known as “sudden” 28 , “big freeze” 37  and “big brake or big separation” 38
singularities, respectively. It is noticeable that in case IV singularity,
however $a(t)$, $\rho(t)$ and $P(t)$ remain finite, the higher derivatives of
$H(t)$ diverge.

So, to continue we try to classify the possible finite-time future
singularities in the context of massive gravity theory with respect to $q$ and
$n$ values as in the following classifications. It is worth noting, the
spacial curvature $\kappa$ in this theory is free to be choused $\pm 1$ or $0$
as studied in 25  but according to those explained in the introduction we just
work with $\kappa=-1$.

Class I: $\textbf{t}\rightarrow\textbf{t}_{s}$ with
$\textbf{n}\in\left(\textbf{0}\textbf{,}\textbf{1}\right)$ and
$\textbf{q}\in\left(\textbf{0}\textbf{,}\textbf{1}\right]$

In this class we obtain

$a(t_{s})\rightarrow
a_{s},~{}~{}~{}\dot{a}(t_{s})\rightarrow+\infty,~{}~{}~{}H(t_{s})\rightarrow+\infty,~{}~{}~{}\ddot{a}(t_{s})\rightarrow+\infty,$
(23)

$\rho(t_{s})\rightarrow\rho_{s}>0,~{}~{}~{}|P(t_{s})|\rightarrow\infty.$ (24)

According to the table $1$, this is “sudden singularity”. To more consider, we
have plotted the evolution of the quantities $a(t)$, $\rho(t)$ and $P(t)$ in
Fig. $1$, for two typical values of $q$. It can be seen that increasing $q$
value results in that $\rho(t)$ and $P(t)$ increase more rapidly for $t$
$\rightarrow$ $t_{s}$.

Figure 1: The scale factor $a(t)$, density $\rho(t)$ and pressure $P(t)$
evolution with two typical values $q=0.2$ and $q=0.8$ in the upper and lower
plots, respectively. We have fixed $a_{s}=2000$, $c_{4}=0.3$, $t_{s}=1000$,
$M_{g}^{2}=1$ and $\kappa=-1$.

Class II: $\textbf{t}\rightarrow\textbf{t}_{s}$ with $\textbf{n}=1$ and
$\textbf{q}\in\left(\textbf{0}\textbf{,}\textbf{1}\right]$

In this class one finds

$a(t_{s})\rightarrow
a_{s},~{}~{}~{}\dot{a}(t_{s})\rightarrow\dot{a}_{S}>0,~{}~{}~{}H(t_{s})\rightarrow
H_{s}>0,~{}~{}~{}\ddot{a}(t_{s})\rightarrow\ddot{a}_{s}\leq 0,$ (25)

$\rho(t_{s})\rightarrow\rho_{s}>0,~{}~{}~{}|P(t_{s})|\rightarrow P_{s}.$ (26)

Clearly, there is no finite-time future singularity because all $a_{s}$,
$\dot{a}_{s}$, $H_{s}$, $\ddot{a}_{s}$, $\rho_{s}$ and $P_{s}$ are finite. It
is noticeable that for $q=1$ we obtain the zero acceleration $\ddot{a}_{s}=0$.
Similar to class I we have plotted $a(t)$, $\rho(t)$ and $P(t)$ evolutions in
Fig. $2$, for two typical values of $q$. We can see that by increasing the
value of $q$ the slopes of the $\rho(t)$ and $P(t)$ for $t\rightarrow t_{s}$
is increasing.

Figure 2: The scale factor $a(t)$, density $\rho(t)$ and pressure $P(t)$
evolution by choosing $n=1$ with two typical values $q=0.2$ and $q=0.8$ in the
upper and lower plots, respectively. We take $a_{s}=2000$, $c_{4}=0.3$,
$t_{s}=1000$, $M_{g}^{2}=1$ and $\kappa=-1$.

Class III: $\textbf{t}\rightarrow\textbf{t}_{s}$ with
$\textbf{n}\in\left(\textbf{1},\textbf{2}\right)$ and
$\textbf{q}\in\left(\textbf{0}\textbf{,}\textbf{1}\right]$

In this class one finds

$a(t_{s})\rightarrow
a_{s},~{}~{}~{}\dot{a}(t_{s})\rightarrow\dot{a}_{S}>0,~{}~{}~{}H(t_{s})\rightarrow
H_{s}>0,~{}~{}~{}\ddot{a}(t_{s})\rightarrow-\infty,$ (27)

$\rho(t_{s})\rightarrow\rho_{s}>0,~{}~{}~{}|P(t_{s})|\rightarrow P_{s}.$ (28)

Apparently, there is no finite-time future singularity because $a_{s}$,
$\rho_{s}$ and $P_{s}$ are finite. As a result, there is no finite-time future
singularity in this class. Similar to previous cases, $a(t)$, $\rho(t)$ and
$P(t)$ evolutions are plotted in Fig. $3$. By increasing $q$ the slopes of
$\rho(t)$ and $P(t)$ is increasing for $t\rightarrow t_{s}$.

Figure 3: The scale factor $a(t)$, density $\rho(t)$ and pressure $P(t)$
evolution with two typical values $q=0.2$ and $q=0.8$ in the upper and lower
plots, respectively. We have fixed $a_{s}=2000$, $c_{4}=0.3$, $t_{s}=1000$,
$M_{g}^{2}=1$ and $\kappa=-1$.

Class IV: $\textbf{t}\rightarrow\textbf{t}_{s}$ with
$\textbf{n}\in\left[\textbf{2},\infty\right)$ and
$\textbf{q}\in\left(\textbf{0}\textbf{,}\textbf{1}\right]$

In this class we face with

$a(t_{s})\rightarrow
a_{s},~{}~{}~{}\dot{a}(t_{s})\rightarrow\dot{a}_{S}>0,~{}~{}~{}H(t_{s})\rightarrow
H_{s}>0,~{}~{}~{}\ddot{a}(t_{s})\rightarrow\ddot{a}_{s}\leq 0,$ (29)

$\rho(t_{s})\rightarrow\rho_{s}>0,~{}~{}~{}|P(t_{s})|\rightarrow P_{s}.$ (30)

Obviously, in this class there is not any singularity and all quantities are
finite. The behavior of the scale factor $a(t)$, density $\rho(t)$ and
pressure $P(t)$ for two typical values of $q$ are plotted in Fig. $4$.

Figure 4: The scale factor $a(t)$, density $\rho(t)$ and pressure $P(t)$
evolution with two typical values $q=0.2$ and $q=0.8$ in the upper and lower
plots, respectively. We have fixed $a_{s}=2000$, $c_{4}=0.3$, $t_{s}=1000$,
$M_{g}^{2}=1$ and $\kappa=-1$.

In the next section, we introduce the massive bigravity theory by extracting
the field equations. Moreover we find the field equations of the most
simplified version of this theory called minimal massive bigravity model to
consider the finite-time future singularities in the cosmological backgrounds.
We show that this simplified massive bigravity model suffers from “sudden” and
“big brake” singularities.

##  IV Friedmann equations in the Massive Bigravity theory

In this section we review the ghost free massive bigravity model and its
equations of motion in the background of FRW universe. The most general action
for two interacting metrics $g_{\mu\nu}$ and $f_{\mu\nu}$ with a non
derivative potential determined by the requirement of the absence of the
Boulware-Deser ghost considering the coupling to matter via the matter
lagrangian $\mathcal{L}_{m}$ 16 , is

$\displaystyle S=-\frac{M_{g}^{2}}{2}\int
dx^{4}\sqrt{-detg}R(g)-\frac{M_{f}^{2}}{2}\int
dx^{4}\sqrt{-detf}R(f)+m^{2}M_{g}^{2}\int
dx^{4}\sqrt{-detg}\sum^{4}_{n=0}\beta_{n}e_{n}\left(g^{-1}f\right)+\int
dx^{4}\sqrt{-detg}\mathcal{L}_{m},$ (31)

that $\beta_{n}$ are free parameters and $e_{n}(X)$’s are elementary symmetric
polynomials of the eigenvalues of $X$:

$\displaystyle e_{0}(X)=$ $\displaystyle
1,~{}~{}e_{1}(X)=[X],~{}~{}e_{2}(X)=\frac{1}{2}\left([X]^{2}-[X^{2}]\right),$
$\displaystyle e_{3}(X)=$
$\displaystyle\frac{1}{6}\left([X]^{3}-3[X][X^{2}]+2[X^{3}]\right),$
$\displaystyle e_{4}(X)=$
$\displaystyle\frac{1}{24}\left([X]^{4}-6[X]^{2}[X^{2}]+3[X^{2}]^{2}+8[X][X^{3}]-6[X^{4}]\right),$
$\displaystyle e_{i}(X)=$ $\displaystyle 0~{}~{}{\rm for}~{}~{}i>4.$ (32)

The square root matrix $\sqrt{g^{-1}f}$ is defined by
$\left(\sqrt{g^{-1}f}\right)^{\mu}~{}_{\rho}\left(\sqrt{g^{-1}f}\right)^{\rho}~{}_{\nu}=g^{\mu\rho}f_{\rho\nu}=X^{\mu}~{}_{\nu}$.

We consider the FRW metric with spatial curvature $\kappa$, for both line
elements as follows

$ds_{g}^{2}=-dt^{2}+a(t)^{2}\left(\frac{dr^{2}}{1-\kappa
r^{2}}+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\varphi^{2}\right),$ (33)
$ds_{f}^{2}=-X(t)^{2}dt^{2}+Y(t)^{2}\left(\frac{dr^{2}}{1-\kappa
r^{2}}+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\varphi^{2}\right).$ (34)

Where $a(t)$ is the cosmic scale factor related to $g_{\mu\nu}$ and $Y(t)$ is
the one related to $f_{\mu\nu}$. Obviously, $X(t)$ (the lapse function of
$f_{\mu\nu}$ metric) is a function of time and note that we are not allowed to
choose $X=1$ nor $X=Y$. According to 19  in which, the Bianchi constraint
(implied by general covariance) has been extracted in details, turns out

$X=\frac{\dot{Y}}{\dot{a}}=\frac{dY}{da},$ (35)

which leads to leaving just two free functions to work with. Using this
relation beside the ansatz (33) and (34) give us the modified Friedmann
equations of the $g_{\mu\nu}$ metric

$\rho=M_{g}^{2}\left(3\left(\frac{\dot{a}}{a}\right)^{2}+\frac{3\kappa}{a^{2}}-m^{2}\left(\beta_{0}+3\beta_{1}\frac{Y}{a}+3\beta_{2}\frac{Y^{2}}{a^{2}}+\beta_{3}\frac{Y^{3}}{a^{3}}\right)\right),$
(36)
$P=M_{g}^{2}\left(-2\frac{\ddot{a}}{a}-\frac{\dot{a}^{2}}{a^{2}}-\frac{\kappa}{a^{2}}+m^{2}\left(\beta_{0}+2\beta_{1}\left(\frac{Y}{a}+X\right)+\beta_{2}\left(\frac{Y^{2}}{a^{2}}+\frac{2YX}{a}\right)+\beta_{3}\frac{Y^{2}X}{a^{2}}\right)\right).$
(37)

It should be noted that here we have the source structure $T^{0}_{0}=-\rho$,
$T^{1}_{1}=P$ because we have assumed a perfect fluid source
$T_{\mu\nu}=\left(\rho+P\right)u_{\mu}u_{\nu}+Pg_{\mu\nu}$ with
$T^{1}_{1}=T^{2}_{2}=T^{3}_{3}$ (consistent with the symmetries of the space
time). As $m^{2}\rightarrow 0$, these $g_{\mu\nu}$ Friedmann equations reduce
to the ordinary Friedmann equations of cosmology but the cosmological
solutions will not always be well defined in this limit. Moreover, we can find
the modified Friedmann equations of the metric $f_{\mu\nu}$ as follows

$-3\left(\frac{\dot{a}}{Y}\right)^{2}-\frac{3\kappa}{Y^{2}}+\frac{m^{2}}{M_{*}^{2}}\left(\beta_{4}+3\beta_{3}\frac{a}{Y}+3\beta_{2}\frac{a^{2}}{Y^{2}}+\beta_{1}\frac{a^{3}}{Y^{3}}\right)=0,$
(38)

$-2\frac{\ddot{a}}{YX}-\left(\frac{\dot{a}}{Y}\right)^{2}-\frac{\kappa}{Y^{2}}+\frac{m^{2}}{M_{*}^{2}}\left(\beta_{4}+\beta_{3}\left(2\frac{a}{Y}+\frac{1}{X}\right)+\beta_{2}\left(\frac{a^{2}}{Y^{2}}+\frac{2a}{YX}\right)+\beta_{1}\frac{a^{2}}{Y^{2}X}\right)=0.$
(39)

In which we have the dimensionless ratio of Planck masses as

$M_{*}^{2}\equiv\frac{M_{f}^{2}}{M_{g}^{2}}.$ (40)

Meanwhile, for simplicity we assume

$\gamma\equiv\frac{Y}{a}.$ (41)

Therefore, equations (36) and (38) can be read as

$\frac{\beta_{3}}{3}\gamma^{3}+\beta_{2}\gamma^{2}+\beta_{1}\gamma+\frac{\beta_{0}}{3}+\frac{\rho}{3m^{2}M_{g}^{2}}-\frac{H^{2}}{m^{2}}-\frac{\kappa}{m^{2}a^{2}}=0,$
(42)

and

$\frac{\beta_{4}}{3M_{*}^{2}}\gamma^{2}+\frac{\beta_{3}}{M_{*}^{2}}\gamma+\frac{\beta_{1}}{3M_{*}^{2}}\frac{1}{\gamma}+\frac{\beta_{2}}{M_{*}^{2}}-\frac{H^{2}}{m^{2}}-\frac{\kappa}{m^{2}a^{2}}=0.$
(43)

Obviously, by using two above equations then eliminating $H^{2}$ we find the
following quadratic equation for $\gamma$

$\frac{\beta_{3}}{3}\gamma^{4}+\left(\beta_{2}-\frac{\beta_{4}}{3M_{*}^{2}}\right)\gamma^{3}+\left(\beta_{1}-\frac{\beta_{3}}{M_{*}^{2}}\right)\gamma^{2}+\left(\frac{\rho}{3m^{2}M_{g}^{2}}+\frac{\beta_{0}}{3}-\frac{\beta_{2}}{M_{*}^{2}}\right)\gamma-\frac{\beta_{1}}{3M_{*}^{2}}=0.$
(44)

###  IV.1 The minimal model with $\beta_{1}=\beta_{3}=0$

Although, setting only $\beta_{3}=0$, in the action (31) eliminates the
highest order interaction term in $\sqrt{g^{-1}f}$, because of the invariance
of this action under the following exchange

$g\leftrightarrow
f,~{}~{}~{}~{}~{}~{}~{}\beta_{n}\leftrightarrow\beta_{4-n},~{}~{}~{}~{}~{}~{}~{}~{}~{}M_{g}\leftrightarrow
M_{f},~{}~{}~{}~{}~{}~{}~{}~{}~{}m^{2}\leftrightarrow\frac{m^{2}M_{g}^{2}}{M_{f}^{2}},$
(45)

we still have a cubic order interaction term in $\sqrt{f^{-1}g}$ which we
eliminate it by choosing $\beta_{1}=0$ to obtain the minimal massive bigravity
model which is the simplest in the class. Actually, in this class we have
$\beta_{2}<\frac{\beta_{4}}{3M_{*}^{2}}$. As a result, for these values
equation (44) gives

$\gamma^{2}=\frac{\frac{\rho}{m^{2}M_{g}^{2}}+\beta_{0}-3\beta_{2}M_{*}^{-2}}{\beta_{4}M_{*}^{-2}-3\beta_{2}}.$
(46)

Referring to equations (38) and (39) with $\beta_{1}=\beta_{3}=0$, we can find
the solution

$\frac{\dot{Y}}{\dot{a}}=X(t)=\frac{\frac{-2a\ddot{a}}{\gamma}+\frac{2m^{2}a^{2}\beta_{2}}{\gamma
M_{*}^{2}}}{\frac{\dot{a}^{2}}{\gamma^{2}}+\frac{\kappa}{\gamma^{2}}-\frac{m^{2}a^{2}\beta_{4}}{M_{*}^{2}}-\frac{m^{2}a^{2}\beta_{2}}{M_{*}^{2}\gamma^{2}}}.$
(47)

This relation helps us to simplify the pressure equation (37) to become able
to plot the functions of $P(t)$ and $\rho(t)$ to analyze their behaviors.

##  V Future cosmological singularities in Massive Bigravity Theory

As section III, we are going to follow that approach again to study the
possible finite-time future singularities in minimal massive bigravity
universe. Considering (19), (20) and (21) we calculate some classifications in
the following.

Class I: $\textbf{t}\rightarrow\textbf{t}_{s}$ with
$\textbf{n}\in\left(\textbf{0}\textbf{,}\textbf{1}\right)$ and
$\textbf{q}\in\left(\textbf{0}\textbf{,}\textbf{1}\right]$

In this class we can write

$a(t_{s})\rightarrow
a_{s},~{}~{}~{}\dot{a}(t_{s})\rightarrow\dot{a}_{s}>0,~{}~{}~{}H(t_{s})\rightarrow
H_{s}>0,~{}~{}~{}\ddot{a}(t_{s})\rightarrow\ddot{a}_{s}\geq 0,$ (48)

$\rho(t_{s})\rightarrow 0,~{}~{}~{}|P(t_{s})|\rightarrow 0.$ (49)

This result represents the case IV singularity as it mentioned in Table $1$,
namely “big brake or big separation”. We have plotted the evolution of the
quantities $a(t)$, $\rho(t)$ and $P(t)$ in Fig. $5$, for two typical values of
$q$. We can see that increasing $q$ value results in that $\rho(t)$ and $P(t)$
increase more rapidly for $t$ $\rightarrow$ $t_{s}$.

Figure 5: The scale factor $a(t)$, density $\rho(t)$ and pressure $P(t)$
evolution with two typical values $q=0.2$ and $q=0.8$ in the upper and lower
plots, respectively. We have fixed $a_{s}=2000$, $M_{*}=M_{g}=m^{2}=1$,
$t_{s}=1000$, $\beta_{0}=9.9$, $\beta_{1}=0$, $\beta_{2}=0.3$, $\beta_{3}=0$,
$\beta_{4}=4$ and $\kappa=-1$.

Class II: $\textbf{t}\rightarrow\textbf{t}_{s}$ with $n=1$ and
$\textbf{q}\in\left(\textbf{0}\textbf{,}\textbf{1}\right]$

In this class we have

$a(t_{s})\rightarrow
a_{s},~{}~{}~{}\dot{a}(t_{s})\rightarrow\dot{a}_{s}>0,~{}~{}~{}H(t_{s})\rightarrow
H_{s}>0,~{}~{}~{}\ddot{a}(t_{s})\rightarrow\ddot{a}_{s}\geq 0,$ (50)

$\rho(t_{s})\rightarrow\rho_{s}<0,~{}~{}~{}|P(t_{s})|\rightarrow P_{s}.$ (51)

This result represents no future singularity. We have plotted the evolution of
the quantities $a(t)$, $\rho(t)$ and $P(t)$ in Fig. $6$, for two typical
values of $q$. We show that increasing $q$ value results in that $\rho(t)$ and
$P(t)$ increase more rapidly for $t$ $\rightarrow$ $t_{s}$.

Figure 6: The scale factor $a(t)$, density $\rho(t)$ and pressure $P(t)$
evolution with two typical values $q=0.2$ and $q=0.8$ in the upper and lower
plots, respectively. We have fixed $a_{s}=2000$, $M_{*}=M_{g}=m^{2}=1$,
$t_{s}=1000$, $\beta_{0}=9.9$, $\beta_{1}=0$, $\beta_{2}=0.3$, $\beta_{3}=0$,
$\beta_{4}=4$ and $\kappa=-1$.

Class III: $\textbf{t}\rightarrow\textbf{t}_{s}$ with
$\textbf{n}\in\left(\textbf{1}\textbf{,}\textbf{2}\right)$ and
$\textbf{q}\in\left(\textbf{0}\textbf{,}\textbf{1}\right]$

This class is identified by

$a(t_{s})\rightarrow
a_{s},~{}~{}~{}\dot{a}(t_{s})\rightarrow\dot{a}_{s}>0,~{}~{}~{}H(t_{s})\rightarrow
H_{s}>0,~{}~{}~{}\ddot{a}(t_{s})\rightarrow\ddot{a}_{s}\geq 0,$ (52)

$\rho(t_{s})\rightarrow\rho_{s}<0,~{}~{}~{}|P(t_{s})|\rightarrow\infty.$ (53)

This result represents “sudden” singularity. The evolution of the quantities
$a(t)$, $\rho(t)$ and $P(t)$ have been plotted in Fig. $7$, for two typical
values of $q$. The figure implies that increasing $q$ value results in that
$\rho$ and $P$ increase more rapidly for $t$ $\rightarrow$ $t_{s}$. In this
singularity we can find the root of the pressure singularity in the relation
(37) in which we can see that the dynamical parameter $\dot{Y}$ is responsible
for the infinity pressure at a finite time. If $\dot{Y}$ goes to infinity then
according to (47) we have $X(t)\rightarrow\infty$. In order to avoid this
singularity we just can again refer to (37) and suppose that
$\beta_{{}_{2}}\gamma$ takes very small value so

$\beta_{{}_{2}}\gamma\frac{\dot{Y}}{\dot{a}}\rightarrow~{}~{}\emph{\rm
finite}.$ (54)

Note that in the limit $\beta_{{}_{2}}\gamma\rightarrow 0$ the cosmological
equations in the minimal massive bigravity approach to the Friedmann
equations, but it seems that in this limit the asymptotic behavior of the
minimal massive bigravity is dissimilar from GR at least in the sudden
cosmological singularity case.

Figure 7: The scale factor $a(t)$, density $\rho(t)$ and pressure $P(t)$
evolution with two typical values $q=0.2$ and $q=0.8$ in the upper and lower
plots, respectively. We have fixed $a_{s}=2000$, $M_{*}=M_{g}=m^{2}=1$,
$t_{s}=1000$, $\beta_{0}=9.9$, $\beta_{1}=0$, $\beta_{2}=0.3$, $\beta_{3}=0$,
$\beta_{4}=4$ and $\kappa=-1$.

Class IV: $\textbf{t}\rightarrow\textbf{t}_{s}$ with
$\textbf{n}\in\left[\textbf{2}\textbf{,}\infty\right)$ and
$\textbf{q}\in\left(\textbf{0}\textbf{,}\textbf{1}\right]$

In this class we face with

$a(t_{s})\rightarrow
a_{s},~{}~{}~{}\dot{a}(t_{s})\rightarrow\dot{a}_{s}>0,~{}~{}~{}H(t_{s})\rightarrow
H_{s}>0,~{}~{}~{}\ddot{a}(t_{s})\rightarrow\ddot{a}_{s}\geq 0,$ (55)

$\rho(t_{s})\rightarrow\rho_{s}<0,~{}~{}~{}|P(t_{s})|\rightarrow P_{s}.$ (56)

This result represents no future singularity. The evolution of the quantities
$a(t)$, $\rho(t)$ and $P(t)$ have been plotted in Fig. $8$, for two typical
values of $q=0.2$ and $q=0.8$. The figure implies that increasing $q$ value
results in that $\rho(t)$ and $P(t)$ increase more rapidly for $t$
$\rightarrow$ $t_{s}$.

Figure 8: The scale factor $a(t)$, density $\rho(t)$ and pressure $P(t)$
evolution with two typical values $q=0.2$ and $q=0.8$ in the upper and lower
plots, respectively. We have set $a_{s}=2000$, $M_{*}=M_{g}=m^{2}=1$,
$t_{s}=1000$, $\beta_{0}=9.9$, $\beta_{1}=0$, $\beta_{2}=0.3$, $\beta_{3}=0$,
$\beta_{4}=4$ and $\kappa=-1$.

##  VI Conclusions

In this paper, we have considered the future cosmological singularities in the
ghost-free massive gravity and massive bigravity models. By regarding the
types of future singularities, we have studied the possible classes of finite-
time future cosmological singularities for sudden, big rip, big freeze and big
brake singularities. In an expanding open Fridmann universe (which was
mentioned in the introduction that people could have found open cosmological
solutions in dRGT model), we have obtained the “sudden singularity” in which
the physical quantities $\dot{a}(t)$, $H(t)$, $\ddot{a}(t)$ and $P(t)$ all go
to infinity and cause such a singularity in a finite-time future accompanying
a finite fluid positive density and we can not eliminate these extracted
singularities by no means of our theory parameters such as graviton mass
because there is no coupling between for instance $m$ and dynamical quantities
like $H(t)$ and $\ddot{a}(t)$, nonetheless the singularities such as sudden
finite-time future singularity that pressure approaches to infinity can
basically be avoidable by taking a realistic equation of state for matter
fields $P(\rho)$ by which the pressure profile is upper bounded by some well-
defined density function. As a result, we can eliminate the pressure
singularity at a finite-time in the future whereas the strong energy condition
holds, see the first reference in 28 .
Also, to complete our study in the finite-time future cosmological
singularities we have focused on the massive bigravity theory and reviewed the
modified Friedmann equations. In particular, for more simplicity we have
considered the minimal massive bigravity model with $\beta_{1}=\beta_{3}=0$.
We have found that there exists the possibility of occurring the sudden and
big brake singularities. The sudden singularity can be eliminated in the
minimal massive bigravity theory by taking $\beta_{{}_{2}}\gamma
X(t)\rightarrow~{}\emph{\rm finite}$, nevertheless the theory approaches to GR
regime.

In order to compare the finite-time future singularities in massive gravity
and massive bigravity theories with cosmological singularities in GR, we can
say that the results of our studding are similar to GR in which the sudden
singularity arises and in fact it was predictable because the modification
terms in the massive theories are just able to moderate the singularity not
able to eliminate it and thus the sudden singularity in massive gravity is
common whit GR. The distinction is just about the big brake singularity and
also the behavior of the expanding scale factor in the comparison of massive
gravity and massive bigravity theories and finally in the massive bigravity
theories sudden singularity can be avoided contrary to massive gravity and GR.

## References

  * (1) S. Perlmutter et al., Astrophys. J. 517 (1999) 565.
  * (2) A. G. Riess et al., Astron. J. 116 (1998) 1009.
  * (3) D. N. Spergel et al., Astrophys. J. Suppl. Ser. 148 (2003) 125.
  * (4) D. J. Einstein et al., Astron. J. 633 (2005) 560.
  * (5) M. Tegmark et al., Phys. Rev. D 69 (2004) 103501.
  * (6) M. Fierz and W. Pauli, Proa. R. Soc. A 173 (1939) 211 .
  * (7) H. Van Dam and M. J. G. Veltman, Nucl. Phys. B 22 (1970) 397;
V. I. Zakharov, JETP Lett. 12 (1970) 312.

  * (8) A. I. Vainshtein, Phys. Lett. B 39. 393 (1972).
  * (9) D. G. Boulware and S. Deser, Phys. Rev. D 6 (1972) 3368.
  * (10) N. Arkani-Hamed, H. Goergi and M. D Schwartz, Ann. Phys. 305(2003) 96;
P. Creminelli, A. Nicolis, M. Papucci, E. Trincherini and J. W. Rombouts,
Phys. Rev. D 72 (2005) 124007.

  * (11) C. de Rham and G. Gabadadze, Phys. Rev. D 82 (2010) 044026;
C. de Rham, G. Gabadadze and A. J. Tolley, Phys. Rev. Lett. 106 (2011) 231101.

  * (12) K. Hinterbichler, Rev. Mod. Phys. 84 (2012) 671.
  * (13) S. F. Hassan and R. A. Rosen. Phys. Rev. Lett. 108 (2012) 041101.
  * (14) G. D, Amico, C. de Rham, S. Dubovsky, G. Gabadadze, D. Pirtskhalava and A. J. Tolley, Phys. Rev. D 84 (2011) 124046.
  * (15) A. E. Gumrukcuoglu, C. Linand and S. Mukohyama, JCAP 11 (2011) 030 .
  * (16) S. F. Hassan and R. A. Rosen. JHEP, 02 (2012) 126 .
  * (17) S. F. Hassan, R. A. Rosen and A. Schmidt-May, JHEP 02 (2012) 026 .
  * (18) D. Comeli, M. Crisostomi, F. Nesti and L. Pilo, JHEP 03 (2012) 067.
  * (19) M. Von Strauss, A. Schmidt-May, J. Enander, E. Mortsell and S. F. Hassan, JCAP 03 (2012) 042.
  * (20) M. Volkov, JHEP, 01 (2012) 035.
  * (21) F. Konnig, and L. Amendola, Phys. Rev. D 90 (2014) 044030.
  * (22) Y. Akrami, T. S. Koivisto and M. Sandstad, JHEP, 03 (2013) 099.
  * (23) K. Koyoma, G. Niz and G. Tasinato, Phys. Rev. Lett. 107 (2011) 131101.
  * (24) G. F. R. Ellis and R. Maartens, Class. Quant. Grav. 21 (2004) 223;
G. F. R. Ellis, J. Murugan and C. G. Tsagas, Class. Quant. Grav. 21 (2004)
233.

  * (25) L. Parisi, N. Radicella and G. Vilasi, Phys. Rev, D 86 (2012) 024035.
  * (26) M. Mousavi and F. Darabi, Nucl. Phys. B 919 (2017) 523.
  * (27) R. R. Caldwell, Phys. Lett. B 545 (2002) 23;
R. R. Caldwell, M. Kamionkowski and N. N. Weinberg, Phys. Rev. Lett. 91 (2003)
071301;
I. Albarran, M. Bouhmadi-López, C. Kiefer, J. Marto and P. V. Moniz, Phys.
Rev. D 94 (2016) 063536;
E. Mörtsell, JCAP 02 (2017) 051;
M. Högås and E. Mörtsell, JCAP 05(2021) 001.

  * (28) J. D. Barrow, Class. Quant. Grav. 21 (2004) L79 ;
J. D. Barrow and C. G. Tsagas, Class. Quant. Grav. 22 (2005) 1563 ;
K. Lake, Class. Quant. Grav. 21 (2004) L129.

  * (29) M. Bouhmadi-Lopez, P. F. Gonzalez-Diaz and P. Martin-Moruno, Phys. Lett. B 659 (2008) 1.
  * (30) Y. Heydarzade, Eur. Phys. J. C 79 (2019) 923 .
  * (31) J. L. Rosa, F. S. N. Lobo and D. Rubiera-Garcia, JCAP 07 (2021) 009.
  * (32) M. P. Dabrowski and T. Denkiewicz, AIP Conference Proceedings 1241 (2010) 561 .
  * (33) A. H. Chamseddine and M. S. Volkov, Phys. Lett. B 704 (2011) 652.
  * (34) K. Koyoma, G. Niz and G. Tasinato, Phys. Rev. D 84 (2011) 064033.
  * (35) J. D. Barrow, Class. Quant. Grav, 21 (2004) 5619.
  * (36) S. Nojiri, S.D. Odinstov and S. Tsujikawa, Phys. Rev. D 71 (2005) 063004 .
  * (37) S. Nojiri and S. D. Odintsov, Phys. Rev. D 70 (2004) 103522.
  * (38) V. Gorini, A. Yu. Kamenshchik, U. Moschella and V. Pasquier, Phys. Rev. D 69 (2004) 123512.