# Improving Neural Cross-Lingual Abstractive Summarization via Employing
Optimal Transport Distance for Knowledge Distillation

Thong Nguyen1, Luu Anh Tuan2 Corresponding author

###### Abstract

Current state-of-the-art cross-lingual summarization models employ multi-task
learning paradigm, which works on a shared vocabulary module and relies on the
self-attention mechanism to attend among tokens in two languages. However,
correlation learned by self-attention is often loose and implicit, inefficient
in capturing crucial cross-lingual representations between languages. The
matter worsens when performing on languages with separate morphological or
structural features, making the cross-lingual alignment more challenging,
resulting in the performance drop. To overcome this problem, we propose a
novel Knowledge-Distillation-based framework for Cross-Lingual Summarization,
seeking to explicitly construct cross-lingual correlation by distilling the
knowledge of the monolingual summarization teacher into the cross-lingual
summarization student. Since the representations of the teacher and the
student lie on two different vector spaces, we further propose a Knowledge
Distillation loss using Sinkhorn Divergence, an Optimal-Transport distance, to
estimate the discrepancy between those teacher and student representations.
Due to the intuitively geometric nature of Sinkhorn Divergence, the student
model can productively learn to align its produced cross-lingual hidden states
with monolingual hidden states, hence leading to a strong correlation between
distant languages. Experiments on cross-lingual summarization datasets in
pairs of distant languages demonstrate that our method outperforms state-of-
the-art models under both high and low-resourced settings.

## Introduction

Cross-Lingual Summarization (CLS) is the task of condensing a document of one
language into its shorter form in the target language. Most of contemporary
works can be classified into two categories, i.e. low-resourced and high-
resourced CLS approaches. In high-resourced scenarios, models are provided
with an enormous number of document /summary pairs on which they can be
trained (Zhu et al. 2019; Cao, Liu, and Wan 2020; Zhu et al. 2020). On the
other hand, in low-resourced settings, those document/summary pairs are
scarce, which restrains the amount of information that a model can learn.
While high-resourced settings are preferred, in reality it is difficult to
attain a sufficient amount of data, especially for less prevalent languages.

Figure 1: Examples of Chinese-English cross-lingual summarization. Here we
present the English output generated by the NCLS model of (Zhu et al. 2019)

Most previous works resolving the issue of little training data concentrate on
multi-task learning framework by utilizing the relationship of Cross-Lingual
Summarization (CLS) with Monolingual Summarization (MLS) or Neural Machine
Translation (NMT). Their approach can be further divided into two groups. The
first group equips their module with two independent decoders, one of them
targets the auxiliary task (MLS or NMT). Nevertheless, since two decoders do
not share their parameters, this approach undermines the model’s ability to
align between two tasks (Bai, Gao, and Huang 2021), making the ancillary and
the main task less relied upon each other. Hence, the trained model might
produce output that does not match up the topic, or miss important spans of
text. In Figure 1, we list two samples of documents with their gold and
generated summaries by the model NCLS of (Zhu et al. 2019). As it can be seen,
their cross-lingual outputs do not include key spans from the summary, e.g., “
_underground_ ”, “ _parking lot_ ” in sample 1, and “ _consumers_ ” in sample
2. In both samples, the content of the cross-lingual summary diverges
significantly from the one of the monolingual gold summary.

The second group decides to employ a single decoder dealing with both CLS and
MLS tasks. To this end, the method concatenates the monolingual to cross-
lingual summary and designate the model to sequentially generate the
monolingual summary, and then the cross-lingual one. Unfortunately,
notwithstanding lessening the computational overhead during training by using
solely one decoder, this method is not efficacious in capturing the connection
between two languages in the output, consequently producing representations
that do not take into account language relationships (Luo et al. 2021). In
that case, the correlation of cross-lingual representations will be
tremendously impacted by the structural and morphological similarity of those
languages (Bjerva et al. 2019). As a result, in case of summarizing the
document from one language to another that possesses distinct morphology and
structure properties, such as from Chinese to English, the decoder might be
prone to underperformance, due to the dearth of language correlation between
two sets of hidden representations in the bilingual vector space (Luo et al.
2021).

To solve the aforementioned problem, we propose a novel Knowledge-Distillation
framework for Cross-Lingual Summarization task. Particularly, our framework
consists of a teacher model targetting Monolingual Summarization, and a
student for Cross-Lingual Summarization. We initiate our procedure by
finetuning the teacher model on monolingual document/summary pairs.
Subsequently, we continue to distill summarization knowledge of the trained
teacher into the student model. Because the hidden vectors of the teacher and
student lie upon two disparate monolingual and cross-lingual spaces,
respectively, we propose a Sinkhorn-Divergence-based Knowledge Distillation
loss, for the distillation process. Whereas multiple distances such as Cosine
Distance or Euclidean Distance demand two sets share the sample size and are
sensitive to outliers (Zimek, Schubert, and Kriegel 2012), Sinkhorn divergence
does not enforce any requirement that relates to the number of samples and is
also robust to noise (Séjourné et al. 2019). Furthermore, compared with other
types of divergences such as KL divergence, the computation of Sinkhorn
divergence does not require two distributions to lie on the same probability
space. This is important because two languages might possess distinct features
that cannot be projected one-to-one, such as the vocabulary set. Consequently,
employing divergences different from Sinkhorn would need additional constraint
to the distillation loss. Lastly, Sinkhorn divergence is able to capture
geometric nature (Feydy et al. 2019) which has been shown to benefit myriad
cross-lingual and multilingual representation learning settings (Huang et al.
2021). We will empirically prove the superiority of Sinkhorn divergence in the
Experiment section.

Since the proposed module perpetuates the one-decoder employment, our
framework is able to explicitly correlate representations from two languages,
thus resolving the issue of two distant languages without demanding any
additional computation overhead. To evaluate the efficacy of our framework, we
proceed to conduct the experiments on myriad datasets containing document
/summary pairs of couples of distant languages, for example, English-to-
Chinese, English-to-Arabic, Japanese-to-English, etc. The empirical results
demonstrate that our model outperforms previous state-of-the-art Cross-Lingual
Summarization approaches. In sum, our contributions are three-fold:

  * •

We propose a Knowledge Distillation framework for Cross-Lingual Summarization
task, which seeks to enhance the summarization performance on distant
languages by aligning the cross-lingual with monolingual summarization,
through distilling the knowledge of monolingual teacher into cross-lingual
student model.

  * •

We propose a novel Knowledge Distillation loss using Optimal-Transport
distance, i.e. Sinkhorn Diveregence, with a view to coping with the spatial
discrepancy formed by the hidden representations produced by teacher and
student model.

  * •

We conducted extensive experiments in both high and low-resourced settings on
multiple Cross-Lingual Summarization datasets that belong to pairs of
morphologically and structurally distant languages, and found that our method
significantly outperforms other baselines in both automatic metrics and by
human evaluation.

## Related Work

### Neural Cross-Lingual Summarization

Due to the advent of Transformer architecture with its self-attention
mechanism, Text Generation has received ample attention from researchers
(Tuan, Shah, and Barzilay 2020; Lyu et al. 2021; Zhang et al. 2021),
especially Document Summarization (Zhang et al. 2020; Nguyen et al. 2021). In
addition to Monolingual Summarization, Neural Cross-Lingual Summarization has
been receiving a tremendous amount of interest, likely due to the burgeoning
need in cross-lingual information processing.

Conventional approaches designate a pipeline in two manners. The first one is
translate-then-summarize, which copes with the task by initially translating
the document into the target language and then performing the summarization
(Wan, Li, and Xiao 2010; Ouyang, Song, and McKeown 2019; Wan 2011; Zhang,
Zhou, and Zong 2016). The second approach is summarize-then-translate, which
firstly summarizes the document and then creates its translated version in the
target language (Lim, Kang, and Lee 2004; Orǎsan and Chiorean 2008; Wan, Li,
and Xiao 2010). Nonetheless, both of these approaches are vulnerable to error
propagation caused by undertaking multiple steps (Zhu et al. 2019).

Recent works apply a general architecture combined with large-scale training
to conduct Cross-Lingual Summarization. The main approach is to utilize the
multi-task framework, in which CLS task benefits from the process of making
use of other tasks such as Monolingual Summarization or Machine Translation
(Zhu et al. 2019). Further approaches design ancillary mechanisms such as
pointer-generator to exploit the translation scheme in the cross-lingual
summary (Zhu et al. 2020). Other work uses a pair of encoders and decoders to
co-operate the cross-lingual alignment with summarization (Cao, Liu, and Wan
2020).

### Optimal Transport in Natural Language Processing

Introduced in 19th century as a method to find the optimal solution to
transport a mass from one place to another destination, researchers have found
its use in a wide variety of scientific fields, such as computational fluid
mechanics (Benamou and Brenier 2000), economics (Carlier, Oberman, and Oudet
2015), physics (Cole et al. 2021), and notably machine learning (Peyré, Cuturi
et al. 2019; Cuturi 2013; Courty et al. 2016; Danila et al. 2006).

Recently, beside Contrastive Learning framework (Nguyen and Luu 2021; Pan et
al. 2021a, b), Optimal Transport has been omnivorously employed in Natural
Language Processing field, as used through Optimal Transport distance, for
instance Word Mover’s Distance (Werner and Laber 2019), to estimate the
necessary quantity of alignment. Its application includes text classification
(Kusner et al. 2015), capturing spatial alignment in word embedding (Alvarez-
Melis and Jaakkola 2018), machine translation (Chen et al. 2019), abstractive
summarization (Chen et al. 2019), etc. Nevertheless, the adaptation of Optimal
Transport distance, especially Sinkhorn divergence, for Neural Cross-Lingual
Summarization task has been attracting limited amount of research effort.

## Background

### Neural Cross-Lingual Summarization

Given a document $X^{L_{1}}=\\{x_{1},x_{2},…,x_{N}\\}$, a monolingual
summarization model’s task is to create a summary
$Y^{L_{1}}=\\{y^{L_{1}}_{1},y^{L_{1}}_{2},…,y^{L_{1}}_{M_{1}}\\}$, where both
$X^{L_{1}}$ and $Y^{L_{1}}$ are in language $L_{1}$. On the contrary, a cross-
lingual summarization model will produce a cross-lingual summary
$Y^{L_{2}}=\\{y^{L_{2}}_{1},y^{L_{2}}_{2},…,y^{L_{2}}_{M_{2}}\\}$ that is in
language $L_{2}$. It is worth noting here that $M_{1}<N$ and $M_{2}<N$.

Analogous to monolingual summarization, current state-of-the-art cross-lingual
summarization methods employ the Transformer-based architecture. Relying
mainly on self-attention mechanism, Transformer-based architecture consists of
an encoder and a decoder. The bidirectional self-attention in the encoder will
extract contextualized representations of the input, which will be fed to the
decoder to generate the output. Due to its generation nature, the decoder will
use unidirectional self-attention to learn the context of previously generated
tokens. During training procedure, the whole framework is updated based upon
the cross-entropy loss as follows

$\mathcal{L}_{\text{CLS}}=-\sum_{t=1}^{M_{2}}\log
P(y^{L_{2}}_{t}|y^{L_{2}}_{<t},X^{L_{1}})$ (1)

### Knowledge Distillation (KD)

Proposed by (Hinton, Vinyals, and Dean 2015), knowledge distillation is a
method to train a model, called the student, by leveraging valuable
information provided by soft targets output by another model, called the
teacher. In particular, the framework initially trains a model on one
designated task to extract useful features. Subsequently, given a dataset
$D=\\{(X_{1},Y_{1}),(X_{2},Y_{2}),…(X_{|D|},Y_{|D|})\\}$, where $|D|$ is the
size of the dataset, the teacher model will generate the output
$H^{T}_{i}=\\{\mathbf{h}^{T}_{1},\mathbf{h}^{T}_{2},…,\mathbf{h}^{T}_{L_{T}}\\}$
for each input $X_{i}$. Dependent on the researchers’ decision, the output
might be hidden representations or final logits. As a consequence, in order to
train the student model, the framework will use a KD loss that discriminates
the output of the student model
$H^{S}_{i}=\\{\mathbf{h}^{S}_{1},\mathbf{h}^{S}_{2},…,\mathbf{h}^{S}_{L_{S}}\\}$
given input $X_{i}$ from the teacher output $H^{T}_{i}$. Eventually, the KD
loss for input $X_{i}$ will possess the form as follows

$\mathcal{L}_{\text{KD}}=dist(H^{T}_{i},H^{S}_{i})$ (2)

where $dist$ is a distance function to estimate the discrepancy of teacher and
student outputs.

The explicated Knowledge Distillation framework has shown its efficiency in a
tremendous number of tasks, such as Neural Machine Translation (Tan et al.
2019; Wang et al. 2021; Li and Li 2021; Sun et al. 2020), Question Answering
(Hu et al. 2018; Arora, Khapra, and Ramaswamy 2019; Yang et al. 2020b), Image
Classification (Yang et al. 2020a; Chen, Chang, and Lee 2018; Fu et al. 2020),
etc. Nonetheless, its application for Neural Cross-Lingual Summarization has
received little interest.

## Methodology

Figure 2: Diagram of Knowledge Distillation Framework for Cross-Lingual
Summarization

To resolve the issue of distant languages, the output representations from two
vector spaces denoting two languages should be indistinguishable, or easily
transported from one space to another. In order to accomplish that goal, we
seek to relate the cross-lingual output of the student model to the
monolingual output of the teacher model, via utilizing Knowledge Distillation
framework and Sinkhorn Divergence calculation. The complete framework is
illustrated in Figure 2.

### Knowledge Distillation Framework for Cross-Lingual Summarization

We inherit the architecture of Transformer model for our module. In
particular, both the teacher and student model uses the encoder-architecture
paradigm combined with two fundamental mechanisms. Firstly, the self-attention
mechanism will attempt to learn the context of the tokens by attending tokens
among each other in the input and output document. Secondly, there is a cross-
attention mechanism to correlate the contextualized representations of the
output tokens to ones of the input tokens.

In our KD framework, we initiate the process by training the teacher model on
monolingual summarization task. In detail, given an input
$X^{L_{1}}=\\{x_{1},x_{2},…,x_{N}\\}$, the teacher model will aim to generate
its monolingual summary
$Y^{L_{1}}=\\{y^{L_{1}}_{1},y^{L_{1}}_{2},…,y^{L_{1}}_{M_{1}}\\}$. Similar to
previous monolingual summarization schemes, our model is trained by maximizing
the likelihood of the groundtruth tokens, which takes the cross-entropy form
as follows

$\mathcal{L}_{\text{MLS}}=-\sum_{t=1}^{M_{1}}\log
P(y^{L_{1}}_{t}|y^{L_{1}}_{<t},X^{L_{1}})$ (3)

After finetuning the teacher model, we progress to train the student model,
which also employs the Transformer architecture. Contrary to the teacher, the
student model’s task is to generate the cross-lingual output
$Y^{L_{2}}=\\{y^{L_{2}}_{1},y^{L_{2}}_{2},…,y^{L_{2}}_{M_{2}}\\}$ in language
${L_{2}}$, given the input document $X^{L_{1}}$ in language ${L_{1}}$. We
update the parameters of the student model by minimizing the objective
function that is formulated as follows

$\mathcal{L}_{\text{CLS}}=-\sum_{t=1}^{M_{2}}\log
P(y^{L_{2}}_{t}|y^{L_{2}}_{<t},X^{L_{1}})$ (4)

With a view to pulling the cross-lingual and monolingual representations
nearer, we implement a KD loss to penalize the large distance of two vector
spaces. Particularly, let
$H^{T}=\\{\mathbf{h}^{T}_{1},\mathbf{h}^{T}_{2},…,\mathbf{h}^{T}_{L_{T}}\\}$
denote the contextualized representations produced by the decoder of the
teacher model, and
$H^{S}=\\{\mathbf{h}^{S}_{1},\mathbf{h}^{S}_{2},…,\mathbf{h}^{S}_{L_{S}}\\}$
denote the representations from the decoder of the student model, we define
our KD loss as follows

$L_{\text{KD}}=dist(H^{T},H^{S})$ (5)

where $dist$ is the Optimal-Transport distance to evaluate the difference of
two representations, which we will delineate in the following section.

### Sinkhorn Divergence for Knowledge Distillation Loss

Due to the dilemma that the hidden representations of the teacher and student
model stay upon two disparate vector spaces (as they represent two different
languages), we will consider the distance of the two spaces as the distance of
two probability measures. To elaborate, we choose to adapt Sinkhorn
divergence, a variant of Optimal Transport distance, to calculate the
aforementioned spatial discrepancy. Let $H^{T}$, $H^{S}$ denote the
representations of the teacher decoder and the student decoder, we encode the
sample measures of them

$\alpha=\sum_{i=1}^{L_{T}}\alpha_{i}\delta_{\mathbf{h}_{i}^{T}},\quad\beta=\sum_{j=1}^{L_{S}}\beta_{j}\delta_{\mathbf{h}_{j}^{S}}$
(6)

where $\alpha$ and $\beta$ are probability distributions that satisfy
$\sum_{i=1}^{L_{T}}\alpha_{i}=1$ and $\sum_{j=1}^{L_{S}}\beta_{j}=1$.

Inspired by (Feydy et al. 2019), we estimate the difference of the
representations through determining the Sinkhorn divergence between them

$dist(H^{T},H^{S})=\text{OT}(\alpha,\beta)-\frac{1}{2}\text{OT}(\alpha,\alpha)-\frac{1}{2}\text{OT}(\beta,\beta)$
(7)

where

$\text{OT}(\alpha,\beta)=\sum_{i=1}^{N}\alpha_{i}f_{i}+\sum_{j=1}^{M}\beta_{j}g_{j}\\\
$ (8)

in which $f_{i},g_{j}$ are estimated by Sinkhorn loop. We thorougly delineate
the loop in Algorithm 1.

Algorithm 1 Sinkhorn loop

1:Probability distributions $\alpha,\beta$, regularization hyperparameter
$\varepsilon$, number of iterations $N_{I}$, log-sum-entropy function
$\text{LSE}_{k=1}^{N}(z_{k})=\log\sum_{k=1}^{N}\exp{(z_{k})}$, distance
function $C(\mathbf{x},\mathbf{y})=||\mathbf{x}-\mathbf{y}||^{2}$

2:for $i=1$ to $N_{I}$ do

3: Compute
$f_{i}=\varepsilon\cdot\text{LSE}_{k=1}^{L_{S}}[\log(\beta_{k})+\frac{1}{\varepsilon}g_{k}-\frac{1}{\varepsilon}C(\mathbf{h}^{T}_{i},\mathbf{h}^{S}_{k})]$

4: Compute
$g_{j}=\varepsilon\cdot\text{LSE}_{k=1}^{L_{T}}[\log(\alpha_{k})+\frac{1}{\varepsilon}f_{k}-\frac{1}{\varepsilon}C(\mathbf{h}^{T}_{k},\mathbf{h}^{S}_{j})]$

5:end for

### Training Objective

We amalgamate the Cross-Lingual Summarization and Knowledge Distillation
objective to obtain the ultimate objective function. Mathematically, for each
input, our training loss is computed as follows

$\mathcal{L}=\mathcal{L}_{\text{CLS}}+\lambda\cdot\mathcal{L}_{\text{KD}}$ (9)

where $\lambda$ is the hyperparameter that controls the influence of the
cross-lingual alignment of two vector spaces.

## Experiments

### Datasets

We evaluate the effectiveness of our methods on En2Zh and Zh2En datasets
processed by (Bai, Gao, and Huang 2021). We also inherit their minimum,
medium, and maximum settings in order to verify the effectiveness of our
method under limited-resourced settings. The sample size under each setting is
depicted in Table 2. Furthermore, to further evaluate the performance of our
method in various languages, we also preprocess datasets of Wikilingua (Ladhak
et al. 2020) in the manner that every sample is converted to a triple of
document, MLS summary, and CLS summary. We choose 4 variants of Wikilingua to
proceed our evaluation, i.e. English to Arabic (En2Ar), English to Japanese
(En2Ja), Japanese to English (Ja2En), and English to Vietnamese (En2Vi). It
should be noted here that (En, Ja), (En, Ar), (En, Zh), and (En, Vi) are all
couples of languages that are distant in terms of structure or morphology. The
statistics of the datasets is demonstrated in Table 1.

Dataset | $l_{\text{Input}}$ | $l_{\text{CLS}}$ | $l_{\text{MLS}}$
---|---|---|---
Zh2En | 105 | 19 | 19
En2Zh | 912 | 97 | 69
En2Ar | 1589 | 227 | 133
En2Ja | 1463 | 212 | 133
Ja2En | 2103 | 133 | 212
En2Vi | 1657 | 175 | 135
Table 1: Statistics of Cross-Lingual Summarization datasets. Scenarios | Minimum | Medium | Maximum | Full-dataset
---|---|---|---|---
Zh2En | 5,000 | 25,000 | 50,000 | 1,693,713
En2Zh | 1,500 | 7,500 | 15,000 | 364,687
Table 2: Dataset sizes of multiple low-resource scenarios for CLS datasets.

### Implementation Details

We initialize the encoder with multilingual BERT (Devlin et al. 2018), whereas
the decoder with Xavier intialization (Glorot and Bengio 2010). The dimensions
of our encoder and decoder hidden states are both 768. We use two seperate
Adam optimizers for encoder and decoder, and the learning rate for encoder and
decoder is 0.002 and 0.2, respectively. The model is trained with the warmup
phase of 25000 steps. We train the model on one Nvidia GeForce A100 GPU that
accumulates gradient every 5 steps. Moreover, we apply Dropout probability of
0.1 to all fully-connected layers in the model. The teacher and student model
shares the architecture and scale of parameters in our Knowledge Distillation
framework. To estimate the Sinkhorn divergence, we employ the entropic
regularization rate $\varepsilon$ of 0.0025 and the iteration length $N_{I}$
of 14. The weight $\lambda$ of KD Loss in Equation 9 is set to 1.

### Baselines

We compare our proposed architecture against the following baselines:

  * •

NCLS (Zhu et al. 2019): a Transformer-based model to conduct CLS.

  * •

NCLS + MS (Zhu et al. 2019): a multi-task framework that leverages an
auxiliary MS decoder to enhance cross-lingual summarization performance.

  * •

TLTran (Bai, Gao, and Huang 2021): a CLS pipeline that firstly performs MLS
and then utilizes a finetuned NMT model to translate the monolingual summary
into the target language.

  * •

MCLAS (Bai, Gao, and Huang 2021): a multi-task framework that sequentially
performs MLS, and CLS which is based upon the MLS result.

### Automatic Evaluation

#### Full-dataset Scenario

The experimental results under the full-dataset scenario are given in Table 3,
4, 5, 6, 7, and 8.

For Zh2En dataset, our method outperforms MCLAS model by ROUGE-1 of 1.3
points, ROUGE-2 of 4.0 points, ROUGE-3 of 0.4 point, and ROUGE-L of 1.2
points. Our model also improves the performance of NCLS model for dataset
En2Zh, with 0.6 point in ROUGE-1, 1.5 points in ROUGE-2, 0.1 point in ROUGE-3,
and 0.8 point in ROUGE-L. For Arabic language, our model achieves the
enhancement compared against NCLS model by 0.1 in ROUGE-1 score, 2.9 in
ROUGE-2 score, 1.6 in ROUGE-3 score, and 5.1 in ROUGE-L score. In En2Ja
dataset, we outperformed previous best method MCLAS by 0.6 point in ROUGE-1,
0.2 point in ROUGE-2, 0.2 point in ROUGE-3, and 0.5 point in ROUGE-L.

Additionally, for the reverse dataset Ja2En, our method significantly achieves
higher performance with the improvement of 1.0 point of ROUGE-1, 0.5 point of
ROUGE-2, 0.4 point of ROUGE-3, and 0.4 point of ROUGE-L, compared with MCLAS
model. Those results substantiate our hypothesis that our framework is able to
enhance the capability of apprehending and summarizing a document into a
summary of another distant language, since English alphabet does not have any
character in common with Japanese, Arabic, and Chinese counterparts.

For En2Vi dataset, our method also obtains notable improvement over other
state-of-the-art methods. As shown in Table 8, our model outperforms MCLAS
model by 0.1 in ROUGE-1, 2.9 in ROUGE-2, 1.6 in ROUGE-3, and 5.1 in ROUGE-L.
This demonstrates that our method is also capable of buttressing the model
capacity in situations where two languages are slightly morphologically or
structurally similar, since Vietnamese and English do share a number of
characters in their alphabets.

Model | R1 | R2 | R3 | RL
---|---|---|---|---
TLTran | 33.64 | 15.58 | 12.02 | 29.74
NCLS | 35.60 | 16.78 | 12.57 | 30.27
NCLS+MS | 34.84 | 16.05 | 12.28 | 29.47
MCLAS | 35.65 | 16.97 | 12.78 | 31.14
Our Model | 36.93 | 20.99 | 13.20 | 32.33
Table 3: Full-dataset Cross-Lingual Summarization results in Zh2En dataset Model | R1 | R2 | R3 | RL
---|---|---|---|---
TLTran | 30.20 | 12.20 | 11.79 | 27.02
NCLS | 44.16 | 24.28 | 17.13 | 30.23
NCLS+MS | 42.68 | 23.51 | 15.62 | 29.24
MCLAS | 42.27 | 24.60 | 16.07 | 30.09
Our Model | 44.75 | 25.76 | 17.20 | 31.05
Table 4: Full-dataset Cross-Lingual Summarization results in En2Zh dataset Model | R1 | R2 | R3 | RL
---|---|---|---|---
NCLS | 36.80 | 17.36 | 10.79 | 27.25
NCLS+MS | 35.53 | 17.01 | 10.33 | 26.36
MCLAS | 36.28 | 17.27 | 10.81 | 27.56
Our Model | 36.89 | 20.28 | 12.40 | 32.38
Table 5: Full-dataset Cross-Lingual Summarization results in En2Ar dataset Model | R1 | R2 | R3 | RL
---|---|---|---|---
NCLS | 29.55 | 15.99 | 10.25 | 23.03
NCLS+MS | 29.42 | 15.83 | 10.12 | 23.00
MCLAS | 29.60 | 16.08 | 10.14 | 33.20
Our Model | 30.21 | 16.27 | 10.46 | 23.90
Table 6: Full-dataset Cross-Lingual Summarization results in En2Ja dataset Model | R1 | R2 | R3 | RL
---|---|---|---|---
NCLS | 32.78 | 12.66 | 6.33 | 26.43
NCLS+MS | 32.50 | 12.02 | 6.15 | 26.41
MCLAS | 33.20 | 12.57 | 6.33 | 27.27
Our Model | 34.21 | 13.08 | 6.70 | 27.63
Table 7: Full-dataset Cross-Lingual Summarization results in Ja2En dataset Model | R1 | R2 | R3 | RL
---|---|---|---|---
NCLS | 36.75 | 16.37 | 8.04 | 28.69
NCLS+MS | 36.28 | 16.14 | 8.03 | 28.61
MCLAS | 36.31 | 15.91 | 7.75 | 28.62
Our Model | 37.38 | 16.20 | 8.09 | 28.97
Table 8: Full-dataset Cross-Lingual Summarization results in En2Vi dataset

#### Low-resource Scenario

We denote results of the experiments conducted under minimum, medium, and
maximum scenarios in Table 9, 10, and 11.

For the minimum setting, our model achieves the improvement over previous
methods. In particular, we outperformed MCLAS model by 1.3 points of ROUGE-1,
0.5 point of ROUGE-2, 0.2 point of ROUGE-3, and 0.3 point of ROUGE-L in Zh2En
dataset. For En2Zh dataset, we obtain an increase of 3.6 points in ROUGE-1,
0.6 point in ROUGE-2, 0.3 point in ROUGE-3, and 1.4 points in ROUGE-L.

Under the medium setting, the performance of our method is also higher than
MCLAS model with 0.1 point in ROUGE-1, 1.1 points in ROUGE-2, 0.5 point in
ROUGE-3, and 3.0 points in ROUGE-L. The improvement is more critical for
dataset En2Zh with an increase of 3.0 in ROUGE-1, 1.9 in ROUGE-2, 0.6 in
ROUGE-3, and 0.5 in ROUGE-L.

Last but not least, in maximum scenario, for dataset Zh2En, our gains compared
against MCLAS model are 0.4 point in ROUGE-1, 0.4 point in ROUGE-2, 0.5 point
in ROUGE-3, and 0.7 point in ROUGE-L. In dataset En2Zh, our improvements are
2.9 points in ROUGE-1, 0.3 point in ROUGE-2, 0.4 point in ROUGE-3, and 0.7
point in ROUGE-L.

Those aforementioned results have shown that our method is also capable of
elevating the Cross-Lingual Summarization performance when the available
training dataset is scarce.

Models | Zh2En | En2Zh
---|---|---
NCLS | 20.93/5.88/2.47/17.58 | 34.14/12.45/4.38/21.20
NCLS+MS | 20.50/5.45/2.22/17.25 | 33.96/12.38/4.36/21.07
MCLAS | 21.03/6.03/2.68/18.16 | 32.03/13.17/4.28/21.17
Our Model | 22.37/6.50/2.91/18.47 | 35.59/13.77/4.57/22.56

Table 9: Minimum Cross-Lingual Summarization Results

Models | Zh2En | En2Zh
---|---|---
NCLS | 26.42/8.90/4.49/22.05 | 35.98/15.88/8.97/23.79
NCLS+MS | 26.86/9.06/4.58/22.47 | 38.95/18.09/9.73/25.39
MCLAS | 27.84/10.41/4.91/24.12 | 37.28/18.10/9.48/25.26
Our Model | 27.97/11.51/5.37/27.16 | 40.30/20.01/10.05/25.79

Table 10: Medium Cross-Lingual Summarization Results

Models | En2Zh | Zh2En
---|---|---
NCLS | 29.05/10.88/6.56/24.32 | 40.18/19.86/10.33/26.52
NCLS+MS | 28.63/10.63/6.24/24.00 | 39.86/19.87/10.23/26.64
MCLAS | 30.73/12.26/6.98/26.51 | 38.35/19.75/10.64/26.41
Our Model | 31.08/12.70/7.45/27.16 | 41.24/20.01/11.00/27.06

Table 11: Maximum Cross-Lingual Summarization Results

### Human Evaluation

Because automatic metrics do not completely betray the quality of the methods,
we conduct further human evaluation for more precise assessment. To fulfil our
objective, we design two tests in order to elicit human judgements in two
manners.

In the first experiment, we present summaries generated by NCLS, MCLAS, our
model, and the gold summary, then asked seven professional English speakers to
indicate the best and worst summaries in terms of informativeness,
faithfulness, topic coherence, and fluency. We randomly sampled 50 summaries
from En2Vi dataset and 50 others from Ja2En dataset. The score of a model will
be estimated as the percentage of times it was denoted as the best minus the
percentage of times it was denoted as the worst.

For the second experiment, we decide to adapt Question Answering (QA) paradigm
to our framework. For each sample, we create two independent questions that
underscore the key information from the input document. Participants would
read and answer each question as best as they could. The score of a system
will be equal to the proportion of questions that the participants answer
correctly.

Fleiss’ Karpa scores of our experiments are shown in Table 12. It is obvious
that the scores prove a strong inter-agreement among the participants.

The experimental results in Table 13 indicate that our model generates
summaries that are conducive to human judgements, and have more likelihood to
preserve important content in the original documents than summaries of other
systems.

Test | Fleiss’ Kappa | Overall Agreement
---|---|---
Preference | 0.57 | 64.95%
QA | 0.64 | 82.15%
Table 12: Fleiss’ Kappa and Overall Agreement percentage of each human evaluation test. Higher score indicates better agreement. Models | Preference Score | QA score
---|---|---
NCLS | -0.123 | 51.11
MCLAS | 0.169 | 59.26
Our Model | 0.498 | 71.85
Gold Summary | 0.642 | 95.52
Table 13: Human evaluation

### Analysis on Distance Methods

We compare our implemented Sinkhorn Divergence with other distance methods.
Particularly, we perform the mean or max-pooling of the teacher and student
hidden representations. Subsequently, we evaluate the teacher and student
discrepancy via Cosine Similarity (CS) or Mean Squared Error (MSE) of two
pooled vectors. We show the numerical results in Table 14. The results
demonstrate the superiority of Sinkhorn Divergence over other approaches. We
hypothesize that those approaches do not efficaciously capture the geometry
nature of cross-lingual output representations.

Distance Methods | R-1 | R-2 | R-3 | R-L
---|---|---|---|---
Mean-CS | 44.20 | 24.54 | 16.96 | 30.27
Mean-MSE | 44.14 | 24.27 | 16.22 | 30.19
Max-CS | 44.29 | 25.65 | 17.07 | 30.82
Max-MSE | 44.23 | 24.61 | 16.44 | 30.21
Our Method | 44.75 | 25.76 | 17.20 | 31.05
Table 14: Results when applying different distance methods in En2Zh dataset
under full-dataset setting.

### Case Study

Figure 1 demonstrates case studies on the summarization results of NCLS,
MCLAS, and our method. It is noticeable that summaries generated by NCLS and
MCLAS missed important terms denoted by monolingual ones, such as “
_consumers_ ”, “ _cost_ ”, “ _underground_ ”, “ _parking lot_ ”, etc.
Consequently, those summaries unfortunately fail to mention the main idea of
the document, particularly the comparison of parking prices in sample 1, and
the announcement of “ _sinopec_ ” in sample 2. In contrast, our outputs cover
almost all of them, closely related the cross-lingual to the monolingual
summary. This shows that attracting CLS and MLS representations towards each
other through Sinkhorn divergence helps CLS model grasp key information, which
is an advantage of our proposed method.

### Impact of Sinkhorn Divergence on Geometric Distance of Cross-Lingual
Representations

We propose to adapt Sinkhorn Divergence to align the cross-lingual decoder
hidden states of the student model with monolingual decoder hidden states of
the teacher model. Nevertheless, whether this geometrically brings two sets of
representations nearer remains a quandary. To further verify the benefit of
leveraging Sinkhorn Divergence, we estimate the distances of those hidden
vectors by using other metrics, i.e. Cosine Similarity and Mean Squared Error.
Particularly, for each input, after getting the decoder to generate the hidden
vectors of the output tokens, we take the average of those vectors and measure
the distance between the mean of the vectors generated by the CLS model (NCLS,
MCLAS, and Our Model) with the mean of the vectors created by the MLS model.
We denote the expected value and standard deviation of each method in Table
15. As it can be obviously seen, employing Sinkhorn Divergence actually pulls
the vectors in the cross-lingual spaces towards one another.

Models | Cosine Similarity | Mean Squared Error
---|---|---
NCLS | 0.165 $\pm$ 0.038 | 19.434 $\pm$ 7.252
MCLAS | 0.064 $\pm$ 0.057 | 17.207 $\pm$ 4.028
Our Model |  0.034 $\pm$ 0.054 |  13.517 $\pm$ 4.013
Table 15: Results when applying different distance methods in Zh2En dataset
under full-dataset setting.

## Conclusion

In this paper, we propose a novel Knowledge Distillation framework to tackle
Neural Cross-Lingual Summarization for morphologically or structurally distant
languages. Our framework trains a monolingual teacher model, and then
finetunes the cross-lingual student model which is distilled knowledge from
the aforementioned teacher. Since the hidden representations of the teacher
and student model lie upon two different lingual spaces, we continually
proposed to adapt Sinkhorn Divergence to efficiently estimate the cross-
lingual discrepancy. Extensive experiments show that our method significantly
outperforms other approaches under both low-resourced and full-dataset
settings.

## References

  * Alvarez-Melis and Jaakkola (2018) Alvarez-Melis, D.; and Jaakkola, T. S. 2018.  Gromov-Wasserstein alignment of word embedding spaces.  _arXiv preprint arXiv:1809.00013_.
  * Arora, Khapra, and Ramaswamy (2019) Arora, S.; Khapra, M. M.; and Ramaswamy, H. G. 2019.  On knowledge distillation from complex networks for response prediction.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 3813–3822.
  * Bai, Gao, and Huang (2021) Bai, Y.; Gao, Y.; and Huang, H. 2021.  Cross-Lingual Abstractive Summarization with Limited Parallel Resources.  _arXiv preprint arXiv:2105.13648_.
  * Benamou and Brenier (2000) Benamou, J.-D.; and Brenier, Y. 2000.  A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem.  _Numerische Mathematik_ , 84(3): 375–393.
  * Bjerva et al. (2019) Bjerva, J.; Östling, R.; Veiga, M. H.; Tiedemann, J.; and Augenstein, I. 2019\.  What do language representations really represent?  _Computational Linguistics_ , 45(2): 381–389.
  * Cao, Liu, and Wan (2020) Cao, Y.; Liu, H.; and Wan, X. 2020.  Jointly Learning to Align and Summarize for Neural Cross-Lingual Summarization.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , 6220–6231.
  * Carlier, Oberman, and Oudet (2015) Carlier, G.; Oberman, A.; and Oudet, E. 2015.  Numerical methods for matching for teams and Wasserstein barycenters.  _ESAIM: Mathematical Modelling and Numerical Analysis_ , 49(6): 1621–1642.
  * Chen et al. (2019) Chen, L.; Zhang, Y.; Zhang, R.; Tao, C.; Gan, Z.; Zhang, H.; Li, B.; Shen, D.; Chen, C.; and Carin, L. 2019.  Improving sequence-to-sequence learning via optimal transport.  _arXiv preprint arXiv:1901.06283_.
  * Chen, Chang, and Lee (2018) Chen, W.-C.; Chang, C.-C.; and Lee, C.-R. 2018.  Knowledge distillation with feature maps for image classification.  In _Asian Conference on Computer Vision_ , 200–215. Springer.
  * Cole et al. (2021) Cole, S.; Eckstein, M.; Friedland, S.; and Życzkowski, K. 2021.  Quantum Optimal Transport.  _arXiv preprint arXiv:2105.06922_.
  * Courty et al. (2016) Courty, N.; Flamary, R.; Tuia, D.; and Rakotomamonjy, A. 2016.  Optimal transport for domain adaptation.  _IEEE transactions on pattern analysis and machine intelligence_ , 39(9): 1853–1865.
  * Cuturi (2013) Cuturi, M. 2013.  Sinkhorn distances: Lightspeed computation of optimal transport.  _Advances in neural information processing systems_ , 26: 2292–2300.
  * Danila et al. (2006) Danila, B.; Yu, Y.; Marsh, J. A.; and Bassler, K. E. 2006.  Optimal transport on complex networks.  _Physical Review E_ , 74(4): 046106.
  * Devlin et al. (2018) Devlin, J.; Chang, M.-W.; Lee, K.; and Toutanova, K. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_.
  * Feydy et al. (2019) Feydy, J.; Séjourné, T.; Vialard, F.-X.; Amari, S.-i.; Trouvé, A.; and Peyré, G. 2019.  Interpolating between optimal transport and MMD using Sinkhorn divergences.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2681–2690. PMLR.
  * Fu et al. (2020) Fu, S.; Li, Z.; Xu, J.; Cheng, M.-M.; Liu, Z.; and Yang, X. 2020.  Interactive knowledge distillation.  _arXiv preprint arXiv:2007.01476_.
  * Glorot and Bengio (2010) Glorot, X.; and Bengio, Y. 2010.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , 249–256. JMLR Workshop and Conference Proceedings.
  * Hinton, Vinyals, and Dean (2015) Hinton, G.; Vinyals, O.; and Dean, J. 2015.  Distilling the knowledge in a neural network.  _arXiv preprint arXiv:1503.02531_.
  * Hu et al. (2018) Hu, M.; Peng, Y.; Wei, F.; Huang, Z.; Li, D.; Yang, N.; and Zhou, M. 2018.  Attention-guided answer distillation for machine reading comprehension.  _arXiv preprint arXiv:1808.07644_.
  * Huang et al. (2021) Huang, K.-H.; Ahmad, W. U.; Peng, N.; and Chang, K.-W. 2021.  Improving Zero-Shot Cross-Lingual Transfer Learning via Robust Training.  _arXiv preprint arXiv:2104.08645_.
  * Kusner et al. (2015) Kusner, M.; Sun, Y.; Kolkin, N.; and Weinberger, K. 2015.  From word embeddings to document distances.  In _International conference on machine learning_ , 957–966. PMLR.
  * Ladhak et al. (2020) Ladhak, F.; Durmus, E.; Cardie, C.; and McKeown, K. 2020.  WikiLingua: A new benchmark dataset for cross-lingual abstractive summarization.  _arXiv preprint arXiv:2010.03093_.
  * Li and Li (2021) Li, Y.; and Li, W. 2021.  Data Distillation for Text Classification.  _arXiv preprint arXiv:2104.08448_.
  * Lim, Kang, and Lee (2004) Lim, J.-M.; Kang, I.-S.; and Lee, J.-H. 2004.  Multi-Document Summarization Using Cross-Language Texts.  In _NTCIR_.
  * Luo et al. (2021) Luo, F.; Wang, W.; Liu, J.; Liu, Y.; Bi, B.; Huang, S.; Huang, F.; and Si, L. 2021\.  VECO: Variable and Flexible Cross-lingual Pre-training for Language Understanding and Generation.
  * Lyu et al. (2021) Lyu, C.; Shang, L.; Graham, Y.; Foster, J.; Jiang, X.; and Liu, Q. 2021.  Improving Unsupervised Question Answering via Summarization-Informed Question Generation.  _arXiv preprint arXiv:2109.07954_.
  * Nguyen and Luu (2021) Nguyen, T.; and Luu, A. T. 2021.  Contrastive Learning for Neural Topic Model.  _Advances in Neural Information Processing Systems_ , 34.
  * Nguyen et al. (2021) Nguyen, T.; Luu, A. T.; Lu, T.; and Quan, T. 2021.  Enriching and controlling global semantics for text summarization.  _arXiv preprint arXiv:2109.10616_.
  * Orǎsan and Chiorean (2008) Orǎsan, C.; and Chiorean, O. A. 2008.  Evaluation of a cross-lingual romanian-english multi-document summariser.
  * Ouyang, Song, and McKeown (2019) Ouyang, J.; Song, B.; and McKeown, K. 2019.  A robust abstractive system for cross-lingual summarization.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 2025–2031.
  * Pan et al. (2021a) Pan, L.; Hang, C.-W.; Sil, A.; Potdar, S.; and Yu, M. 2021a.  Improved Text Classification via Contrastive Adversarial Training.  _arXiv preprint arXiv:2107.10137_.
  * Pan et al. (2021b) Pan, X.; Wang, M.; Wu, L.; and Li, L. 2021b.  Contrastive learning for many-to-many multilingual neural machine translation.  _arXiv preprint arXiv:2105.09501_.
  * Peyré, Cuturi et al. (2019) Peyré, G.; Cuturi, M.; et al. 2019.  Computational optimal transport: With applications to data science.  _Foundations and Trends® in Machine Learning_ , 11(5-6): 355–607.
  * Séjourné et al. (2019) Séjourné, T.; Feydy, J.; Vialard, F.-X.; Trouvé, A.; and Peyré, G. 2019.  Sinkhorn divergences for unbalanced optimal transport.  _arXiv preprint arXiv:1910.12958_.
  * Sun et al. (2020) Sun, H.; Wang, R.; Chen, K.; Utiyama, M.; Sumita, E.; and Zhao, T. 2020.  Knowledge distillation for multilingual unsupervised neural machine translation.  _arXiv preprint arXiv:2004.10171_.
  * Tan et al. (2019) Tan, X.; Ren, Y.; He, D.; Qin, T.; Zhao, Z.; and Liu, T.-Y. 2019.  Multilingual neural machine translation with knowledge distillation.  _arXiv preprint arXiv:1902.10461_.
  * Tuan, Shah, and Barzilay (2020) Tuan, L. A.; Shah, D.; and Barzilay, R. 2020.  Capturing greater context for question generation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, 9065–9072.
  * Wan (2011) Wan, X. 2011.  Using bilingual information for cross-language document summarization.  In _Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies_ , 1546–1555.
  * Wan, Li, and Xiao (2010) Wan, X.; Li, H.; and Xiao, J. 2010.  Cross-language document summarization based on machine translation quality prediction.  In _Proceedings of the 48th Annual Meeting of the Association for Computational Linguistics_ , 917–926.
  * Wang et al. (2021) Wang, F.; Yan, J.; Meng, F.; and Zhou, J. 2021.  Selective Knowledge Distillation for Neural Machine Translation.  _arXiv preprint arXiv:2105.12967_.
  * Werner and Laber (2019) Werner, M.; and Laber, E. 2019.  Speeding up Word Mover’s Distance and its variants via properties of distances between embeddings.  _arXiv preprint arXiv:1912.00509_.
  * Yang et al. (2020a) Yang, J.; Martinez, B.; Bulat, A.; and Tzimiropoulos, G. 2020a.  Knowledge distillation via adaptive instance normalization.  _arXiv preprint arXiv:2003.04289_.
  * Yang et al. (2020b) Yang, Z.; Shou, L.; Gong, M.; Lin, W.; and Jiang, D. 2020b.  Model compression with two-stage multi-teacher knowledge distillation for web question answering system.  In _Proceedings of the 13th International Conference on Web Search and Data Mining_ , 690–698.
  * Zhang et al. (2021) Zhang, A.; Wu, K.; Wang, L.; Li, Z.; Xiao, X.; Wu, H.; Zhang, M.; and Wang, H. 2021\.  Data Augmentation with Hierarchical SQL-to-Question Generation for Cross-domain Text-to-SQL Parsing.  _arXiv preprint arXiv:2103.02227_.
  * Zhang et al. (2020) Zhang, J.; Zhao, Y.; Saleh, M.; and Liu, P. 2020.  Pegasus: Pre-training with extracted gap-sentences for abstractive summarization.  In _International Conference on Machine Learning_ , 11328–11339. PMLR.
  * Zhang, Zhou, and Zong (2016) Zhang, J.; Zhou, Y.; and Zong, C. 2016.  Abstractive cross-language summarization via translation model enhanced predicate argument structure fusing.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 24(10): 1842–1853.
  * Zhu et al. (2019) Zhu, J.; Wang, Q.; Wang, Y.; Zhou, Y.; Zhang, J.; Wang, S.; and Zong, C. 2019.  NCLS: Neural cross-lingual summarization.  _arXiv preprint arXiv:1909.00156_.
  * Zhu et al. (2020) Zhu, J.; Zhou, Y.; Zhang, J.; and Zong, C. 2020.  Attend, translate and summarize: An efficient method for neural cross-lingual summarization.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , 1309–1321.
  * Zimek, Schubert, and Kriegel (2012) Zimek, A.; Schubert, E.; and Kriegel, H.-P. 2012.  A survey on unsupervised outlier detection in high-dimensional numerical data.  _Statistical Analysis and Data Mining: The ASA Data Science Journal_ , 5(5): 363–387.


# Testing $f(R)$ gravity with the simulated data of gravitational waves from
the Einstein Telescope

Yu Pan1, Yuan He1, Jing-Zhao Qi2∗, Jin Li3, Shuo Cao4†, Tonghua Liu4, and Jun
Wang 5 1\. School of Science, Chongqing University of Posts and
Telecommunications, Chongqing 400065, China;
2\. Department of Physics, College of Sciences, Northeastern University,
Shenyang 110004, China<EMAIL_ADDRESS>
3\. Department of Physics, Chongqing University, Chongqing 400030, China;
4\. Department of Astronomy, Beijing Normal University, Beijing 199875, China;
<EMAIL_ADDRESS>
5\. School of Physics and Astronomy, Yunnan University, Kunming 650091, China

###### Abstract

In this paper we analyze the implications of gravitational waves (GWs) as
standard sirens on the modified gravity models by using the third-generation
gravitational wave detector, i.e., the Einstein Telescope. Two viable models
in $f(R)$ theories within the Palatini formalism are considered in our
analysis ($f_{1}(\mathcal{R})=\mathcal{R}-\frac{\beta}{\mathcal{R}^{n}}$ and
$f_{2}(\mathcal{R})=\mathcal{R}+\alpha\ln{\mathcal{R}}-\beta$), with the
combination of simulated GW data and the latest electromagnetic (EM)
observational data (including the recently released Pantheon type Ia
supernovae sample, the cosmic chronometer data, and baryon acoustic
oscillation distance measurements). Our analysis reveals that the standard
sirens GWs, which provide an independent and complementary alternative to
current experiments, could effectively eliminate the degeneracies among
parameters in the two modified gravity models. In addition, we thoroughly
investigate the nature of geometrical dark energy in the modified gravity
theories with the assistance of $Om(z)$ and statefinder diagnostic analysis.
The present analysis makes it clear-cut that the simplest cosmological
constant model is still the most preferred by the current data. However, the
combination of future naturally improved GW data most recent EM observations
will reveal the consistency or acknowledge the tension between the
$\Lambda$CDM model and modified gravity theories.

###### Subject headings:

modified gravity, gravitational wave

##  1\. Introduction

The measurement of cosmic acceleration is crucial for our understanding of the
nature of the Universe. A number of scenarios have been proposed to explain
this remarkable discovery, which fall into two general categories. In the
framework of one approach, an exotic component called dark energy is
introduced to act the driving force behind the cosmic acceleration (Cao et
al., 2011; Cao & Zhu, 2014). Until now, all of the independent astrophysical
observations are well consistent with the existence of vacuum energy density
or a non-vanishing cosmological constant $\Lambda$ (Cao & Zhu, 2012; Cao,
Covone & Zhu, 2012; Cao et al., 2012; Pan et al., 2015; Cao et al., 2015a, b,
2017a, 2017b; Li et al., 2017; Risaliti & Lusso, 2018; Liu et al., 2019; Ma et
al., 2019). However, such simplest cosmological constant model is still
embarrassed by the coincidence problem and fine-tuning problem, which
triggered great efforts to understand the evolution of the Universe at early
time (Abbott et al., 2017; Qi et al., 2019a; Cao et al., 2019a, 2020; Zheng et
al., 2020). On the other hand, such an accelerated expansion is possible with
the modification of the gravitational sector in Einstein-Hilbert action (Li et
al., 2007; Carroll et al., 2005; Capozziello et al., 2005; De Felice &
Tsujikawa, 2010; Sotiriou & Faraon, 2010). Following this direction, a lot of
extended theories of gravity have been proposed in the past decades, i.e.,
Gauss-Bonnet gravity (Nojiri & Odintsov, 2005), $f(T)$ gravity (Bengochea &
Ferraro, 2009; Cai et al., 2015; Qi et al., 2016), and $f(R)$ gravity (De
Felice & Tsujikawa, 2010; Sotiriou & Faraon, 2010; Wang et al., 2018).
Especially, the $f(R)$ theories of gravity (in which $f(R)$ is a function of
the Ricci scalar $R$) are widely employed in modern cosmology, benefit from
its advantages of describing the large scale structure distribution of the
universe (Dos Santos et al., 2016; Voivodic et al., 2017; Wang & Wang, 2013)
and explaining the late-time cosmic acceleration. Note that there are mainly
two different variational approaches providing various interesting mechanisms,
i.e., the metric formalism and the Palatini formalism (Li et al., 2007; Santos
et al., 2008; Capozziello & Francaviglia, 2008). In the former formalism, by
assuming the connections to be the Christoffel symbols defined in terms of the
metric, one can get the generalized field equations by varying the action with
respect to the metric, which leads to the fourth-order equations difficult to
deal with in practice. Furthermore, it was found in the recent studies that
some models within such framework may not only exhibit violent instabilities
in the weak gravity regime of matter (Dolgov & Kawasaki, 2003; Nojiri &
Odintsov, 2003; Wang et al., 2010), but also fail to produce a standard
matter-dominated era followed by an accelerated expansion (Amendola et al.,
2006a, b). In this paper, we will mainly investigate the constraints on the
Palatini formalism model ($R=g^{\mu\nu}R_{\mu\nu}(\Gamma)$), in the framework
of which the metric and connection are treated as independent dynamical
variables in the action (Vollick, 2003; De Felice & Tsujikawa, 2010; Sotiriou
& Faraon, 2010). The advantage of the Palatini approach is that the resultant
equations are of second order, which are different from the ones derived in
the metric approach but more concordant with the field equations in other
branches of physics (Vollick, 2003; De Felice & Tsujikawa, 2010; Sotiriou &
Faraon, 2010). Specifically, Palatini gravity have been studied in the context
of dark energy models, with the late-time constraints from various
astrophysical observations (Bekov et al., 2020; Baghram & Rahvar, 2020; Aoki &
Shimada, 2020; Leanizbarrutia et al., 2017; Bohmer et al., 2013; Capozziello
et al., 2013; Harko et al., 20212; Gu et al., 2018; Rosa et al., 2020, 2017;
Borowiec et al., 2012). Meanwhile, it is interesting to note that some models
could exhibit very promising behaviors, i.e., they not only correctly fit
different observational data, but also well reproduce the cosmological
evolution with early and late-time acceleration (Yang & Chen, 2009; Fay &
Tavakol, 2007; Seokcheon, 2007). We refer to Sotiriou & Faraon (2010) for a
review on the $f(R)$ gravity.

From observational point of view, it is important to investigate whether these
modified gravity theories are indeed compatible with different kinds of
observational data. Currently, such analysis has been performed in the EM
domain, focusing on different types of standard cosmological probes: type Ia
supernovae (SNe Ia), baryon acoustic oscillations (BAO), cosmic microwave
background (CMB), and intermediate-luminosity radio quasars (ILQSO)
(Amarzguioui et al., 2006; Fay & Tavakol, 2007; de la Cruz-Dombriz et al.,
2016; Song et al., 2007; Xu et al., 2018; Santos et al., 2008). For instance,
the first attempt to derive constraints on the $f(R)$ model
$f(R)=R-\frac{\beta}{R^{n}}$ was presented in Amarzguioui et al. (2006), which
was then extended in (Fay & Tavakol, 2007) by employing the first year Super-
Nova Legacy Survey (SNLS) data, the baryon acoustic oscillation peak in the
SDSS luminous red galaxy sample and the CMB shift parameter. It was found that
in subsequent analysis that such $f(R)$ model could produce the sequence of
radiation-dominated, matter-dominated, and accelerating periods without the
inclusion of dark energy (Santos et al., 2008). However, it is necessary and
important to have other complementary cosmological probes in the GW domain,
especially the inspiraling and merging compact binaries consisting of neutron
stars (NSs) and black holes (BHs) (Schutz, 1986). It is well known that the
first direct detection of gravitational waves (GWs) by the LIGO/Virgo
collaboration (Abbott et al., 2016) has opened the era of GW astronomy. As a
promising high-redshift complementary tool to SNe Ia, the greatest advantage
of GWs lies in the fact that the distance calibration of such standard siren
is independent of any other distance ladders (Taylor et al., 2012). Using the
waveform signal to directly measure the luminosity distance $D_{L}$ to the GW
sources, the possible cosmological application of these standard sirens has
been extensively discussed in the literature (Holz & Hughes, 2005; Dalal et
al., 2006; Zhao et al., 2011; Cai & Yang, 2017; Wei, 2019; Qi et al., 2019b;
Liao, 2019; Cao et al., 2019b; Wu et al., 2020; Liu et al., 2020). More
interestingly, it was originally proposed (D’Agostino & Nunes, 2019) that
standard sirens from binary neutron star mergers could provide strong
observational bounds on modified gravity theories (scalar-tensor theories).
Specially, the final results showed that the $f(R)$ gravity (i.e., the Hu-
Sawicki $f(R)$ gravity model) can be tested at very high precision (95%
confidence level), based on the third generation ground-based GW detector
(Einstein telescope). Such conclusion was further supported by recent
measurements of Hubble constant in the context of modified gravity theories
($f(R)$ and $f(T)$ models), based on the combined geometrical EM data sets
obtained in a model-independent way (D’Agostino & Nunes, 2020). In this paper,
we will mainly investigate the constraints on two viable models in $f(R)$
theories within the Palatini formalism
($f_{1}(\mathcal{R})=\mathcal{R}-\frac{\beta}{\mathcal{R}^{n}}$ and
$f_{2}(\mathcal{R})=\mathcal{R}+\alpha\ln{\mathcal{R}}-\beta$), focusing on an
updated sample of GW events based on the third-generation GW ground-based
detector, Einstein Telescope (ET). The latest electromagnetic (EM)
observational data, including the recently released Pantheon type Ia
supernovae sample, the cosmic chronometers sample and the recent baryon
acoustic oscillation measurements are also included in our analysis, in order
to achieve a reasonable and compelling constraints on the modified gravity
models in both the electromagnetic (EM) and gravitational wave (GW) window.

This paper is organized as follows. In Sec. 2, we review the $f(R)$ theories
of gravity, mainly with two variable models within the Palatini formalism. The
observational data (GW+EM) and the corresponding constraint results are
introduced in Sec. 3 and 4, respectively. In Sec. 5, we present some model
diagnostics which results in distinction between $f(R)$ models and
$\Lambda$CDM. Finally, we discuss our results and present our main conclusions
in Sec. 6.

##  2\. The $f(R)$ theory in Palatini formalism

We briefly summarize the basic features of $f(R)$ theories of gravity in the
Palatini formalism (Sotiriou & Faraon, 2010; Will, 1981). The modified
Einstein-Hilbert action is given by

$S=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}f(\mathcal{R})+S_{m}\,,$ (1)

where $G$ is the gravitational constant, $g$ is the determinant of the metric
tensor, and $S_{m}$ is the standard action for the matter fields. In the
Palatini formalism, the metric and the connection are completely independent
variables. Therefore, the definition of the Ricci scalar can be expressed as
$\mathcal{R}=g^{\alpha\beta}\mathcal{R}_{\alpha\beta}(\tilde{\Gamma}_{\mu\nu}^{\rho})$,
where $\tilde{\Gamma}_{\mu\nu}^{\rho}$ is the affine connection, which is
different from the Levi-Civita connection $\Gamma_{\mu\nu}^{\rho}$ in the
metric formalism.

Varying the action with respect to the metric and the connection, the field
equations can be obtained as

$f^{\prime}\mathcal{R}_{(\mu\nu)}-\frac{f}{2}g_{\mu\nu}=8{\pi}GT_{\mu\nu}\,,$
(2) $\tilde{\nabla}_{\alpha}[f^{\prime}(\mathcal{R})\sqrt{-g}g^{\mu\nu}]=0\,,$
(3)

with $f^{\prime}\equiv df/d\mathcal{R}$ and $f^{\prime\prime}\equiv
d^{2}f/d\mathcal{R}^{2}$. Here
$T_{\mu\nu}\equiv-\frac{2}{\sqrt{-g}}\frac{\delta S_{m}}{\delta g^{\mu\nu}}$
denotes the energy-momentum tensor of matter, and $\tilde{\nabla}_{\alpha}$
represents the covariant derivative corresponding to the affine connection
$\tilde{\Gamma}_{\mu\nu}^{\rho}$. It is obvious that Eq. (2) will reduce to
the Einstein’s field equations when $f(\mathcal{R})=\mathcal{R}$, while the
Einstein’s field equations in $\Lambda$CDM model can be reproduced when
$f(\mathcal{R})=\mathcal{R}-2\Lambda$.

By introducing a metric conformal
$h_{\mu\nu}=f^{\prime}(\mathcal{R})g_{\mu\nu}$ to Eq. (3), one could derive
the definition of the usual Levi-Civita connection in terms of the new metric
$h_{\mu\nu}$:

$\tilde{\Gamma}^{\lambda}_{\
\mu\nu}=\frac{h^{\lambda\sigma}}{2}(h_{\nu\sigma,\mu}+h_{\mu\sigma,\nu}-h_{\mu\nu,\sigma}).$
(4)

Then the relation between the generalized affine connection
$\tilde{\Gamma}^{\lambda}_{\ \mu\nu}$ and the metric $g_{\mu\nu}$ can be
obtained as

$\tilde{\Gamma}^{\lambda}_{\ \mu\nu}=\Gamma^{\lambda}_{\
\mu\nu}+\frac{1}{2f^{\prime}}[2\delta^{\lambda}_{\
(\mu}\partial_{\nu)}f^{\prime}-g^{\lambda\tau}g_{\mu\nu}\partial_{\tau}f^{\prime}].$
(5)

Now the relationship between the generalized Ricci tensor
$\mathcal{R}_{\mu\nu}$ and the original Ricci tensor $R_{\mu\nu}$ is given by

$\mathcal{R}_{\mu\nu}=R_{\mu\nu}-\frac{3}{2}\frac{\nabla_{\mu}f^{\prime}\nabla_{\nu}f^{\prime}}{f^{\prime
2}}+\frac{\nabla_{\mu}\nabla_{\nu}f^{\prime}}{f^{\prime}}+\frac{1}{2}g_{\mu\nu}\frac{\nabla^{\mu}\nabla_{\nu}f^{\prime}}{f^{\prime}},$
(6)

where

$\mathcal{R}_{\mu\nu}\equiv\tilde{\Gamma}^{\alpha}_{\mu\nu,\alpha}-\tilde{\Gamma}^{\alpha}_{\mu\alpha,\nu}+\tilde{\Gamma}^{\alpha}_{\alpha\lambda}\tilde{\Gamma}^{\lambda}_{\mu\nu}-\tilde{\Gamma}^{\alpha}_{\mu\lambda}\tilde{\Gamma}^{\lambda}_{\alpha\nu}\,.$
(7)

In the framework of the flat Friedman-Robertson-Walker (FRW) metric

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)(\mathrm{d}x^{2}+\mathrm{d}y^{2}+\mathrm{d}z^{2})\,$
(8)

we take the energy-momentum tensor of matter as a perfect fluid

$T_{\mu\nu}=(\rho+p)U_{\mu}U_{\nu}+pg_{\mu\nu},$ (9)

where $\rho$ is the energy density, $p$ is the pressure and the four-velocity
$U_{\mu}$ satisfies $U_{\mu}U^{\mu}=-1$ and $U^{\mu}U_{\mu;\nu}=0$. The
Friedmann equation of the $f(R)$ model can be written as

$6(H+\frac{1}{2}\frac{\dot{f}^{\prime}}{f^{\prime}})^{2}=\frac{3f-\mathcal{R}f^{\prime}}{f^{\prime}},$
(10)

where $H\equiv\dot{a}/a$ is the Hubble parameter, and the overdot denotes the
derivative with respect to the cosmic time $t$. Based on the trace of Eq. (2),
i.e., $\mathcal{R}f^{\prime}(\mathcal{R})-2f(\mathcal{R})=\kappa T$ and the
conservation equation, i.e., $\dot{\rho}_{\rm m}+3H\rho_{\rm m}=0$, we obtain

$\dot{\mathcal{R}}=\frac{3H\kappa\rho_{\rm
m}}{\mathcal{R}f^{\prime\prime}-f^{\prime}}.$ (11)

Substituting Eq. (11) into Eq. (10), one can numerically compute the Friedman
equation as

$H^{2}=\frac{3f-\mathcal{R}f^{\prime}}{6f^{\prime}\eta^{2}},$ (12)

where

$\eta=1-\frac{3}{2}\frac{f^{\prime\prime}}{f^{\prime}}\frac{\mathcal{R}f^{\prime}-2f}{\mathcal{R}f^{\prime\prime}-f^{\prime}}.$
(13)

According to the definition of the angular diameter distance at redshift $z$,
one also has

$\displaystyle D_{\rm A}(z)$ $\displaystyle=$
$\displaystyle\frac{1}{1+z}\int_{0}^{z}\frac{\mathrm{d}z}{H(z)}$
$\displaystyle=$
$\displaystyle\frac{1}{3}(2f-\mathcal{R}f^{\prime})^{-\frac{1}{3}}\int^{\mathcal{R}_{z}}_{\mathcal{R}_{0}}{\frac{\mathcal{R}f^{\prime\prime}-f^{\prime}}{(2f-\mathcal{R}f^{\prime})^{\frac{2}{3}}}\frac{d\mathcal{R}}{H(\mathcal{R})}}.$

Following the procedure proposed in Fay & Tavakol (2007), the generalized
Friedmann equation can be written in terms of redshift $z=a_{0}/a-1$ and the
matter density parameter $\Omega_{m}\equiv\kappa\rho_{m0}/(3H_{0}^{2})$ as

$\frac{H^{2}}{H_{0}^{2}}=\frac{3\Omega_{m0}(1+z)^{3}+f/H_{0}^{2}}{6f^{\prime}\xi^{2}}\,,$
(15)

where

$\xi=1+\frac{9}{2}\,\frac{f^{\prime\prime}}{f^{\prime}}\,\frac{H_{0}^{2}\Omega_{m0}(1+z)^{3}}{Rf^{\prime\prime}-f^{\prime}}.$
(16)

The trace of Eq. (2) leads to

$\mathcal{R}f^{\prime}-2f=-3H_{0}^{2}\Omega_{m0}(1+z)^{3}\,,$ (17)

based on which one may find that Eq. (15) will reduce to General Relativity
for the Einstein-Hilbert Lagrangean case ($f(\mathcal{R})=\mathcal{R}$).
Moreover, it can be clearly seen that by taking
$f(\mathcal{R})=\mathcal{R}-\beta/\mathcal{R}^{n}$, Eq. (17) evaluated at
$z=0$ imposes the following relation

$\beta=\frac{\mathcal{R}_{0}^{n+1}}{n+2}\,\left(1-\frac{3\Omega_{m}H^{2}_{0}}{\mathcal{R}_{0}}\right)\,,$
(18)

where the present value of the Ricci scalar $\mathcal{R}_{0}$ is determined
from the algebraic equation resulting from Eq. (17) and (15) for $z=0$.
Therefore, in the framework of $f(R)$ parametrization in the form of
$f(\mathcal{R})=\mathcal{R}-\beta/\mathcal{R}^{n}$ and
$f(\mathcal{R})=\mathcal{R}+\alpha{ln\mathcal{R}}-\beta$ (Fay & Tavakol,
2007), the constraint results on the $f(R)$ model can be shown in
the($n,\Omega_{m}$) and ($\alpha,\Omega_{m}$) plane.

Now it is necessary to comment on the theoretical applicability of the
parameterizations used in this paper. For the former form of
$f(\mathcal{R})=\mathcal{R}-\beta/\mathcal{R}^{n}$, Capozziello et al. (2003,
2004) have extensively discussed the possibility of using the models in the
metric approach to generate a late-time acceleration of the Universe. Their
results showed that for $n>0$ the metric approach was unable to give rise to a
standard matter-dominated era followed by a cosmic acceleration (Amendola et
al., 2006a). However, such theories using the Palatini approach could provide
an effective solution to the above difficulties, as was proved in the
subsequent analysis of (Fay & Tavakol, 2007). Similarly, for the latter form
of $f(\mathcal{R})=\mathcal{R}+\alpha{ln\mathcal{R}}-\beta$, the numerical
results have also proved the capability of the Palatini approach to produce
the sequence of radiation-dominated, matter-dominated, and de Sitter periods
(Fay & Tavakol, 2007). Therefore, in this paper we shall focus on a variable
$f(R)$ theories recently put forward within the Palatini formalism, the
cosmological dynamics of which is still not fully understood up to now.
Despite these advantages, some concerns have been expressed on several
problems of Palatini $f(R)$ gravity, due to the differential structure of its
field equations (Sotiriou & Faraon, 2010), as well as its possible difficulty
in passing the solar system tests and providing the correct Newtonian limit
(Flanagan, 2004a, b). From theoretical point of view, such issue could be
potentially fixed by the inclusion of extra terms quadratic in the Ricci
and/or Riemann tensor, which is hard to be rigorously accounted in context of
cosmological studies like in this paper.

##  3\. GW simulations and EM observations

In this subsection we briefly introduce the method of simulating GW events
from Einstein Telescope. Compared with the current advanced ground-based
detectors (i.e., the advanced LIGO and Virgo detectors), such third-generation
GW detector composed of three collocated underground detectors is designed to
be ten times more sensitive, especially in the frequency range of $1-10^{4}$
Hz 111The Einstein Telescope Project, https://www.et-gw.eu/et/.

We focus on the GW signals generated by the coalescence of binary systems with
component masses $m_{1}$ and $m_{2}$, defining the total mass $M=m_{1}+m_{2}$,
the symmetric mass ratio $\eta=m_{1}m_{2}/M^{2}$, and the chirp mass
$\mathcal{M}_{c}=M\eta^{3/5}$. From observational point of view, the observed
chirp mass in the observer frame can be written as $\mathcal{M}_{c,\rm
obs}=(1+z)\mathcal{M}_{c,\rm phys}$. Note that the GW amplitude depends on the
so-called chirp mass and the luminosity distance $D_{L}$, while the former can
be measured from the GW signal’s phasing. Therefore, the chirping GW signals
from inspiraling compact binary stars (NS and BH) can provide an absolute
measure of the luminosity distance. Following the previous analysis of Zhao et
al. (2011); Sathyaprakash et al. (2010), the mass of each black hole and
neutron star is assumed to be uniformly distributed in the range [3,10]
$M_{\bigodot}$ and [1,2] $M_{\bigodot}$. The ratio between NS-NS and BH-NS
systems is taken to be $0.03$, as is predicted by the Advanced LIGO-Virgo
network (Abadie, 2010).

GW detectors based on the ET could measure the strain (or the time domain
waveform $h(t)$) of GWs, the Fourier transform $\mathcal{H}(f)$ of which could
be derived by applying the stationary phase approximation

$\displaystyle\mathcal{H}(f)=\mathcal{A}f^{-7/6}\exp[i(2\pi
ft_{0}-\pi/4+2\psi(f/2)-\varphi_{(2.0)})],$ (19)

where the Fourier amplitude $\mathcal{A}$ is given by

$\displaystyle\mathcal{A}=$
$\displaystyle~{}~{}\frac{1}{D_{L}}\sqrt{F_{+}^{2}(1+\cos^{2}(\iota))^{2}+4F_{\times}^{2}\cos^{2}(\iota)}$
$\displaystyle~{}~{}\times\sqrt{5\pi/96}\pi^{-7/6}\mathcal{M}_{c}^{5/6},$ (20)

See Zhao et al. (2011) for the definition of the epoch of the merger $t_{0}$,
the functions of $\psi$ and $\varphi_{(2.0)}$, the angle of inclination of the
binary’s orbital angular momentum $\iota$, and the two interferometers’
antenna pattern functions ($F_{+}$, $F_{\times}$) that depend on the position
of GW source, as well as the location and orientation of the detector.

The uncertainty of luminosity distance extracted from the GW signals can be
divided into two parts. The first source of uncertainty is quantified by the
instrumental uncertainty with the form of

$\displaystyle\sigma_{D_{L}}^{\rm
inst}\simeq\sqrt{\left\langle\frac{\partial\mathcal{H}}{\partial
D_{L}},\frac{\partial\mathcal{H}}{\partial D_{L}}\right\rangle^{-1}},$ (21)

Assuming that this parameter is uncorrelated with any other GW parameters, we
add a factor of 2 in the instrumental error to take the effect of inclination
into consideration (Li, 2015)

$\displaystyle\sigma_{D_{L}}^{\rm inst}\simeq\frac{2D_{L}}{\rho}.$ (22)

Based on the noise power spectral density $S_{h}(f)$ (PSD) of ET given in Zhao
et al. (2011), the Signal-to-Noise ratio (SNR) of the network of three
independent interferometers can be calculated as

$\rho=\sqrt{\sum\limits_{i=1}^{3}\left\langle\mathcal{H}^{(i)},\mathcal{H}^{(i)}\right\rangle}.$
(23)

Here the inner product is defined as

$\left\langle{a,b}\right\rangle=4\int_{f_{\rm lower}}^{f_{\rm
upper}}\frac{\tilde{a}(f)\tilde{b}^{\ast}(f)+\tilde{a}^{\ast}(f)\tilde{b}(f)}{2}\frac{df}{S_{h}(f)},$
(24)

where the lower limit in frequency of ET is $f_{\rm lower}=1$Hz, the upper
cutoff frequency is set as $f_{\rm upper}=2/(6^{3/2}2\pi M_{\rm obs})$ and
$M_{\rm obs}=(1+z)M_{\rm phys}$ is the observed total mass. Note that a GW
detection is confirmed if a network SNR is $\rho>8.0$, the SNR threshold
currently used by LIGO/Virgo network (Cai & Yang, 2017). The second source of
uncertainty is generated by gravitational lensing, which is always modeled as
$\sigma_{D_{L}}^{lens}/D_{L}=0.05z$ (Sathyaprakash et al., 2010; Li, 2015; Cai
& Yang, 2017). Thus, the total uncertainty of $D_{L}$ is taken to be

$\displaystyle\sigma_{D_{L}}$ $\displaystyle~{}~{}=\sqrt{(\sigma_{D_{L}}^{\rm
inst})^{2}+(\sigma_{D_{L}}^{\rm lens})^{2}}$
$\displaystyle~{}~{}=\sqrt{\left(\frac{2D_{L}}{\rho}\right)^{2}+(0.05zD_{L})^{2}}.$
(25)

Thus, the total uncertainty on luminosity distance is

$\displaystyle\sigma_{D_{L}}^{2}=(\sigma_{D_{L}}^{\rm
inst})^{2}+(\sigma_{D_{L}}^{\rm lens})^{2}.$ (26)

Focusing on the expected rates of BNS and BHNS detections per year for the ET,
only a small fraction of binary mergers may have the observation of detectable
EM counterparts like short and intense bursts of $\gamma$-rays (SGRB). In this
work we follow the recent analysis of Cai & Yang (2017) and assume that ET can
detect up to 1000 GW events in a 10 year observation (with measurable source
redshift). As for the redshift distribution of the observable GW sources, we
follow the function taking into account evolution and stellar synthesis
(Schneider et al., 2001; Cutler & Holz, 2009; Sathyaprakash et al., 2010),
which has been widely applied in gravitational wave cosmology (Qi et al.,
2019b). The fiducial cosmological model we choose is the concordant model
based on the most recent Plank results (Ade et al., 2016). The GW mock data
covering the redshift range $z\sim 5$ are shown in Fig. 1, with the luminosity
distance measurements obeying the Gaussian distribution
$D_{L}^{mean}=\mathcal{N}(D_{L}^{\rm fid},\sigma_{D_{L}})$.

Figure 1.— The 1000 simulated GW events from the third-generation
gravitational wave detector, i.e., the Einstein Telescope.

In this analysis, we also include the latest observational data derived in the
electromagnetic (EM) window, including the recently released Pantheon type Ia
supernovae sample, the baryon acoustic oscillation (BAO) observations, and
Hubble parameter measurements obtained through the cosmological-model-
independent method using the so-called cosmic chronometers (Jimenez & Loeb,
2002; Jimenez et al., 2003). The type Ia supernovae sample used in our
analysis include 1048 SNe Ia observed by the Pan-STARRS1 (PS1) Medium Deep
Survey, covering the redshift range of $0.01<z<2.3$ (Scolnic et al., 2018).
For the BAO data, we use the six measurement of distance ratio
$\frac{d_{A}(z*)}{D_{V}(z_{BAO})}$ in the redshift range $0.106\leq z\leq
0.73$, from the 6-degree Field Galaxy Survey (6dFGS), the WiggleZ galaxy
survey, and the Sloan Digital Sky Survey (SDSS-DR7) catalog combined with
galaxies from 2dFGRS (Beutler et al., 2011; Ross et al., 2015; Alam et al.,
2017; Zhao et al., 2018; Des Bourboux et al., 2017). Here $D_{V}$ is the
dilation scale

$D_{V}(z)=\left(d_{A}(z)^{2}\frac{z}{H(z)}\right)^{1/3},$ (27)

and $d_{A}$ is the co-moving angular diameter distance (Xu et al., 2018)

$d_{A}(z)=\int^{z}_{0}\frac{dz^{\prime}}{H(z^{\prime})}.$ (28)

The specific BAO data and their corresponding fitting procedures are
summarized in the source papers (Simon et al., 2005; Stern et al., 2010;
Moresco et al., 2012; Moresco, 2015; Moresco et al., 2016; Zhang et al.,
2014). For the Hubble parameter data, 30 measurements of $H(z)$ covering the
redshift range $0.07\leq z\leq 1.965$ have been obtained, based on the
differential age evolution of passively evolving early-type galaxies (Moresco,
2015; Moresco et al., 2016). The possible cosmological application of these
cosmic chronometers has been extensively discussed in the literature
(Ratsimbazafy et al., 2017), such as the reconstruction of cosmological
distances (Wu et al., 2020; Zheng et al., 2021) and the spatial curvature of
the universe (Liu et al., 2020) in a cosmological model-independent way. In
particular, Nunes, Pan & Saridakis (2017) have recently studied the advantage
of utilizing the Hubble parameter data to impose improved constraints on the
viable and most used $f(R)$ gravity models, which encourages us to improve and
develop it further in this analysis.

Figure 2.— The likelihood distributions of $\Omega_{m}$, $n$ and $H_{0}$ in the $f_{1}(\mathcal{R})$ model, from the simulated GW data and its combination with the current EM observations. Figure 3.— The likelihood distributions of $\Omega_{m}$, $\alpha$ and $H_{0}$ in the $f_{2}(\mathcal{R})$ model, from the simulated GW data and its combination with the current EM observations. Data |  $\Omega_{m}$ | $H_{0}$ | $n$
---|---|---|---
EM+GW |  $0.285\pm 0.019(1\sigma)\pm{0.038}\sigma)$ |  $67.94\pm 0.36(1\sigma)\pm 0.70(2\sigma)$ |  $-0.055\pm 0.086(1\sigma)^{+0.17}_{-0.16}(2\sigma)$
EM |  $0.303^{+0.024}_{-0.028}(1\sigma)^{+0.054}_{-0.049}(2\sigma)$ |  $64.6\pm 2.7(1\sigma)^{+5.3}_{-5.2}(2\sigma)$ |  $0.04^{+0.11}_{-0.14}(1\sigma)^{+0.25}_{-0.24}(2\sigma)$
GW |  $0.303^{+0.073}_{-0.050}(1\sigma)^{+0.12}_{-0.13}(2\sigma)$ |  $67.91\pm 0.42(1\sigma)\pm 0.82(2\sigma)$ |  $0.01^{+0.27}_{-0.24}(1\sigma)^{+0.50}_{-0.53}(2\sigma)$
Data |  $\Omega_{m}$ | $H_{0}$ | $\alpha$
EM+GW |  $0.293^{+0.017}_{-0.015}(1\sigma)^{+0.035}_{-0.046}(2\sigma)$ |  $68.03^{+0.29}_{-0.33}(1\sigma)^{+0.69}_{-0.62}(2\sigma)$ |  $-0.06\pm 0.35(1\sigma)^{+0.70}_{-0.97}(2\sigma)$
EM |  $0.310^{+0.023}_{-0.029}(1\sigma)^{+0.056}_{-0.048}(2\sigma)$ |  $64.6\pm 2.6(1\sigma)^{+5.2}_{-5.1}(2\sigma)$ |  $0.35^{+0.31}_{-0.59}(1\sigma)^{+1.2}_{-0.93}(2\sigma)$
GW |  $0.347\pm 0.063(1\sigma)^{+0.13}_{-0.12}(2\sigma)$ |  $67.99\pm 0.42(1\sigma)^{+0.83}_{-0.81}(2\sigma)$ |  $0.7\pm 1.1(1\sigma)^{+2.1}_{-2.0}(2\sigma)$
Table 1The marginalized $1\sigma$ uncertainties of the parameters
$\Omega_{m}$, $H_{0}$ and $n$ ($\alpha$) for the two viable models in $f(R)$
theories.

Figure 4.— The evolution of $O_{m}(z)$ versus the redshift $z$ for
$f_{1}(\mathcal{R})$ and $f_{2}(\mathcal{R})$ model from the simulated EM+GW
data, with the $1\sigma$ uncertainty denoted by blue shades. The black dashed
lines represent the standard $\Lambda$CDM model.

Figure 5.— The evolution of $O_{m}(z)$ versus the redshift $z$ for
$f_{1}(\mathcal{R})$ and $f_{2}(\mathcal{R})$ model from the EM+GW data and EM
data, with the $1\sigma$ uncertainty denoted by blue shades(EM+GW), with the
$1\sigma$ uncertainty denoted by oronge shades(EM), The black dashed lines
represent the standard $\Lambda$CDM model.

##  4\. Constraints and model diagnostics

###  4.1. Observational constraints on the $f(R)$ models

In this section we will investigate observational bounds on two $f(R)$ models
from a statistical analysis involving three classes of data combinations,
i.e., GW, EM (SNe+BAO+$H(z)$), and GW+EM (GW+SNe+BAO+$H(z)$). The Monte Carlo
Markov Chain (MCMC) method is applied to perform the minimum likelihood method
of $\chi^{2}$ fit (Lewis & Bridle, 2002), with the fitting results shown in
Table 1. The marginalized probability distribution of each parameter
($\Omega_{m}$, $n$ and $H_{0}$ for $f_{1}(\mathcal{R})$; $\Omega_{m}$,
$\alpha$ and $H_{0}$ for $f_{2}(\mathcal{R})$) and the marginalized 2D
confidence contours are presented in Fig. 2-3.

Let us start from the first $f(R)$ theory within the Palatini formalism:
$f_{1}(\mathcal{R})=\mathcal{R}-\frac{\beta}{\mathcal{R}^{n}}$. For the
simulated GW data, the best-fit values of the parameters in the
$f_{1}(\mathcal{R})$ model are $\Omega_{m}=0.303^{+0.073}_{-0.050}$,
$H_{0}=67.91\pm 0.42$ km/s/Mpc, and $n=0.01^{+0.27}_{-0.24}$. In order to have
a good comparison, the fitting results derived from the EM observations
(SNe+BAO+$H(z)$) are also shown in Fig. 2, with the best-fit parameters of
$\Omega_{m}=0.303^{+0.024}_{-0.028}$, $H_{0}=64.6\pm 2.7$km/s/Mpc, and
$n=0.04^{+0.11}_{-0.14}$. One should note that, although GWs alone seem not
able to provide the constraints as good as the current EM observations
(SNe+BAO+$H(z)$), the degeneracy between different $f(\mathcal{R})$ parameters
obtained from the standard sirens is different from the statistical standard
probes in the EM window. Such tendency could also be clearly seen from the
comparison between the plots obtained with the joint GW+EM analysis, i.e.,
with the combined standard probe data sets in the GW and EM data sets, the
best-fit value for the parameters are $\Omega_{m}=0.285\pm 0.019$,
$H_{0}=67.94\pm 0.36$km/s/Mpc, and $-0.055\pm 0.086$. Therefore, the shift in
the best-fitted parameters (with remarkably reduced allowed region)
illustrates how the combination of the future GW observations can be used to
improve the model parameters in $f_{1}(\mathcal{R})$ cosmology. This has been
noted by the previous analysis studying the constraint ability of the
gravitational wave (GW) as the standard siren on the cosmological parameters
(Cai & Yang, 2017).

The best-fit values of the parameters along with their $1\sigma$ uncertainties
from three different data combinations are explicitly presented in Table 1.
Note also that in the modified $f_{1}(\mathcal{R})$ gravity, the best-fit
value for the matter density parameter and the Hubble parameter, i.e.,
$\Omega_{m}=0.285\pm 0.019$ and $H_{0}=67.94\pm 0.36$km/s/Mpc, are well
consistent with the current estimates in the framework of $\Lambda$CDM
cosmology given by recent Planck data release (Aghanim et al., 2018). More
importantly, the deviation parameter $n$ from three combinations of data
analysis (EM+GW, EM and GW) indicates that $\Lambda$CDM model ($n=0$) is still
included within 68.3% confidence level. However, we can see that the best
value of $n=-0.055$ seems to be slightly smaller than $0$, which suggests that
there still exists some possibility that $\Lambda$CDM may not be the best
cosmological model preferred by the current observations. Such conclusion is
marginally different from the previous results obtained in Santos et al.
(2008); Amarzguioui et al. (2006); Fay & Tavakol (2007), using the supernova
Gold and the SNLS data sets respectively.

Next we turn to study the ability of the standard siren to infer the
parameters in the $f_{2}(\mathcal{R})$ model:
$f_{2}(\mathcal{R})=\mathcal{R}+\alpha\ln{\mathcal{R}}-\beta$. The best fits
takes place at $\Omega_{m}=0.347\pm 0.063$, $H_{0}=67.99\pm 0.42$ km/s/Mpc and
$\alpha=0.7\pm 1.1$ with the simulated GW data. The constraining power of the
standard sirens in breaking degeneracy between model parameters is slightly
existence, as can be seen from the marginalized 1$\sigma$ and 2$\sigma$
contours of each parameter in Fig. 3. By fitting the $f_{2}(\mathcal{R})$
model to the GW+EM observations, we obtain
$\Omega_{m}=0.293^{+0.017}_{-0.015}$, $H_{0}=68.03^{+0.29}_{-0.33}$ km/s/Mpc
and $\alpha=-0.06\pm 0.35$. In Table 1 we also summarize the best-fit values
for the three combined data sets, respectively. In broad terms, the estimated
values of $\Omega_{m}$ and $H_{0}$ are basically consistent with each other at
the 68.3% confidence level. Moreover, compared with the case in the
$f_{1}(\mathcal{R})$ model, the largest difference happens on the constraint
of $\alpha$: the deviation from $\Lambda$CDM also tends to be slightly smaller
than $0$, although the concordance cosmological scenario is still included
within 68.3% confidence level. Therefore, there still exists some possibility
that $\Lambda$CDM may not the best cosmological model preferred by the current
and future high-redshift GW+EM observations. The constraining ability of
standard siren data can be quantified by comparing the results at 1$\sigma$ C.
L. based on GW+EM analysis, i.e., we obtain the error bar smaller than that
from EM data alone, when the $f_{2}(\mathcal{R})$ model is considered.

Figure 6.— The plot shows the behavior of the $f_{1}(\mathcal{R})$ model and
$f_{2}(\mathcal{R})$ model in the $(r,q)$ plane. The red and blue lines
represent the evolution of $(r,q)$ within the 1$\sigma$ error range.

Figure 7.— The same as Fig. 6, but in the $(r,s)$ plane.

###  4.2. Model diagnostics

Now two diagnostics analysis (the $Om(z)$ and the statefinder diagnostics)
will be performed in this section, acting as alternative probes capable of
discriminating between $f(R)$ cosmology and $\Lambda$CDM model.

On the one hand, the expansion rates at different redshifts, or the Hubble
parameters $H(z)$, has opened a new chapter of using the so-called $Om(z)$
diagnostics to discriminate different cosmological models as well as
$\Lambda$CDM model (Sahni et al., 2008; Ding et al., 2015; Qi et al., 2018;
Zheng et al., 2016, 2019). In this method, only a single time derivative of
$a(t)$ is used and a new diagnostics is defined as

$Om(z)=\frac{E^{2}(z)-1}{(1+z)^{3}-1}$ (29)

where $E(z)=\frac{H(z)}{H_{0}}$ represents the dimensionless expansion rate.
It is obvious that in the flat $\Lambda$CDM model, the $Om(z)$ evaluated at
any redshift is always equal to the present mass density parameter
$\Omega_{m}$. Therefore, such diagnostic can be used as a cosmological probe
to directly illustrate the difference between $\Lambda$CDM and other
cosmological models.

Based on the best-fit model parameters and their $1\sigma$ uncertainties
derived from the simulated EM+GW data, we have different values of $Om(z)$ as
the function of redshift for the $\Lambda$CDM and two $f(R)$ models, which are
presented in Fig. 4. It is obvious that the $O_{m}(z)$ curve of the
$f_{1}(\mathcal{R})$ model will coincide with $\Lambda$CDM at high redshifts
($z>4$), which indicates that the $O_{m}(z)$ for the $f_{1}(\mathcal{R})$
model cannot be distinguished $\Lambda$CDM in the early universe. However, the
$f_{1}(\mathcal{R})$ cosmology begins to deviate from $\Lambda$CDM in the
redshift range of $z<4$, while the deviation between $f_{2}(\mathcal{R})$
model and $\Lambda$CDM takes place at the redshift of $z>1.5$. In Fig. 5 we
also compare the $O_{m}(z)$ function based on the fits from EM+GW and EM data.
One may clearly see that in the framework of two $f(R)$ models, the
combination of future naturally improved GW data most recent EM observations
will reveal the consistency or acknowledge the tension between the
$\Lambda$CDM model and modified gravity theories.

On the other hand, in order to more effectively differentiate and compare
between different cosmological models, we turn to a discussion on the so-
called statefinder diagnostic. For a specific cosmological model, the Hubble
parameter $H$ and the deceleration parameter $q$ are respectively related as

$q(z)=-\frac{\ddot{a}}{aH^{2}}=\frac{E^{\prime}(z)}{E(z)}(1+z)-1$ (30)

In the analysis of statefinde diagnostic (Sahni et al., 2008), the new
dimensionless parameters ${r,s}$ are constructed from the scale factor $a(t)$
and its time derivatives similar to the geometrical parameters $H(z)$ and
$q(z)$:

$\displaystyle r(z)$ $\displaystyle=$ $\displaystyle
q(z)(1+2q(z))+q^{\prime}(z)(1+z)$ (31) $\displaystyle s(z)$ $\displaystyle=$
$\displaystyle\frac{1}{3}(r(z)-1)/(q(z)-1/2)$

where $q^{\prime}(z)\equiv\frac{dq(z)}{dz}$. Any obtained model can be
compared with the standard cosmological constant model, from which one may see
how a model approaches or deviates from $\Lambda$CDM in the $s-r$ plane and
$q-r$ plane. The evolution of starefinder pair $(s,r)$ in the two $f(R)$
models is shown in Fig. 7. On the one hand, it is obvious that the statefinder
pair exhibits the similar behaviors in the framework of two $f(R)$ models,
i.e., originating from the same point of $(r,q)=(1,0.5)$, evolving along
different trajectory, and finally converging on the same point of ($1,-1$). On
the other hand, compared with the deceleration parameter $q$, the parameter
$r$ exhibiting apparent fluctuations is more effective in discriminating
different cosmological models. The trajectories of the two $f(R)$ models in
the $r-q$ plane are also presented in Fig. 6, in which the point $(r,s)=(1,0)$
corresponding to the $\Lambda$CDM universe is also added for comparison. It is
noteworthy that, although the corresponding value for the two $f(R)$ models
significantly deviates from $\Lambda$CDM at the present epoch, both of two
modified gravity models discussed in this analysis is resembling to the
standard cosmological model in the future and ultimately freezing to it. Such
a tendency is more obvious when the $1\sigma$ uncertainties of the $f(R)$
parameters are taken into consideration.

##  5\. Conclusions

In this paper we have investigated the constraint ability of GW events on the
modified gravity models using the simulated data from the third-generation
gravitational wave detector, the Einstein Telescope. We focus on two $f(R)$
models which introduces a perturbation of the Ricci scalar $R$ in the
Einstein-Hilbert action, which can be considered as an important candidate
models alternative to dark energy. More specifically, in the framework of two
viable models in $f(R)$ theories within the Palatini formalism
($f_{1}(\mathcal{R})=\mathcal{R}-\frac{\beta}{\mathcal{R}^{n}}$ and
$f_{2}(\mathcal{R})=\mathcal{R}+\alpha\ln{\mathcal{R}}-\beta$), our results
show that the sensitivity achieved by the ET detector or a similar third
generation interferometer is enough to improve the current estimates on the
free parameter within $f(R)$ gravity. In practise more GW events may be needed
in order to achieve the accuracy of EM data. More importantly, the standard
sirens GWs could effectively eliminate the degeneracies among parameters in
the two modified gravity models, in joint analysis with the latest
electromagnetic (EM) observational data (including the recently released
Pantheon type Ia supernovae sample, Hubble parameter data from cosmic
chronometers, and baryon acoustic oscillation distance measurements).

In addition, we thoroughly investigate the nature of geometrical dark energy
in the modified gravity theories with the assistance of $Om(z)$ and
statefinder diagnostic analysis. The present analysis makes it clear that the
simplest cosmological constant model is still the most preferred by the
current data. However, the combination of future naturally improved GW data
most recent EM observations will reveal the consistency or acknowledge the
tension between the $\Lambda$CDM model and modified gravity theories.
Therefore, in the case of detections with the ET, our forecast analysis
indicates that it is possible to provide an independent and complementary
alternative to current experiments, once that these events begin to be
detected. One should note that, dedicated observations of the sky position of
each host galaxy (i.e., with negligible measurement error) would be necessary,
in order to obtain significant improvements in the distance, and hence in the
constraints on modified gravity theories. Although the GW distance posterior
changes slowly over the sky and therefore is not sensitive to the precise
location of the counterpart, obtaining such measurements for a sample of
different types of GW events would require substantial follow-up efforts
(Zhang et al., 2020).

## Acknowledgments

We thank Prof. Z.-H. Zhu for useful discussions. This work was supported by
National Key R&D Program of China No. 2017YFA0402600; the National Natural
Science Foundation of China under Grants Nos. 12021003, 11690023, and
11633001; Chongqing Municipal Science and Technology Commission Fund
(cstc2018jcyjAX0192); Beijing Talents Fund of Organization Department of
Beijing Municipal Committee of the CPC; the Strategic Priority Research
Program of the Chinese Academy of Sciences, Grant No. XDB23000000; the
Interdiscipline Research Funds of Beijing Normal University; and the Opening
Project of Key Laboratory of Computational Astrophysics, National Astronomical
Observatories, Chinese Academy of Sciences. J. Wang was supported by the
National Natural Science Foundation of China under Grant No.11647079, Science
and Technology Department of Yunnan Province - Yunnan University Joint Funding
(2019FY003005), Donglu Youth Teacher Plan of Yunnan University and the Key
Laboratory of Astroparticle Physics of Yunnan Province.

## References

  * Abadie (2010) Abadie, J., et al. 2010, Nucl. Instrum. Methods Phys. Res., 624, 223
  * Abbott et al. (2016) Abbott, B. P., et al. 2016, PRL, 116, 061102
  * Abbott et al. (2017) Abbott, B. P., et al. 2017, Nature, 551, 85
  * Ade et al. (2016) Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * Aghanim et al. (2018) Aghanim, N., Akrami, Y., Ashdown, M., et al. 2018, arXiv:1807.06209
  * Alam et al. (2017) Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617
  * Amendola et al. (2006a) Amendola, L., et al. 2006a, arXiv:0603703
  * Amendola et al. (2006b) Amendola, L., et al. 2006b, arXiv:0612180
  * Amarzguioui et al. (2006) Amarzguioui, M., et al. 2006, A&A, 454, 707
  * Aoki & Shimada (2020) Aoki, K. & Shimada, K. 2020, PRD, 98, 044038
  * Baghram & Rahvar (2020) Baghram, S. & Rahvar, S. 2020, PRD, 80, 124049
  * Bekov et al. (2020) Bekov, S., Myrzakulov, K.,Myrzakulov, R. & Gomez, D. S. C. 2020, arXiv:2010.12360v2 [gr-qc]
  * Bengochea & Ferraro (2009) Bengochea, G. R. & Ferraro, R. 2009, PRD, 79, 124019
  * Beutler et al. (2011) Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017
  * Bohmer et al. (2013) Bohmer, C. G., Lobo, F. S. N. & Tamanini, N. 2013, PRD, 88, 104019
  * Borowiec et al. (2012) Borowiec, A., Kamionka, M., Kurek, A. & Szydlowski, M. 2012, JCAP, 2, 027
  * Cai et al. (2015) Cai, Y.-F., et al. 2015, arXiv:1511.07586
  * Cai & Yang (2017) Cai, R.-G. & Yang, T. 2017, PRD, 95, 044024
  * Cao et al. (2011) Cao, S., Zhu, Z.-H., & Zhao, R. 2011, PRD, 84, 023005
  * Cao & Zhu (2012) Cao, S., & Zhu, Z.-H. 2012, A&A, 538, A43
  * Cao, Covone & Zhu (2012) Cao, S., Covone, G., & Zhu, Z.-H. 2012, ApJ, 755, 31
  * Cao et al. (2012) Cao, S., Pan, Y., Biesiada, M., Godlowski, W., & Zhu Z.-H. 2012, JCAP, 03, 016
  * Cao & Zhu (2014) Cao, S., & Zhu, Z.-H. 2014, PRD, 90, 083006
  * Cao et al. (2015a) Cao, S., et al. 2015a, AJ, 149, 3
  * Cao et al. (2015b) Cao, S., Biesiada, M., Gavazzi, R., Piórkowska, A., & Zhu, Z.-H. 2015b, ApJ, 806, 185
  * Cao et al. (2017a) Cao, S., et al. 2017, JCAP, 02, 012
  * Cao et al. (2017b) Cao, S., et al. 2017, A&A, 606, A15
  * Cao et al. (2019a) Cao, S., Qi, J. Z., Biesiada, M., et al. 2019a, PDU, 24, 100274
  * Cao et al. (2019b) Cao, S., et al. 2019b, NatSR, 9, 11608
  * Cao et al. (2020) Cao, S., Qi, J. Z., Biesiada, M., Liu, T. H., Zhu, Z.-H. 2020, ApJL, 888, L25
  * Capozziello et al. (2003) Capozziello, S., et al. 2003, PRD, 12, 1969
  * Capozziello et al. (2004) Capozziello, S., et al. 2004, PRD, 70, 043528
  * Capozziello et al. (2013) Capozziello, S., Harko, T., Koivisto, T. S., Lobo, F. S. N. & Olmo, G. J. 2013, JCAP, 4, 011
  * Capozziello et al. (2005) Capozziello, V. F. S., et al. 2005, PRD, 71, 043503
  * Capozziello & Francaviglia (2008) Capozziello, S. & Francaviglia, M. 2008, Gen. Relativ. Gravit, 40, 357
  * Carroll et al. (2005) Carroll, S. M., et al. 2005, PRD, 71, 063513
  * Chen (2019) Chen, L., Huang, Q.-G., & Wang, K. 2019, JCAP, 02, 028
  * Cutler & Holz (2009) Cutler, C. & Holz, D. E. 2009, PRD, 80, 104009
  * D’Agostino & Nunes (2019) D’Agostino, R. & Nunes, R. C. 2019, PRD, 100, 044041
  * D’Agostino & Nunes (2020) D’Agostino, R. & Nunes, R. C. 2020, PRD, 101, 103505
  * Dalal et al. (2006) Dalal, N., et al. 2006, PRD, 74, 063006
  * De Felice & Tsujikawa (2010) De Felice, A. & Tsujikawa, S. 2010, Living Rev. Relativ, 13, 3
  * de la Cruz-Dombriz et al. (2016) de la Cruz-Dombriz, A., et al. 2016, PRD, 93, 084016
  * Des Bourboux et al. (2017) Des Bourboux, H. D. M., Le Goffff, J.-M., Blomqvist, M., et al. 2017, A&A, 608, A130
  * Ding et al. (2015) Ding, X., Biesiada, M., Cao, S., Li, Z., & Zhu, Z.-H. 2015, ApJL, 803, L22
  * Dolgov & Kawasaki (2003) Dolgov, A. D. & Kawasaki, M. 2003, PLB, 573, 1
  * Dos Santos et al. (2016) Dos Santos, M. V., et al. 2016, A&A, 587, A132
  * Fay & Tavakol (2007) Fay, S. & Tavakol, R. 2007, PRD, 75, 063509
  * Flanagan (2004a) Flanagan, E. E. 2004a, PRL, 92, 071101
  * Flanagan (2004b) Flanagan, E. E. 2004b, CQG, 21, 417
  * Gu et al. (2018) Gu, B. M., Liu, Y. X. & Zhong, Y. 2018, PRD, 98, 024027
  * Jimenez & Loeb (2002) Jimenez, R. & Loeb, A. 2002, ApJ, 573, 37
  * Jimenez et al. (2003) Jimenez, R., Verde, L., Treu, T. & Stern, D. 2003, ApJ, 593, 622
  * Harko et al. (20212) Harko, T., Koivisto, T. S., Lobo, F. S. N. & Olmo, G. J. 2012, PRD, 85, 084016
  * Holz & Hughes (2005) Holz, D. E. & Hughes, S. A. 2005, ApJ, 629, 15
  * Leanizbarrutia et al. (2017) Leanizbarrutia, I., Lobo, F. S. N. & Saez-Gomez, D. 2017, PRD, 95, 084046
  * Lewis & Bridle (2002) Lewis, A. & Bridle, S. 2002, PRD, 66, 103511
  * Li et al. (2007) Li, B., Chan, K. C., & Chu, M.-C. 2007, PRD, 76, 024002
  * Li (2015) Li, T. G. F., Extracting Physics from Gravitational Waves (Springer Theses, New York (2015)
  * Li et al. (2017) Li, X. L., et al. 2017, EPJC, 77, 677
  * Liao (2019) Liao, K. 2019, ApJ, 885,70
  * Liu et al. (2019) Liu, T. H., et al. 2019, ApJ, 886, 94
  * Liu et al. (2020) Liu, T. H., et al. 2020, ApJ, 899, 71
  * Liu et al. (2020) Liu, Y. T., et al. 2020, ApJ, 901, 129
  * Ma et al. (2019) Ma, Y. B., et al. 2019, EPJC, 79, 121
  * Moresco (2015) Moresco, M. 2015, MNRAS, 450, L16
  * Moresco et al. (2016) Moresco, M., et al. 2016, JCAP, 5, 014
  * Moresco et al. (2012) Moresco, M., Verde, L., Pozzetti, L., Jimenez, R. & Cimatti, A. 2012, JCAP, 7, 053
  * Nojiri & Odintsov (2003) Nojiri, S. & Odintsov, S. D. 2003, PRD, 68, 123512
  * Nojiri & Odintsov (2005) Nojiri, S. & Odintsov, S. D. 2005, PLB, 631, 1
  * Nunes, Pan & Saridakis (2017) Nunes, R. C., Pan, S., & Saridakis, E. N. 2017, JCAP, 01, 005
  * Pan et al. (2015) Pan, Y., et al. 2015, ApJ, 808, 78
  * Qi et al. (2016) Qi, J.-Z., et al. 2016, 16, 002
  * Qi et al. (2018) Qi, J.-Z., et al. 2018, RAA, 18, 66
  * Qi et al. (2019a) Qi, J.-Z., et al. 2019a, MNRAS, 483, 1
  * Qi et al. (2019b) Qi, J.-Z., et al. 2019b, PRD, 99, 6
  * Ratsimbazafy et al. (2017) Ratsimbazafy, A. L., et al. 2017, MNRAS, 467, 3239
  * Risaliti & Lusso (2018) Risaliti, G., & Lusso, E. 2019, Nature Astronomy, 3, 272
  * Ross et al. (2015) Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS, 449, 835
  * Rosa et al. (2020) Rosa, J. L., Carloni, S. & Lemos, J. P. S. 2020, PRD, 101, 104056
  * Rosa et al. (2017) Rosa, J. L., Carloni, S., Lemos, J. P. S. & Lobo, F. S. N.. 2017, PRD, 95, 124035
  * Sahni et al. (2008) Sahni, V., Shafieloo, A., & Starobinsky, A. A. 2008, PRD, 78, 103502
  * Santos et al. (2008) Santos, J., et al. 2008, PLB, 669, 14
  * Sathyaprakash et al. (2010) Sathyaprakash, B. S., et al. 2010, CQG, 27, 215006
  * Schneider et al. (2001) Schneider, R., et al. 2001, MNRAS, 324, 797
  * Schutz (1986) Schutz, B. F. 1986, Nature, 323, 310
  * Scolnic et al. (2018) Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101
  * Seokcheon (2007) Seokcheon L. 2007, arXiv:0710.2395
  * Simon et al. (2005) Simon, J., Verde, L. & Jimenez, R. 2005, PRD, 71, 123001
  * Song et al. (2007) Song, Y.-S., Peiris, H., & Hu, W. 2007, PRD, 76, 063517
  * Sotiriou & Faraon (2010) Sotiriou, T. P. & Faraoni, V. 2010, Reviews of Modern Physics, 82, 451
  * Stern et al. (2010) Stern, D., Jimenez, R., Verde, L., Kamionkowski, M. & Stanford, S. A. 2010, JCAP, 2, 008
  * Taylor et al. (2012) Taylor, S. R., et al. 2012, PRD, 85, 023535
  * Voivodic et al. (2017) Voivodic, R., et al. 2017, PRD, 95, 024018
  * Vollick (2003) Vollick, D. N. 2003, PRD, 68, 063510
  * Wang et al. (2010) Wang, J., et al. 2010, PLB, 689, 133
  * Wang et al. (2018) Wang, J., et al. 2018, PDU, 19, 60
  * Wang & Wang (2013) Wang, J. & Wang, H. 2013, PLB, 724, 5
  * Wei (2019) Wei, J.-J. 2019, ApJ, 876, 1
  * Will (1981) Will, C. M. 1981, Theory and Experiment in Gravitational Physics (Cambridge University Press, New York)
  * Wu et al. (2020) Wu, Y., et al. 2020, ApJ, 888, 113
  * Xu et al. (2018) Xu, T., et al. 2018, JCAP, 06, 042
  * Yang & Chen (2009) Yang, X. J. & Chen, D. M. 2009, MNRAS, 394, 3
  * Zhang et al. (2014) Zhang, C. et al. 2014, RAA, 14, 1221
  * Zhang et al. (2020) Zhang, S. X., et al. 2020, IJMPD, 2020, 2050105
  * Zhao et al. (2011) Zhao, W., et al. 2011, PRD, 83, 023005
  * Zhao et al. (2018) Zhao, G.-B., Wang, Y., Saito, S., et al. 2018, MNRAS, 482, 3
  * Zheng et al. (2016) Zheng, X.-G., Ding, X., Biesiada, M., Cao, S., & Zhu, Z.-H. 2016, ApJ, 825, 17
  * Zheng et al. (2019) Zheng, X.-G., et al. 2019, EPJC, 79, 637
  * Zheng et al. (2020) Zheng, X.-G., et al. 2020, ApJ, 892, 103
  * Zheng et al. (2021) Zheng, X.-G., et al. 2021, SCMPA, in press [arXiv:2012.14607]