# Nowruz at SemEval-2022 Task 7: Tackling Cloze Tests with Transformers and
Ordinal Regression

Mohammadmahdi Nouriborji
NODET, Iran11footnotemark: 1
<EMAIL_ADDRESS>
Omid Rohanian
Department of Engineering Science
University of Oxford
Oxford, UK
<EMAIL_ADDRESS>
&David Clifton
Department of Engineering Science
University of Oxford
Oxford, UK
<EMAIL_ADDRESS>

###### Abstract

This paper outlines the system using which team Nowruz participated in SemEval
2022 Task 7 “Identifying Plausible Clarifications of Implicit and
Underspecified Phrases” for both subtasks A and B (Roth et al., 2022). Using a
pre-trained transformer as a backbone, the model targeted the task of multi-
task classification and ranking in the context of finding the best fillers for
a cloze task related to instructional texts on the website Wikihow.

The system employed a combination of two ordinal regression components to
tackle this task in a multi-task learning scenario. According to the official
leaderboard of the shared task, this system was ranked 5th in the ranking and
7th in the classification subtasks out of 21 participating teams. With
additional experiments, the models have since been further optimised. The code
used in the experiments is going to be freely available at
https://github.com/mohammadmahdinoori/Nowruz-at-SemEval-2022-Task-7.

Nowruz at SemEval-2022 Task 7: Tackling Cloze Tests with Transformers and
Ordinal Regression

Mohammadmahdi Nouriborji NODET, Iran††footnotemark: †
<EMAIL_ADDRESS>

Omid Rohanian Department of Engineering Science University of Oxford Oxford,
UK<EMAIL_ADDRESS>David Clifton Department of Engineering Science
University of Oxford Oxford, UK<EMAIL_ADDRESS>

00footnotetext: National Organization For Development of Exceptional Talents

##  1 Introduction

Oxford dictionary defines cloze test as “a test of readability or
comprehension in which a person is required to supply words which have been
deliberately omitted from a passage” (Oxford University Press, 2022). In the
context of NLP, a cloze format task is one in which the context is one or more
sentences with masked spans and the model is expected to predict a suitable
filler for each span. Cloze-format datasets have become popular in NLP
recently as they are relatively easy to create automatically and provide high
quality resources for model training (Rogers et al., 2021).

SemEval 2022 task 7 is framed as a cloze task in which the goal is to rank or
classify fillers within a given context based on their suitability. The texts
are taken from actual articles on an instructional website and the masked
spans are placed at locations of edits made by users. We participated in this
shared task in both ranking and classification parts and developed a
transformer-based model that utilises both classification and regression
components at the top layer. The code and the data used in these experiments
are publicly available.

##  2 Related Work

Cloze tasks have been a subject of interest in Natural Language Understanding
(NLU) in recent years, especially within the context of reading comprehension,
story understanding, and summarisation (Deutsch and Roth, 2019; Sharma et al.,
2018; Mostafazadeh et al., 2016). There is evidence that cloze tasks can be
used to effectively pretrain or finetune language models in order to perform
few shot learning (Schick and Schütze, 2021; Liu et al., 2021).

WikiHow is a community-edited open domain repository that hosts how-to
articles on a variety of different subjects. It is possible to track edits
made by users and compare different versions. There have been some recent
computational works exploring this resource, including modelling of revision
requirements (Bhat et al., 2020), and the effect of edits on fluency (Anthonio
and Roth, 2020) and vagueness (Debnath and Roth, 2021).

###  2.1 Masked Language Modelling

Masked Language Modeling (MLM) is a pre-training task which is widely used in
transformer-based models. This task forms a self-supervised cloze test by
randomly removing some of the tokens of the sentence which will be further
filled by the model during training. Prominent transformers including BERT
(Devlin et al., 2019) , RoBERTa (Liu et al., 2019), DeBERTa (He et al., 2020)
and T5 (Raffel et al., 2019) are trained using MLM as an auxiliary objective.
Given the successes of transformers in most areas of machine learning and NLP,
it is standard practice to fine-tune them for various down-stream tasks.

###  2.2 Ordinal Regression

Ordinal Regression (also known as Ordinal Classification) is a type of
classification in which labels have order with respect to each other. Examples
of ordinal regression tasks include age estimation (Niu et al., 2016),
assessment of damage (Ci et al., 2019), and monocular depth estimation Fu et
al. (2018). In ordinal regression, performance of the model is sensitive to
the order of the predictions with regards to the labels. For instance, in the
age estimation task, the error of the model should be higher when it
incorrectly predicts the age of $30$ as $10$, as opposed to when it predicts
the age of $30$ as $20$. Ordinal Regression is commonly done by breaking the
multi-class classification task into several binary classification subtasks
within a multi-task learning scenario. The output of these binary
classification subtasks should be rank-consistent to achieve good performance.

##  3 System Description

Figure 1 shows the overall architecture of our model. A pre-trained
transformer sits at the base of the architecture. Different variations of this
transformer are tried and reported in Sec 5. These include BERT (Devlin et
al., 2019), RoBERTa (Liu et al., 2019), DeBERTa (v1 & 3) (He et al., 2020,
2021), and T5111Since T5 is originally an encoder-decoder language model, we
only use its pre-trained encoder for our experiments Raffel et al. (2019) in
different configurations.

At the bottom of the network, there are two components, one for regression and
the other for classification. At training time, both of these are trained in
tandem using ordinal regression and a combined loss. The loss function ensures
that the ranks of the labels are kept consistent across these task (Sec. 3.6).

Figure 1: Overall architecture of the model. Note how a pooled representation
is produced using the last token of the filler and the CLS token

###  3.1 Representing the Filler

To classify each filler, it is first placed in the location of the blank in
the text to form a full context. The pre-trained transformer is equipped with
a word-piece tokeniser that breaks down the filler to subword units and also
uses a [CLS] token to represent the entire context 222T5 does not use CLS for
context representation. In that case we just classify the last word-piece with
no additional concatenation. To build a pooled representation for the filler,
the contextualised representations for the last word piece and the entire
context (i.e. representation of [CLS]) are concatenated and passed on to the
next layer.

###  3.2 Multi-task learning with Ordinal Regression

Once a combined contextualised representation is obtained from concatenation
of [CLS] and the filler, we address both subtasks of the shared task in a
single multi-task learning architecture. The objective of this model is to
predict both the class and the suitability score for each provided filler. The
concatenated representation is fed to a fully connected feed-forward layer
followed by a GELU activation function (Hendrycks and Gimpel, 2016). This
layer projects the representation to a lower dimensional space. Subsequently,
the output is passed on to two separate classification and regression heads.

###  3.3 Decomposing the Problem into Binary Classification tasks

As mentioned in Sec. 3, our model uses a special loss function named coral to
perform ordinal regression. This loss function and the training procedure it
requires are explained in Cao et al. (2020). Coral layers break a
$\mathrm{K}$-class classification problem into $\mathrm{K}-1$ binary
classification tasks as part of a multi-task learning scenario. For a coral
layer with $\mathrm{K}-1$ units 333A $\mathrm{T}$-unit coral layer, is
comprised of $\mathrm{T}$ binary classification units which share the same
weights but have different biases., the value of loss is constructed from the
sum of $\mathrm{K}-1$ separate binary cross entropy losses belonging to each
unit.

To train a coral layer, it is necessary to first transform the original labels
to sets of binary labels. This is the step where the notion of order is
introduced into the model. Given $\mathrm{K}$ classes, we convert each label
to a collection of binary labels as follows:

$\displaystyle f(y,k)=\begin{cases}1\quad&\text{if}\,k<y\\\
0\quad&\text{if}\,k\geq y\\\ \end{cases}$ (1) $\displaystyle
Y^{(i)}_{ordinal}=\\{f(y^{(i)},k)|k\in\mathbb{W}\land k<\mathrm{K}-1\\}$

where $y^{(i)}$ is the original label for the $i$th training instance and
$0\leq y^{(i)}\leq\mathrm{K}-1$.

###  3.4 Classification and Regression Heads

At the top end of the architecture there are two components for classification
and regression. The classification head is a coral layer with two units which
is used to address the classification subtask where we have three labels,
namely, Implausible, Neutral, and Plausible. Note that the labels in the
classification task have inherent order and for the model to be trained
effectively, it is important that the training objective penalises
misclassifications based on this underlying assumption. For instance, the
error of the model should be more when it predicts an Implausible sample as
Plausible compared to when it predicts an Implausible sample as Neutral. This
is the motivation to use ordinal regression losses like coral.

The other head is assigned to the regression subtask and rates each filler for
suitability. We converted continuous scores in the $1-5$ range to labels with
discrete values in the $0-4$ range by either rounding or flooring the scores.
For example, the original scores of $\\{1.333,1.75,2.5,3.75,4.25\\}$ are
mapped to $\\{0,1,2,3,3\\}$ by rounding and $\\{0,0,1,2,3\\}$ by flooring. We
can frame the regression task as ordinal classification with five labels by
binning the values. A $4$-unit coral layer is used to perform this
classification. These heads are jointly trained using a combined loss.

###  3.5 Constructing Labels and Ranks from Heads

Since coral layers differ from normal dense layers, their output can not be
directly converted to labels. Furthermore, the regression subtask is also
framed as an additional classification task in our methodology. However, since
the purpose is to report continuous scores rather than discrete labels, a
unique conversion is necessary for the output of the regression head.

For the classification head, the goal is to output discrete labels in range of
$0-2$. Given the output of the classification head for one sample $\hat{C}$,
the conversion to the label is defined as follows:

$\displaystyle f(y)=\begin{cases}1\quad&\text{if}\,y>0.5\\\
0\quad&\text{if}\,y\leq 0.5\end{cases}$ (2) $\displaystyle
c=\sum_{k=1}^{2}f(\sigma(\hat{C}_{k}))$

where $c$ is the final label and $\sigma$ is the sigmoid function.

For the regression head, the goal is to output continuous scores in range of
$1-5$. Given the output of the regression head for one sample $\hat{R}$, the
conversion to the continuous scores in range of $1-5$ is defined as follows:

$\displaystyle r=(\sum_{k=1}^{4}\sigma(\hat{R}_{k}))+1$ (3)

###  3.6 Computation of Loss

Since we have three labels for the classification task, ordinal labels would
be sets of binary values with the length of two. Given $C^{(i)}$ as the
converted classification label for the $i_{th}$ training sample and the
$\hat{C}^{(i)}$ as the output of the classification head for the $i_{th}$
training sample, we can define the classification loss as:

$\displaystyle l^{(i)}_{c}=\sum_{k=1}^{2}-C^{(i)}_{k}\log(\hat{C}^{(i)}_{k})$
(4) $\displaystyle-(1-C^{(i)}_{k})\log(1-\hat{C}^{(i)}_{k})$

For the regression task there are $5$ different labels, and accordingly,
ordinal labels would be sets of binary values with the length of four. Given
$R^{(i)}$ as the converted regression label for the $i_{th}$ training sample
and the $\hat{R}^{(i)}$ as the output of the regression head for the $i_{th}$
training sample we can define the regression loss in a similar way:

$\displaystyle l^{(i)}_{r}=\sum_{k=1}^{4}-R^{(i)}_{k}\log(\hat{R}^{(i)}_{k})$
(5) $\displaystyle-(1-R^{(i)}_{k})\log(1-\hat{R}^{(i)}_{k})$

Since we are aiming to perform both classification and regression, a joint
loss is needed to combine losses from the two heads. For a given training
batch of size $n$:

$\displaystyle
L_{total}=\frac{1}{n}\sum_{i=1}^{n}\lambda_{c}l^{(i)}_{c}+\lambda_{r}l^{(i)}_{r}$
(6)

where $\lambda_{c}$ is the weight associated with the classification loss and
the $\lambda_{r}$ is associated with regression.

##  4 Dataset and Experimental Setup

The shared task is based on Anthonio et al. (2020), where instructional texts
from the website Wikihow 444A wiki-style online collection of how-to articles
accessible at https://www.wikihow.com are used to create training instances
for a cloze task. The deletions in the dataset are based on actual edits made
by online users and they are assumed to represent certain types of
clarifications to make a point more clear or disambiguate a sentence. Based on
the actual edits, alternative fillers are automatically extracted and added to
build a number of possible fillers. These texts along with the fillers were
later annotated by humans and given plausibility scores from $1$ to $5$. There
is also a separate type of annotation available in which there are $3$ labels
with discrete plausibility values. The shared task was organised in two
separate tracks of classification and regression, depending on what kind of
annotation was used for modelling the task. Table 2 & 3 show the basic
statistics of the dataset in the shared task.

The dataset for this task consists of six features, named Resolved Pattern,
Article title, Section header, Previous context, Sentence, and Follow-up
context along with five different fillers for each sample. Resolved Pattern is
one of the four following categories: IMPLICIT REFERENCE, ADDED COMPOUND,
METONYMIC REFERENCE, and FUSED HEAD which indicates the relationship of the
fillers with the context. Article title is the name of the article from which
the paragraph is selected. Section header is the section from which the
article is selected. Previous context is a few sentences before the sentence
that contains the filler. Sentence is the sentence that contains the filler.

We obtained our best results when we used a custom formatting using which we
can feed all the features to the model as textual input. Table 1 is an example
of how we represent each training instance. Note how the ‘Text’ feature is
constructed by concatenating previous and follow-up contexts with the target
sentence.

Example input while using all features
---
Resolved pattern: ADDED COMPOUND
Section header: Following a Basic Routine
Article title: How to Get Rid of Peeling Skin
Text: (…) 6. Never tear away loose skin. (…) 7. Protect your skin from
sunlight. Exposure to direct sunlight can weaken your skin further and
complicate the [Filler] problem. This is true regardless of whether your skin
is peeling due to a sunburn or due to dryness.
Table 1:  An example of how each input is formatted for training and inference. Each identifier followed by ":" represents a feature of the dataset represented in textual format. Variation | implausible | neutral | plausible | total
---|---|---|---|---
train | 5474 | 7162 | 7339 | 19975
| (27.40%) | (35.85%) | (36.75%) | (100%)
dev | 982 | 602 | 916 | 2500
| (39.28%) | (24.08%) | (36.64%) | (100%)
Table 2:  Statistics of the data for subtask A Variation | 1 | 2 | 3 | 4 | 5 | total
---|---|---|---|---|---|---
train F | 2254 | 4123 | 6259 | 5321 | 2018 | 19975
train R | 1053 | 4421 | 4034 | 8441 | 2026 | 19975
dev F | 645 | 458 | 481 | 596 | 320 | 2500
dev R | 386 | 596 | 376 | 639 | 503 | 2500
Table 3:  Statistics of the data for subtask B. R and F represent rounded and
floored scores, respectively

##  5 Results and Discussion

The official test results for both subtasks are presented in Table 5. These
are the best results that we have obtained on the test set prior to the end of
the competition. On the official leaderboard of the shared task, our system
was ranked $5$th and $7$th for the ranking and classification tasks,
respectively.

We have since performed extensive analyses on the effect of hyperparameters on
different variations of our models, and since we do not have access to the
true labels in the test set, we report our best results on the dev set. As can
be seen in Table 4, we have produced our best results using DeBERTa-V3large.
We have noticed two important factors that influence the final performance of
the models. The first factor is batch size. Our best results were obtained on
lower batch sizes of $8$ and $16$. The choice of rounded or floored numbers
for subtask B is also a significant factor. The reason for this is that the
distribution of labels changes depending on the normalisation method used.

Backbone Transformer |  Model |  Batch Size |  Scores |  Accuracy (subtask A) |  Spearman’s rank (subtask B)
---|---|---|---|---|---
|  Bertbase |  8 |  R |  57.80% |  0.6020
|  Bertbase |  8 |  F |  57.48% |  0.6038
|  Bertbase |  16 |  R |  56.68% |  0.5950
BERT |  Bertbase |  16 |  F |  57.92% |  0.5960
|  Bertlarge |  8 |  R |  57.60% |  0.6323
|  Bertlarge |  8 |  F |  59.12% |  0.6341
|  RoBERTabase |  8 |  R |  59.16% |  0.6610
|  RoBERTabase |  8 |  F |  58.68% |  0.6582
|  RoBERTabase |  16 |  R |  57.32% |  0.6743
RoBERTa |  RoBERTabase |  16 |  F |  59.44% |  0.6695
|  RoBERTalarge |  8 |  R |  58.68% |  0.6799
|  RoBERTalarge |  8 |  F |  60.20% |  0.6928
|  DeBERTa-V3base |  8 |  R |  61.72% |  0.7194
|  DeBERTa-V3base |  8 |  F |  62.64% |  0.7260
|  DeBERTa-V3base |  16 |  R |  63.12% |  0.7287
DeBERTa-V3 |  DeBERTa-V3base |  16 |  F |  63.84% |  0.7326
|  DeBERTa-V3large |  8 |  R |  64.12% |  0.7406
|  DeBERTa-V3large |  8 |  F |  64.12% |  0.7411
|  T5-Encoderbase |  8 |  R |  58.88% |  0.6598
|  T5-Encoderbase |  8 |  F |  59.48% |  0.6573
|  T5-Encoderbase |  16 |  R |  59.24% |  0.6526
T5 |  T5-Encoderbase |  16 |  F |  59.08% |  0.6568
|  T5large |  8 |  R |  61.76% |  0.6949
|  T5large |  8 |  F |  62.24% |  0.6907
Table 4:  Results on dev set for both subtasks, R and F represent rounded and floored scores respectively (as mentioned in Sec. 3.4) Model |  Batch Size |  Scores |  Accuracy (subtask A) |  Spearman’s rank (subtask B)
---|---|---|---|---
RoBERTalarge |  8 |  F |  61.00% |  0.6700
DeBERTa-V1large |  8 |  F |  61.00% |  0.6900
T5-Encoderlarge |  16 |  F |  62.40% |  0.7070
Table 5:  Official results on leaderboard for both subtasks

##  6 Conclusion

In this work we developed a set of transformer-based models powered with
ordinal regression to tackle an NLP cloze task as part of the SemEval 2022
shared task 7. The goal was to assign suitability scores or labels to several
different provided fillers given each context and masked span.

Using a combined architecture based on ordinal regression that used training
labels from both subtasks, we developed and trained models with a multi-task
learning objective. The proposed system was ranked $5$th out of 21 teams in
the ranking track of the shared task. In the subsequent analyses in the post
evaluation phase, we have showed the effectiveness of this architecture in
addressing this task. We compared different variations of our models and
explored the effects of hyperparameters on model performance. The code and
analyses are going to be publicly available.

## References

  * Anthonio et al. (2020) Talita Anthonio, Irshad Bhat, and Michael Roth. 2020.  wikiHowToImprove: A resource and analyses on edits in instructional texts.  In _Proceedings of the 12th Language Resources and Evaluation Conference_ , pages 5721–5729, Marseille, France. European Language Resources Association.
  * Anthonio and Roth (2020) Talita Anthonio and Michael Roth. 2020.  What can we learn from noun substitutions in revision histories? In _Proceedings of the 28th International Conference on Computational Linguistics_ , pages 1359–1370, Barcelona, Spain (Online). International Committee on Computational Linguistics.
  * Bhat et al. (2020) Irshad Bhat, Talita Anthonio, and Michael Roth. 2020.  Towards modeling revision requirements in wikiHow instructions.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 8407–8414, Online. Association for Computational Linguistics.
  * Cao et al. (2020) Wenzhi Cao, Vahid Mirjalili, and Sebastian Raschka. 2020.  Rank consistent ordinal regression for neural networks with application to age estimation.  _Pattern Recognition Letters_ , 140:325–331.
  * Ci et al. (2019) Tianyu Ci, Zhen Liu, and Ying Wang. 2019.  Assessment of the degree of building damage caused by disaster using convolutional neural networks in combination with ordinal regression.  _Remote Sensing_ , 11(23):2858.
  * Debnath and Roth (2021) Alok Debnath and Michael Roth. 2021.  A computational analysis of vagueness in revisions of instructional texts.  In _Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Student Research Workshop_ , pages 30–35, Online. Association for Computational Linguistics.
  * Deutsch and Roth (2019) Daniel Deutsch and Dan Roth. 2019.  Summary cloze: A new task for content selection in topic-focused summarization.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 3720–3729, Hong Kong, China. Association for Computational Linguistics.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  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)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Fu et al. (2018) Huan Fu, Mingming Gong, Chaohui Wang, Kayhan Batmanghelich, and Dacheng Tao. 2018\.  Deep ordinal regression network for monocular depth estimation.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2002–2011.
  * He et al. (2021) Pengcheng He, Jianfeng Gao, and Weizhu Chen. 2021.  Debertav3: Improving deberta using electra-style pre-training with gradient-disentangled embedding sharing.
  * He et al. (2020) Pengcheng He, Xiaodong Liu, Jianfeng Gao, and Weizhu Chen. 2020.  Deberta: Decoding-enhanced bert with disentangled attention.  In _International Conference on Learning Representations_.
  * Hendrycks and Gimpel (2016) Dan Hendrycks and Kevin Gimpel. 2016.  Bridging nonlinearities and stochastic regularizers with gaussian error linear units.  _CoRR_ , abs/1606.08415.
  * Liu et al. (2021) Pengfei Liu, Weizhe Yuan, Jinlan Fu, Zhengbao Jiang, Hiroaki Hayashi, and Graham Neubig. 2021.  Pre-train, prompt, and predict: A systematic survey of prompting methods in natural language processing.  _arXiv preprint arXiv:2107.13586_.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.  _ArXiv_ , abs/1907.11692.
  * Mostafazadeh et al. (2016) Nasrin Mostafazadeh, Nathanael Chambers, Xiaodong He, Devi Parikh, Dhruv Batra, Lucy Vanderwende, Pushmeet Kohli, and James Allen. 2016.  A corpus and cloze evaluation for deeper understanding of commonsense stories.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 839–849, San Diego, California. Association for Computational Linguistics.
  * Niu et al. (2016) Zhenxing Niu, Mo Zhou, Le Wang, Xinbo Gao, and Gang Hua. 2016.  Ordinal regression with multiple output cnn for age estimation.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_.
  * Oxford University Press (2022) Oxford University Press. 2022.  cloze, n.
  * Raffel et al. (2019) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J Liu. 2019.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _arXiv preprint arXiv:1910.10683_.
  * Rogers et al. (2021) Anna Rogers, Matt Gardner, and Isabelle Augenstein. 2021.  Qa dataset explosion: A taxonomy of nlp resources for question answering and reading comprehension.  _arXiv preprint arXiv:2107.12708_.
  * Roth et al. (2022) Michael Roth, Talita Anthonio, and Anna Sauer. 2022.  SemEval-2022 Task 7: Identifying plausible clarifications of implicit and underspecified phrases in instructional texts.  In _Proceedings of the 16th International Workshop on Semantic Evaluation (SemEval-2022)_. Association for Computational Linguistics.
  * Schick and Schütze (2021) Timo Schick and Hinrich Schütze. 2021.  Exploiting cloze-questions for few-shot text classification and natural language inference.  In _Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Main Volume_ , pages 255–269, Online. Association for Computational Linguistics.
  * Sharma et al. (2018) Rishi Sharma, James Allen, Omid Bakhshandeh, and Nasrin Mostafazadeh. 2018.  Tackling the story ending biases in the story cloze test.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , pages 752–757, Melbourne, Australia. Association for Computational Linguistics.

##  Appendix A Appendix

In this section, we are going to share the details of the hyperparameters used
for the fine-tuning of our models and the final training procedure used for
the submission.

###  A.1 Hyperparameters

We fine-tuned all models for $5$ epochs while keeping the embedding layers of
all models frozen. We used the AdamW as our optimizer with a cosine learning
rate schedular from the Hugging Face library and a weight decay of
$0.00123974$ and an initial learning rate of $1.90323e-05$. Also, the
$\lambda_{l}$ and $\lambda_{r}$ (as mentioned in 3.6) are set to $0.5$ in all
of the reported experiments as shown in Table 6.

Hyperparamters |
---|---
Epochs |  5
Optimizer |  AdamW
Learaning Rate Scheduler |  Cosine
Initial Learning Rate |  1.90323e-05
Weight Decay |  0.00123974
$\lambda_{l}$ |  0.5
$\lambda_{r}$ |  0.5
Table 6: Details of the Hyperparamters used for fine-tuning

###  A.2 Training Procedure for Submission

Once we found the best hyperparameters and models using Dev dataset, We used a
combination of Training and Dev data to train our final models for submission.
With this approach, we have been able to achieve a $3\%$ improvement on
accuracy and up to $0.05$ improvement on Spearman’s rank correlation.
Additionally, we found that combining Train and Dev data is noticeably less
effective in the smaller models such as Bertbase or RoBERTabase, however, this
is a better strategy when it comes to larger models such as RoBERTalarge or
T5large. Based on this, we suppose that with further tuning. DeBERTa-V3large
can potentially surpass the state-of-the-art if it has access to the
combination of Train and Dev data.


# Multiwavelength Study of the Quiescent States of Six Brightest Flat Spectrum
Radio Quasars detected by Fermi-LAT

Abhradeep Roy1, S. R. Patel2, A. Sarkar1, A. Chatterjee1, V. R. Chitnis1
1Department of High Energy Physics, Tata Institute of Fundamental Research,
Mumbai-400005, India
2Deutsches Elektronen-Synchrotron, Platanenallee 6, D-15738 Zeuthen, Germany
E-mail<EMAIL_ADDRESS>
###### Abstract

The regular monitoring of flat-spectrum radio quasars (FSRQs) in $\gamma$-rays
by Fermi-LAT since past 12 years indicated six sources who exhibited extreme
$\gamma$-ray outbursts crossing daily flux of 10-5 photons cm-2 s-1. We
obtained nearly-simultaneous multi-wavelength data of these sources in radio
to $\gamma$-ray waveband from OVRO, Steward Observatory, SMARTS, Swift-UVOT,
Swift-XRT and Fermi-LAT. The time-averaged broadband Spectral Energy
Distributions (SEDs) of these sources in quiescent states were studied to get
an idea about the underlying baseline radiation processes. We modeled the SEDs
using one-zone leptonic synchrotron and inverse-Compton emission scenario from
broken power-law electron energy distribution inside a spherical plasma blob,
relativistically moving down a conical jet. The model takes into account
inverse-Compton scattering of externally and locally originated seed photons
in the jet. The big blue bumps visible in quiescent state SEDs helped to
estimate the accretion disk luminosities and central black hole masses. We
found a correlation between the magnetic field inside the emission region and
the ratio of emission region distance to disk luminosity, which implies that
the magnetic field decreases with an increase in emission region distance and
decrease in disk luminosity, suggesting a disk-jet connection. The high-energy
index of the electron distribution was also found to be correlated with
observed $\gamma$-ray luminosity as $\gamma$-rays are produced by high energy
particles. In most cases, kinetic power carried by electrons can account for
jet radiation power as jets become radiatively inefficient during quiescent
states.

###### keywords:

galaxies: active – galaxies: jets – galaxies: individual (3C 273, 3C 279, 3C
454.3, CTA 102, PKS 1510-089, PKS B1222+216) – radiation mechanisms: non-
thermal

††pubyear: 2021††pagerange: Multiwavelength Study of the Quiescent States of
Six Brightest Flat Spectrum Radio Quasars detected by Fermi-LAT–References

##  1 Introduction

Table 1:  Details of studied blazars and selected epochs for quiescent state
study.

Name | R.A. | dec. | z | Tstart | Tstop | $L_{disk}$ | $M_{BH}$
---|---|---|---|---|---|---|---
| (deg) | (deg) |  | DD-MM-YYYY (MJD) | DD-MM-YYYY (MJD) | ($10^{46}$ erg/s) | ($10^{8}$ $M_{\sun}$)
CTA 102 | 338.158 | 11.728 | 1.037 | 04-09-2011 (55808) | 18-10-2011 (55852) | 4.5 a | 8.5 b
3C 279 | 194.045 | -5.786 | 0.536 | 14-01-2010 (55210) | 28-06-2010 (55375) | 0.1 c | 7.9 d
3C 273 | 187.266 | 2.051 | 0.158 | 22-06-2012 (56100) | 11-08-2012 (56150) | 4.8 a | 21-30 e
3C 454.3 | 343.493 | 16.149 | 0.859 | 19-05-2011 (55700) | 30-09-2012 (56200) | 6.75 f | 40 g
PKS 1510-089 | 228.210 | -9.106 | 0.361 | 08-02-2011 (55600) | 19-05-2011 (55700) | 0.5 h | 13 g
PKS B1222+216 | 186.226 | 21.382 | 0.432 | 28-04-2015 (57140) | 07-06-2015 (57180) | 3.5 i | 6 i

  * a

Ghisellini et al. (2010);

  * b

Zamaninasab et al. (2014);

  * c

Paliya et al. (2015);

  * d

Nilsson et al. (2009);

  * e

Paltani & Türler (2005);

  * f

Bonnoli et al. (2011);

  * g

Gu et al. (2001);

  * h

Nalewajko et al. (2012);

  * i

Farina et al. (2012)

Blazars form a subclass of Active Galactic Nuclei (AGN) whose relativistic
particle outflows (jets) are aligned close to our line of sight (Urry &
Padovani, 1995). Blazars are characterized by their significantly variable
emission at all wavelengths across the electromagnetic spectrum ranging from
radio frequencies to very high $\gamma$-ray energies. Blazars can be divided
into two subclasses- BL Lacertae like object (BL Lac) and Flat Spectrum Radio
Quasar (FSRQ). BL Lacs have almost featureless optical-UV spectra. Sometimes
very weak narrow emission lines are present. FSRQs show stronger broad
emission lines in optical spectra. The Fermi Large Area Telescope (LAT) has
detected 3137 blazars above 4$\sigma$ significance in 50 MeV to 1 TeV energy
range (4FGL catalogue), out of which 22 BL Lacs and 43 FSRQs were identified
in other wavebands also (Abdollahi et al., 2020).

Blazars display flux variability on timescales of the order minutes to years
when observed at $\gamma$-ray energies by Fermi-LAT and ground-based Cherenkov
telescopes (Aleksić et al., 2011b; Shukla et al., 2018). From light travel-
time argument, short variability timescales imply radiation from compact
emission regions. Blazars are known for exhibiting occasional outbursts or
flares. But, they spend most of their lifetime in the quiescent state.
Studying these states allow us to probe the most common emission processes at
work, estimating baseline parameters which can be useful in studying high flux
activities. The direct thermal emission from accretion disks of FSRQs in
quiescent states are generally visible in optical-UV waveband, but in other
wavebands the jet emission dominates. Thus the study of quiescent states may
help us to understand the underlying disk-jet connection and jet-empowering
mechanisms.

We selected the six FSRQs (CTA 102, 3C 273, 3C 279, 3C 454.3, PKS 1510-089 &
PKS B1222+216) from the Fermi-LAT monitored source list, also studied in Meyer
et al. (2019), which are known for exhibiting brightest $\gamma$-ray flares
with average daily fluxes over $10^{-5}$ cm-2s-1 within 1$\sigma$ statistical
uncertainties above 100 MeV (Meyer et al., 2019). Detailed multiwavelength
studies on brightest $\gamma$-ray flares of all these sources have been done
by various authors (CTA 102, Gasparyan et al. (2018); 3C 279, Bottacini et al.
(2016), Larionov et al. (2020), Pittori et al. (2018); 3C 273, Rani et al.
(2013), Esposito et al. (2015); 3C 454.3, Diltz & Böttcher (2016); PKS
1510-089, Aleksić et al. (2014); PKS B1222+216, Ackermann et al. (2014),
Bhattacharya et al. (2020)). Except 3C 273, for the other five sources
comparative study on quiescent states with flares have been done (Zacharias et
al., 2017; Hayashida et al., 2012; Boettcher & Baring, 2019; Acciari et al.,
2018). In this work, we carried out a detailed study on quiescent state
spectral energy distributions of these six brightest blazars.

Broadband emission from blazars produces a typical double-humped SED extended
from radio to $\gamma$-ray with low energy peak at IR to X-ray and the high
energy peak at the $\gamma$-ray band. The blazar SEDs are often explained
using leptonic and hadronic models which predict that the synchrotron emission
from electrons rotating in a magnetic field of the jet produces the low energy
peak. It is generally assumed that radiation is emitted by accelerated charged
particles in plasma blobs moving down the jet at relativistic speeds. In
leptonic processes, the high energy peak is produced due to inverse-Compton
(IC) scattering of less energetic seed photons by the relativistic electrons
present in the jet. Seed photons for IC scattering in leptonic model can be
the synchrotron photons themselves (Synchrotron Self-Compton) or photons
entering the emission region from accretion disk, broad-line region (BLR) and
dusty torus region (External Compton). The spectral shape of the $\gamma$-rays
generated by IC scattering is simply related to the energy distribution of the
scatterer particles if the Comtponization lies in the Thomson regime depending
on the energy of target seed photons. For Comptonization in the Klein-Nishina
regime, the relation between $\gamma$-ray spectral shape and energy
distribution of scatterer becomes more complicated. So, $\gamma$-ray spectral
shape can act as an indicator of underlying particle acceleration mechanisms
in the jet (Singh et al., 2020). Earlier studies suggest that underlying
particle acceleration processes in blazar jet can be explained using various
mechanisms like Fermi first order and second-order processes, magnetic
reconnection, acceleration at recollimation shocks etc (Asano & Hayashida,
2020; Joshi et al., 2012; Hervet et al., 2016).

In past, various kinds of emission scenarios have been applied to model the
observed SEDs of these six sources in different states. In many past works,
the blazar SEDs have been described using a hadronic model including photo-
hadronic interactions and pair production process (Spanier & Weidinger, 2012;
Aleksić et al., 2011a; Tchernin et al., 2013). In hadronic models, the high
energy hump in SED is explained as synchrotron emission from massive hadrons
present in the jet, photo-hadronic interactions, pion decay, $e^{-}/e^{+}$
pair production and cascade emission. These hadronic models require very high
magnetic fields ($\sim$100 G) and much harder proton spectrum, which seems to
be unnatural (Bottacini et al., 2016). Some studies invoked multiple emission
regions in order to get more satisfactory fit (Boettcher & Baring, 2019;
Berger et al., 2011; Prince et al., 2019; Tavecchio et al., 2012). Time-
dependant leptonic and hadronic models were used by some authors to calculate
the radiating particle distribution solving the Fokker-Planck equation
considering various acceleration processes and these models are helpful to
understand the evolution of blazar SED during significant flux variations
(Asano & Hayashida, 2018; Zheng et al., 2013; Das et al., 2020; Saito et al.,
2015).

A generalised continuous jet model was applied by Potter & Cotter (2012,
2013a, 2013b, 2013c) to explain blazar jet emission in a quiescent state. They
have assumed a conical jet structure with a magnetic field dominated parabolic
base region where particles get accelerated as suggested from VLBI image of M
87 jet. This model considers emission from the entire jet and so it can
reproduce the observed radio emission when applied on different blazars’ data.
High energy emission has been explained by inverse-Compton scattering of seed
photons coming from cosmic microwave background and narrow-line region.

We used a leptonic one-zone synchrotron and inverse-Compton emission model to
explain observed SEDs during quiescent source states. We imposed a strict
conical jet model to constrain the emission region distance from the central
black hole and found that the low $\gamma$-ray flux in the quiescent state can
be explained by radiation from softer distribution of high energy particles
inside large emission regions placed at distances estimated from a conical jet
model. In section 2 of this paper details of the instruments used and
corresponding data analysis procedures are briefly discussed. Details of the
SED model used for explaining observed data and the results obtained by
modelling are mentioned in section 3. Our inferences on physical processes
inside the jet and correlations between different model parameters and
observables are discussed in section 4. We listed the final conclusions in
section 5.

Figure 1: Fermi-LAT weekly aperture photometry lightcurves of all sources. The
selected quiescent state epochs are shown within vertical dashed blue lines.

##  2 Observations and Data Analysis

Some details of the selected sources and their quiescent epochs (Tstart–Tstop)
are listed in Table 1. Publicly available radio lightcurves from Owens Valley
Radio Observatory (OVRO), optical and infrared (IR) data from SPOL-CCD of
Steward Observatory and Small and Moderate Aperture Research Telescope System
(SMARTS) were used in this work. Optical-UV and X-ray data were obtained from
UV/Optical Telescope (UVOT) and X-ray Telescope (XRT) onboard Neil Gehrels
Swift observatory respectively and analysed. High energy $\gamma$-ray data
were taken from Fermi-Large Area Telescope (LAT) and analysed. In this
section, analysis procedures are discussed. Already published flux measurement
data of all six sources were obtained from the archive of the Space Science
Data Center (SSDC), a facility of the Italian Space Agency
(ASI)111https://www.ssdc.asi.it/.

###  2.1 Radio Observations from OVRO

Except PKSB 1222+216, all other sources are regularly monitored at 15 GHz by
40m telescope of OVRO, which is a part of Fermi monitoring programme. Details
of calibration, observations and data reduction technique of OVRO is given in
Richards et al. (2011). Each source is observed twice a week with a minimum
flux density of about 4 mJy and typical 3% uncertainties. Average flux density
of each source within the epochs of our interest were computed in SED unit
from OVRO lightcurves. Lightcurves of 3C 279, 3C 273 and 3C 454.3 are publicly
available in the OVRO website222https://www.astro.caltech.edu/ovroblazars/.
CTA 102 and PKS 1510-089 lightcurves were kindly provided by OVRO
collaboration. Observations of 3C 454.3 at 22 GHz and 37 GHz were obtained
from 22 m radio telescope (Nesterov et al., 2000) of the Crimean Astrophysical
Observatory (CrAO) and 14 m radio telescope (Teraesranta et al., 1998) of
Aalto University Metsähovi Radio Observatory respectively.

###  2.2 Optical-NIR Obserations from SPOL-CCD of Steward Observatory and
SMARTS

SPOL-CCD imaging/Spectropolarimeter at Steward Observatory of the University
of Arizona is also a part of Fermi multiwavelength support program and
regularly monitors all the six blazars. Publicly available optical V and
R-band photometric, polarimetric data and spectra were obtained from SPOL
archive333http://james.as.arizona.edu/~psmith/Fermi/. Details about these
telescopes, instrumentation and data analysis procedures are discussed in
Smith et al. (2009). Except CTA 102 and PKS B1222+216, observations of other
sources within selected epochs in optical (B, V, R) and near-infrared (J, K)
wavebands were obtained from the SMARTS
archive444http://www.astro.yale.edu/smarts/glast/home.php. Details of SMARTS
telescope design, detectors and data analysis procedures are given in Bonning
et al. (2012) and Buxton et al. (2012). Observed optical-IR fluxes were
corrected for Galactic reddening and extinction using an online
tool555https://irsa.ipac.caltech.edu/applications/DUST/ contributed by IPAC at
Caltech.

###  2.3 Optical-UV Observations from Swift-UVOT

Optical-UV and X-ray data collected by space telescopes on-board Neil Gehrels
Swift-mission were downloaded from NASA HEASARC
archive666https://heasarc.gsfc.nasa.gov/docs/cgro/db-
perl/W3Browse/w3browse.pl. Optical-UV data were obtained from Swift-UVOT
(Roming et al., 2005) which operates in imaging mode. It consists of three UV
filters UM2, UW1, UW2 and three optical filters V, B, U. Swift-UVOT and XRT
data were analysed using tools bundled in HEASOFT (v 6.26) package. Data from
all observations within our selected epochs were integrated using uvotimsum
tool. A circular source region of radius 10″was selected around the source
position in the integrated images from each filter, and background was
extracted from a circular region of radius 50″centred at the source-free
region. Flux magnitudes were obtained using uvotsource tool and then corrected
for galactic extinction of E(B-V) according to Schlegel et al. (1998) for all
filters using a web tool777http://ned.ipac.caltech.edu/forms/calculator.html.
Then the corrected observed magnitudes from all filters were converted into
the fluxes in SED unit using zero-point magnitudes (Poole et al., 2008).

###  2.4 X-Ray Observations from Swift-XRT

Swift-XRT (Burrows et al., 2005) is a grazing incidence telescope which
focusses soft X-rays in 0.3–8 keV energy range on CCD. XRT data were processed
with XRTDAS (v3.5.0) software available within the HEASOFT package (v6.26).
Data from all sources except PKS 1510-089 were recorded in photon counting
mode. PKS 1510-089 data were recorded in windowed timing mode. Event files
were cleaned and calibrated using xrtpipeline(v0.13.5) and source spectra were
obtained using xrtproducts (v0.4.2). Source photons were extracted forming a
circular region of 20 pixel-radius centred at RA and dec of the corresponding
source. Background region was constructed as a circle having a radius of 40
pixels beside and completely detached from the source region. All XRT spectra
within our time intervals of interest were combined using addspec (v1.3.0) and
then grouped to ensure a minimum of 20 counts in each energy bin using grppha
(v3.1.0). Standard auxiliary response files and response matrices were used
during spectral fitting. Combined XRT spectra were fitted in xspec (v12.0.1f)
using power-law and log-parabola (whichever fits better) with the line of
sight absorption in interstellar gas in terms of neutral hydrogen column
density (Kalberla et al., 2005). The power-law model used is given as,

$\frac{dN}{dE}=kE^{-\Gamma}$ (1)

where $\Gamma$ is the spectral index. The log-parabola model is given as,

$\frac{dN}{dE}=k\left(\frac{E}{E_{b}}\right)^{-\alpha-\beta\log\left(E/E_{b}\right)}$
(2)

where $\alpha$ is the spectral index of at $E_{b}$ . During fitting $E_{b}$
was kept fixed at 1 keV. $\beta$ is curvature parameter.

Table 2: 15 GHz average radio fluxes calculated from OVRO lightcurves. Source | CTA 102 | 3C 279 | 3C 273 | 3C 454.3 | PKS 1510-089 | PKS B1222+216
---|---|---|---|---|---|---
flux (Jy) | 2.67$\pm$0.05 | 11.5$\pm$0.2 | 28.9$\pm$0.3 | 11.3$\pm$0.2 | 1.96$\pm$0.04 | –
Table 3: Reddening corrected average SPOL-CCD (R, V-filter) and SMARTS (K, J,
R, V, B-filter) fluxes in SED unit.

Source / Filters $\rightarrow$ | Unit | K | J | R | V | B | R | V
---|---|---|---|---|---|---|---|---
$\downarrow$ | (erg/cm2/s) |  |  |  |  |  | (SPOL) | (SPOL)
CTA 102 | $10^{-12}$ | – | – | – | – | – | 3.33$\pm$0.04 | 3.52$\pm$0.02
3C 279 | $10^{-12}$ | 7.345$\pm$0.005 | 3.647$\pm$0.003 | 2.286$\pm$0.002 | 1.983$\pm$0.002 | 1.601$\pm$0.003 | 1.924$\pm$0.002 | 1.853$\pm$0.002
3C 273 | $10^{-10}$ | 1.255$\pm$0.001 | 0.974$\pm$0.001 | 1.276$\pm$0.001 | 1.766$\pm$0.001 | 1.987$\pm$0.001 | 1.33$\pm$0.02 | 1.61$\pm$0.03
3C 454.3 | $10^{-12}$ | 3.5$\pm$0.4 | 5.3$\pm$0.2 | 5.9$\pm$0.2 | 7.3$\pm$0.1 | 6.7$\pm$0.1 | 6.049$\pm$0.007 | 6.78$\pm$0.01
PKS 1510-089 | $10^{-12}$ | 6.40$\pm$0.02 | 4.64$\pm$0.01 | 5.03$\pm$0.01 | 5.21$\pm$0.01 | 5.86$\pm$0.01 | 4.32$\pm$0.01 | 4.79$\pm$0.02
PKS B1222+216 | $10^{-11}$ | – | – | – | – | – | 1.25$\pm$0.06 | 1.5$\pm$0.1

  * Note: All the SPOL-CCD and SMARTS observations within the selected epochs were averaged and corrected for reddening to obtain the fluxes mentioned above.

Table 4: Summary of the Swift-UVOT data analysis. The fluxes of six UVOT
filters are reported. Same observation IDs listed here were used for both UVOT
and XRT analysis.

Source / Filters $\rightarrow$ |  | V | B | U | W1 | M2 | W2
---|---|---|---|---|---|---|---
$\downarrow$ | Unit |  |  |  |  |  |
CTA 102 |  $10^{-12}$ erg cm-2 s-1 | 4.1$\pm$0.1 | 4.03$\pm$0.09 | 5.09$\pm$ 0.08 | 5.8$\pm$0.1 | 6.6$\pm$0.1 | 5.10$\pm$0.09
3C 279 |  $10^{-12}$ erg cm-2 s-1 | 2.1$\pm$0.1 | 1.86$\pm$0.07 | 1.38$\pm$0.04 | 1.17$\pm$0.03 | 1.21$\pm$0.04 | 1.11$\pm$0.03
3C 273 |  $10^{-10}$ erg cm-2 s-1 | 1.71$\pm$0.04 | 1.98$\pm$0.06 | – | 3.0$\pm$0.1 | 3.5$\pm$0.1 | 3.5$\pm$0.1
3C 454.3 |  $10^{-12}$ erg cm-2 s-1 | 9.2$\pm$0.5 | 8.4$\pm$0.3 | 9.3$\pm$0.3 | 8.8$\pm$0.3 | 9.9$\pm$0.4 | 7.7$\pm$0.2
PKS 1510-089 |  $10^{-12}$ erg cm-2 s-1 | 6.2$\pm$0.5 | 7.0$\pm$0.3 | 7.5$\pm$0.3 | 5.8$\pm$0.3 | 7.1$\pm$0.3 | 6.7$\pm$0.2
PKS B1222+216 |  $10^{-11}$ erg cm-2 s-1 | – | – | – | – | – | 2.46$\pm$0.04

  * Note: Swift observation IDs used are listed below as range of starting observation ID to last observation ID within selected epochs.
CTA 102: (00091094001–00091094018); 3C 279: (00035019021–00035019037);
3C 273: (00035017114–00035017120); 3C 454.3: (00035030204–00035030224);
PKS 1510-089: (00031173064–00031173071); PKS B1222+216:
(00092193001–00092193007).

Table 5: Results from spectral fit of Swift-XRT data. For each of the flux states power law index ($\Gamma$), log parabola index ($\alpha$), curvature parameter ($\beta$), normalization factor of differential spectrum($k$), observed flux in 0.3 to 8.0 keV band (F${}_{0.3-8\ \text{keV}}$) and reduced chi square values ($\chi^{2}_{r}$) are reported.  Source | nH | model | $\Gamma$ | $\alpha$ | $\beta$ | $k$ | F${}_{0.3-8\ \text{keV}}$ | $\chi^{2}_{r}$
---|---|---|---|---|---|---|---|---
| (1020 cm-2) |  |  |  |  | (ph cm-2 s-1 keV-1) | (erg cm-2 s-1) |
CTA 102 | 4.81 | log parabola | – | 1.03$\pm$0.05 | 0.28$\pm$0.08 | (4.05$\pm$0.09)$\times 10^{-4}$ | 3.81$\times 10^{-12}$ | 1.12
3C 279 | 2.24 | log parabola | – | 1.47$\pm$0.03 | 0.27$\pm$0.07 | (1.70$\pm$0.03)$\times 10^{-3}$ | 1.07$\times 10^{-11}$ | 1.05
3C 273 | 1.69 | power law | 1.44$\pm$0.01 | – | – | (8.56$\pm$0.06)$\times 10^{-3}$ | 6.60$\times 10^{-11}$ | 1.24
3C 454.3 | 6.78 | power law | 1.36$\pm$0.05 | – | – | (6.1$\pm$0.3)$\times 10^{-4}$ | 5.18$\times 10^{-12}$ | 0.98
PKS 1510-089 | 7.13 | power law | 1.41$\pm$0.04 | – | – | (8.6$\pm$0.3)$\times 10^{-4}$ | 6.85$\times 10^{-12}$ | 1.07
PKS B1222+216 | 1.72 | power law | 1.49$\pm$0.07 | – | – | (4.2$\pm$0.2)$\times 10^{-4}$ | 3.11$\times 10^{-12}$ | 1.52
Table 6: Results from spectral fit of Fermi-LAT data. For all sources the normalization parameter ($N_{0}$), spectral index ($\alpha$), curvature parameter of the spectrum ($\beta$), integrated flux in 0.1–300 GeV (F${}_{0.1-300\ \text{GeV}}$) and test statistic value of unbinned likelihood analysis (TS) are listed.  Source | $N_{0}$ | $\alpha$ | $\beta$ | F${}_{0.1-300\ \text{GeV}}$ | TS
---|---|---|---|---|---
| (10-10 ph cm-2 s-1 MeV-1) |  |  | (10-7 ph cm2 s-1) |
CTA 102 | 1.0$\pm$0.1 | 2.5$\pm$0.1 | 0.10$\pm$0.08 | 2.1$\pm$0.2 | 305.4
3C 279 | 0.84$\pm$0.05 | 2.46$\pm$0.05 | 0.11$\pm$0.04 | 2.0$\pm$0.1 | 1109.89
3C 273 | 2.0$\pm$0.3 | 2.9$\pm$0.2 | 0.2$\pm$0.2 | 1.9$\pm$0.2 | 158.73
3C 454.3 | 0.26$\pm$0.01 | 2.54$\pm$0.05 | 0.08$\pm$0.04 | 0.99$\pm$0.05 | 1091.34
PKS 1510-089 | 0.500$\pm$0.007 | 2.41$\pm$0.01 | 0.031$\pm$0.006 | 4.17$\pm$0.09 | 1909.65
PKS B1222+216 | 0.59$\pm$0.07 | 2.6$\pm$0.1 | 0 | 1.3$\pm$0.2 | 143.16

###  2.5 $\gamma$-Ray Observations from Fermi-LAT

Fermi Large Area Telescope (LAT) is a pair production space telescope (Atwood
et al., 2009). LAT has field of view $\sim$2.3 Sr and it covers 30 MeV to 1
TeV energy range. Data were analysed using standard software package
fermitools-v1.2.1 provided by Fermi-LAT collaboration and user-contributed
enrico python script (Sanchez & Deil, 2013). Data were collected in the energy
range 0.1–300 GeV from Fermi-LAT data
archive888https://fermi.gsfc.nasa.gov/ssc/data/access/. A circular region of
interest (ROI) having a radius of 15∘ centred at the source was chosen for
event reconstruction from the events belonging to SOURCE class. Events having
zenith angle less than 95∘ were selected to get rid of $\gamma$-ray
contribution from the Earth’s albedo. Good time intervals were selected using
a filter “DATA_QUAL$>$0"&& “LAT_CONFIG==1". The galactic diffused emission
component gll_iem_v07.fits and an isotropic background emission model
iso_P8R3_SOURCE_V2_v1.txt were used as background models. With an instrumental
response function P8R3_SOURCE_V2, unbinned maximum likelihood analysis was
carried out to obtain source spectrum. All the sources lying within ROI+10∘
radius around the source according to fourth Fermi-LAT catalogue (4FGL) were
included in the XML file. All parameters except the scaling factors were
allowed to vary during the fitting process for sources within 5∘ from source
position. The source-spectra were modelled using log-parabola as mentioned in
4FGL catalogue. The flux determination and spectral fitting were carried out
by the likelihood analysis method using GTLIKE tool. Likelihood analysis was
done iteratively by removing all sources having significances less than
1$\sigma$ after each fitting process. The entire energy range of each source
was divided into few bins for obtaining flux points in SED unit. For flux
points having test statistics less than 9 (i.e. $<$3$\sigma$ significance),
flux upperlimits were estimated at 95% confidence level using profile
likelihood method. Fermi unfiltered aperture photometry light-curves were
obtained from their
website999https://fermi.gsfc.nasa.gov/ssc/data/access/lat/msl_lc.

##  3 Results

Fermi-LAT weekly aperture photometry lightcurves were used to select quiescent
states for all the sources (Figure 1). The quiescent states were selected such
that there are a bunch of Swift-XRT observations present during these epochs
so that we can make nearly simultaneous SEDs. It was also checked whether the
sources were in quiescent states in radio, optical-UV, X-ray and $\gamma$-ray
wavebands altogether. Average radio fluxes obtained from OVRO are listed in
Table 2. Average fluxes in optical-IR waveband obtained from Steward
Observatory and SMARTS are shown in Table 3. Results of spectral analysis for
combined UVOT and XRT observations from all sources within selected epochs are
listed in Table 4 and Table 5 respectively. Spectral fit results of Fermi-LAT
data are listed in Table 6. We have made time-averaged broadband spectral
energy distributions of all the six FSRQs during their quiescent states and
modeled them with one-zone leptonic scenario using numerical code ‘Jetset’
provided by Andrea Tramacere (Massaro et al., 2006; Tramacere et al., 2009,
2011; Tramacere, 2020). The model parameters for all sources were studied and
compared in detail.

###  3.1 The SED model

SEDs of the six FSRQs were fitted with leptonic one-zone synchrotron and
inverse Compton model. The inverse Compton process takes into account the
respective contributions from synchrotron photons generated inside the jet
(synchrotron self-Compton or SSC), external photon field coming directly from
the accretion disk, and reprocessed disk photon field coming from BLR and
dusty torus (external Compton or EC).

In this model, the source of broadband emission is assumed to be a spherical
plasma blob of radius ‘$R$’ located at a distance ‘$d$’ from the central
supermassive black hole of mass ‘$M_{BH}$’. The blob is moving
relativistically down the jet with a bulk Lorentz factor ‘$\Gamma$’. Magnetic
field (B) inside the blob is assumed to be same and isotropic everywhere.

A non-thermal population of electrons having energy distribution of broken
power-law shape was considered.

$n\left(\gamma\right)=\left\\{\begin{array}[]{cc}k\gamma^{-p}&\hskip
14.22636pt\gamma\leq\gamma_{break}\\\
k\gamma_{break}^{(p_{1}-p)}\gamma^{-p_{1}}&\hskip
14.22636pt\gamma>\gamma_{break}\\\ \end{array}\right.$ (3)

This electron population gets cooled by synchrotron emission due to
interaction with the magnetic field inside emission blob and generates the low
energy hump in SED. The synchrotron photons get Compton upscattered due to
collision with the synchrotron emitting relativistic electron population
(SSC). Photons emitted by the accretion disk can directly enter the emission
region or can get reprocessed from BLR and dusty torus and enter the emission
blob and get Compton upscattered (EC). Thus, SSC and EC process produces high
energy hump in the SED. As the emission blob moves down the jet, the radiation
gets Doppler boosted in observer’s frame by a factor $\delta$ along our line
of sight given as, $\delta=[\Gamma\left(1-\beta cos\theta\right)]^{-1}$, where
$\Gamma$ is the bulk Lorentz factor of the emission blob and $\theta$ is the
angle between the jet axis and our line of sight. These values for all sources
were taken from literature. The radius of the emission region was constrained
using the light travel-time argument as,

$R\leq\frac{c\delta t_{var}}{2\left(1+z\right)}$ (4)

where $t_{var}$ is the variability timescale. Distance of the emission region
from the central engine, ‘$d$’, was constrained using conical jet model so
that the emission blob fills the entire jet cross-section.

Using accretion disk luminosity ‘$L_{d}$’, BLR distance ($R_{BLR}$) and dusty
torus distance ($R_{DT}$) from central engine were calculated by applying
scaling relations derived from reverberation mapping technique (Ghisellini et
al., 2010).

$R_{BLR}=10^{17}\times(L_{d}/10^{45})^{1/2}$ (5) $R_{DT}=2.5\times
10^{18}\times(L_{d}/10^{45})^{1/2}$ (6)

Assuming that BLR is a thin spherical shell of ionized gases, the inner and
outer radii of BLR were selected around $R_{BLR}$ such that, $R_{BLR}$ =
($R_{BLR}^{in}$+$R_{BLR}^{out}$)/2 and ($R_{BLR}^{out}$-$R_{BLR}^{in}$) =
$2\times 10^{16}$ cm. The temperature of the dusty torus was kept at typical
1000 K for all sources. Multi-temperature blackbody type accretion disk model
was used where the temperature of a portion on the disk depends on its
distance from the core as,

$T^{4}\left(r\right)=\frac{3R_{S}L_{d}}{16\epsilon\pi\sigma_{SB}r^{3}}\left(1-\sqrt{\frac{3R_{S}}{r}}\right)$
(7)

where $\epsilon$ is the accretion efficiency which was set to 0.08 and
$\sigma_{SB}$ is the Stefan-Boltzmann constant. The accretion disk was assumed
to be extended from $3R_{S}$ distance from central engine to a radius of
$500R_{S}$. $\tau_{BLR}$ and $\tau_{DT}$ represents the fraction of accretion
disk emission intercepted and reprocessed by BLR and dusty torus respectively.
Photon fields from BLR and dusty torus get Doppler boosted by a factor
$\sim$$\Gamma^{2}$ in the blob-comoving frame till the blob is inside BLR and
torus. This model first calculates the energy densities and luminosities in
the blob comoving frame. Then, luminosities are converted into flux in
observer frame by calculating luminosity distance ($D_{L}$) from the given
redshift using cosmological model with $\Omega_{\Lambda}=0.685$,
$\Omega_{M}=0.315$ and Hubble’s constant H0= 67.3 km s-1 Mpc-1 (Planck
Collaboration et al., 2014).

###  3.2 SED modelling approach

(i) Radio data do not constrain the emission model. Radio flux point did not
get fitted by synchrotron radiation component because this model takes into
account self-absorption of synchrotron emission at low frequencies below
$\sim$$10^{12}$ Hz by synchrotron emitting electrons in the compact emission
region (synchrotron self-absorption). This implies that other extended regions
in the jet significantly contribute to radio emission where cross-section for
synchrotron self-absorption process is less.
(ii) Optical-UV flux points form a small bump in the SED of most FSRQs. This
bump is interpreted as direct thermal emission from the accretion disk.
Varying disk luminosity ($L_{d}$) and $M_{BH}$ within 10% around the values
quoted in literature (see Table 1) the height and peak position of the thermal
radiation component can be fitted. Good UVOT observations are required for
this process.
(iii) Available data do not give enough coverage to infer peak frequencies of
synchrotron and inverse Compton radiation component. Due to the lack of
simultaneous mm-IR observations, it is not possible to constrain synchrotron
peak position.
(iv) The high energy slope of the electron energy distribution ($p_{1}$) was
constrained using the slope of Fermi-LAT flux points. Value of $p_{1}$ greater
than 3 implies that the peak of the synchrotron and inverse Compton radiation
is produced by the electrons at the break of their energy distribution (i.e.
$\gamma_{peak}=\gamma_{break}$). Sometimes, for low $\gamma_{peak}$ values,
the synchrotron peaks are produced below the self-absorption frequencies. In
these cases, self-absorption frequencies become peak synchrotron frequencies.
(v) The jet viewing angle ($\theta$) and bulk Lorentz factor ($\Gamma$) were
quoted from literature.
(vi) The magnetic field ($B$) and electron density ($N_{e}$) were constrained
by fitting the synchrotron component to the IR observations from SMARTS. But
in the absence of IR data in some cases, low energy end of the X-ray spectrum
is fitted with significant contribution from SSC component (e.g. CTA 102, PKS
1510-089) which help to constrain magnetic field. For blazars having low disk
luminosities (e.g. 3C 279), synchrotron emission dominates over thermal disk
emission in the optical-UV region and thus $B$ and $N_{e}$ can be constrained.
Sometimes even when a big blue bump is visible, the shape of the optical-UV
spectrum ensures some contribution from synchrotron emission along with
thermal disk emission (e.g. CTA 102, 3C 454.3, PKS B1222+216).
(vii) Other parameters like $p$, $\gamma_{min}$ and $\gamma_{max}$ were varied
within feasible values to fit the SEDs by eye estimation due to sparse data
sampling.

Figure 2: Fitted SEDs of the six brightest blazars in quiescent state. Black
points represent analysed flux measurements within the selected epochs and
pink points represent archived data from earlier studies. Yellow, blue,
orange, cyan, black, green and magenta dashed lines indicate synchrotron, disk
thermal, dusty-torus thermal, SSC, ECBLR, ECDT and ECDisk emission component
respectively. The solid red line represent total broadband emission. Table 7:
The input parameters of the model used to reproduce the observed SEDs.

Sr. | Parameters | CTA 102 | 3C 279 | 3C 273 | 3C 454.3 | PKS | PKS
---|---|---|---|---|---|---|---
No. |  |  |  |  |  | 1510-089 | 1222+216
1. |  $R$ (1016 cm) | 4.9 | 6.0 | 2.6 | 1.3 | 2.0 | 6.0
2. |  $\theta$ (degree) | 3.7 a | 2.4 a | 3.0 b | 1.3 a | 2.2 * | 2.5 d
3. | $\Gamma$ | 15.5 a | 20.9 a | 9.0 b | 19.9 a | 20.0 c | 23.0 *
4. | $\delta$ | 15.5 * | 23.67 * | 14.7 * | 33.0 * | 25.0 c | 23.0 d
5. |  $d$ (1017 cm) | 7.5 | 15.0 | 2.0 | 2.5 | 4.0 | 14.0
6. | B (G) | 2.3 | 0.9 | 12.0 | 8.0 | 0.6 | 1.35
7. | $\gamma_{min}$ | 1 | 1.1 | 1.7 | 1.4 | 3 | 1
8. | $\gamma_{max}$ | 4000 | 4500 | 1200 | 5000 | 15000 | 15000
9. | $\gamma_{break}$ | 72 | 340 | 140 | 25 | 350 | 480
10. |  $N_{e}$ (cm-3) | 4200 | 3000 | 5000 | 3000 | 10000 | 1550
11. | $p$ | 1.9 | 2.3 | 2.75 | 2.0 | 2.3 | 2.3
12. | $p_{1}$ | 3.3 | 4.0 | 4.2 | 3.54 | 3.8 | 4.2
|  | Parameters external to the jet |  |
13. |  $L_{d}$ (1046 erg s-1) | 5.0 | 0.1 | 4.8 | 6.75 | 0.5 | 3.6
14. | MBH (M☉) | 8.5E+8 | 7.9E+8 | 3.0E+9 | 4.0E+9 | 1.3E+9 | 6E+8
15. |  $R_{BLR}^{in}$ (1017 cm) | 7.0 | 0.9 | 6.9 | 8.0 | 2.1 | 5.9
16. |  $R_{BLR}^{out}$ (1017 cm) | 7.2 | 1.1 | 7.1 | 8.2 | 2.3 | 6.1
17. | $\tau_{BLR}$ | 0.15 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1
18. | TDT (K) | 1000 | 1000 | 1000 | 1000 | 1000 | 1000
29. |  $R_{DT}$(1018 cm) | 16.7 | 2.5 | 17.3 | 20.5 | 5.6 | 15
20. | $\tau_{DT}$ | 0.1 | 0.39 | 0.1 | 0.1 | 0.1 | 0.1

  * a

Hovatta et al. (2009);

  * b

Paliya et al. (2017);

  * c

Acciari et al. (2018);

  * d

Kushwaha et al. (2014);

  * *

Computed using, $\delta=[\Gamma\left(1-\beta cos\theta\right)]^{-1}$

Table 8: Energy densities, powers and mass accretion rates calculated from the parameters of SED modelling. Sr. | Parameters | CTA 102 | 3C 279 | 3C 273 | 3C 454.3 | PKS | PKS
---|---|---|---|---|---|---|---
No. |  |  |  |  |  | 1510-089 | 1222+216
1. |  $U_{e}$ (erg/cc) | 1.46e-02 | 9.92e-03 | 1.58e-02 | 1.26e-02 | 8.52e-02 | 4.75e-03
2. |  $U_{B}$ (erg/cc) | 4.07e-01 | 3.22e-02 | 5.73 | 2.55 | 1.43e-02 | 7.25e-02
3. |  $\eta$ ($U_{e}$/$U_{B}$) | 0.036 | 0.308 | 0.003 | 0.005 | 5.958 | 0.065
4. |  $P_{rad}$ (erg/s) | 3.32e+45 | 2.70e+44 | 7.63e+43 | 1.42e+44 | 8.65e+43 | 1.28e+44
5. |  $P_{kin}$ (erg/s) | 4.74e+46 | 7.30e+46 | 3.35e+46 | 1.61e+46 | 2.41e+46 | 5.56e+46
6. |  $P_{e}$ (erg/s) | 7.92e+44 | 1.47e+45 | 8.17e+43 | 7.92e+43 | 1.28e+45 | 8.53e+44
7. |  $P_{p}$ (erg/s) | 2.45e+46 | 6.68e+46 | 3.88e+45 | 2.84e+45 | 2.26e+46 | 4.18e+46
8. |  $P_{B}$ (erg/s) | 2.21e+46 | 1.30e+46 | 2.95e+46 | 1.61e+46 | 2.15e+44 | 8.74e+44
9. |  $P_{jet}$ (erg/s) | 4.74e+46 | 5.63e+46 | 3.35e+46 | 1.90e+46 | 2.41e+46 | 5.03e+46
10. |  $L_{Edd}$ (erg/s) | 1.07e+47 | 9.95e+46 | 3.78e+47 | 5.04e+47 | 1.64e+47 | 7.56e+46
11. | $L_{d}/L_{Edd}$ | 0.47 | 0.01 | 0.13 | 0.13 | 0.03 | 0.48
12. |  $\dot{M}_{accr}$ ($M_{\sun}/yr$) | 10.95 | 0.22 | 10.51 | 14.78 | 1.10 | 7.88
13. |  $\dot{M}_{out}$ ($M_{\sun}/yr$) | 0.028 | 0.056 | 0.007 | 0.003 | 0.020 | 0.032
14. | $\dot{M}_{out}/\dot{M}_{accr}$ | 0.003 | 0.256 | 0.001 | 0.0002 | 0.018 | 0.004
15. | $\dot{M}_{accr}/\dot{M}_{Edd}$ | 5.84 | 0.13 | 1.59 | 1.67 | 0.38 | 5.95

###  3.3 SED fit results

SEDs of six sources fitted with a one-zone leptonic model are shown in Figure
2 and the sets of fit parameters are listed in Table 7. The trends shown by
the archived data helped to model the SEDs. When all UVOT filters were not
available, the underlying trend of archived data helped us to decide whether
optical-UV were forming a part of the disk thermal emission bump. These trends
also helped in modelling the height and shape of the synchrotron peak. In some
earlier studies of FSRQs, the X-ray spectra were fitted by SSC component and
$\gamma$-ray spectra were fitted with EC component (Yan et al., 2015;
Sahakyan, 2020; Pacciani et al., 2014; Hayashida et al., 2012; Paliya et al.,
2015). Though hard X-ray data were not available in our case, the underlying
trends in archived data implied that soft X-ray and $\gamma$-ray spectra form
the rising and decaying part of a single broad EC component. For this purpose,
the $\gamma_{min}$ values were kept low (within 1–10). Except for 3C 279 and
PKS B1222+216, the single broad EC components could not describe the low
energy ends of the soft X-ray spectra. The contributions from EC emissions of
Doppler de-boosted disk radiation ($EC_{disk}$) were required in 3C 273 and 3C
454.3 SED, and SSC contribution was required in CTA 102 and PKS 1510-089 SED.
Synchrotron emission flux from 3C 273 is the highest and that of 3C 454.3 is
the lowest. PKS B1222+216 has the lowest EC flux and 3C 273 has the highest.

Derived black hole masses ($M_{BH}$) were found to be in the range
$(7-40)\times 10^{8}M_{\sun}$ and Disk luminosities were found to be in the
range $(1-68)\times 10^{45}$ erg s-1 for our selected sources. All these
values are consistant with values mentioned in literature. In the case of 3C
279, the optical-UV data did not indicate a dominance of disk thermal emission
over synchrotron radiation. So, in this case, $L_{d}$ and $M_{BH}$ values were
quoted from literature. Derived disk luminosities are in between 1 and 48 per
cent of the Eddington luminosity ($L_{Edd}$).

For four out of six sources, inverse Compton emission was found to be
dominated by $EC_{BLR}$ component, though the emission region is inside BLR in
only two sources (3C 273 and 3C 454.3). $\tau_{BLR}$ values were found to be
10% (except 15% for CTA 102) and $\tau_{DT}$ values are 10% (except 38% for 3C
279).

Magnetic fields were found in range ($0.6-12$) G. Electron densities of
($10^{3}-10^{4}$) cm-3 in the emission region indicate the presence of very
powerful jets in FSRQs. Low energy slopes of the particle distribution ($p$)
were found to be greater than 2 for four sources except CTA 102 and 3C 454.3.
Values of $p_{1}$ were found to be greater than 3 and break energies were
found in the range ($30-500$).

###  3.4 Power estimation

Relativistic jet contains power in form of radiation ($P_{rad}$), electrons
($P_{e}$), protons ($P_{p}$) and magnetic field ($P_{B}$). Powers are
estimated mainly using,

$P_{i}\simeq\pi R^{2}\Gamma^{2}\beta cU_{i}$ (8)

where $U_{i}$ is the energy density of the $i^{th}$ power carrying component
in the comoving frame. Thus, power carried by the relativistic electron
population is given as,

$P_{e}\simeq\pi R^{2}\Gamma^{2}\beta\langle\gamma\rangle N_{e}m_{e}c^{3}$ (9)

and power carried by ‘cold’ protons (as they are massive they do not get
accelerated enough to move relativistically) are estimated as,

$P_{p}\simeq\pi R^{2}\Gamma^{2}\beta N_{p}m_{p}c^{3}$ (10)

where $\langle\gamma\rangle$ is the average Lorentz factor of relativistic
electrons, $m_{e}$ and $m_{p}$ are rest masses of electron and proton
respectively, $N_{e}$ and $N_{p}$ are electron and proton density respectively
inside the emission blob. $N_{e}$ goes as input parameter in the model. We
assumed presence of one proton per ten electrons inside blob (Ghisellini et
al., 2014) to calculate $P_{p}$. The power carried by magnetic field as
Poynting flux is given as,

$P_{B}\simeq\frac{1}{8}R^{2}\Gamma^{2}\beta cB^{2}$ (11)

Total jet power ($P_{jet}$) is the sum of powers carried by electrons, protons
and magnetic field. Power radiated by the jet ($P_{rad}$) can be restated
using radiation energy density in comoving frame ($U_{r}=L^{\prime}/(4\pi
R^{2}c)$) as,

$P_{rad}=L^{\prime}\frac{\Gamma^{2}}{4}=L\frac{\Gamma^{2}}{4\delta^{4}}$ (12)

where $L^{\prime}$ is the total non-thermal radiation luminosity in comoving
frame and $L$ is the total luminosity in observer’s frame. Four out of six
FSRQs have jet powers in quiescent state greater than or almost equal to disk
luminosities, which indicates that FSRQs have very powerful jets (Figure 3).
Calculated powers are listed in Table 8. Total jet powers of all sources were
found to be less than the Eddington luminosities.

###  3.5 Estimation of mass outflow rate

Mass accretion rate can be estimated from disk luminosity as,

$\dot{M}_{accr}=\frac{L_{d}}{\epsilon c^{2}}$ (13)

Mass outflow rate can be estimated as,

$\dot{M}_{out}=\frac{P_{jet}}{\Gamma c^{2}}$ (14)

Except 3C 279, all other sources require more than 1 $M_{\sun}$ yr-1 mass
accretion rate. But the mass outflow rate is $\sim 0.003-0.3$ times of
accretion rate. This can be understood as $\epsilon/\Gamma<<1$ and $P_{jet}$
and $L_{d}$ have same orders of magnitude.

Figure 3: Disk luminosity vs jet power plot. Figure 4: Dependence of break
energy on difference of emission region distance and BLR distance. Figure 5:
Dependence of magnetic field (B) in jet on central black hole mass ($M_{BH}$).
Figure 6: Dependence of magnetic field on distance of emission region ($d$),
black hole mass ($M_{BH}$) and disk luminosity ($L_{d}$) fitted with
bootstrapped linear regression. Blue lines are bootstrapped regression lines
and the black dotted line is the average of them. Figure 7: Dependence of
intrinsic $\gamma$-luminosity (0.1 – 300 GeV) on high energy slope of particle
energy distibution ($p_{1}$) fitted with bootstrapped linear regression. Blue
lines are bootstrapped regression lines and the black dotted line is the
average of them.

###  3.6 Comparison with earlier studies

In this subsection the obtained results are compared with different results
published earlier. For reference see Table 6, Table 7 and Table 8.
(i) CTA 102: Gasparyan et al. (2018) studied the extreme outburst of 2016-17
and also the pre-flare low state. The low state and flaring states’ Fermi-LAT
spectral indices are mentioned to be 2.39 and 2.01–1.81 respectively. They
reported $P_{e}$$\sim$$10^{45}$ erg/s and $P_{B}$$\sim$$10^{46}$ erg/s in the
low state, and $P_{e}$$\sim$$(10^{45}-10^{47})$ erg/s and
$P_{B}$$\sim$$(10^{42}-10^{39})$ erg/s in the flaring states. In our quiescent
state analysis, we obtained steeper Fermi-LAT spectral index (2.5), lower
$P_{e}$ ($\sim$$10^{44}$ erg/s) and similar $P_{B}$ ($\sim$$10^{46}$ erg/s)
values as expected.
(ii) 3C 279: Comparative study of low state, flares and post-flare state SED
has been carried out by Pittori et al. (2018). The Fermi-LAT spectral indices
during two low states and two flares are (1.62, 1.77) and (1.53, 1.47)
respectively. We obtained a steeper Fermi index of 2.46. Pittori et al. (2018)
reported $P_{e}$$\sim$$10^{44}$ erg/s, $P_{B}$$\sim$$10^{44}$ erg/s and
$P_{p}$$\sim 10^{46}$ across all states, and $P_{rad}$$\sim$$10^{44}$,
$10^{46}$ erg/s during low and high states respectively. We obtained higher
$P_{e}$ and $P_{B}$ ($\sim$$10^{45}$ erg/s), similar $P_{p}$ but lower
$P_{rad}$ ($\sim$$10^{44}$ erg/s) in quiescent state.
(iii) 3C 273: Rani et al. (2013) did a comparative study of Fermi-LAT spectrum
during a low state and five different flare states. They reported spectral
index of 2.86 during low state and 2.5–2.8 during the flare states. Broadband
SED analysis of quiescent state is not yet done for this source. We got
similar index of 2.9 during the selected quiescent state.
(iv) 3C 454.3: Diltz & Böttcher (2016) did broadband SED modeling of an
average state using data from Abdo et al. (2010) and reported $P_{e}=6.1\times
10^{45}$ and $P_{B}=1.18\times 10^{45}$ erg/s, which are greater and less than
values obtained by us respectively.
(v) PKS 1510-089: Acciari et al. (2018) did exclusive study of low-state
broadband SED with MAGIC observations and reported spectral index of 2.56.
Whereas we obtained spectral index of 2.41. They did not mention resulting
powers carried by jet. So, we used their SED models to calculate jet powers
for comparison purpose and obtained $P_{e}$, $P_{B}$ and $P_{rad}$ of the
order of $10^{45}$, $10^{44}$ and $10^{44}$ erg/s respectively, which support
our results.
(vi) PKS B1222+216: Bhattacharya et al. (2020) carried out detailed comparison
between flaring and quiescent states and reported Fermi spectral indices of
2.5 during quiescence and 1.7–2.3 during flares, which is comparable to
spectral index obtained in this work (2.6). They also did not mention
resulting jet powers from SED modelings. Kushwaha et al. (2014) reported total
jet power of $2\times 10^{46}$ erg/s during the brightest flare, which is less
than what we obtained.
Thus, during quiescent states we get steeper Fermi-LAT spectra as compared to
flare states. Calculated jet power budgets during quiescent state in this work
are comparable to the values mentioned in literature. But, there is no clear
trend of variation of $P_{e}$ and $P_{B}$ with change in flux. Their changes
depend on the type of the process used to describe the flux variation (eg.
propagation of shock-front, geometric effects, magnetic reconnection etc).

##  4 Discussion

(i) Black hole mass and disk luminosity: It is worth noticing that in a
quiescent state, thermal emission from accretion disk dominates over
synchrotron radiation in the optical-UV band for all sources (except 3C 279).
The trend shown by archival data underlying the actual analysed data in the
optical-UV band indicates the presence of a big blue bump (BBB) which has a
thermal origin. Thus, fitting the optical-UV band using thermal emission from
multi-temperature blackbody type accretion disk helps to estimate disk
luminosity and the central supermassive black hole mass. According to blazar
sequence, generally synchrotron emission in FSRQs peak in the infrared band
and contribute less in the optical-UV band. Thus, thermal emission from disk
becomes visible in the form of BBB in SEDs. As the selected FSRQs are the
brightest $\gamma$-ray emitters, they contain black holes with very large
masses ($\ga 8\times 10^{8}M_{\sun}$) as expected. They also have very high
mass accretion rate of (0.1–6) in Eddington units ($>$1$M_{\sun}$ yr-1).

(ii) Emission region size: The results listed in the previous section lead to
the discussion on how the study of quiescent states of FSRQs is important.
Though FSRQs are known for showing highly variable emission in all wavebands
and occasional outbursts, they remain mostly in quiescent state during their
lifetime. Previous studies on FSRQ flares have reported significant fast
variability timescales of less than a day. Even timescales of the order of
minutes have been reported in $\gamma$-ray waveband (Shukla et al., 2018;
Ackermann et al., 2016; Nalewajko et al., 2012; Aleksić et al., 2011c). But,
from the Fermi-lightcurves (Figure 1) it can be seen that FSRQs show lower
flux variability during quiescent states. The emission region sizes ($R$) of
the order of 1016 cm used for modeling quiescent SEDs (Table 7) correspond to
variability timescales of 1 to 5 days, which is larger than minute-hour
variability timescales observed for these sources during flares. As expected,
this implies slower variability in quiescence for all sources which
corresponds to larger emission regions according to light travel-time
argument.

(iii) Emission mechanisms: According to conical jet model, larger emission
regions were to be placed at larger distances from the central engines. Thus,
emission regions for CTA 102, 3C 279, PKS 1510-089 and PKS B1222+216 were
found to be outside BLR. Emission regions for 3C 273 and 3C 454.3 were found
to be inside BLR. When $d<R_{BLR}$, seed photon contribution from BLR is
maximum. High energy $\gamma$-ray data is then fitted by dominant $EC_{BLR}$
component with negligible contribution from $EC_{DT}$ component. When the
emission region is placed outside BLR, $EC_{DT}$ component becomes important
in describing $\gamma$-ray emission. In case of CTA 102 and PKS 1510-089 the
emission regions were close to the outer edge of the BLR and far away from the
dusty torus. That is why, despite being Doppler de-boosted in blob rest-frame,
the EC components from BLR photon field ($EC_{BLR}$) dominate over EC
component from torus photon field ($EC_{DT}$). In PKS B1222+216 and 3C 279,
the emission regions are further away from BLR which indicate the dominance of
$EC_{DT}$ over $EC_{BLR}$ as non-thermal high energy emission process.
$\tau_{BLR}$ values for all sources were kept at 10% except for CTA 102
($\tau_{BLR}$ = 15%), which is quite reasonable. A bit higher $\tau_{DT}$
value was required for 3C 279 to fit the SED. As dusty torus is a huge
doughnut-shaped cloud surrounding the AGN and in 3C 279 it is quite close to
the central engine due to less radiation pressure generated by low disk
luminosity, we may assume that it covers larger solid angle to intercept a
higher fraction of disk luminosity. Having the emission region placed outside
BLR imply reduced photon field inside the emission region which makes the
cooling process weaker. Thus, it is expected to have higher break energy in
electron distribution where the difference between $d$ and $R_{BLR}$ is
higher, i.e. the cooling process is weaker. Figure 4 shows that CTA 102, PKS
1510-089, PKS B1222+216 and 3C 279 follow this trend. Due to high variability
in flare states, the emission region, being smaller in size, generally placed
inside BLR. Thus, Doppler-boosted high $EC_{BLR}$ component is produced which
can explain bright $\gamma$-ray flux.

(iv) Magnetic fields: According to Ghisellini et al. (2010), magnetic field in
FSRQ jet is (1–10) G on an average. But we found that in quiescent state two
sources (3C 279 and PKS 1510-089) have magnetic fields less than 1 G. From
Figure 5, it seems that FSRQs containing more massive black holes have a
higher magnetic field in the jet. But, this does not tell the whole picture.
As the source of the magnetic field is the mass accretion process, we can
expect a correlation between disk luminosity and magnetic field. The magnetic
field inside the emission region should decrease with an increase in emission
region distance from the central engine. Thus, all these effects together can
explain the nice correlation (Pearson Correlation = -0.83) between ‘$B$’ and
‘$\log_{10}(d/(R_{S}*L_{d}))$’ shown in Figure 6 as, $B$ = (-133 $\pm$ 38) -
(3.2$\pm$0.9) $\log_{10}(d/(R_{S}*L_{d}))$.

(v) Particle energy distribution: FSRQs have very powerful jets containing a
highly dense population of charged particles. So, these particles get cooled
very efficiently due to collision with each other. Because of this short
cooling time-scale, injected particle distribution cannot get accelerated up
to very high energy and cools down rapidly before escaping the emission
region. So, $\gamma_{min}$ value for all sources were kept less than 10. Use
of low $\gamma_{min}$ value yields very wide $EC_{BLR}$ or $EC_{DT}$ component
which can simultaneously describe both X-ray and $\gamma$-ray data in flare
state (Aleksić et al., 2011a). But, in quiescent state SEDs, a significant
contribution of SSC or $EC_{disk}$ component is required depending on the
distance of emission region from the central engine to fit low energy X-ray
data along with $EC_{BLR}$ or $EC_{DT}$. Electron densities ($N_{e}$) found
are quite high ($\sim$$10^{3}$–$10^{4}$ cm-3) which imply powerful jets. We
used electron energy distribution of broken power-law shape. This kind of
particle distribution can be produced by underlying Fermi first and second-
order particle acceleration processes (Lewis et al., 2019). The slopes of both
synchrotron and EC component above peak frequencies are dependent on the high
energy index of particle energy distribution ($p_{1}$). For all sources,
‘$p_{1}$’ value was found to be more than 3, which implies that electrons
having energy around the break in energy distribution radiate the maximum
amount of energy. Higher $p_{1}$ value implies less number of radiating
electrons at higher energies which in turn result in lower $\gamma$-ray
luminosity. As expected, the $\gamma$-ray luminosity observed by Fermi-LAT was
found to be decreasing (Pearson Correlation = -0.85) with $p_{1}$ as, $p_{1}$
= (24$\pm$4) - (0.4 $\pm$ 0.1)$\log_{10}L_{\gamma}$ (Figure 7). We found
$\gamma_{max}$ values of the order of $\sim$($10^{3}$–$10^{4}$). To describe
flare SEDs, higher $\gamma_{max}$ and $\gamma_{break}$ values and harder
particle spectrum are required (Bottacini et al., 2016; Pacciani et al., 2014;
Ahnen et al., 2017).

(vi) Power budget: Radiated power of only two sources (CTA 102 and 3C 454.3)
are higher than the powers carried by electrons. Considering one proton per
ten electrons, kinetic power carried by protons ($P_{p}$) is found to be much
higher than $P_{e}$ in all sources. Even in quiescent states, Poynting flux
carries a significant amount of power ($P_{B}\ga P_{e}$). Thus, even though
power carried by electrons is not always sufficient to account for the
observed radiation power, proton kinetic power and Poynting flux power are
always individually sufficient in a quiescent state. Ghisellini et al. (2010)
reported that, for FSRQs in average state $P_{rad}$ is greater than $P_{e}$.
But, in the quiescent state, we can see that most of the brightest FSRQs have
sufficient kinetic power carried by electrons to account for observed
radiation. Jets can be powered by mass accretion process and the black hole
spin. To explain the $P_{jet}\geq L_{d}$ scenario, Jolley et al. (2009)
proposed that a large fraction of accretion power is spent to power the jet
and a smaller fraction ($\epsilon\sim 0.08$) produces disk luminosity.
Ghisellini et al. (2010) mentioned that we can have $P_{jet}\geq L_{d}$ if
there is an efficient way to extract energy from the spinning central black
hole. According to the Blandford-Znajek (BZ) mechanism of extracting the
energy of a spinning black hole (Blandford & Znajek, 1977), $P_{jet}$ is
proportional to $L_{d}$ where the value of proportionality constant can be
greater than unity.

(vii) Application: FSRQs mostly remain in a state of low activity during their
lifetime. So, the quiescent state SED models represent the underlying
consistent emission process in the FSRQs. The high activity states can be
considered as perturbations of some physical parameters in the jet environment
or emergence of any additional emission region. Thus, for modeling brighter
and more variable states of the FSRQs selected in this work, emission from a
smaller region should be added to the quiescent state SED model while fitting
the SED data points. Morever, time-dependant leptonic models can be used to
describe the evolution of flares from the quiescent state (Diltz & Böttcher,
2014). The quiescent SED model parameters can be used as the starting point of
the time-evolution process. The observed cross-band correlations and
associated time-lags can be explained as perturbations of some of those
initial parameters.

##  5 Conclusion

After extensive study of quiescent state spectral energy distributions of six
brightest FSRQ detected by Fermi-LAT, we have arrived at the following
conclusions.
(i) For most of the FSRQs in quiescence, low synchrotron emission is dominated
by thermal emission from accretion disk in optical-UV waveband. Thus, disk
luminosity and black hole mass can be estimated using Swift-UVOT data.
Estimated black hole masses and disk luminosities are consistent with values
mentioned in literature (Table 1 and Table 7).
(ii) Slow variability of $>$ 1-day timescale in the quiescent state indicates
a large emission region size.
(iii) Following conical jet structure and assuming that the emission regions
fill the whole jet cross-section, the large emission regions were found to be
mostly outside BLR. Thus, low $\gamma$-ray flux in quiescent state can be
described by Doppler-de-boosted (in blob rest frame) weak $EC_{BLR}$ component
and Doppler-boosted (in blob rest frame) weak $EC_{DT}$ component (as
$R_{DT}\simeq 25R_{BLR}$) depending on emission region distance from BLR.
(iv) In most of the quiescent FSRQ SED models, a small bump is clearly visible
in the IR waveband due to the contribution of thermal emission from dusty
torus. If nearly simultaneous mm-IR observations were available, it would have
been possible to constrain torus temperature and the fraction of disk
luminosity intercepted by the dusty torus.
(v) In the quiescent state, blazar jets become radiatively inefficient. In
most of the brightest FSRQs, the kinetic power carried by electrons can
account for observed radiation. But still, the jets remain powerful enough to
radiate more than disk luminosity.
(vi) We could not constrain the synchrotron peak positions properly in the SED
models due to lack of simultaneous mm-IR observations. So we can roughly say,
in quiescent states of CTA 102, 3C 279 and 3C 454.3 synchrotron and IC
emission peak at lower frequencies than those in flare states.
(vii) Magnetic field inside the emission region has an expected correlated
with ‘$d/(R_{S}*L_{d})$’ for these six sources during quiescent state.
(viii) The high energy indices of particle spectrum ($p_{1}$) of these six
sources are found to be correlated with quiescent state $\gamma$-ray
luminosities ($L_{\gamma}$) observed by Fermi-LAT in 0.1–300 GeV energy range.
For FSRQs, Fermi-LAT data form the decaying part of the high energy hump in
SED, i.e. $L_{\gamma}$ represents the luminosity emitted by particles with
higher energies. Higher $p_{1}$ value implies lower number of high energy
particles (Equation 3), which in turn reduces $L_{\gamma}$.

## Acknowledgements

We would like to thank the anonymous referee whose valuable suggestions truly
helped to improve our work. This research has made use of data from the OVRO
40-m monitoring program (Richards, J. L. et al. 2011, ApJS, 194, 29) which is
supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and
NSF grant AST-0808050 and AST-1109911. We are grateful to Sebastian Kiehlmann
for providing us non-public OVRO lightcurves of some sources. Data from the
Steward Observatory spectropolarimetric monitoring project were used. This
program is supported by Fermi Guest Investigator grants NNX08AW56G,
NNX09AU10G, NNX12AO93G, and NNX15AU81G. This research has made use of the XRT
Data Analysis Software (XRTDAS) developed under the responsibility of the
Space Science Data Center (SSDC) managed by Italian Space Agency (ASI). Data
obtained at the Metsähovi Radio Observatory operated by the Aalto University
were also used in this publication. We have used Fermi-LAT data, obtained from
the Fermi Science Support Center, provided by NASA’s Goddard Space Flight
Center (GSFC). The data and analysis software were obtained from NASA’s High
Energy Astrophysics Science Archive Research Center (HEASARC), a service of
GSFC. We used a community-developed Python package named Enrico to make Fermi-
LAT data analysis easier and more convenient (Sanchez & Deil, 2013). Finally,
We acknowledge the support of the Department of Atomic Energy, Government of
India, under Project Identification No. RTI 4002.

##  6 DATA AVAILABILITY STATEMENT

OVRO data is available at: https://www.astro.caltech.edu/ovroblazars/.
Steward Observatory SPOL-CCD data is available at:
http://james.as.arizona.edu/~psmith/Fermi/.
SMARTS data is available at: http://www.astro.yale.edu/smarts/glast/home.php.
Swift-mission data is available at:
https://heasarc.gsfc.nasa.gov/docs/archive.html.
Fermi-LAT data is available at: https://fermi.gsfc.nasa.gov/ssc/data/access/.
SSDC (ASI) archive: https://www.ssdc.asi.it/.
Swift-mission data analysis software package is available at:
https://heasarc.gsfc.nasa.gov/lheasoft/download.html.
Fermi-LAT data analysis software is available at:
https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/.
Software for Broadband SED modeling is available at:
https://github.com/andreatramacere/jetset.

## References

  * Abdo et al. (2010) Abdo A. A., et al., 2010, ApJ, 716, 30
  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Acciari et al. (2018) Acciari V., et al., 2018, Astron. Astrophys., 619, A159
  * Ackermann et al. (2014) Ackermann M., et al., 2014, ApJ, 786, 157
  * Ackermann et al. (2016) Ackermann M., et al., 2016, The Astrophysical Journal, 824, L20
  * Ahnen et al. (2017) Ahnen M. L., et al., 2017, A&A, 603, A29
  * Aleksić et al. (2011a) Aleksić J., et al., 2011a, A&A, 530, A4
  * Aleksić et al. (2011b) Aleksić J., et al., 2011b, ApJ, 730, L8
  * Aleksić et al. (2011c) Aleksić J., et al., 2011c, The Astrophysical Journal, 730, L8
  * Aleksić et al. (2014) Aleksić J., et al., 2014, A&A, 569, A46
  * Asano & Hayashida (2018) Asano K., Hayashida M., 2018, ApJ, 861, 31
  * Asano & Hayashida (2020) Asano K., Hayashida M., 2020, ApJ, 890, 56
  * Atwood et al. (2009) Atwood W. B., et al., 2009, ApJ, 697, 1071
  * Berger et al. (2011) Berger K., Becerra González J., Bonnoli G., Maraschi L., Tavecchio F., Domínguez A., Lindfors E., Stamerra A., 2011, arXiv e-prints, p. arXiv:1110.6371
  * Bhattacharya et al. (2020) Bhattacharya D., Mohana A K., Bhattacharyya S., Bhatt N., Stalin C. S., 2020, Monthly Notices of the Royal Astronomical Society, 500, 1127
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Boettcher & Baring (2019) Boettcher M., Baring M. G., 2019, arXiv e-prints, p. arXiv:1903.12381
  * Bonning et al. (2012) Bonning E., et al., 2012, ApJ, 756, 13
  * Bonnoli et al. (2011) Bonnoli G., Ghisellini G., Foschini L., Tavecchio F., Ghirlanda G., 2011, MNRAS, 410, 368
  * Bottacini et al. (2016) Bottacini E., Böttcher M., Pian E., Collmar W., 2016, ApJ, 832, 17
  * Burrows et al. (2005) Burrows D. N., et al., 2005, Space Sci. Rev., 120, 165
  * Buxton et al. (2012) Buxton M. M., Bailyn C. D., Capelo H. L., Chatterjee R., Dinçer T., Kalemci E., Tomsick J. A., 2012, AJ, 143, 130
  * Das et al. (2020) Das A. K., Prince R., Gupta N., 2020, ApJS, 248, 8
  * Diltz & Böttcher (2014) Diltz C., Böttcher M., 2014, Journal of High Energy Astrophysics, 1, 63
  * Diltz & Böttcher (2016) Diltz C., Böttcher M., 2016, ApJ, 826, 54
  * Esposito et al. (2015) Esposito V., Walter R., Jean P., Tramacere A., Türler M., Lähteenmäki A., Tornikoski M., 2015, A&A, 576, A122
  * Farina et al. (2012) Farina E. P., Decarli R., Falomo R., Treves A., Raiteri C. M., 2012, MNRAS, 424, 393
  * Gasparyan et al. (2018) Gasparyan S., Sahakyan N., Baghmanyan V., Zargaryan D., 2018, ApJ, 863, 114
  * Ghisellini et al. (2010) Ghisellini G., Tavecchio F., Foschini L., Ghirlanda G., Maraschi L., Celotti A., 2010, Monthly Notices of the Royal Astronomical Society, 402, 497
  * Ghisellini et al. (2014) Ghisellini G., Tavecchio F., Maraschi L., Celotti A., Sbarrato T., 2014, Nature, 515, 376
  * Gu et al. (2001) Gu M., Cao X., Jiang D. R., 2001, MNRAS, 327, 1111
  * Hayashida et al. (2012) Hayashida M., et al., 2012, ApJ, 754, 114
  * Hervet et al. (2016) Hervet O., Boisson C., Sol H., 2016, A&A, 592, A22
  * Hovatta et al. (2009) Hovatta T., Valtaoja E., Tornikoski M., Lähteenmäki A., 2009, A&A, 498, 723
  * Jolley et al. (2009) Jolley E. J. D., Kuncic Z., Bicknell G. V., Wagner S., 2009, Monthly Notices of the Royal Astronomical Society, 400, 1521
  * Joshi et al. (2012) Joshi M., et al., 2012, arXiv e-prints, p. arXiv:1206.6147
  * Kalberla et al. (2005) Kalberla P. M. W., Burton W. B., Hartmann D., Arnal E. M., Bajaja E., Morras R., Pöppel W. G. L., 2005, A&A, 440, 775
  * Kushwaha et al. (2014) Kushwaha P., Singh K. P., Sahayanathan S., 2014, ApJ, 796, 61
  * Larionov et al. (2020) Larionov V. M., et al., 2020, Monthly Notices of the Royal Astronomical Society, 492, 3829
  * Lewis et al. (2019) Lewis T., Finke J., Becker P. A., 2019, in High Energy Phenomena in Relativistic Outflows VII. p. 75 (arXiv:1912.10260)
  * Massaro et al. (2006) Massaro E., Tramacere A., Perri M., Giommi P., Tosti G., 2006, A&A, 448, 861
  * Meyer et al. (2019) Meyer M., Scargle J. D., Blandford R. D., 2019, ApJ, 877, 39
  * Nalewajko et al. (2012) Nalewajko K., Sikora M., Madejski G. M., Exter K., Szostek A., Szczerba R., Kidger M. R., Lorente R., 2012, ApJ, 760, 69
  * Nesterov et al. (2000) Nesterov N. S., Volvach A. E., Strepka I. D., 2000, Astronomy Letters, 26, 204
  * Nilsson et al. (2009) Nilsson K., Pursimo T., Villforth C., Lindfors E., Takalo L. O., 2009, A&A, 505, 601
  * Pacciani et al. (2014) Pacciani L., Tavecchio F., Donnarumma I., Stamerra A., Carrasco L., Recillas E., Porras A., Uemura M., 2014, ApJ, 790, 45
  * Paliya et al. (2015) Paliya V. S., Sahayanathan S., Stalin C. S., 2015, ApJ, 803, 15
  * Paliya et al. (2017) Paliya V. S., Marcotulli L., Ajello M., Joshi M., Sahayanathan S., Rao A. R., Hartmann D., 2017, The Astrophysical Journal, 851, 33
  * Paltani & Türler (2005) Paltani S., Türler M., 2005, A&A, 435, 811
  * Pittori et al. (2018) Pittori C., et al., 2018, The Astrophysical Journal, 856, 99
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Poole et al. (2008) Poole T. S., et al., 2008, MNRAS, 383, 627
  * Potter & Cotter (2012) Potter W. J., Cotter G., 2012, MNRAS, 423, 756
  * Potter & Cotter (2013a) Potter W. J., Cotter G., 2013a, MNRAS, 429, 1189
  * Potter & Cotter (2013b) Potter W. J., Cotter G., 2013b, MNRAS, 431, 1840
  * Potter & Cotter (2013c) Potter W. J., Cotter G., 2013c, MNRAS, 436, 304
  * Prince et al. (2019) Prince R., Gupta N., Nalewajko K., 2019, ApJ, 883, 137
  * Rani et al. (2013) Rani B., Lott B., Krichbaum T. P., Fuhrmann L., Zensus J. A., 2013, A&A, 557, A71
  * Richards et al. (2011) Richards J. L., et al., 2011, ApJS, 194, 29
  * Roming et al. (2005) Roming P. W. A., et al., 2005, Space Sci. Rev., 120, 95
  * Sahakyan (2020) Sahakyan N., 2020, A&A, 635, A25
  * Saito et al. (2015) Saito S., Stawarz Ł., Tanaka Y. T., Takahashi T., Sikora M., Moderski R., 2015, ApJ, 809, 171
  * Sanchez & Deil (2013) Sanchez D. A., Deil C., 2013, in International Cosmic Ray Conference. p. 2784 (arXiv:1307.4534)
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525
  * Shukla et al. (2018) Shukla A., et al., 2018, ApJ, 854, L26
  * Singh et al. (2020) Singh K. K., Meintjes P. J., Ramamonjisoa F. A., 2020, Ap&SS, 365, 33
  * Smith et al. (2009) Smith P. S., Montiel E., Rightley S., Turner J., Schmidt G. D., Jannuzi B. T., 2009, arXiv e-prints, p. arXiv:0912.3621
  * Spanier & Weidinger (2012) Spanier F., Weidinger M., 2012, in International Journal of Modern Physics Conference Series. pp 293–298, doi:10.1142/S2010194512004734
  * Tavecchio et al. (2012) Tavecchio F., Roncadelli M., Galanti G., Bonnoli G., 2012, Phys. Rev. D, 86, 085036
  * Tchernin et al. (2013) Tchernin C., Aguilar J. A., Neronov A., Montaruli T., 2013, A&A, 555, A70
  * Teraesranta et al. (1998) Teraesranta H., et al., 1998, A&AS, 132, 305
  * Tramacere (2020) Tramacere A., 2020, JetSeT: Numerical modeling and SED fitting tool for relativistic jets (ascl:2009.001)
  * Tramacere et al. (2009) Tramacere A., Giommi P., Perri M., Verrecchia F., Tosti G., 2009, A&A, 501, 879
  * Tramacere et al. (2011) Tramacere A., Massaro E., Taylor A. M., 2011, ApJ, 739, 66
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, PASP, 107, 803
  * Yan et al. (2015) Yan D., Zhang L., Zhang S.-N., 2015, MNRAS, 454, 1310
  * Zacharias et al. (2017) Zacharias M., Böttcher M., Jankowsky F., Lenain J. P., Wagner S. J., Wierzcholska A., 2017, ApJ, 851, 72
  * Zamaninasab et al. (2014) Zamaninasab M., Clausen-Brown E., Savolainen T., Tchekhovskoy A., 2014, Nature, 510, 126
  * Zheng et al. (2013) Zheng Y. G., Zhang L., Huang B. R., Kang S. J., 2013, MNRAS, 431, 2356