# Text Augmentation in a Multi-Task View

Jason Wei${}^{\includegraphics[height=5.25pt]{images/lit.png}}$ Chengyu
Huang${}^{\includegraphics[height=5.25pt]{images/world.png}}$ Shiqi
Xu${}^{\includegraphics[height=5.25pt]{images/lit.png}}$ Soroush
Vosoughi${}^{\Springtree}$
ProtagoLabs
International Monetary Fund
Dartmouth College
<EMAIL_ADDRESS>
###### Abstract

Traditional data augmentation aims to increase the coverage of the input
distribution by generating augmented examples that strongly resemble original
samples in an online fashion where augmented examples dominate training.

In this paper, we propose an alternative perspective—a multi-task view (MTV)
of data augmentation—in which the primary task trains on original examples and
the auxiliary task trains on augmented examples. In MTV data augmentation,
both original and augmented samples are weighted substantively during
training, relaxing the constraint that augmented examples must resemble
original data and thereby allowing us to apply stronger levels of
augmentation.

In empirical experiments using four common data augmentation techniques on
three benchmark text classification datasets, we find that the MTV leads to
higher and more robust performance improvements than traditional augmentation.

##  1 Introduction

Most data augmentation techniques aim to generate augmented examples for
training that are similar to original data. In computer vision, operations
such as flipping, cropping, and color jittering are both widely used and
highly effective—it is self-evident that augmented examples closely resemble
original data, and so we generate augmented data in an online fashion during
each minibatch such that no original, unmodified examples are seen during
training Krizhevsky et al. (2012); Zagoruyko and Komodakis (2016); Huang et
al. (2017).

In language, on the other hand, even slight modifications can cause
significant semantic changes, and so it is not always clear whether augmented
examples resemble original data. Despite this uncertainty, many augmentation
techniques in NLP still generate examples stochastically and ignore original
data Zhang et al. (2015); Sennrich et al. (2016a); Xie et al. (2017); Li et
al. (2017); Kobayashi (2018); Wang et al. (2018). When it is unclear whether
augmented examples resemble original data—as is often the case—is it wise to
neglect the original training data?

Traditional Data Augmentation
---
• Intuition: Increase coverage of input distribution by using augmented
examples for training. • Guideline: Augmented examples should be similar to
original data. • Training: Dominated by augmented examples that are generated
stochastically.
Multi-Task View (MTV) of Data Augmentation
• Intuition: Auxiliary task of classifying augmented examples acts as
regularization for the primary task of classifying original examples. •
Guideline: It might be a good idea for augmented samples to resemble original
data, but they can be anything that boosts performance. • Training: Both
original and augmented data receive substantive weighting during training.
Table 1:  Summary of traditional data augmentation versus MTV data
augmentation.

Our paper questions this practice by proposing to include original data during
training. Specifically, we make two contributions:

  1. 1.

We propose a multi-task view of data augmentation (MTV data augmentation),
which trains on both original and augmented examples and therefore allows us
to relax the constraint that augmented examples must resemble original data.
The MTV facilitates augmentation using a higher strength parameter.

  2. 2.

We show empirically that four common data augmentation techniques provide
higher and more robust performance gains using the MTV compared with
traditional augmentation.

##  2 Traditional Data Augmentation111We closely follow the intuition and
notation of Wang et al. Wang et al. (2018)

Situation. During regular training, the canonical maximum likelihood objective
minimizes the cost of the original training set $J_{\textrm{O}}$:

$J_{\textrm{O}}(\theta)=\mathbb{E}_{x,y\sim\boldsymbol{\hat{\textrm{p}}}(X,Y)}\left[-\log\boldsymbol{\textrm{p}}_{\theta}(y|x)\right],$

where $\boldsymbol{\hat{\textrm{p}}}(X,Y)$ is the empirical distribution of
training pairs $x,y$ and $\boldsymbol{\textrm{p}}_{\theta}(y|x)$ is the
parameterized model that we aim to learn (e.g., a neural network). As
$\boldsymbol{\hat{\textrm{p}}}(X,Y)$ is typically the observed data, it will
likely have some mismatch with the true data distribution
$\boldsymbol{\textrm{p}}(X,Y)$. When the mismatch is dramatic—for instance,
when $\boldsymbol{\hat{\textrm{p}}}(X,Y)$ does not sufficiently cover the
training space—model performance will likely suffer.

Remedy. In practice, we often use data augmentation to mitigate the inadequacy
of $\boldsymbol{\hat{\textrm{p}}}(X,Y)$ by providing additional training data.
We generate an augmented distribution
$\boldsymbol{\textrm{q}}(\hat{X},\hat{Y})$ and now minimize the cost of this
augmented training set $J_{\textrm{aug}}$:

$J_{\textrm{aug}}(\theta)=\mathbb{E}_{x,y\sim\boldsymbol{\textrm{q}}(\hat{X},\hat{Y})}\left[-\log\boldsymbol{\textrm{p}}_{\theta}(y|x)\right].$

As we now optimize solely on $\boldsymbol{\textrm{q}}(\hat{X},\hat{Y})$, our
goal is to find $(\hat{x},\hat{y})$ pairs that are likely to fall in the true
distribution p. Assuming the smoothness of p, similar $(x,y)$ pairs will have
similar probabilities, and therefore if an augmented example is more similar
to an observed example, it is more likely to be sampled under the true
distribution. In other words, good augmented examples resemble the observed
data, and we aim to find them. Conversely, if an augmented example diverges
too far from any observed data, it is likely invalid and thus harmful for
training; we don’t want to train on these examples.

The majority of prior work follows this framework of augmented examples
resembling real data. As popular techniques, semantic noising substitutes
tokens with synonyms Wang and Yang (2015); Zhang et al. (2015); Li et al.
(2017); Pervasive Dropout randomly removes words from the input sequence
Sennrich et al. (2016a); and SwitchOut (for machine translation) replaces some
words in both source and target sentences with other words from their
corresponding vocabularies Wang et al. (2018).

Moreover, most of these techniques perform augmentation on every training
example in an online fashion, implicitly assuming that augmented examples so
closely resemble original data that directly training on original examples is
not even worth considering. As we shall see in the next section, adding in
these original examples during training might actually be a worthwhile idea.

##  3 MTV Data Augmentation

Multi-task optimization jointly trains on a primary task and one or more
auxiliary tasks—the intuition is that requiring an algorithm to also learn an
auxiliary task can act as better regularization than penalizing all complexity
uniformly. Prior work has found that multi-task models work particularly well
when the tasks are similar, but can also improve performance even on unrelated
tasks Paredes et al. (2012); Hajiramezanali et al. (2018).

We propose a multi-task view of data augmentation that has a primary task that
optimizes regular training on original examples and an auxiliary task that
optimizes training on augmented data. This MTV jointly optimizes the primary
and auxiliary task(s) using a weighted cost function so that both original and
augmented data receive substantial weight during training:

$J(\theta)=\gamma_{\textrm{O}}\cdot
J_{\textrm{O}}(\theta)+\gamma_{\textrm{aug}}\cdot J_{\textrm{aug}}(\theta)\ ,$

where $\gamma_{\textrm{O}}$ is the weight of original data and
$\gamma_{\textrm{aug}}$ is the weight of augmented data, and
$\gamma_{\textrm{O}}+\gamma_{\textrm{aug}}=1$. In this context, observe that
vanilla training uses $\gamma_{\textrm{O}}\hskip 1.70717pt$=$\hskip
1.70717pt1$ and $\gamma_{\textrm{aug}}\hskip 1.70717pt$=$\hskip 1.70717pt0$,
and traditional data augmentation uses $\gamma_{\textrm{O}}\hskip
1.70717pt$=$\hskip 1.70717pt0$ and $\gamma_{\textrm{aug}}\hskip
1.70717pt$=$\hskip 1.70717pt1$.

The MTV gives us an important freedom that is not offered by the traditional
data augmentation framework. Since traditional data augmentation only trains
on augmented examples, performance suffers detrimentally when augmented data
differs too much from the true distribution—therefore, most studies aim to
generate augmented examples that resemble original data. MTV data
augmentation, however, jointly trains on both original and augmented data,
thereby allowing us to relax the constraint that original and augmented
examples come from the same distribution. In fact, accepting that the original
and augmented distributions might differ or could even be unrelated—as work in
multi-task learning has done Paredes et al. (2012); Hajiramezanali et al.
(2018); Rai and Daumé (2010)—liberates us to apply stronger levels of data
augmentation, which, as we will demonstrate in the next section, leads to
higher and more robust performance.

##  4 Experiments

This section compares multi-task view augmentation to traditional augmentation
for various datasets and augmentation techniques.

###  4.1 Experimental Setup

Datasets. We conduct experiments on three text classification tasks often used
as benchmarks Kim (2014): (1) Stanford Sentiment Treebank (SST2) Socher et al.
(2013) of movie reviews classified as positive/negative, (2)
subjectivity/objectivity dataset (SUBJ) Pang and Lee (2004), where sentences
are classified as either subjective or objective, and (3) question type
dataset (TREC) Li and Roth (2002), in which questions ask for either a
description, entity, abbreviation, human, location, or number.

Models and Experimental Procedures. For text classification, we use BERT
Devlin et al. (2019) (bert-base-uncased from HuggingFace) to extract features
by averaging the last hidden states of the input tokens. To reduce the number
of model hyperparameters and save computation time, we classify these features
using a linear SVM trained for 1000 epochs.333This setup is not state-of-the-
art but allows for experiments to be performed on CPU. Since training data
size depends on the amount of augmented data, we adjust the number of training
epochs so that all models receive the same number of updates. All experiments
are run for five random seeds. Our baseline models without data augmentation
achieved 84.5%, 93.1%, and 83.9% accuracy respectively on the SST2, SUBJ, and
TREC tasks.

Augmentation Techniques. In this paper, we experiment with four simple and
common data augmentation techniques studied in Wei and Zou (2019): (1) Token
Substitution Zhang et al. (2015) replaces words with WordNet Miller (1995)
synonyms; (2) Pervasive Dropout Sennrich et al. (2016a) applies word-level
dropout; (3) Token Injection Wei and Zou (2019) insert a synonym of a random
token in the sequence into a random position in that sequence; (4) Positional
Shuffling Wei and Zou (2019) randomly chooses two tokens and swaps their
positions. For all four techniques, a parameter $\alpha$ indicates
augmentation strength by dictating how many perturbations are performed. For a
given $\alpha$, we perform $n\hskip 1.70717pt$=$\hskip 1.70717pt\alpha$ $l$
perturbations, where $l$ is the sequence length.

Aug. Technique | MTV | Best $\alpha$ | Avg. Boost | ($\Delta_{\textrm{MTV}}$)
---|---|---|---|---
Token Substitution | ✗ | 0.05 | 1.3% | -
✓ | 0.3 | 2.1% | (+0.8%)
Pervasive Dropout | ✗ | 0.1 | 1.8% | -
✓ | 0.4 | 2.5% | (+0.7%)
Token Injection | ✗ | 0.05 | 0.7% | -
✓ | 0.5 | 2.2% | (+1.5%)
Positional Shuffling | ✗ | 0.05 | 1.4% | -
✓ | 0.4 | 2.5% | (+1.1%)
Table 2:  Average performance boost on three text classification tasks for
four augmentation techniques using the best-performing augmentation strength
from $\alpha\in\\{0.05,0.1,0.2,0.3,0.4,0.5\\}$. Traditional data augmentation
works best at low $\alpha$, whereas MTV data augmentation provides the
strongest performance for high $\alpha$. $\Delta_{\textrm{MTV}}$ indicates
additional boost from using the MTV compared with traditional augmentation.

###  4.2 Stronger augmentation for more gains

Table 2 summarizes results for data augmentation in the MTV using
$\gamma_{\textrm{O}}\hskip 1.70717pt$=$\hskip
1.70717pt\gamma_{\textrm{aug}}\hskip 1.70717pt$=$\hskip 1.70717pt0.5$ compared
with traditional augmentation for the best-performing augmentation strength
from $\alpha\in\\{0.05,0.1,0.2,0.3,$ $0.4,0.5\\}$. In the traditional
framework, pervasive dropout had the strongest performance boost of 1.8% using
$\alpha\hskip 1.70717pt$=$\hskip 1.70717pt0.1$. The MTV, however, allowed for
stronger augmentation (i.e., $\alpha\geq 0.3$) that resulted in all four
techniques to achieving boosts of more than 2.0%.

Perhaps strikingly, token injection and positional shuffling, which are less
intuitive and not as commonly used as token substitution and pervasive
dropout, achieve the strongest gains ($>1.0\%$) from using the MTV. One
potential reason for this is that, compared with token substitution and
pervasive dropout, token injection and positional shuffling are non-
destructive in that they do not remove any of the original words, and so the
nature of examples augmented at high $\alpha$ could be more conducive for the
MTV.

A: $-4$$-2$$0$$2$$4$Boost (%)SUBJ $-4$$-2$$0$$2$$4$SST2 $-6$$-3$$0$$3$$6$TREC
$-4$$-2$$0$$2$$4$AverageMTVTraditional
B: $-4$$-2$$0$$2$$4$Boost (%) $-4$$-2$$0$$2$$4$ $-6$$-3$$0$$3$$6$
$-4$$-2$$0$$2$$4$
C: $-4$$-2$$0$$2$$4$Boost (%) $-4$$-2$$0$$2$$4$ $-6$$-3$$0$$3$$6$
$-4$$-2$$0$$2$$4$
D: $0$$0.2$$0.4$$-4$$-2$$0$$2$$4$Strength $\alpha$Boost (%)
$0$$0.2$$0.4$$-4$$-2$$0$$2$$4$Strength $\alpha$
$0$$0.2$$0.4$$-6$$-3$$0$$3$$6$Strength $\alpha$
$0$$0.2$$0.4$$-4$$-2$$0$$2$$4$Strength $\alpha$

Figure 1:  Performance boosts on three tasks using the traditional and multi-
task view (MTV) frameworks for four data augmentation techniques: Token
Substitution Zhang et al. (2015) (A), Pervasive Dropout Sennrich et al.
(2016a) (B), Token Injection Wei and Zou (2019) (C), Positional Shuffling Wei
and Zou (2019) (D). In the traditional framework, improvements are largest
when augmentation strength $\alpha$ is small, with performance deteriorating
for large $\alpha$. The MTV, on the other hand, jointly optimizes for both
original and augmented data, leveraging higher $\alpha$ to provide higher and
more robust performance gains.

###  4.3 More-robust gains at high $\alpha$

When using data augmentation with high $\alpha$, high levels of noising are
employed and augmented data are therefore more likely to diverge from their
original examples. Figure 1 takes a closer look at how performance is affected
by varying $\alpha$. Whereas traditional augmentation often negatively
affected performance at high $\alpha$, the multi-task view, which jointly
optimizes the original distribution, had robust performance gains at high
augmentation strengths.

Figure 2:  Performance boost (%) with varying augmentation strengths $\alpha$
and weights of original data $\gamma_{\textrm{O}}$ during training.
Traditional data augmentation (yellow solid box in lower left) uses modest
augmentation strength ($\alpha\hskip 1.70717pt$=$\hskip 1.70717pt0.05,1$) with
no original examples for training ($\gamma_{\textrm{O}}\hskip
1.70717pt$=$\hskip 1.70717pt0$). The MTV data augmentation approach (green
dashed box) suggests substantive weighting of original examples (e.g.,
$\gamma_{\textrm{O}}\hskip 1.70717pt$=$\hskip 1.70717pt0.5$) which allows for
much stronger augmentation (e.g., $\alpha\geq 0.3$).

###  4.4 Choosing $\gamma_{\textrm{O}}$ and $\gamma_{\textrm{aug}}$ weighting

As our experiments so far have used the MTV with balanced weighting of
original and augmented data ($\gamma_{\textrm{O}}\hskip 1.70717pt$=$\hskip
1.70717pt\gamma_{\textrm{aug}}\hskip 1.70717pt$=$\hskip 1.70717pt0.5$), in
this section we explore different weightings of $\gamma_{\textrm{O}}$ and
$\gamma_{\textrm{aug}}$. Figure 2 shows these results averaged over all three
datasets and all four augmentation techniques. Traditional data augmentation,
which uses modest augmentation strength (e.g., $\alpha\in\\{0.05,0.1\\}$) and
does not train on original data ($\gamma_{\textrm{O}}\hskip 1.70717pt$=$\hskip
1.70717pt0.0$), achieves reasonable performance gains. As expected, when
stronger augmentations were applied (e.g., $\alpha\geq 0.4$), training with
only augmented data hurts performance. When training on both augmented and
original data, however, performance improved with stronger augmentation and
remained robust for varying augmentation strengths $0.2\leq\alpha\leq 0.5$ and
original data weights $0.3\leq\gamma_{\textrm{O}}\leq 0.7$.

##  5 Further Related Work

Prior work on data augmentation, to our knowledge, generally follows the
traditional data augmentation framework. In addition to the methods mentioned
in $\S\ref{sec:traditional}$, Xie et al. Xie et al. (2017) replaced words with
samples from the unigram frequency distribution; Yu et al. Yu et al. (2018)
translated English sentences to French and back to English (backtranslation);
and Kobayashi Kobayashi (2018) replaced words with other words based on a
language model. All these methods could potentially be formulated in the MTV.

Some prior work has also drawn connections between seeing data augmentation as
multiple tasks. Similar to how we optimize augmented data as a separate task,
Meyerson and Miikkulainen (2018) created fake tasks by using multiple distinct
decoders to train a shared structure to solve the same problem in different
ways. In machine translation, Sennrich et al. (2016b) used monolingual
training examples as parallel examples with an empty source side, noting that
their setup could be seen as multi-task learning with the tasks as translation
with known sources and language modeling with unknown sources. Compared with
these papers that create multiple tasks in very specialized scenarios, the
multi-task view that we have presented here can be used for any type of text
data augmentation.

To be clear, our study is not the first to mix original and augmented data in
training. For instance, Wang and Yang (2015) use a ratio of 1:5 original to
augmented examples, but this weight of original data is much smaller than the
$0.3\leq\gamma_{\textrm{O}}\leq 0.7$ that we advocate for. Sennrich et al.
(2016b) also include original data when training with back-translation
augmentation, but the given ratios of original and augmented data they use
appear to dictated by the speed of their back-translation models rather than
an intentionally-motivated design choice. We see our work as the first to
explicitly formulate the MTV, advocate for a joint optimization function, and
comprehensively explore its implications on common text augmentation
techniques.

As a limitation, our study has focused on label-preserving augmentation
techniques, and our line of reasoning may not apply when augmentation
techniques intentionally change the label. Moreover, we have only studied text
classification with simple models using task-agnostic augmentation techniques.
Future work in this direction could experiment with larger-scale models or
study task-specific augmentation.

##  6 Conclusions

We have proposed a multi-task view that gives both original and augmented
examples substantial weight during training, contrasting prior work that
performs stochastic data augmentation and ignores original training data. For
four common augmentation techniques, we found experimentally that this
alternative view allows for stronger levels of augmentation, which in turn
leads to better and more robust performance than traditional augmentation. We
hope our paper inspires future work using text data augmentation to think more
explicitly about how much augmented examples resemble original data and
consider substantive weighting of original data when using data augmentation
to improve model performance.

To close, we leave the enthusiastic reader with one last thought. Most
existing text data augmentation techniques have obediently followed the
paradigm from computer vision of generating augmented examples that are
similar to the original data. Who’s to say that’s how data augmentation ought
to work in NLP? In this paper, we’ve shown how to search for relative freedom
from this constraint, simply by taking a different view of the underlying
assumptions. Now, a bigger question arises on the horizon—what new text
augmentation techniques are unlocked when augmented data are not forced to
resemble the original?

## References

  * 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.
  * Hajiramezanali et al. (2018) Ehsan Hajiramezanali, Siamak Zamani Dadaneh, Alireza Karbalayghareh, Mingyuan Zhou, and Xiaoning Qian. 2018.  Bayesian multi-domain learning for cancer subtype discovery from next-generation sequencing count data.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NIPS’18, page 9133–9142, Red Hook, NY, USA. Curran Associates Inc.
  * Huang et al. (2017) Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger. 2017.  Densely connected convolutional networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 4700–4708.
  * Kim (2014) Yoon Kim. 2014.  Convolutional neural networks for sentence classification.  In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1746–1751, Doha, Qatar. Association for Computational Linguistics.
  * Kobayashi (2018) Sosuke Kobayashi. 2018.  Contextual augmentation: Data augmentation by words with paradigmatic relations.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 2 (Short Papers)_ , pages 452–457, New Orleans, Louisiana. Association for Computational Linguistics.
  * Krizhevsky et al. (2012) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. 2012.  Imagenet classification with deep convolutional neural networks.  In _Advances in neural information processing systems_ , pages 1097–1105.
  * Li and Roth (2002) Xin Li and Dan Roth. 2002.  Learning question classifiers.  In _Proceedings of the 19th International Conference on Computational Linguistics - Volume 1_ , COLING ’02, pages 1–7, Stroudsburg, PA, USA. Association for Computational Linguistics.
  * Li et al. (2017) Yitong Li, Trevor Cohn, and Timothy Baldwin. 2017.  Robust training under linguistic adversity.  In _Proceedings of the 15th Conference of the European Chapter of the Association for Computational Linguistics: Volume 2, Short Papers_ , pages 21–27, Valencia, Spain. Association for Computational Linguistics.
  * Meyerson and Miikkulainen (2018) Elliot Meyerson and Risto Miikkulainen. 2018.  Pseudo-task augmentation: From deep multitask learning to intratask sharing-and back.  In _Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 3508–3517. PMLR.
  * Miller (1995) George A Miller. 1995.  Wordnet: a lexical database for english.  _Communications of the ACM_ , 38(11):39–41.
  * Pang and Lee (2004) Bo Pang and Lillian Lee. 2004.  A sentimental education: Sentiment analysis using subjectivity summarization based on minimum cuts.  In _Proceedings of the 42Nd Annual Meeting on Association for Computational Linguistics_ , ACL ’04, Stroudsburg, PA, USA. Association for Computational Linguistics.
  * Paredes et al. (2012) Bernardino Romera Paredes, Andreas Argyriou, Nadia Berthouze, and Massimiliano Pontil. 2012.  Exploiting unrelated tasks in multi-task learning.  In _Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics_ , volume 22 of _Proceedings of Machine Learning Research_ , pages 951–959, La Palma, Canary Islands. PMLR.
  * Rai and Daumé (2010) Piyush Rai and Hal Daumé. 2010.  Infinite predictor subspace models for multitask learning.  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , volume 9 of _Proceedings of Machine Learning Research_ , pages 613–620, Chia Laguna Resort, Sardinia, Italy. PMLR.
  * Sennrich et al. (2016a) Rico Sennrich, Barry Haddow, and Alexandra Birch. 2016a.  Edinburgh neural machine translation systems for WMT 16.  In _Proceedings of the First Conference on Machine Translation: Volume 2, Shared Task Papers_ , pages 371–376, Berlin, Germany. Association for Computational Linguistics.
  * Sennrich et al. (2016b) Rico Sennrich, Barry Haddow, and Alexandra Birch. 2016b.  Improving neural machine translation models with monolingual data.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 86–96, Berlin, Germany. Association for Computational Linguistics.
  * Socher et al. (2013) Richard Socher, John Bauer, Christopher D. Manning, and Andrew Y. Ng. 2013.  Parsing with compositional vector grammars.  In _Proceedings of the 51st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 455–465, Sofia, Bulgaria. Association for Computational Linguistics.
  * Wang and Yang (2015) William Yang Wang and Diyi Yang. 2015.  That’s so annoying!!!: A lexical and frame-semantic embedding based data augmentation approach to automatic categorization of annoying behaviors using #petpeeve tweets.  In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing_ , pages 2557–2563, Lisbon, Portugal. Association for Computational Linguistics.
  * Wang et al. (2018) Xinyi Wang, Hieu Pham, Zihang Dai, and Graham Neubig. 2018.  SwitchOut: an efficient data augmentation algorithm for neural machine translation.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 856–861, Brussels, Belgium. Association for Computational Linguistics.
  * Wei and Zou (2019) Jason Wei and Kai Zou. 2019.  EDA: Easy data augmentation techniques for boosting performance on text classification tasks.  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 6382–6388, Hong Kong, China. Association for Computational Linguistics.
  * Xie et al. (2017) Ziang Xie, Sida I. Wang, Jiwei Li, Daniel Lévy, Aiming Nie, Dan Jurafsky, and Andrew Y. Ng. 2017.  Data noising as smoothing in neural network language models.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net.
  * Yu et al. (2018) Adams Wei Yu, David Dohan, Minh-Thang Luong, Rui Zhao, Kai Chen, Mohammad Norouzi, and Quoc V. Le. 2018.  Qanet: Combining local convolution with global self-attention for reading comprehension.  In _6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings_. OpenReview.net.
  * Zagoruyko and Komodakis (2016) Sergey Zagoruyko and Nikos Komodakis. 2016.  Wide residual networks.  In _Proceedings of the British Machine Vision Conference (BMVC)_ , pages 87.1–87.12. BMVA Press.
  * Zhang et al. (2015) Xiang Zhang, Junbo Zhao, and Yann LeCun. 2015.  Character-level convolutional networks for text classification.  In _Proceedings of the 28th International Conference on Neural Information Processing Systems - Volume 1_ , NIPS’15, page 649–657, Cambridge, MA, USA. MIT Press.


11institutetext: Facultad de Ciencias Astronómicas y Geofísicas, Universidad
Nacional de La Plata, Paseo del Bosque, B1900FWA La Plata, Argentina
22institutetext: Instituto Argentino de Radioastronomía (CCT La Plata,
CONICET; CICPBA; UNLP), C.C.5, (1894) Villa Elisa, Buenos Aires, Argentina
33institutetext: Physics and Astronomy Department Galileo Galilei, University
of Padova, Vicolo dell’Osservatorio 3, I-35122, Padova, Italy 44institutetext:
INFN - Padova, Via Marzolo 8, I–35131 Padova, Italy 55institutetext:
Departamento de Física (EPS), Universidad de Jaén, Campus Las Lagunillas s/n,
A3, 23071 Jaén, Spain

#  NuSTAR and XMM-Newton observations of the binary 4FGL J1405.1-6119

A $\gamma$-ray emitting microquasar?

Enzo A. Saavedra 11 Federico A. Fogantini 22 Gastón J. Escobar 3344 Gustavo E.
Romero 1122
Jorge A. Combi 112255 Estefania Marcel 11

(Received; accepted)

###### Abstract

Context. 4FGL J1405.1-6119 is a high-mass $\gamma$-ray-emitting binary that
has been studied at several wavelengths. The nature of this type of binary is
still under debate, with three possible scenarios usually invoked to explain
the origin of the $\gamma$-ray emission: collisions between the winds of a
rapidly rotating neutron star and its companion, collisions between the winds
of two massive stars, and nonthermal emission from the jet of a microquasar.

Aims. We analyzed two pairs of simultaneous NuSTAR and XMM-Newton observations
to investigate the origin of the radio, X-ray, and $\gamma$-ray emissions.

Methods. We extracted light curves between 0.5 and 78 keV from two different
epochs, wich we call Epoch 1 and Epoch 2. We then extracted and analyzed the
associated spectra to gain insight into the characteristics of the emission in
each epoch. To explain these observations, along with the overall spectral
energy distribution, we developed a model of a microquasar jet. This allowed
us to make some inferences about the origin of the observed emission and to
discuss the nature of the system.

Results. A power-law model combined with the inclusion of a blackbody
accurately characterizes the X-ray spectrum. The power-law index
($E^{-\Gamma}$) was found to be $\sim$ 1.7 for Epoch 1 and $\sim$ 1.4 for
Epoch 2. Furthermore, the associated blackbody temperature was $\sim$ 1 keV
and with a modeled emitting region of size ${\lesssim}\leavevmode\nobreak\ 16$
km. The scenario we propose to explain the observations involves a parabolic,
mildly relativistic, lepto-hadronic jet. This jet has a compact acceleration
region that injects a hard spectrum of relativistic particles. The dominant
nonthermal emission processes include synchrotron radiation of electrons,
inverse Compton scattering of photons from the stellar radiation field, and
the decay of neutral pions resulting from inelastic proton-proton collisions
within the bulk matter of the jet. These estimates are in accordance with the
values of a super-Eddington lepto-hadronic jet scenario. The compact object
could be either a black hole or a neutron star with a weak magnetic field.
Most of the X-ray emission from the disk could be absorbed by the dense wind
that is ejected from the same disk.

Conclusions. We conclude that the binary 4FGL J1405.1-6119 could be a
supercritical microquasar similar to SS433.

###### Key Words.:

X-rays: binaries – gamma-rays: stars – stars: individual (CXOU
J053600.0-673507, 4FGL J1405.1-6119)

##  1 Introduction

Binary sources containing neutron stars (NSs) or black holes (BHs) dominate
the Galactic X-ray emission above 2 keV (see, e.g., Grimm et al., 2002). These
systems are called X-ray binaries and are usually divided into two major
classes, high-mass X-ray binaries and low-mass X-ray binaries, according to
the mass of the donor star (mass $\gtrsim 8\leavevmode\nobreak\ M_{\odot}$ for
the former and $\lesssim 8\leavevmode\nobreak\ M_{\odot}$ for the latter).

Within the high-mass X-ray binary class, three types of non-transient systems
can emit $\gamma$-ray radiation (Chernyakova et al., 2019; Chernyakova &
Malyshev, 2020):

1\. Colliding-wind binaries involve the interaction of two massive stars,
whose nonrelativistic winds collide and produce $\gamma$-ray emission.
Prominent examples of such systems are $η$-Carinae, WR11, and Apep. In
colliding-wind binaries, intense shocks occur in the wind collision region,
leading to the formation of a very hot plasma (¿ $10^{6}$ K). In addition,
these systems have the ability to accelerate relativistic particles (Eichler &
Usov, 1993; Benaglia & Romero, 2003), which classifies them as particle-
accelerating colliding-wind binaries (e.g., De Becker & Raucq, 2013; del
Palacio et al., 2023, and references therein).

Figure 1: Background-corrected light curves of 4FGL J1405.1-6119 observed by
NuSTAR (FPMA+B; red) and XMM-Newton (EPIC pn+MOS; blue) with a binning of
350s, starting at 58711.705573 MJD. The first section of the light curve
corresponds to the observation from August 18 (Epoch 1), while the second
light curve is associated with the observation from August 25 (Epoch 2).

2\. Gamma-ray binaries are characterized by the presence of a young,
magnetized, and rapidly rotating NS that emits a relativistic wind. This wind
collides with the nonrelativistic wind of the companion OB star, and through
this interaction they emit nonthermal radiation that peaks at high energies
($E>100$ MeV). These systems are thought to represent a short phase in the
evolution of massive binaries, which comes after the birth of the compact
object (CO) and is followed by the X-ray binary phase, in which the CO
accretes matter from its companion (e.g., Dubus et al., 2017; Saavedra et al.,
2023, and references therein). In the latter phase, the nonthermal emission of
the system peaks in X-rays. Examples of $\gamma$-ray binaries are LS
I+61$^∘$303, PSR J1259-63, PSR J2032+4127, LS 5039, 1FGL J1018.6-5856, and
HESS J0632+057.

3\. Microquasars (MQs) differ from the previous two categories in that their
emission of $\gamma$-rays does not come from wind collisions. Instead, it
comes from the jets ejected by the CO (BH or NS) and their interaction with
the environment. Examples of MQs are Cyg X-1, Cyg X-3, and SS433

In total, there are nine known non-transient $\gamma$-ray-emitting binaries in
our Galaxy (e.g., Corbet et al., 2019; Dubus et al., 2017; Chernyakova &
Malyshev, 2020, and references therein). Five of these systems contain an
O-type star: LS 5039, LMC P3, 1FGL J1018.6-5856, HESS J1832-093, and 4FGL
J1405.1-6119. The remaining four contain a Be star: HESS J0632+057, LS
I+61$^∘$303, PSR B1259-63, and PSR J2032+4127. Radio pulsations were observed
in three of these systems – PSR J2032+4127 (Camilo et al. 2009), PSR B1259-63
(Moldón et al. 2014), and LS I+61$^∘$303 (Weng et al. 2022) – evidence in
favor of the presence of a NS. Many more galactic $\gamma$-ray binaries are
expected to exist. We recall that the nature of many of the $\gamma$-ray
sources detected by telescopes such as Fermi-LAT is not known, and some of the
sources may be $\gamma$-ray binaries. In fact, the number of detectable
Galactic $\gamma$-ray binaries is estimated to be $\sim 100$ (Dubus et al.,
2017).

4FGL J1405.1-6119 (also known as 3FGL J1405.4-6119) is a high-mass
$\gamma$-ray-emitting binary first studied by Corbet et al. (2019). Using
Fermi and Swift/XRT data, Corbet et al. (2019) found a strong modulation of
$13.713\leavevmode\nobreak\ \pm\leavevmode\nobreak\ 0.002$ days associated
with the system’s orbital period. The companion is classified as an O6.5 III
star with a mass of about $25-35$ M⊙ (Mahy et al., 2015). The absence of
partial and total eclipses suggests that this system has a low inclination
($60\leavevmode\nobreak\ \degr$).

To explain the origin of its $\gamma$-ray emission, Corbet et al. (2019)
proposed that 4FGL J1405.1-6119 may be a $\gamma$-ray binary. This hypothesis
draws on analogies with other similar systems, such as 1FGL J1018.6-5856 and
LMC P3. Xingxing et al. (2020) modeled the GeV time behavior of 4FGL
J1405.1-6119 in the $\gamma$-ray binary scenario, assuming a binary consisting
of a young pulsar and an O-type main sequence star. Conversely, the radio (5.5
GHz and 9 GHz) and X-ray (0.2–10 keV) luminosities show a positive correlation
(Corbet et al., 2019), as expected in MQs (Falcke et al., 2004).

We were able to perform a detailed temporal and spectral study of 4FGL
J1405.1-6119 over a broad energy range by analyzing simultaneous XMM-Newton
and NuSTAR observations. In this paper we present our results and conclusions
about this source. In Sect. 2 we present the X-ray observations and the
corresponding tools used for their analysis. Our main results are presented in
Sect. 3. A jet model for the source is introduced and discussed in Sect. 4. We
discuss our results and present our conclusions in Sect. 5.

##  2 Observation and Data Analysis

###  2.1 XMM-Newton data

The XMM-Newton observatory is equipped with an optical instrument and two
X-ray instruments: the Optical Monitor, which is mounted on the mirror support
platform and provides coverage between 170 nm and 650 nm of the central 17
arcminute square region; the European Photon Imaging Camera (EPIC); and the
Reflecting Grating Spectrometers (RGS). The EPIC instrument comprises three
detectors – the pn camera (Strüder et al., 2001) and two MOS cameras (Turner
et al., 2001) – which are most sensitive in the $0.3-10$ keV energy range. The
RGS instrument comprises two high-resolution spectrographic detectors
sensitive in the energy range $0.3-2$ keV.

XMM-Newton observed 4FGL J1405.1-6119 on August 17, 2019, with an exposure
time of 32 ks (ObsID 0852020101) and on August 24, 2019, with an exposure time
of 44 ks (ObsID 0852020201). In both observations, the MOS cameras were in
large window mode, and the pn camera was in timing mode.

We reduced the XMM-Newton data using Science Analysis System (SAS) v20.0 and
the latest available calibration in early 2022. To process the observation
data files, we used the EPPROC and EMPROC tasks. We selected circular regions
with radii of 18 arcsec and 36 arcsec for the source and the background,
respectively, with the latter away from any source contamination. We then
filtered the raw events lists, removing the high-energy, single-pattern
particle backgrounds and thus creating cleaned event lists for each camera and
observation. The resulting exposure times after background filtering are 19 ks
(59% of total) for the first observation, and 44 ks (100% of total) for the
second observation. We studied the presence of pile-up with the EPATPLOT task
and did not find any deviation of the data from the expected models; thus, no
excision radii were applied. We barycentered each cleaned event list with the
barycen task in order to perform precise timing studies. EPIC light curves
were summed using the LCMATH task, with proper scaling factors for the
different source photons collecting areas. We extracted and grouped spectra
with a minimum of 25 counts per bin in the $0.5-12$ keV energy band.

###  2.2 NuSTAR data

The NuSTAR X-ray observatory was sent into orbit in the year 2012 and is
notable for its exceptional sensitivity at hard X-rays. It is equipped with
two X-ray grazing incidence telescopes, designated FPMA and FPMB, which are
arranged in parallel and contain 2x2 solid-state CdZnTe detectors each. NuSTAR
can operate in the energy range of 3-79 keV and can achieve an angular
resolution of 18 arcsec as reported in Harrison et al. (2013).

NuSTAR observed 4FGL J1405.1-6119 on August 16, 2019 (58711.7031 MJD - ObsID
30502015002), with an exposure of time of ${\sim}61$ ks and on August 24, 2019
(58719.3511 MJD - ObsID 30502015004) with an exposure time of ${\sim}86$ ks.
We processed NuSTAR data using NuSTARDAS-v.2.1.2 from HEASoft v.6.30 and CALDB
(v.20211020) calibration files. We took source events that accumulated within
a circular region of 85 arcsec around the focal point. The chosen radius
encloses $\sim 90$% of the point spread function. We took a circular source-
free region with a radius of 160 arcsec to obtain the background events within
the same CCD.

Figure 2: NuSTAR and XMM-Newton folded light curve using 48 phase bins, an
orbital period of 13.713 days, and with 56498.7 MJD as the reference epoch
(Corbet et al., 2019). A sine function fit is shown in orange (see the main
text for details). The observed orbital modulation is similar to that shown by
Corbet et al. (2019) using Swift/XRT data.

We used the nupipeline task to create level 2 data products, with SAA
parameters saacalc=1 saamode=strict tentacle=yes to filter for high-energy
particle background, obtained from SAA filtering
reports111https://nustarsoc.caltech.edu/NuSTAR_Public/NuSTAROperationSite/SAA_Filtering/SAA_Filter.php.
We extracted light curves and spectra with the nuproducts task. We obtained
the barycenter-corrected light curves using barycorr task with
nuCclock20100101v136 clock correction file. We used celestial coordinates
${\alpha}=211.2472^{\circ}$ and ${\delta}=+61.3234^{\circ}$ for the
barycentric correction with JPL-DE200 as the Solar System ephemeris. Finally,
we subtracted the background from each detector’s light curve. Then we used
the LCMATH task to create FPMA+B light curves. We extracted and re-binned
spectra with a minimum of 25 counts per bin in the $3-78$ keV energy band.

We used the XSPEC v12.12.1 package (Arnaud, 1996) to model XMM-Newton and
NuSTAR spectra, with parameters uncertainties reported at the 90% confidence
level.

##  3 Results

###  3.1 Analysis of the light curves

Figure 1 shows the background-corrected light curves obtained from XMM-Newton
(0.5-10 keV, top panel) and NuSTAR (3–78 keV, bottom panel) missions, with a
bin time of 350 s. The observation conducted on August 17 is labeled as Epoch
1, while the observation on August 24 is labeled as Epoch 2. The long-term
exposures of NuSTAR, both on the order of ${\sim}$1.5 d, show a hard X-ray
flux modulation of ${\sim}1$ d seen in both epochs. The shorter but continuous
exposures of XMM-Newton do not capture this behavior. Instead, it captures
very short (on the order of some ks) changes in soft X-ray flux, as seen in
Epoch 1.

Figure 2 shows the orbital flux modulation associated with each observation
and mission. Epoch 1 occurred within the orbital phases of 0.93–1.08, while
Epoch 2 occurred within the orbital phases of 0.37–0.48. The observed orbital
behavior of XMM-Newton and NuSTAR data is very similar to that reported by
Corbet et al. (2019) using Swift/XRT data.

We used a sinusoidal model with Gaussian measurement errors to visualize the
orbital modulation through $\sim$104 simulations (Buchner, 2021).
Specifically, the following was used:

$y=A\,\sin\left(2\pi\left[\frac{t}{P}+t_{0}\right]\right)+B+\epsilon,$ (1)

where $\epsilon$ $\sim$ Normal(0, $\sigma$), that is, a normal distribution
with a mean of zero and a standard deviation of $\sigma$. We obtain the
following values: $A=0.56\pm 0.3$, $P=1.002\pm 0.002$, $t_{0}=0.623\pm 0.001$
and $B=1.13\pm 0.01$. The fitted model is shown in Fig. 2.

From phase $\sim$ 0.93, the flux starts to decay, and from phase $\sim$ 0.37
the source has the maximum local emission. This modulation is anticorrelated
with the modulation obtained from the Fermi data in the energy range 200 MeV
to 500 GeV (Corbet et al., 2019).

We employed spectral timing routines provided by Stingray software
(Huppenkothen et al., 2019) to conduct a comprehensive search for any
potential pulsation linked to the X-ray source. Light curves from both
telescopes across various energy ranges do not show any significant pulsations
above noise on the 0.1–100 mHz frequency range.

Figure 3: Spectral modeling results corresponding to Epoch 1 (left column) and
Epoch 2 (right column) derived from simultaneous XMM-Newton and NuSTAR data
(top panels). $\chi^{2}$ residuals correspond to an absorbed power-law (middle
panel) and absorbed power-law with a black-body component (bottom panel).

###  3.2 Spectral analysis

We simultaneously modeled source and background spectra extracted from all
five detectors (pn, MOS1, MOS2, FPMA and FPMB). We introduced calibration
constants in our models to account for disparities in effective areas between
instruments. The pn constant was fixed to unity, while the remaining
calibration constants were permitted to vary: $C_{\rm MOS1}$ = 1.10$\pm$0.16,
$C_{\rm MOS2}$ = 1.05$\pm$0.15, $C_{\rm FPMA}$=1.40$\pm$0.16, and $C_{\rm
FPMB}$=1.35$\pm$0.16. Each epoch was modeled separately.

The NuSTAR background showed significant activity at energies above 20 keV. As
a result, we limited our analysis to the energy range 3–20 keV for both
epochs. In the case of the XMM-Newton background, it was significant at
energies below 2 keV. Therefore, we focused our analysis on the 2–10 keV
energy range for both epochs. Consequently, the total spectrum used for each
epoch included the 2–20 keV energy range.

The interstellar absorption was modeled using the Tübingen-Boulder model
(tbabs), with solar abundances defined by Wilms et al. (2000) and effective
cross sections of Verner et al. (1996). Several continuum models were used to
fit the time-averaged spectra, including combinations of power-law variants
such as powerlaw, highecut+powerlaw, cutoffpl, and bknpow, and thermal models
such as apec, bbodyrad, and diskbb. After several fits, we selected the best
fits: a power law (tbabs powerlaw, hereafter Model 1) and a combination of a
power law and a blackbody (tbabs(powerlaw + bbodyrad), hereafter Model 2). No
emission or absorption lines above the continuum were observed in the spectra
of all epochs.

Model 1 yielded a $\chi^{2}$ of 221 with 191 degrees of freedom
($\chi^{2}_{\nu}=1.16$) for Epoch 1, while Epoch 2 resulted in
$184/194{\sim}0.95$. On the other hand, if we apply Model 2, we obtain
$216/189{\sim}1.14$ for Epoch 1 and $181/192{\sim}0.94$ for Epoch 2.

| Epoch 1 | Epoch 2
---|---|---
Parameters | Model 1 | Model 2 | Model 1 | Model 2
${\rm N_{H}}$ [1022 cm-2] | 7.0${}_{-1}^{+1.5}$ | 5$\pm$3 | 8.0${}_{-1}^{+1.5}$ | $9\pm 3$
$\Gamma$ | 1.6$\pm$0.2 | 1.15$\pm$0.10 | 1.65$\pm$0.15 | 1.2$\pm$0.5
Norm${}_{\rm\leavevmode\nobreak\ PL}$ [$10^{-4}$] | 1.7${}^{+0.5}_{-0.4}$ | 0.5${}^{+1.5}_{-0.5}$ | 1.00$\pm$0.04 | 0.02$\pm$0.01
kT [keV] | - | 1.0$\pm$0.3 | - | 0.8$\pm$0.2
$R_{\rm bb}$ [km] | - | 3${}_{-2}^{+8}$ | - | 6${}_{-3}^{+10}$
$L_{x}$ [1033 erg s-1] | $7.5\pm 0.3$ | $4.7\pm 0.3$
$\chi^{2}/{\rm dof}$ | 221/191 | 216/189 | 184/194 | 181/192

222† Parameter is fixed. Parameter confidence intervals are reported within
90% significance. Unabsorbed luminosities, $L_{x}$, are calculated between 2
and 20 keV assuming a distance of 7.7 kpc.

Table 1: Best-fit parameters of XMM-Newton +NuSTAR time-averaged spectral
modeling of 4FGL J1405.1-6119 with const*tbabs(powerlaw) (Model 1) and
const*tbabs(powerlaw+bbodyrad) (Model 2).

To assess the significance level of the blackbody components, we performed
spectral simulations using the same observational data using fake-it command
from XSPEC. Each fake spectra was constructed from arrays of randomly sampled
parameters using the simpars command. By generating a cumulative distribution
function of F-values from the simulated data, we determined the minimum
significance level of detection by comparing it with the F-value derived from
the real data, as described in Hurkett et al. (2008).

Each F-value, for both the real and simulated data, is calculated as
$F=(\nu_{0}/\delta\nu)\times(\delta\chi^{2}/\chi^{2}_{1})$, where the
subscripts ${0,1}$ correspond to the null hypothesis (powerlaw) and the tested
hypothesis (powerlaw+bbody). Each hypothesis is associated with a total
$\chi^{2}$ and $\nu$ degrees of freedom. The significance level is determined
by computing the corresponding $p$ value, which is the number of simulated
spectra with $F$ values greater than the $F$ value derived from fitting the
actual data. The uncertainty of this quantity can be calculated using
$\sqrt{p(1-p)/N{\rm s}}$, where $N{\rm s}$ is the total number of simulated
spectra.

We ran $\sim$ 105 simulations for both epochs and found that the blackbody
component is significant at $\sim$ 2.6$\sigma$ level for Epoch 1 and $\sim$
2.8$\sigma$ for Epoch 2.

From our analysis we conclude that the spectrum was nonthermal dominated
during the observed period, possibly of synchrotron origin. This implies the
presence of particles with TeV energies for the typical magnetic field
strengths in this type of system. In the next section we explore this
hypothesis in more detail.

In Fig. 3 we present the time averaged spectra and residuals associated with
Model 1 and Model 2 of Epoch 1 (left panel) and Epoch 2 (right panel), while
the corresponding best-fit parameters and uncertainties are detailed in Table
1. The powerlaw normalization component is expressed in units of photons keV-1
cm-2 s-1 at 1 keV. The bbodyrad normalization is equal to $R_{\rm
km}^{2}/D_{10}^{2}$, where $R_{\rm km}$ is the source radius in km and
$D_{10}$ is the distance to the blackbody source in units of 10 kpc. At a
distance of 6.7 kpc, the size of the emitting region ranged from 1.5 to 5.4 km
during Epoch 1 and from 2.9 to 8.9 km during Epoch 2. Alternatively, at a
distance of 8.7 kpc, the size of the emitting region ranged from 2 to 7 km for
Epoch 1 and from 4 to 11.6 km for Epoch 2. In Table 1, we report the values
assuming a distance of 7.7 kpc.

To compute the unabsorbed flux, we used the convolution model cflux. Assuming
a distance of 7.7 kpc (Corbet et al., 2019), the unabsorbed 2–20 keV
luminosity ranges between $7.2-7.8\leavevmode\nobreak\ \times 10^{33}$ erg s-1
for Epoch 1 and $4.4-5\leavevmode\nobreak\ \times 10^{33}$ erg s-1 for Epoch
2.

##  4 Model

###  4.1 Jet nonthermal radiation

We adopted the hypothesis that a jet is present in the $\gamma$-ray-emitting
binary 4FGL J1405.1-6119, and tried to evaluate whether it can adequately
explain the nonthermal spectral energy distribution (SED) of the system. The
radiative jet model used follows that presented in detail in Escobar et al.
(2022), which in turn is based on Romero & Vila (2008). In the following, the
model setup is outlined.

The scenario consists of a CO accreting material from the companion star with
an accretion power $L_{\mathrm{acc}}$, which can be expressed in terms of the
Eddington luminosity as

$\displaystyle L_{\mathrm{acc}}=q\leavevmode\nobreak\ L_{\mathrm{Edd}}\approx
q\leavevmode\nobreak\ 1.3\times
10^{38}\left(\frac{M}{M_{\odot}}\right)\leavevmode\nobreak\
\mathrm{erg}\leavevmode\nobreak\ \mathrm{s}^{-1},$ (2)

where $q$ is a constant that represents the accretion efficiency in terms of
the Eddington limit. Coupled with the inner accretion disk we assume the
presence of a lepto-hadronic jet of kinetic luminosity, $L_{\mathrm{jet}}$.
This jet power relates to the accretion power through

$\displaystyle L_{\mathrm{jet}}=q_{\mathrm{jet}}L_{\mathrm{acc}},$ (3)

where $q_{\mathrm{jet}}$ is another constant indicating the fraction of the
accretion power that is transferred to the jet. The parameters $q$ and
$q_{\mathrm{jet}}$ define the accretion regime of the MQ. To compute the SED
and the normalizations with the emission power, we chose to use
$L_{\mathrm{jet}}$ directly as the parameter instead of a combination of $q$
and $q_{\mathrm{jet}}$; we defer discussion of the interpretation of the
accretion regime to Sect. 5.

The jet propagates with a bulk velocity $v_{\mathrm{jet}}$, corresponding to a
bulk Lorentz factor $\Gamma_{\mathrm{jet}}$. A fraction $q_{\mathrm{rel}}$ of
this power is converted into relativistic particles by an acceleration
mechanism. The relativistic proton and electron luminosities, $L_{\mathrm{p}}$
and $L_{\mathrm{e}}$, are distributed according to the power ratio
$a=L_{\mathrm{p}}/L_{\mathrm{e}}$.

We use cylindrical coordinates, with the coordinate $z$ along the jet axis and
the origin at the CO. The jet is started at a distance $z_{0}$ from the CO.
The region in which the particles are accelerated extends from
$z_{\mathrm{min}}$ to $z_{\mathrm{max}}$, and its shape is described by

$\displaystyle r(z)=r_{0}\left(\frac{z}{z_{0}}\right)^{\varepsilon},$ (4)

Figure 4: Acceleration, escape, and cooling rates of electrons (top row) and
protons (bottom row). Each column shows the aforementioned rates calculated at
$z_{\rm min}$ (left column), at the logarithmic midpoint (middle column), and
at $z_{\rm max}$ (right column). Figure 5: Nonthermal SED derived from our jet
model. We considered the following radiative processes: for protons,
synchrotron emission and proton-proton collisions with the cold protons of the
jet bulk and the companion’s wind; for electrons, synchrotron emission and
inverse Compton scattering off the radiation field of the companion. The
figure also shows luminosities derived from XMM-Newton +NuSTAR data (our work)
as well as ATCA and Fermi data taken from Corbet et al. (2019) and the 4FGL-
DR2 catalog (Abdollahi et al., 2020; Ballet et al., 2020), respectively. All
the references are in the figure.

where $r_{0}$ is the radius of the jet at its base and $0<\varepsilon\leq 1$
describes its geometry. We note that the degree of collimation of the jet
increases with decreasing values of $\varepsilon$, with the extreme value
$\varepsilon=1$ representing a conical shape. The magnetic field along the
jet, $B$, decreases with $z$ following a power law of index $m$, namely,
$B(z)=B_{0}(z/z_{0})^{-m}$, where $1\leq m\leq 2$ (e.g., Krolik, 1999). The
value of $B_{0}$ is obtained by assuming equipartition between magnetic and
kinetic energy at the jet base.

We parameterize the injection function of energy to relativistic particles as
a power law with an exponential cutoff,

$\displaystyle
Q_{i}=Q_{i,0}E_{i}^{-p}\exp\left(\frac{E}{E_{i,\mathrm{max}}}\right),$ (5)

where $i=\mathrm{e},\mathrm{p}$ accounts for electrons and protons,
respectively; $Q_{i,0}$ is obtained normalizing the injection function with
the total power of each particle population, $L_{i}$; and $E_{i,\mathrm{max}}$
is the maximum reachable energy, achieved when the acceleration rate equals
that of energy losses.

Relativistic particles are accelerated at a rate of

$\displaystyle t^{-1}_{\mathrm{acc}}=\frac{\eta ceB}{E_{i}},$ (6)

where $\eta\leq 1$ is the acceleration efficiency, $c$ is the speed of light,
and $e$ is the elementary electric charge. On the other hand, relativistic
particles lose energy via both radiative and non-radiative mechanisms. The
latter are adiabatic and escape losses for both proton and electron
populations. Regarding the radiative mechanisms, we consider synchrotron
emission for both types of particles. In the case of protons, we also
considered photons from the decay of neutral pions, which are the product of
inelastic collisions between relativistic and cold protons; the latter consist
mainly of protons in the jet bulk and those in the stellar wind of the
companion star. In the case of electrons, we also computed the inverse Compton
scattering of photons from the stellar radiation.

We assume that the relativistic particle populations reach a steady state.
Their spectral densities are obtained by solving the steady state transport
equation, taking injection, escape, and continuous losses into account. For a
discussion of the general form of the transport equation, we refer the reader
to Ginzburg & Syrovatskii (1964).

The complete formulae for calculating of radiation processes in this work can
be found in Blumenthal & Gould (1970), Begelman et al. (1990), Atoyan & Dermer
(2003), Kafexhiu et al. (2014), and references therein.

###  4.2 Spectral energy distribution of 4FGL J1405.1-6119

To fit the SED of the source with our model, we used the radio observations of
Corbet et al. (2019), the X-ray flux obtained in this work, and the
$\gamma$-ray flux from the 4FGL-DR2 catalog (Abdollahi et al., 2020; Ballet et
al., 2020).

In Fig. 4 we show the acceleration, escape, and cooling rates of the
relativistic particles at different locations in the jet. We find that for the
electron population, the losses are dominated by synchrotron cooling
throughout the emitting region. For the protons, both adiabatic and proton-
proton collisions with bulk protons are the dominant mechanisms of energy loss
in the region close to the base; the escape rate competes with them in regions
further from the base and dominates the losses toward the end of the emitting
region.

In Fig. 5 we show the derived nonthermal SED of the jet, $L_{\gamma}$, which
includes all the above radiative processes, where

$\displaystyle L_{\gamma}=E_{\gamma}^{2}\frac{dN}{dE_{\gamma}dt},$ (7)

and $dN$ is the number of photons with energies between
$(E_{\gamma},E_{\gamma}+dE_{\gamma})$ emitted during a time $dt$. The assumed
and derived parameters of the model, with uncertainties reported at the 90%
confidence level for the free parameters, are listed in Table 2. The X-ray
data points shown in Fig. 5 correspond to Epoch 1. Since there is no
significant change in the X-ray flux between the two epochs (see Fig. 1), the
same set of parameters also fits the observations of Epoch 2.

To obtain the reported model parameters, we ran a first set of ${\sim}100$
simulations333We use the term ”simulation” to mean the computation of the
total SED with a defined set of parameters, i.e., a running of our model for a
specific implementation of the parameters. covering a wide range of values,
and decided which of these remained fixed, apart from those resulting from
observed properties of the system. Then we ran another set of $120$
simulations considering the variation of all free parameters (i.e., using $40$
simulations for each parameter variation at a time; $m$, $p$, $\varepsilon$,
and $\eta$), and chose the model with the minimum $\chi^{2}/{\rm dof}$. The
set of values of the free parameter space was covered with a uniform grid for
$p$, $\varepsilon$, $m$, and $\log\eta$. To estimate the free parameter
errors, we computed the chi-squared for each simulation and chose the
$\chi_{\mathrm{min}}^{2}+2.706$ contour for each free parameter, which
represents the $90\leavevmode\nobreak\ \%$ credible region (e.g., Frodesen et
al., 1979). The errors are reported in Table 2. In the particular case of $m$
($\eta$) the upper (lower) bound on the error comes from restricting the
possible values to those reported in the fourth column of Table 2, while for
the case of $\varepsilon$ all the values within the bounds fall below the
aforementioned contour. We obtained a minimum chi-squared goodness of fit of
$\chi^{2}/$d.o.f$=1.06$, with $7$ d.o.f.

We show that the broad spectrum of 4FGL J1405.1-6119 can be explained by the
nonthermal emission associated with a mildly relativistic
($\Gamma_{\mathrm{jet}}=1.9$), lepto-hadronic model of a MQ jet. In
particular, this model represents an approximately parabolic jet
($\varepsilon\approx 0.56$), with a compact emitting region of size $\approx
1.0\times 10^{12}\leavevmode\nobreak\ \mathrm{cm}$ (jet extension could be
orders of magnitude larger), and a magnetic induction field at its base of
$B_{0}\approx 2.8\times 10^{7}\leavevmode\nobreak\ \mathrm{G}$. Relativistic
protons and electrons share the total power in the relativistic populations,
$L_{\mathrm{rel}}=q_{\mathrm{rel}}L_{\mathrm{jet}}=10^{37}\leavevmode\nobreak\
\mathrm{erg}\leavevmode\nobreak\ \mathrm{s}^{-1}$, according to an assumed
hadron-to-lepton ratio of $a\approx 0.11$, and are accelerated through a low-
efficiency mechanism for which $\eta\approx 1.0\times 10^{-4}$. The particle
injection function shows a hard spectrum with spectral index $p=1.98$.

Parameter | Symbol | Adopted/Derived Value | Typical Values | Units
---|---|---|---|---
Luminosity | $L_{\mathrm{jet}}$ | $10^{38}$ | $10^{36}-10^{39}$ | $\mathrm{erg\,s^{-1}}$
Launching distance | $z_{0}$ | $1.0\times 10^{8}$ | $\gtrsim 10^{8}$ | $\mathrm{cm}$
Magnetic field at $z_{0}$ | $B_{0}$ | $2.8\times 10^{7}$ |  $\sim 10^{5}$ – $10^{7}$ | $\mathrm{G}$
Base of acceleration region | $z_{\mathrm{a}}$ | $3.0\times 10^{9}$ | $\gtrsim 10^{8}$ | $\mathrm{cm}$
Top of acceleration region | $z_{\mathrm{t}}$ | $1.0\times 10^{12}$ | up to jet extent | $\mathrm{cm}$
Base half-opening angle† | $\theta_{\mathrm{jet}}$ | $0.3$ | $\lesssim 10$ | ${\mathrm{deg}}$
Viewing angle | $\theta$ | $60$ | $0-180$ | ${\mathrm{deg}}$
Bulk Lorentz factor | $\Gamma_{\mathrm{jet}}$ | $1.9$ | $\sim 1-5$ |
Relativistic power fraction | $q_{\mathrm{rel}}$ | $0.1$ | $\lesssim 0.1$ |
Proton-electron power ratio | $a$ | $0.11$ | $0-100$ |
Magnetic power-law index | $m$ | $1.95\begin{subarray}{c}+0.05\\\ -0.08\end{subarray}$ | $1-2$ |
Acceleration efficiency log | $\log\eta$ | $-4.01\begin{subarray}{c}+0.26\\\ -0.99\end{subarray}$ | $-5\leq\log\eta\leq-1$ |
Injection spectral index | $p$ | $1.98\begin{subarray}{c}+0.35\\\ -0.29\end{subarray}$ | $1.5-2.5$ |
Geometric index | $\varepsilon$ | $0.56\begin{subarray}{c}+0.44\\\ -0.46\end{subarray}$ | $0.1\leq\varepsilon\leq 1$ |

444${}^{\dagger}\,\tan\theta_{\mathrm{jet}}=r_{0}/z_{0}$.

Table 2: Jet model parameters.

In all cases, the values of the parameters are consistent with their commonly
assumed range of values for MQ jets (see, for example, Romero & Vila, 2008;
Vila et al., 2012; Pepe et al., 2015; Escobar et al., 2021). The SED is
computed assuming a viewing angle of $\theta=30^{\circ}$ (i.e., the angle
between the jet axis and the line of sight). With the Lorentz factor from
Table 2, and for lower viewing angles, the same data could be explained with a
less powerful jet ($\sim 0.07$ times the one reported here), while higher
angles would favor a scenario with more powerful jets, up to a factor of $\sim
50$. Instead, maintaining fixed all the values of the parameters but the
viewing angle, the model still manages to explain the observations with at
most a variation of $\approx 10^{\circ}$. As we can see, knowing the viewing
angle is a crucial factor in determining the accretion regime of the emitter.

##  5 Discussion and conclusions

In this work we have analyzed simultaneous XMM-Newton (2–10 keV) and NuSTAR
(3–20 keV) observations of the $\gamma$-ray-emitting binary 4FGL J1405.1-6119.
The difference between the two observed epochs is about six days. We found no
significant pulsation above statistical noise in the XMM-Newton and NuSTAR
light curves. This non-detection leaves open the question of the nature of the
CO.

We fitted the time-averaged spectra of both epochs with empirical models,
which yielded continuum parameters consistent with $\gamma$-ray-emitting
binaries (Kretschmar et al., 2019; An et al., 2015; Yoneda et al., 2021) and
MQs (Natalucci et al., 2014; Soria et al., 2020; Hirsch et al., 2020; Rodi et
al., 2021; Saavedra et al., 2022). For Epoch 1, a power-law model adequately
characterized the spectrum. For Epoch 2, however, the addition of a blackbody
component significantly improved the fit to the observed data. The power-law
index was found to be ${\sim}1.7$ for Epoch 1 and ${\sim}1.4$ for Epoch 2.
Furthermore, the associated blackbody temperature for Epoch 2 was ${\sim}0.8$
keV, corresponding to a compact region of radius
${\lesssim}\leavevmode\nobreak\ 16$ km tentatively associated with the inner
region of an accretion disk.

Corbet et al. (2019) present a different explanation for the emission of 4FGL
J1405.1-6119. They set out a scenario involving the collision of winds between
a rapidly rotating NS and its stellar companion. In addition, Xingxing et al.
(2020) modeled this scenario to explain the GeV emission from 4FGL
J1405.1-6119. In the colliding-wind scenario, the interaction between the
stellar and pulsar winds gives rise to shocked regions that are characterized
by a spiral shape due to Coriolis forces (see Molina & Bosch-Ramon, 2020).
Within this framework, the high-energy emission in a $\gamma$-ray binary
system can be attributed to the up-scattering of photons from the stellar
radiation field due to inverse Compton scattering by the relativistic
electrons and positrons present in the shocked fluid. Synchrotron emission
contributes to a lesser extent, especially in the energy range around
${\sim}1$ MeV (Molina & Bosch-Ramon, 2020).

To explain the new X-ray data and archival multiwavelength observations (ATCA
and Fermi), we implemented a lepto-hadronic jet model under the hypothesis
that 4FGL J1405.1-6119 is a MQ. The parameters associated with the properties
of the relativistic particle populations in the jet and the nonthermal
emission are consistent with the values commonly assumed for this type of
source (see, for example, Romero & Vila, 2008; Vila et al., 2012; Pepe et al.,
2015; Escobar et al., 2021, 2022). The scenario consists of a parabolic jet
with a compact acceleration region, where the high-energy emission is produced
by a hard spectrum of relativistic particles driven by a low-efficiency
acceleration mechanism.

In this scenario, the $\gamma$-ray emission is produced via two processes.
First, inverse Compton scattering produces nonthermal radiation in the energy
range from about 100 keV to 10 GeV. Second, the decay of neutral pions from
proton-proton collisions produces energetic photons at energies above 10 GeV.
At lower energies, the synchrotron emission of electrons completely dominates,
allowing us to accurately model radio and X-ray data. We thus propose that a
lepto-hadronic jet model, which includes both leptonic and hadronic processes,
may be sufficient to explain the multiwavelength emission of 4FGL
J1405.1-6119.

There are essentially two interpretations of the results we obtained with our
jet model. If the accretion regime is sub-Eddington, the CO should be a BH of
at least ${\sim}10\leavevmode\nobreak\ M_{\odot}$. This scenario, like the
case of Cygnus X-1, would require an X-ray luminosity on the order of $\sim
10^{37}\leavevmode\nobreak\ \mathrm{erg}\leavevmode\nobreak\ \mathrm{s}^{-1}$
in the low-hard state (Di Salvo et al., 2001; Makishima et al., 2008), which
is in apparent contradiction with the observations. On the other hand, if the
source is in a super-Eddington accretion regime, this would imply that the CO
consists of a BH of a few solar masses or a NS with a weak magnetic field. In
this case, the X-ray emission from the disk is absorbed by the photosphere of
the wind ejected by the same disk (Abaroa et al., 2023). This is more
consistent with what is observed in the X-ray emission. The supercritical
source may also have an equatorial radio component or equatorial lines with
velocities of $10^{3}-10^{4}$ km s-1, as in the case of SS 433, which could be
observed in the future (Fabrika et al., 2021; Abaroa et al., 2023).

We note that leptonic jet models can also be adopted to explain observations
of other MQs, as in the cases of Cyg X-3 (Zdziarski et al., 2012) and Cyg X-1
(Zdziarski et al., 2014). In our case, however, this type of model would not
fully explain the observations for two main reasons. First, there is no clear
evidence for disk or coronal emission. The X-ray observations could then be
explained by a dominant synchrotron emission, which at these energies hides
the radiation from the disk and/or the corona (see, for example, Fig. 11 of
Bosch-Ramon et al., 2006). On the other hand, a relativistic proton component
seems necessary to explain the observed slope change in the high-energy part
of the Fermi data. In the lepto-hadronic picture, this part of the spectrum is
dominated by the $\gamma$-ray emission, which originates in proton-proton
collisions through neutral pion decays.

As shown in Fig. 2, the flux is higher in Epoch 1 than in Epoch 2. This
behavior is also consistent with the presence of a thermal component in the
spectra of Epoch 2, which can be explained by the contribution of a thermal
inner disk and a reduced nonthermal jet component. In the case of a moderate
or low viewing angle, some of the emission from the disk can escape from the
central funnel in the wind (Abaroa et al., 2023). As for Epoch 1, the
contribution to the total flux may be completely dominated by the jet.

The method we implemented for estimating parameters and uncertainties,
although not very robust, allowed us to explain the multiwavelength behavior
of 4FGL J1405.1-6119, favoring the MQ scenario. A more detailed analysis using
a Markov chain Monte Carlo method would improve these estimates. This method
would also allow a model comparison test to be included (e.g., to compute the
odds ratio between leptonic and lepto-hadronic models, or between the MQ and
colliding-wind scenarios). This analysis is beyond the scope of the current
work, but we will present it in a future paper.

It is important to collect more observational data to evaluate the pros and
cons of the pulsar-star collision wind and MQ scenarios for 4FGL J1405.1-6119.
These additional data should be used to focus on timing analysis, especially
with respect to the orbital period. By analyzing these data, we can gain a
more complete understanding of the nature of this and other $\gamma$-ray-
emitting binary systems.

###### Acknowledgements.

We thank the anonymous reviewer for their valuable comments on this
manuscript. We extend our gratitude to Fiona A. Harrison and Brian W.
Grefenstette for their assistance with the instruments’ technical aspects. We
thank Sergio Campana for his suggestions that helped improve this work. EAS,
FAF and JAC acknowledge support by PIP 0113 (CONICET). FAF is fellow of
CONICET. JAC is CONICET researchers. JAC is a María Zambrano researcher fellow
funded by the European Union -NextGenerationEU- (UJAR02MZ). This work received
financial support from PICT-2017-2865 (ANPCyT). JAC was also supported by
grant PID2019-105510GB-C32/AEI/10.13039/501100011033 from the Agencia Estatal
de Investigación of the Spanish Ministerio de Ciencia, Innovación y
Universidades, and by Consejería de Economía, Innovación, Ciencia y Empleo of
Junta de Andalucía as research group FQM-322, as well as FEDER funds. GER
acknowledges financial support from the State Agency for Research of the
Spanish Ministry of Science and Innovation under grants
PID2019-105510GB-C31AEI/10.13039/501100011033/ and
PID2022-136828NB-C41/AEI/10.13039/501100011033/, and by “ERDF A way of making
Europe”, by the ”European Union”, and through the ”Unit of Excellence María de
Maeztu 2020-2023” award to the Institute of Cosmos Sciences
(CEX2019-000918-M). Additional support came from PIP 0554 (CONICET). GJE
acknowledges financial support from the European Research Council for the ERC
Consolidator grant DEMOBLACK, under contract no. 770017.

## References

  * Abaroa et al. (2023) Abaroa, L., Romero, G. E., & Sotomayor, P. 2023, A&A, 671, A9
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33
  * An et al. (2015) An, H., Bellm, E., Bhalerao, V., et al. 2015, ApJ, 806, 166
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Atoyan & Dermer (2003) Atoyan, A. M. & Dermer, C. D. 2003, ApJ, 586, 79
  * Ballet et al. (2020) Ballet, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020, arXiv e-prints, arXiv:2005.11208
  * Begelman et al. (1990) Begelman, M. C., Rudak, B., & Sikora, M. 1990, ApJ, 362, 38
  * Benaglia & Romero (2003) Benaglia, P. & Romero, G. E. 2003, A&A, 399, 1121
  * Blumenthal & Gould (1970) Blumenthal, G. R. & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237
  * Bosch-Ramon et al. (2006) Bosch-Ramon, V., Romero, G. E., & Paredes, J. M. 2006, A&A, 447, 263
  * Buchner (2021) Buchner, J. 2021, The Journal of Open Source Software, 6, 3001
  * Camilo et al. (2009) Camilo, F., Ray, P. S., Ransom, S. M., et al. 2009, ApJ, 705, 1
  * Chernyakova & Malyshev (2020) Chernyakova, M. & Malyshev, D. 2020, in Multifrequency Behaviour of High Energy Cosmic Sources - XIII. 3-8 June 2019. Palermo, 45
  * Chernyakova et al. (2019) Chernyakova, M., Malyshev, D., Paizis, A., et al. 2019, A&A, 631, A177
  * Corbet et al. (2019) Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2019, ApJ, 884, 93
  * De Becker & Raucq (2013) De Becker, M. & Raucq, F. 2013, A&A, 558, A28
  * del Palacio et al. (2023) del Palacio, S., García, F., De Becker, M., et al. 2023, A&A, 672, A109
  * Di Salvo et al. (2001) Di Salvo, T., Done, C., Życki, P. T., Burderi, L., & Robba, N. R. 2001, ApJ, 547, 1024
  * Dubus et al. (2017) Dubus, G., Guillard, N., Petrucci, P.-O., & Martin, P. 2017, A&A, 608, A59
  * Eichler & Usov (1993) Eichler, D. & Usov, V. 1993, ApJ, 402, 271
  * Escobar et al. (2021) Escobar, G. J., Pellizza, L. J., & Romero, G. E. 2021, A&A, 650, A136
  * Escobar et al. (2022) Escobar, G. J., Pellizza, L. J., & Romero, G. E. 2022, A&A, 665, A145
  * Fabrika et al. (2021) Fabrika, S. N., Atapin, K. E., Vinokurov, A. S., & Sholukhova, O. N. 2021, Astrophysical Bulletin, 76, 6
  * Falcke et al. (2004) Falcke, H., Körding, E., & Markoff, S. 2004, A&A, 414, 895
  * Frodesen et al. (1979) Frodesen, A. G., Skjeggestad, O., & Tofte, H. 1979, Probability and Statistics in Particle Physics (Bergen, Norway: Universitetsforlaget)
  * Ginzburg & Syrovatskii (1964) Ginzburg, V. L. & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays
  * Grimm et al. (2002) Grimm, H. J., Gilfanov, M., & Sunyaev, R. 2002, A&A, 391, 923
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103
  * Hirsch et al. (2020) Hirsch, M., Pottschmidt, K., Smith, D. M., et al. 2020, A&A, 636, A51
  * Huppenkothen et al. (2019) Huppenkothen, D., Bachetti, M., Stevens, A. L., et al. 2019, ApJ, 881, 39
  * Hurkett et al. (2008) Hurkett, C. P., Vaughan, S., Osborne, J. P., et al. 2008, ApJ, 679, 587
  * Kafexhiu et al. (2014) Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90, 123014
  * Kretschmar et al. (2019) Kretschmar, P., Fürst, F., Sidoli, L., et al. 2019, New A Rev., 86, 101546
  * Krolik (1999) Krolik, J. H. 1999, Active galactic nuclei : from the central black hole to the galactic environment
  * Mahy et al. (2015) Mahy, L., Rauw, G., De Becker, M., Eenens, P., & Flores, C. A. 2015, A&A, 577, A23
  * Makishima et al. (2008) Makishima, K., Takahashi, H., Yamada, S., et al. 2008, PASJ, 60, 585
  * Moldón et al. (2014) Moldón, J., Ribó, M., Paredes, J. M., & Johnston, S. 2014, in International Journal of Modern Physics Conference Series, Vol. 28, International Journal of Modern Physics Conference Series, 1460173
  * Molina & Bosch-Ramon (2020) Molina, E. & Bosch-Ramon, V. 2020, A&A, 641, A84
  * Natalucci et al. (2014) Natalucci, L., Tomsick, J. A., Bazzano, A., et al. 2014, ApJ, 780, 63
  * Pepe et al. (2015) Pepe, C., Vila, G. S., & Romero, G. E. 2015, A&A, 584, A95
  * Rodi et al. (2021) Rodi, J., Tramacere, A., Onori, F., et al. 2021, ApJ, 910, 21
  * Romero & Vila (2008) Romero, G. E. & Vila, G. S. 2008, A&A, 485, 623
  * Saavedra et al. (2023) Saavedra, E. A., Romero, G. E., Bosch-Ramon, V., & Kefala, E. 2023, MNRAS, 525, 1848
  * Saavedra et al. (2022) Saavedra, E. A., Sotomayor Checa, P., Vieyro, F. L., & Romero, G. E. 2022, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 63, 274
  * Soria et al. (2020) Soria, R., Blair, W. P., Long, K. S., Russell, T. D., & Winkler, P. F. 2020, ApJ, 888, 103
  * Strüder et al. (2001) Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18
  * Turner et al. (2001) Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27
  * Verner et al. (1996) Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487
  * Vila et al. (2012) Vila, G. S., Romero, G. E., & Casco, N. A. 2012, A&A, 538, A97
  * Weng et al. (2022) Weng, S.-S., Qian, L., Wang, B.-J., et al. 2022, Nature Astronomy, 6, 698
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914
  * Xingxing et al. (2020) Xingxing, H., Jumpei, T., & Qingwen, T. 2020, MNRAS, 494, 3699
  * Yoneda et al. (2021) Yoneda, H., Khangulyan, D., Enoto, T., et al. 2021, ApJ, 917, 90
  * Zdziarski et al. (2014) Zdziarski, A. A., Pjanka, P., Sikora, M., & Stawarz, Ł. 2014, MNRAS, 442, 3243
  * Zdziarski et al. (2012) Zdziarski, A. A., Sikora, M., Dubus, G., et al. 2012, MNRAS, 421, 2956