# Root Cause Analysis In Microservice Using Neural Granger Causal Discovery

Cheng-Ming Lin, Ching Chang, Wei-Yao Wang, Kuang-Da Wang, Wen-Chih Peng

###### Abstract

In recent years, microservices have gained widespread adoption in IT
operations due to their scalability, maintenance, and flexibility. However, it
becomes challenging for site reliability engineers (SREs) to pinpoint the root
cause due to the complex relationships in microservices when facing system
malfunctions. Previous research employed structured learning methods (e.g.,
PC-algorithm) to establish causal relationships and derive root causes from
causal graphs. Nevertheless, they ignored the temporal order of time series
data and failed to leverage the rich information inherent in the temporal
relationships. For instance, in cases where there is a sudden spike in CPU
utilization, it can lead to an increase in latency for other microservices.
However, in this scenario, the anomaly in CPU utilization occurs before the
latency increase, rather than simultaneously. As a result, the PC-algorithm
fails to capture such characteristics. To address these challenges, we propose
RUN, a novel approach for root cause analysis using neural Granger causal
discovery with contrastive learning. RUN enhances the backbone encoder by
integrating contextual information from time series, and leverages a time
series forecasting model to conduct neural Granger causal discovery. In
addition, RUN incorporates Pagerank with a personalization vector to
efficiently recommend the top-k root causes. Extensive experiments conducted
on the synthetic and real-world microservice-based datasets demonstrate that
RUN noticeably outperforms the state-of-the-art root cause analysis methods.
Moreover, we provide an analysis scenario for the sock-shop case to showcase
the practicality and efficacy of RUN in microservice-based applications. Our
code is publicly available at https://github.com/zmlin1998/RUN.

##  1 Introduction

Root cause analysis plays a crucial role in numerous domains such as cloud
system operations (Soldani, Forti, and Brogi 2022), manufacturing processes
(Oliveira, Miguéis, and de Moura Borges 2023), or telecommunications networks
(Chen et al. 2022). For instance, during the manufacturing process of wafers,
if the thickness of a wafer is found to deviate from the standard, root cause
analysis can be employed to identify the specific manufacturing step that
caused the abnormal thickness. By applying root cause analysis to sensor data
with complex relationships, these situations can effectively pinpoint the
underlying causes when a system malfunction occurs.

In recent years, as the companies experience continued growth, the system
operations have scaled up and grown increasingly complex. Therefore, these
organizations have opted to migrate from the so-called monolithic architecture
to microservice architecture (Liu et al. 2021). Microservice architecture
offers benefits such as better scalability, easier maintenance, and greater
flexibility. Each service can be independently scaled based on demand,
enabling efficient resource utilization. Developers can focus on individual
services, making it simpler to debug, test, and deploy changes. Despite the
numerous benefits of microservices, when an anomaly arises in one service, the
interdependencies among services can create a domino effect, resulting in
subsequent issues and ultimately leading to system failure (Zhang et al.
2022). In such scenarios, in-depth analysis becomes imperative to identify the
culprit of the anomaly and mitigate the problem effectively.

Figure 1: An example of causal structure discovery-based techniques for RCA.

However, in microservice monitoring systems, only the operational values of
the system are recorded, without documenting the relationships between them.
Hence, researchers have recently employed causal structure discovery-based
techniques for Root Cause Analysis (RCA) in cloud applications (Wang et al.
2023), aiming to identify the underlying causes of anomalies. Figure 1
illustrates a flow of RCA using the causal structure discovery-based approach.
When the anomaly detection alarm is triggered by an anomalous Key Performance
Indicator (KPI), engineers initially designate that particular KPI as the
trigger point. Subsequently, they aim to identify the underlying root cause of
this trigger point. To achieve this, they construct a causal graph that
establishes relationships between different KPIs, allowing them to precisely
pinpoint the culprit of the anomaly using the insights provided by the causal
graph.

On the other hand, the Granger causality (Granger 1969) analysis is another
widely recognized approach used to assess whether a set of time series $x$ is
a causal factor for another set of time series $y$. Granger causality has
gained attention and is widely acknowledged for its advantages of
interpretability and compatibility with the emergence of deep learning (Tank
et al. 2022). Previous works (Nauta, Bucur, and Seifert 2019) utilized time
series forecasting models and interpreted model parameters as the
relationships between variables to establish causal relationships. However,
prior approaches in neural Granger causal discovery have not effectively
leveraged the contextual information inherent in temporal data.

The previous approach (Yue et al. 2022) generated positive and negative
samples at instance-wise and temporal dimensions for generating fine-grained
representations for any granularity. However, we notice that real-world time
series often exhibit multiple-periodicity. For instance, the temperature tends
to reach its peak at noon every day, there is an influx of people in the city
during weekends, and electricity consumption significantly increases during
the summer every year. Therefore, this approach might inadvertently treat
timestamps with the same periodicity as negative pairs, even though they
actually hold similar significance.

To tackle these challenges, we propose a novel framework RUN, which employs a
self-supervised neural Granger causal discovery approach for conducting root
cause analysis. To capture contextual information in time series, we employ a
self-supervised learning scheme that ensures timestamps with diverse contexts,
but the same timestamps are learned to exhibit proximity, where we only treat
the identical timestamps as positive pairs for contrastive learning. This way
aims to prevent the inclusion of erroneous information stemming from negative
pairs. After enhancing the backbone encoder to acquire improved
representations, we utilize neural Granger causal discovery to explore the
causal relationships among variables. Subsequently, we proceed to identify the
root causes of the trigger point within the acquired causal graph using
Pagerank with a personalized vector. Comprehensive experiments conducted on
both synthetic and real-world microservice-based datasets conclusively
demonstrate that our proposed framework outperforms state-of-the-art root
cause analysis methods.

Our main contributions are summarized as follows:

  * •

We propose a self-supervised neural Granger causal discovery-based root cause
analysis framework RUN, which captures contextual information in temporal
data, and then leverages the time series forecasting model to construct a
causal graph between multivariate time series. RUN employs Pagerank on the
derived causal graph to identify the root cause of the trigger point.

  * •

We introduce an innovative self-supervised learning method for time series
data that exclusively treats identical timestamps with distinct contextual
information as positive pairs. This approach mitigates the problem of
misidentifying timestamps with similar periodicity as negative pairs, thereby
preventing the separation of representations for timestamps sharing the same
periodicity.

  * •

RUN outperforms existing state-of-the-art methods on synthetic datasets and
data generated from a real-world application. We further conduct an ablation
study to validate the effectiveness of our contrastive learning method.

##  2 Related Work

Figure 2: Illustrations of contextual information from the same timestamps but
with different contexts. The red window and green window respectively
represent two distinct types of contextual information. The same timestamp
with different contexts should be close. Figure 3: An illustrated issue of
negative pair selection. Figure 4: Overview of our proposed framework, RUN,
consisting of three stages: 1) Maximizing the positive pair to capture the
contextual information; 2) Neural Granger causal discovery to derive the
causal graph from multivariate time series; and 3) The diagnosis stage infers
the root cause from the obtained causal graph.

###  2.1 Root Cause Analysis in Microservices

Within large organizations, the adoption of microservices has become the
prevailing architecture due to its advantages in terms of enhanced
scalability, simplified maintenance, and increased flexibility. When the
anomaly detection alarm is activated by an anomalous KPI, engineers often have
to invest a significant amount of time to pinpoint the underlying cause behind
this anomaly. Consequently, they engage in the development of automated root
cause analysis systems to mitigate the time spent on the investigation.
$\epsilon$-diagnosis (Shan et al. 2019) considers the time series of KPIs and
computes the similarity between them during normal periods and abnormal
periods, which serves as an anomaly indicator. AutoMAP (Ma et al. 2020)
utilizes the PC-algorithm to construct a causality graph based on the KPIs. It
further conducts a random walk to determine the candidate’s root causes. RCD
(Ikram et al. 2022) employs the hierarchical and localized learning algorithm.
It deems the failures as an intervention, and adopts a divide-and-conquer
scheme to divide the variables into multiple subsets for conducting $\Psi$-PC
to find interventional targets from each subset.

The candidate root causes from all subsets are then combined, and the same
process is iteratively applied until further splitting of candidate root
causes is not possible. Nonetheless, these methods rely on similarity or the
PC-algorithm, both of which overlook the essential role of temporal dependency
in root cause analysis. For instance, when a sudden spike in CPU utilization
occurs, it can lead to increased latency in other microservices, effectively
encompassing the temporal dimension of the information. Therefore, we aim to
develop a neural Granger causal discovery method that can leverage the
temporal dependency to identify the root causes more accurately.

###  2.2 Neural Granger Causal Discovery

As deep Neural Networks (NNs) continue to advance rapidly, researchers have
started using Recurrent Neural Networks or other Temporal Convolutional
Networks to infer nonlinear Granger causality. MSNGC (Fan et al. 2023)
extracts diverse causal information from the data, considering various delay
ranges, and effectively integrates them using learned attention weights. Tank
et al. (2022) applies structured component-wise multilayer perceptrons (MLPs)
combined with sparsity-inducing penalties on the weights. CUTS (Cheng et al.
2023) utilizes EM-Style joint causal graph learning and missing data
imputation for irregular temporal data. Nevertheless, previous works have not
effectively leveraged contextual information in time series data. As a result,
we present an innovative approach for time series contrastive learning that
enables to the capture of contextual information from identical timestamps but
with different contexts, as depicted in Figure 2.

###  2.3 Contrastive Learning

Contrastive learning has been widely used in various domains, such as natural
language processing (NLP) (Gao, Yao, and Chen 2021), computer vision (CV)
(Chen and He 2021), and time series (Yue et al. 2022). TS2Vec (Yue et al.
2022) samples positive pairs which are the same timestamp with different
contexts, and negative pairs which are different timestamps. However, we
observe that real-world time series frequently demonstrate multiple-
periodicity. With the advancement of contrastive learning, (Grill et al. 2020;
Chen and He 2021) argue that it can perform well even without negative pairs
because they aim to prevent the inclusion of incorrect negative pairs. Hence,
we propose an innovative approach that exclusively leverages identical
timestamps as positive pairs for self-supervised learning. This scheme is
designed to prevent the incorporation of information from negative pairs. In
Figure 3, we can observe that the time series data exhibits periodic patterns.
Therefore, if we were to consider the two circled red points on the graph as
negative pairs, it would lead to incorporating incorrect information from
them.

##  3 Problem Formulation

The monitoring system will monitor the operational status of the microservice
system through KPIs, and in typical scenarios, we pinpoint the root cause when
an anomalous KPI triggers the anomaly detection system. We term the anomalous
KPI as a trigger point and subsequently identify the underlying root causes
behind the trigger point. These KPIs belong to multivariate time series data.
In the microservice architecture, identifying the underlying causes for the
anomaly can be formalized as follows. Multivariate time series X consists of
$N$ features, $\textbf{X}=[X_{1},X_{2},...,X_{N}]$. We collect the
corresponding time series with a specified time period $T$, each time series
can be denoted as $X_{i}=[x^{1}_{i},x^{2}_{i},\ldots,x^{T}_{i}]$,
$\textbf{X}\in\mathbb{R}^{N\times T}$. We utilize the neural Granger causal
discovery method to derive the causal graph $\hat{\textbf{G}}=\\{V,E\\}$,
where nodes $V$ denotes each feature $X_{i}$ and edges $E$ denotes that
feature $X_{i}$ Granger causes feature $X_{j}$. Our goal is to identify the
root cause $X_{culprit}$ which leads to the trigger point within the causal
graph $\hat{\textbf{G}}$.

##  4 Methodology

Our proposed framework is shown in Figure 4, which consists of three stages:
the pre-training stage, the neural Granger causal discovery stage, and the
diagnosis stage. In the pre-training stage, we enhance the backbone encoder to
generate informative time series representations that incorporate contextual
information. This enhancement is achieved by maximizing the agreement among
instances with the same timestamp but with different contexts. In the neural
Granger causal discovery stage, we utilize the time series forecasting model
to exploit the causal graph and prune the spurious edge to get the final
Directed Acyclic Graphs (DAGs). In the diagnosis stage, we apply a random walk
algorithm called Pagerank with a personalization vector which recommends the
most likely root causes based on the causal graph.

###  4.1 Pre-training Stage

We utilize DLinear (Zeng et al. 2023) as our backbone encoder for time series
forecasting. In order to enhance our backbone encoder, we design contrastive
learning without relying on negative pairs. We consider that negative pairs
are not suitable for time series because they exhibit periodic patterns as
depicted in Figure 3. When selecting different timestamps as negative pairs,
it is possible that they belong to the same periodicity. If we choose the
incorrect negative pairs, it may lead to separating embeddings that should be
similar. Therefore, we learn the context information by considering instances
with the same timestamp but with different contexts as positive pairs.

First, we utilize random cropping following (Yue et al. 2022) by randomly
sampling two overlapping time segments $y_{1}$ and $y_{2}$,

$y_{1}=[a_{1},b_{1}],y_{2}=[a_{2},b_{2}],0<a_{1}\leq a_{2}\leq b_{1}\leq
b_{2}\leq T.$ (1)

The two views $y_{1}$ and $y_{2}$ are processed by our backbone encoder $f$
and a projector $g$. To learn the contextual information without negative
pairs, we utilize the same schema of SimSiam (Chen and He 2021). We denote the
two representations as $z_{1}=g(f(y_{1}))$ and $p_{2}=f(y_{2})$, and maximize
the similarity between both sides. Here, we define our loss function as
follows:

$L_{con}=-\frac{1}{2}(cos(z_{1},stopgrad(p_{2}))+cos(stopgrad(p_{1}),z_{2})).$
(2)

###  4.2 Neural Granger Causal Discovery Stage

RUN represents the causal graph as a Causal Attention Matrix
$\textbf{G}=\\{0\leq\alpha_{ij}\leq 1,\ \forall i,j\in[1,2,...,N]\\}$, where
the element $\alpha_{ij}$ demonstrates the attention score from $X_{i}$ to
$X_{j}$. Each element $\alpha_{ij}$ serves as a trainable parameter,
contributing to identifying causal relationships among time series. We define
the final causal graph $\bm{\widetilde{G}}=\\{g_{ij}\in[0,1],\forall
i,j\in[0,1]\\}$. If the attention score $\alpha_{ij}$ is beyond a certain
threshold, $g_{ij}$ will be 1, indicating that time series $X_{i}$ Granger
causes $X_{j}$. The proposed neural Granger causal discovery comprises two
components: time series forecasting and causal graph discovery.

#### Time Series Forecasting

In time series forecasting, we propose $N$ independent neural network
$f_{\theta_{i}}$ for time series $X_{i}$ to fit its data generation function
(Cheng et al. 2023). We define the input for our model by applying a sliding
window with a size of $w$ across all the historical time series X and
multiplying with corresponding attention score $\alpha_{ni}$ and the forecast
$\hat{x}_{i}$ at timestamp $t$ is the output of the neural network
$f_{\theta_{i}}$

$\displaystyle\hat{x}^{t}_{i}=f_{\theta_{i}}(\textbf{X}\odot
G)=f_{\theta_{i}}(x^{t-w:t-1}_{1}\odot\alpha_{1i},...,x^{t-w:t-1}_{N}\odot\alpha_{Ni}),$
(3) $\displaystyle i=1,2,...,N.$

An overview of time series forecasting is illustrated in 5. Each independent
neural network has the same architecture but a different target time series.
For a neural network $f_{\theta_{i}}$, the goal is to forecast corresponding
time series $X_{i}$ by minimizing the following loss function

$L_{MSE}=-\frac{1}{N-w}\sum_{t=w+1}^{T}(x^{t}_{i}-\hat{x}^{t}_{i})^{2}$ (4)

Figure 5: Overview of time series forecasting. There are $N$ independent
neural networks for each time series $i$ to predict their causal
relationships.

#### Causal Graph Discovery

After training the networks, the learnable Causal Attention Matrix $G$ can be
used to interpret the causal relationships among time series. When
$\alpha_{ij}$ exceeds a certain threshold $H$, we can infer that time series
$X_{i}$ Granger causes $X_{j}$ and add an edge from $X_{i}$ to $X_{j}$ in the
causal graph $\bm{\widetilde{G}}$, that is

$\bm{\widetilde{G}}_{ij}=\begin{cases}1,&\text{if}\ \alpha_{ij}>H\\\
0,&\text{else}\end{cases}$ (5)

Based on the definition of a causal graph, cycles are not allowed in the
graph. Therefore, we need to incorporate a pruning strategy to remove spurious
edges. We calculate the similarity of each edge, which involves computing the
Pearson correlation between two connected nodes. Subsequently, we eliminate
the edge with the lowest similarity iteratively until $\bm{\widetilde{G}}$
transforms into the final causal graph $\hat{\textbf{G}}$.

###  4.3 Diagnosis Stage

In the diagnosis stage, we follow GrootRank (Wang et al. 2021) to apply
PageRank with node weight personalization to calculate the root cause ranking.
Wang et al. (2021) observes that dangling nodes are more likely to be the root
cause. Hence, we customize the personalization vector as $P_{d}$ and $P_{n}$,
where $P_{d}$ is the personalization score for dangling nodes and $P_{n}$ is
for the remaining nodes. In the case of tied rankings, we calculate the access
distance from the trigger point to resolve the tie. We calculate the access
distance (AD) differently from GrootRank

$\bm{AD}=\begin{cases}D,&\text{distance from trigger point to the variable}\\\
0,&\text{if any “access" is not reachable}\end{cases}$ (6)

In GrootRank, the authors set the distance of unreachable nodes as infinity
and consider that shorter access distance will more likely be the root cause.
However, this may contradict the observed phenomena which suggest that
dangling nodes are more likely to be the root cause. Hence, we set the
distance of the unreachable node as $0$ and consider the bigger access
distance will more likely be the root cause. Finally, RUN outputs the top-k
root causes based on the root cause ranking.

##  5 Experiments

###  5.1 Experiment Settings

#### Dataset

As no publicly available real-world dataset for root cause analysis is
accessible due to data confidentiality, we test on a synthetic dataset and a
test bed utilizing an actual microservice-based application.

  * •

Synthetic data: To generate synthetic data, we follow RCD Ikram et al. (2022)
to randomly generate DAG and generate conditional probability tables (CPT) for
the normal period dataset. Then, we inject failures by randomly choosing one
node as the root cause and regenerating CPT for it. Hence, we can get the
anomalous period dataset and combine them into one case. Finally, we also
utilize an anomaly detection system to identify the trigger point. The data
are generated with node counts of 10, 20, 30, 40, and 50, spanning over 2,000
timestamps.

  * •

Sock-shop (Daniel Holbach 2022): The framework of sock-shop encompasses a
total of 13 microservices, each developed using distinct technologies. These
microservices are deployed on individual virtual machines or containers.
Communication among them occurs through HTTP-based API requests. Additionally,
these microservices offer a substantial volume of statistical data in the form
of various metrics such as CPU and memory utilization. Within the sock-shop
dataset, ten distinct categories of root causes are present, each comprising
five instances.

#### Baselines

We compare our performance with the following baseline approaches:

  * •

$\epsilon$-Diagnosis (Shan et al. 2019): $\epsilon$-Diagnosis analyzes the
time series of KPIs and calculates the similarity between them during both
normal and abnormal periods, which serves as an indicator for anomalies. If
the similarity of a particular KPI between normal and anomalous time periods
falls below a certain confidence threshold, indicating significant changes in
that KPI before and after the anomaly occurrence, it can be considered a
candidate root cause.

  * •

AutoMAP (Ma et al. 2020): AutoMAP creates a weighted causal graph utilizing an
adapted PC algorithm. In this graph, the weight of an edge between two
microservices reflects the degree of dependence of a performance metric. The
process concludes by navigating nodes through a random walk algorithm and
identifying the root cause via the correlation metric from the weighted causal
graph.

  * •

RCD (Ikram et al. 2022): RCD treats failures as interventions on the nodes
representing root causes, and utilizes a localized hierarchical learning
algorithm to identify these root causes of failure. It also employs a divide-
and-conquer scheme to decrease the time required for inferring the root cause
from the entire graph.

  * •

$\Psi$-PC: $\Psi$-PC is a specific instance of RCD. $\Psi$-PC will acquire the
entire causal graph, which might not be essential for root cause
identification. Additionally, localizing learning would demand extra time.

  * •

CausalRCA (Xin, Chen, and Zhao 2023): CausalRCA uses a gradient-based causal
structure learning method to generate weighted causal graphs and a root cause
inference method to localize root cause metrics.

Dataset | # of Time Series | # of Timestamps
---|---|---
Synthetic data | 10 | 2000
20
30
40
50
Sock-shop | 38 | 600
Table 1: Statistics of the synthetic and sock-shop datasets.

CPU hog Memory leak $\epsilon$-Diag. AMAP $\Psi$-PC* RCD* CausalRCA RUN
$\epsilon$-Diag. AMAP $\Psi$-PC* RCD* CausalRCA RUN HR@1 Carts 0.20 0.05 0.02
0.00 0.00 0.20 0.00 0.14 0.02 0.02 0.00 0.00 Catalogue 0.20 0.22 0.02 0.08
0.00 0.40 0.40 0.00 0.10 0.06 0.00 0.00 Orders 0.20 0.01 0.06 0.12 0.80 0.00
0.00 0.26 0.02 0.04 0.00 0.20 Payment 0.60 0.07 0.06 0.04 0.00 0.40 0.20 0.01
0.02 0.02 0.00 0.00 User 0.20 0.29 0.06 0.12 0.20 0.20 0.40 0.08 0.20 0.14
0.00 0.40 Avg. 0.28 0.13 0.04 0.07 0.04 0.40 0.20 0.10 0.07 0.06 0.00 0.12
HR@3 Carts 0.40 0.23 0.04 0.04 0.00 0.40 0.00 0.20 0.02 0.06 0.00 0.20
Catalogue 0.20 0.22 0.04 0.12 0.00 0.60 0.40 0.00 0.06 0.12 0.00 0.40 Orders
0.20 0.36 0.06 0.30 0.00 0.80 0.00 0.37 0.06 0.14 0.00 0.20 Payment 0.60 0.07
0.06 0.12 0.00 0.40 0.20 0.05 0.02 0.16 0.00 0.40 User 0.20 0.53 0.04 0.12
0.40 0.20 0.40 0.27 0.16 0.28 0.20 0.40 Avg. 0.32 0.28 0.05 0.16 0.08 0.48
0.20 0.14 0.06 0.15 0.04 0.32 HR@5 Carts 0.40 0.24 0.08 0.22 0.20 0.40 0.00
0.25 0.02 0.04 0.20 0.20 Catalogue 0.20 0.23 0.06 0.08 0.40 0.80 0.40 0.00
0.08 0.12 0.00 0.60 Orders 0.20 0.51 0.06 0.36 0.00 0.80 0.00 0.49 0.08 0.08
0.00 0.20 Payment 0.60 0.07 0.06 0.06 0.00 0.40 0.20 0.11 0.04 0.10 0.20 0.40
User 0.20 0.55 0.06 0.34 0.40 0.20 0.40 0.30 0.20 0.28 0.20 0.40 Avg. 0.32
0.32 0.06 0.21 0.20 0.52 0.20 0.18 0.08 0.12 0.12 0.36

Table 2: HR@k on sock-shop data. Bold indicates the best performance and
underline represents the second best performance. We note that $\Psi$-PC* and
RCD* indicate different results compared with the original paper, where we
directly reran the experiments from the official codes.

#### Evaluation Metrics

We evaluate our solution by top-k hit ratio (HR@$k$) and mean reciprocal rank
(MRR).

  * •

HR@$k$ represents the probability of getting the correct root cause from the
top-k outputs.

  * •

MRR sums up the reciprocal of the rank of the root cause. If the root cause is
not included in the output, its rank can be deemed as infinity, resulting in a
score of zero.

#### Implementation Details

We implement our method on a machine with AMD EPYC 7302 16-Core CPU, NVIDIA
RTX A5000 graphics cards. In the time series forecasting stage, the window
size $w$ is set to $32$. We use the Adam (Kingma and Ba 2015) optimizer and
set the learning rate as $0.001$ and the batch size as $128$. In the causal
graph discovery stage, the threshold $H$ is set to $0.5$. The training epochs
of the pre-training and fine-tuning stages are set to 50. In the diagnosis
stage, we set the value of the personalization vector $P_{d}$ as 1 and $P_{n}$
as 0.5, similar to (Wang et al. 2021). The k in $HR@k$ is set to 1, 3, and 5.

CPU hog Memory leak $\epsilon$-Diag. RUN $\epsilon$-Diag. RUN MRR Carts 0.273
0.299 0.036 0.158 Catalogue 0.203 0.537 0.412 0.207 Orders 0.212 0.825 0.017
0.267 Payment 0.613 0.494 0.437 0.164 User 0.221 0.264 0.231 0.443 Avg. 0.304
0.484 0.227 0.248

Table 3: MRR on sock-shop data.

###  5.2 Overall Performance

#### Sock-shop

Table 2 demonstrates the HR@k of root cause analysis methods between RUN and
baselines, which demonstrates that our proposed model achieves at least a 63%
improvement in CPU hog, and 80% improvement in Memory leak compared with the
best-performing baseline. Considering baselines, $\epsilon$-diagnosis
outperforms other baselines because the root cause in the sock-shop data
exhibits distinct behaviors during normal and anomaly periods. AutoMAP’s
performance is suboptimal in some cases due to the fact that it constructs the
causal graph utilizing PC-algorithm, which, however, overlooks the temporal
dependency. Although we adhere to the settings used in RCD, $\Psi$-PC and RCD
yield varying root causes in each run and perform worse on the sock-shop
dataset. CausalRCA neglects temporal information within time series data,
leading to suboptimal performance in sock-shop. The comparison of RUN and
causal discovery-based approaches reveals the significance of temporal
dependency. RUN leverages positive pairs which capture the contextual
information, and applies neural Granger causal discovery to consider the
temporal dimension. This underscores the importance of considering temporal
dependency, particularly in time series data within microservice systems.

As $\Psi$-PC and RCD cannot output the rank of the root cause, they are
omitted for comparing performance in terms of MRR. Table 3 illustrates the MRR
of root cause analysis methods between RUN and $\epsilon$-diagnosis. The table
clearly illustrates that even though our performance in HR@1 for Memory leaks
is the second best, it is still noticeable that the root cause maintains a
relatively high ranking. Furthermore, we observe that RUN outperforms
$\epsilon$-diagnosis significantly, underscoring the impressive performance of
RUN in identifying CPU hog incidents.

(a) HR@1

(b) HR@3

(c) HR@5

Figure 6: HR@k on a synthetic dataset Figure 7: Ablation study removing the
pre-training stage and incorporating the negative pairs in contrastive
learning.

#### Synthetic Data

Figure 6 shows the HR@k of root cause analysis methods between RUN and
baselines. However, $\epsilon$-diagnosis falls short due to the fact that
synthetic data lack the distinctive behaviors during both normal and anomaly
periods that are characteristic of microservice-based data when anomalies
occur. Since the synthetic dataset is generated using the function of
generateCPT in PyAgrum, it might not fully exhibit the characteristics of time
series data. Both $\Psi$-PC and RCD, on the other hand, learn causal graphs
using the PC-algorithm from these synthetic datasets. As a result, they
exhibit better performance on synthetic data when compared to the sock-shop
data. However, our model is specifically designed for performing root cause
analysis within microservice systems consisting of abundant time series data,
making it less suitable for application to the synthetic dataset. Nonetheless,
our approach still demonstrates competitive performance with both RCD and
$\Psi$-PC.

###  5.3 Ablation Study

To dissect the contributions arising from different designs of RUN, we
conducted an ablation study by omitting the pre-training stage and
incorporating the consideration of negative pairs in contrastive learning on
sock-shop data for its microservice-based architecture. The results are
presented in Figure 7. As expected, removing the pre-training stage
significantly results in a notable reduction in the performance of RUN. This
indicates that when performing root cause analysis via neural Granger causal
discovery, incorporating contextual information can effectively elevate
performance. Nevertheless, excluding the incorporation of negative pairs in
contrastive learning does not lead to a substantial decrease in performance.
We theorize that this finding can be ascribed to the limited temporal extent
of the sock-shop data, potentially constraining the capacity to identify the
multi-periodicity within the temporal domain.

###  5.4 Case Study: Neural Granger Causal Discovery

In order to investigate the efficacy of neural Granger causal discovery, we
sample a causal graph in which each node has a pathway to the root cause. In
Figure 8, the yellow node represents the trigger point and the red node
represents the root cause. Figure 8 represents causal relationships that are
associated with temporal order based on Granger causality. To identify the
root cause, each edge is reversed from the original causal graph for PageRank
analysis. According to the observation from Wang et al. (2021), the root cause
is a dangling node and impacts the trigger point. Our approach validates the
observation and successfully identifies the root cause from the causal graph.
It also illustrates the neural Granger causal discovery-based approach in
explaining how to identify the root cause.

Figure 8: Construct the causal graph from the multivariate time series.

##  6 Conclusion

In this paper, we propose RUN for identifying the root causes in
microservices. Existing works apply PC-algorithm on time series data which
leads to overlooking the temporal dimension of data. Distinct from these
works, our proposed method is able to capture the temporal dependency. Based
on neural Granger causal discovery architecture, our model learns the graph
structure using a time series forecasting model. However, previous works on
neural Granger causal discovery did not leverage the inherent contextual
information in time series data. Hence, contrastive learning without negative
pairs is proposed to leverage the context of time series to enhance
representations of time series. The quantitative analysis carried out on both
synthetic and sock-shop datasets showcases the effectiveness of our proposed
approach compared to state-of-the-art baselines. For our future work, we
intend to enhance the scalability of our model to accommodate datasets of
varying sizes. This improvement is crucial as it will allow us to not only
build upon the experiments conducted on the sock-shop dataset but also extend
our research to larger-scale and more representative datasets. These larger
datasets will provide a more comprehensive understanding of our model’s
performance in real-world scenarios, further enriching the depth and
applicability of our findings.

## References

  * Chen et al. (2022) Chen, B.; Yu, L.; Luo, W.; Wu, C.; Li, M.; Tan, H.; Huang, J.; and Wan, Z. 2022\.  Hybrid tree model for root cause analysis of wireless network fault localization.  _Web Intell._ , 20(3): 213–223.
  * Chen and He (2021) Chen, X.; and He, K. 2021.  Exploring Simple Siamese Representation Learning.  In _CVPR_ , 15750–15758. Computer Vision Foundation / IEEE.
  * Cheng et al. (2023) Cheng, Y.; Yang, R.; Xiao, T.; Li, Z.; Suo, J.; He, K.; and Dai, Q. 2023.  CUTS: Neural Causal Discovery from Irregular Time-Series Data.  In _ICLR_. OpenReview.net.
  * Daniel Holbach (2022) Daniel Holbach. 2022.  Sock-shop, a microservice demo application.  GitHub repository.
  * Fan et al. (2023) Fan, C.; Wang, Y.; Zhang, Y.; and Ouyang, W. 2023.  Interpretable Multi-Scale Neural Network for Granger Causality Discovery.  In _ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 1–5. IEEE.
  * Gao, Yao, and Chen (2021) Gao, T.; Yao, X.; and Chen, D. 2021.  SimCSE: Simple Contrastive Learning of Sentence Embeddings.  In _EMNLP (1)_ , 6894–6910. Association for Computational Linguistics.
  * Granger (1969) Granger, C. W. 1969.  Investigating causal relations by econometric models and cross-spectral methods.  _Econometrica: journal of the Econometric Society_ , 424–438.
  * Grill et al. (2020) Grill, J.; Strub, F.; Altché, F.; Tallec, C.; Richemond, P. H.; Buchatskaya, E.; Doersch, C.; Pires, B. Á.; Guo, Z.; Azar, M. G.; Piot, B.; Kavukcuoglu, K.; Munos, R.; and Valko, M. 2020.  Bootstrap Your Own Latent - A New Approach to Self-Supervised Learning.  In _NeurIPS_.
  * Ikram et al. (2022) Ikram, A.; Chakraborty, S.; Mitra, S.; Saini, S. K.; Bagchi, S.; and Kocaoglu, M. 2022.  Root Cause Analysis of Failures in Microservices through Causal Discovery.  In _NeurIPS_.
  * Kingma and Ba (2015) Kingma, D. P.; and Ba, J. 2015.  Adam: A Method for Stochastic Optimization.  In _ICLR (Poster)_.
  * Liu et al. (2021) Liu, D.; He, C.; Peng, X.; Lin, F.; Zhang, C.; Gong, S.; Li, Z.; Ou, J.; and Wu, Z. 2021.  MicroHECL: High-Efficient Root Cause Localization in Large-Scale Microservice Systems.  In _ICSE (SEIP)_ , 338–347. IEEE.
  * Ma et al. (2020) Ma, M.; Xu, J.; Wang, Y.; Chen, P.; Zhang, Z.; and Wang, P. 2020.  AutoMAP: Diagnose Your Microservice-based Web Applications Automatically.  In _WWW_ , 246–258. ACM / IW3C2.
  * Nauta, Bucur, and Seifert (2019) Nauta, M.; Bucur, D.; and Seifert, C. 2019.  Causal Discovery with Attention-Based Convolutional Neural Networks.  _Mach. Learn. Knowl. Extr._ , 1(1): 312–340.
  * Oliveira, Miguéis, and de Moura Borges (2023) Oliveira, E. E.; Miguéis, V. L.; and de Moura Borges, J. L. C. 2023.  Automatic root cause analysis in manufacturing: an overview & conceptualization.  _J. Intell. Manuf._ , 34(5): 2061–2078.
  * Shan et al. (2019) Shan, H.; Chen, Y.; Liu, H.; Zhang, Y.; Xiao, X.; He, X.; Li, M.; and Ding, W. 2019\.  ?-Diagnosis: Unsupervised and Real-time Diagnosis of Small- window Long-tail Latency in Large-scale Microservice Platforms.  In _WWW_ , 3215–3222. ACM.
  * Soldani, Forti, and Brogi (2022) Soldani, J.; Forti, S.; and Brogi, A. 2022.  Failure Root Cause Analysis for Microservices, Explained.  In _DAIS_ , volume 13272 of _Lecture Notes in Computer Science_ , 74–91. Springer.
  * Tank et al. (2022) Tank, A.; Covert, I.; Foti, N. J.; Shojaie, A.; and Fox, E. B. 2022.  Neural Granger Causality.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 44(8): 4267–4279.
  * Wang et al. (2023) Wang, D.; Chen, Z.; Ni, J.; Tong, L.; Wang, Z.; Fu, Y.; and Chen, H. 2023.  Hierarchical Graph Neural Networks for Causal Discovery and Root Cause Localization.  _CoRR_ , abs/2302.01987.
  * Wang et al. (2021) Wang, H.; Wu, Z.; Jiang, H.; Huang, Y.; Wang, J.; Köprü, S.; and Xie, T. 2021.  Groot: An Event-graph-based Approach for Root Cause Analysis in Industrial Settings.  In _ASE_ , 419–429. IEEE.
  * Xin, Chen, and Zhao (2023) Xin, R.; Chen, P.; and Zhao, Z. 2023.  CausalRCA: Causal inference based precise fine-grained root cause localization for microservice applications.  _J. Syst. Softw._ , 203: 111724.
  * Yue et al. (2022) Yue, Z.; Wang, Y.; Duan, J.; Yang, T.; Huang, C.; Tong, Y.; and Xu, B. 2022.  TS2Vec: Towards Universal Representation of Time Series.  In _AAAI_ , 8980–8987. AAAI Press.
  * Zeng et al. (2023) Zeng, A.; Chen, M.; Zhang, L.; and Xu, Q. 2023.  Are Transformers Effective for Time Series Forecasting?  In _AAAI_ , 11121–11128. AAAI Press.
  * Zhang et al. (2022) Zhang, C.; Peng, X.; Sha, C.; Zhang, K.; Fu, Z.; Wu, X.; Lin, Q.; and Zhang, D. 2022\.  DeepTraLog: Trace-Log Combined Microservice Anomaly Detection through Graph-based Deep Learning.  In _ICSE_ , 623–634. ACM.


# Multiwavelength investigation of the candidate Galactic PeVatron MGRO
J1908+06

S. Crestan1,2, A. Giuliani2, S. Mereghetti2, L. Sidoli2, F. Pintore3, N. La
Palombara2
1 Università degli Studi dell’Insubria, Via Valleggio 11, 22100 Como, Italy
2 INAF – IASF Milano, Via A. Corti 12, 20133 Milano, Italy
3 INAF – IASF Palermo, Via U. La Malfa 153, 90146 Palermo, Italy
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The candidate PeVatron MGRO J1908+06, which shows a hard spectrum beyond 100
TeV, is one of the most peculiar $\gamma$-ray sources in the Galactic plane.
Its complex morphology and some possible counterparts spatially related with
the VHE emission region, preclude to distinguish between a hadronic or
leptonic nature of the $\gamma$-ray emission. In this paper we illustrate a
new multiwavelength analysis of MGRO J1908+06, with the aim to shed light on
its nature and the origin of its ultra high-energy emission. We performed an
analysis of the 12CO and 13CO molecular line emission demonstrating the
presence of dense molecular clouds spatially correlated with the source
region. We also analyzed 12-years of Fermi-LAT data between 10 GeV and 1 TeV
finding a counterpart with a hard spectrum ($\Gamma\sim 1.6$). Our reanalysis
of XMM-Newton data allowed us to put a more stringent constraint on the X-ray
flux from this source. We demonstrate that a single accelerator cannot explain
the whole set of multiwavelength data, regardless of whether it accelerates
protons or electrons, but a 2-zone model is needed to explain the emission
from MGRO J1908+06. The VHE emission seems most likely the superposition of a
TeV PWN powered by PSR J1907+0602, in the southern part, and of the
interaction between the supernova remnant G40.5-0.5 and the molecular clouds
towards the northern region.

###### keywords:

ISM: individual object: MGRO J1908+06 – ISM: clouds – ISM: cosmic rays – ISM:
supernova remnants

††pubyear: 2020††pagerange: Multiwavelength investigation of the candidate
Galactic PeVatron MGRO J1908+06–References

##  1 Introduction

Cosmic rays (CRs) are high-energy atomic nuclei (90% protons) moving through
space at nearly the speed of light. The local CR proton spectrum is well
described by a power-law up to the so-called ‘knee’ around 1 PeV ($=10^{15}$
eV), indicating the existence of powerful proton accelerators (‘PeVatrons’)
residing in our Galaxy. However, despite decades of efforts, no specific
Galactic source has been securely identified as a proton PeVatron, with the
possible exception of the Galactic Centre (HESS Collaboration et al., 2016; H.
E. S. S. Collaboration et al., 2018; MAGIC Collaboration et al., 2020). The
association of a $\gamma$-ray source with a PeVatron can be indicated by the
absence of an exponential cut-off in the $\gamma$-ray spectrum below 100 TeV,
as expected when CRs are accelerated up to PeV energies and produce
$\gamma$-rays interacting with the ambient material.

The source MGRO J1908+06 is one of the best Galactic PeVatron candidates,
thanks to its hard spectrum reaching energies above 100 TeV as observed by
HAWC (Abeysekara et al., 2020), and with no evidence of a cut-off. It was
discovered by the MILAGRO collaboration (Abdo et al., 2007) and later
confirmed with the HESS atmospheric Cherenkov telescope (Aharonian et al.,
2009). These authors reported the source detection above 300 GeV, with a large
angular size ($\sigma=0.34$°) and a hard spectrum with a photon index of
$2.1$. MGRO J1908+06 was observed at very high energy also with VERITAS (Aliu
et al., 2014), which revealed an extended morphology ($\sigma=0.44$°) with
three peaks of emission and a spectrum with photon index of $2.2$.
Observations of this sky region with HAWC revealed emission up to 100 TeV
spatially consistent with the VERITAS error box of MGRO J1908+06, indicating
the presence of very energetic particles (Abeysekara et al., 2020). Due to its
complex spatial structure, difficult to study with the limited angular
resolution of current instruments, the origin of the $\gamma$-ray emission
from MGRO J1908+06 is still uncertain. As discussed below, a few counterparts
(the supernova remnant (SNR) G40.5-0.5, dense molecular clouds illuminated by
cosmic rays from the nearby SNR, and the pulsar PSR J1907+0602) are compatible
with the $\gamma$-ray source error box ($\sim 0.5$°), preventing a secure
identification of this extreme accelerator and making it difficult to
distinguish between hadronic or leptonic interpretations of its emission. This
source has been studied also at radio wavelength by Duvidovich et al. (2020).
These authors reported the presence of molecular material toward the SNR.
In this work we investigate the molecular cloud environment of MGRO J1908+06
using data from molecular line radio surveys, modelling its multiwavelength
spectral energy distribution using published TeV spectra and our re-analysis
of Fermi-LAT and XMM-Newton data.

Figure 1: VLA Galactic Plane Survey at 1.4 GHz (Stil et al., 2006) map
covering the entire source MGRO J1908+06. The white solid lines are the VHE
contours taken from VERITAS (Aliu et al., 2014), where the significance level
ranges from 3 to 5.2. The yellow and cyan crosses indicate the positions of
PSR J1907+0631 and PSR J1907+0602, respectively. The thin dashed magenta line
indicates the SNR G40.5-0.5 shell, while the red, green and orange dashed
lines mark the VERITAS emission lobes to which we refer in the following.

##  2 Possible counterparts of MGRO J1908+06

Several possible counterparts are present in the region of MGRO J1908+06: in
Fig. 1 we show the continuum emission at 1.4 GHz obtained in the VLA Galactic
Plane Survey (VGPS - (Stil et al., 2006)) with superimposed the TeV contours
obtained by VERITAS (Aliu et al., 2014). The radio data show the shell-like
SNR G40.05-0.5, which is much brighter in the northern region. The positions
of the two pulsars PSR J1907+0602 and PSR J1907+0631 are also indicated.

G40.05-0.5 is a middle aged supernova remnant, with estimated age between 20
and 40 kyr (Downes et al., 1980) and uncertain distance. The $\Sigma-D$
relation gives a distance between 5.5 and 8.5 kpc (Downes et al., 1980). A
more accurate estimate can be obtained from the HI absorption spectrum (see
Ranasinghe & Leahy 2017 for details) but, as reported by Duvidovich et al.
(2020), for SNR G40.5-0.5 this method resulted in a very noisy spectrum,
probably due to the fact that the neutral gas is patchy. The relatively young
and energetic pulsar PSR J1907+0631 (characteristic age $\tau$=11 kyr, spin-
down luminosity $\sim$5$\times 10^{35}$ erg s-1) lies close to the centre of
the SNR (Lyne et al., 2017). Its dispersion measure (DM) implies a distance of
7.9 kpc (Cordes & Lazio, 2002), compatible with the range estimated for
G40.05-0.5 and suggesting an association between these two objects. Duvidovich
et al. (2020) showed that, in principle, PSR J1907+0631 could power the whole
TeV source, since this would require a conversion efficiency from rotational
energy to $\gamma$-rays of about 3%, in line with the efficiency $\leq 10\%$
of other known TeV sources associated to pulsar wind nebulae (Gallant, 2007).
However, this hypothesis is disfavoured by the pulsar position slightly
outside the TeV contours and significantly offset from the centroid of the
$\gamma$-ray emission.

The $\gamma$-ray loud pulsar PSR J1907+0602, discovered with Fermi-LAT (Abdo
et al., 2010), is located in the southern part of MGRO J1908+06, slightly
offset from the peak of the $\gamma$-ray excess counts. This pulsar has a
characteristic age of 19.5 kyr and a spin-down luminosity of $\sim$3$\times
10^{36}$ erg s-1. The source distance was estimated to be 3.2 kpc (Abdo et
al., 2010), as derived from the DM of $\sim 82$ pc cm-3 with the electron
distribution model of Cordes & Lazio (2002).

Another pulsar (PSR J1905+0600, not marked in Fig. 1) lies in this region, but
its distance of $\sim$18 kpc and large characteristic age of 6 Mys (Hobbs et
al., 2004) exclude it as a possible counterpart.

In conclusion, the two main candidates for the $\gamma$-ray emission from MGRO
J1908+06 are the SNR G40.5-0.5, near the northern border of the TeV source,
and PSR J1907+0602, lying in the southern region.

##  3 Data analysis and results

###  3.1 CO Analysis

We have investigated the distribution of the CO gas in the environment of MGRO
J1908+06 in order to identify molecular material in spatial correlation with
the SNR and the $\gamma$-ray emission. Any such association would be relevant
for hadronic models to explain the TeV emission and it would also provide some
information on the source distance. We used the molecular line emission
extracted from the FOREST Unbiased Galactic Plane Imaging (FUGIN)
survey111Available at http://jvo.nao.ac.jp/portal/. This project aims at
investigating the distribution, kinematics, and physical properties of both
diffuse gas and dense molecular clouds in the Galaxy by observing
simultaneously the 12CO, 13CO, and 18CO J=1-0 lines. This survey achieves the
highest angular resolution to date ($\sim$20″) for the Galactic plane, making
it possible to find dense clumps located at farther distances than those seen
in previous surveys.

Figure 2: The 12CO (top) and 13CO (bottom) summed spectra in the region of
MGRO J1908+06. The velocity interval between the two dashed lines (58–78 km
s-1) represents the bulk of the emission, while the red zone marks the
velocity range between 58 and 62 km s-1 (shown in Fig. 3) that is the velocity
range considered for the molecular cloud analysis (section 3.1).

We recovered the spectra in brightness temperature $T_{B}$ as a function of
the local standard of rest velocity ($V_{LSR}$) for the whole region
corresponding to the $3\sigma$ contours of the TeV emission, both in 12CO and
13CO. As shown in Fig. 2, the bulk of the emission is concentrated between 50
and 80 km s-1.

We plot in Fig. 3 the 12CO and 13CO molecular line emission integrated from 58
to 62 km s-1. The contours presented in the figure are those of the VERITAS
TeV emission (Aliu et al., 2014) and of the SNR G40.5-0.5 at 1.4 GHz from the
VGPS. We denote the three maxima of $\gamma$-ray emission as lobes A, B, and C
(see Fig. 1). The maps of Fig. 3 show that lobe A overlaps with CO emission,
lobe B partially overlaps with CO emissions, while no obvious molecular clouds
association is seen for lobe C.

Figure 3: Maps of 12CO (left) and 13CO (right) emission in the MGRO J1908+06
region integrated between 58–62 km s-1. The white solid lines are the same as
in Fig. 1, while the green contours are the continuum emission from SNR
G40.5-0.5 at 1.4 GHz.

We concentrate on the molecular cloud in the 58–62 km s-1 velocity interval,
as it overlaps both the A-B lobes and the southern border of the SNR. We
obtain the distance of the cloud using the Galaxy rotation curve from Clemens
(1985), with $R_{\sun}$= 8.5 kpc and $v_{\sun}=220$ km s-1. The first Galactic
quadrant presents distance ambiguity for positive radial velocities, so
adopting 60 km s-1, we obtain near and far distances of 3.0 and 9.4 kpc,
respectively.

To study the properties of the molecular gas, and in particular to estimate
their density, we use the dendrogram technique (Rosolowsky et al., 2008). A
dendrogram is a topological representation of the significant local maxima in
N-dimensional intensity data and the way these local maxima are connected
along contours (or isosurfaces) of constant intensity. A local maximum, by
definition, has a small region around it containing no data greater than its
value and, hence, a distinct isosurface containing only that local maximum can
be drawn. The local maxima determines the top level of the dendrogram, which
we refer to as the “leaves” , defined as the set of isosurfaces that contain a
single local maximum. We identify and characterize molecular clouds in the CO
data cube between 58–62 km s-1 using astrodendro222Available at
http://www.dendrograms.org/. This python algorithm efficiently constructs a
dendrogram representation of all the emission in the selected region. The
minimum value to consider (any value lower than this will not be considered in
the dendrogram) is set as the “detection level”, namely 5$~{}\sigma_{T}$,
where $\sigma_{T}$ is the median RMS noise level in the dataset, so that only
significant values are included in the dendrogram (T${}_{\rm min}=3$ K).
Another consideration is about how significant a leaf has to be in order to be
considered an independent entity. The significance of a leaf is measured from
the difference between its peak flux and the value at which it is being merged
into the tree. This parameter is set to 1 $\sigma_{T}$, which means that any
leaf that is locally less than 1 $\sigma_{T}$ high is combined with its
neighboring leaf (or branch) and is no longer considered as a separate entity.

Once an index of structures in the data has been produced by the algorithm, it
can be used to catalog the properties of each structure, such as integrated
intensity, centroid position, spatial position angle, spatial extent, and
spectral line-width.

We estimate the luminosity based on the zeroth moment, i.e., the sum of the
intensity, and then translate the moments into estimates of physical
quantities. For these calculations, we consider the pixels in a cloud mask
$\mathcal{M}$, i.e. only the pixels belonging to a single cloud identified by
the segmentation algorithm. We measure the luminosity of each cloud as:

$L_{\rm CO}=A_{\rm pix}\Delta v\sum_{i}T_{i}$ (1)

where $A_{\rm pix}$ is the projected physical area of a cube pixel in pc2,
$\Delta v$ = 4 km s-1 is the channel width, and $T_{i}$ is the brightness of
the cube pixels measured in K in the cloud mask $\mathcal{M}$. We convert from
luminosity to mass, scaling the extrapolated luminosity through the CO-to-H2
conversion factor, $\alpha_{\rm CO}$.

$M_{\rm CO}=L_{\rm CO}\alpha_{\rm CO}$ (2)

where we take $\alpha_{\rm CO}$= 4.35 $M_{\sun}$ pc-2 (km s-1 K)-1 at solar
metallicity (Bolatto et al., 2013). To measure cloud radii we convert from the
deconvolved major and minor sizes, $\sigma_{\rm maj}$ and $\sigma_{\rm min}$,
to a cloud radius measurement using:

$R=\eta\sqrt{\sigma_{\rm maj}\sigma_{\rm min}}$ (3)

The factor $\eta$ depends on the light or mass distribution within the cloud.
We adopt $\eta$=1.91 following Rosolowsky & Leroy (2006). Our model
approximates the cloud as a spherically symmetric object so that R also
characterizes the object in three dimensions. Therefore, we do not apply any
inclination corrections to R. The resulting mean cloud density is $\sim 180$
particles cm-3 assuming a distance of 3 kpc, while it is $\sim 60$ particles
cm-3 assuming 9 kpc.

###  3.2 Fermi-LAT data analysis

We analyzed 12 years of Fermi-LAT data, obtained from 2008-09-01 to
2020-12-16, exploiting the Pass 8 data processing (P8R3) with the public
fermitools (v2.0.0) and fermipy packages (v1.0.0). We selected the Pass 8
‘source’ class and ‘front+back’ type events coming from zenith angles smaller
than 90° and from a circular region of interest (ROI) with radius of 10°
centered at R.A. = 286.97° and Dec. = 6.03° (J2000). The instrument response
function version P8R3-SOURCE-V3 was used. We selected only the events in the
10 GeV–1 TeV energy range, to avoid the contribution from PSR J1907+0602 (see
fig. 4 of Abdo et al. 2010). We included in the background model all the
sources from the 4FGL catalog within the ROI, as well as the Galactic (gll-
iem-v07.fits) and the isotropic (P8R3-SOURCE-V3-v1) diffuse components.

We performed a binned analysis with five bins per energy decade and spatial
pixel size of 0.05° . In the maximum likelihood fitting, the normalization
parameter of all the sources within 3° of the ROI centre, as well as the
diffuse emission components, were left free to vary. Instead the parameters of
all the other sources at more than 3° were fixed to the values given in the
4FGL catalog (Abdollahi et al., 2020). To describe the spatial morphology of
MGRO J1908+06, we used the VERITAS emission region at 3$\sigma$ level (i.e.
the outermost contour in Fig. 1 ), while for the spectral model we assumed a
power law with photon index $\Gamma=1.6$. This leads to a detection
significance $\sqrt{TS}\sim 6$ in the energy band considered. The $\gamma$-ray
flux was obtained by binning the $\gamma$-ray data in the range from 10 to
1000 GeV into four energy intervals, and performing a binned likelihood
analysis in each energy bin. The resulting Fermi-LAT spectral energy
distribution is plotted in Fig. 4.

###  3.3 X-ray Analysis

To study the X-ray emission in the vicinity of PSR J1907+0602 we used a 52 ks
long observation carried out on 2010 April 26 with the XMM-Newton satellite.
We analyzed the data of the EPIC-MOS instrument that was operated in full
frame imaging mode and with the medium thickness optical filter. We excluded
time intervals with high background, resulting in net exposure times of 36 and
38 ks for the two MOS cameras.

Using the Extended Source Analysis Software
(ESAS333https://heasarc.gsfc.nasa.gov/docs/xmm/esas/cookbook/xmm-esas.html),
we extracted the spectra from a circular region of 5 arcmin radius centered at
the position of PSR J1907+0602 (excluding a circle of 30″radius around the
source) and from a concentric annular region with radii 5 and 12.5 arcmin. The
latter was used to estimate the X-ray background (which in this sky region is
dominated by the Galactic Ridge diffuse emission). Comparison of the two
spectra showed no evidence for diffuse emission associated with PSR
J1907+0602, with an upper limit (at 95% c.l.) of 1.2$\times 10^{-15}$ erg cm-2
s-1 arcmin-2 on the surface brightness in the 1–10 keV energy range.

##  4 Origin of the $\gamma$-ray emission

Emission at TeV energies indicates the presence of ultra-relativistic
particles which, in principle, can produce it through Inverse Compton (IC)
scattering of the CMB, IR and/or star-light seed photons by electrons, or
through the decay of neutral pions resulting from proton-proton (and/or other
nuclei) interactions. In this Section, we first explore the possibility that a
single mechanism is responsible for the emission from the whole trilobed
region in either the leptonic or hadronic scenario.

We then consider the possibility of a two-zone model, in which both components
(hadronic and leptonic) are present. This scenario can be originated by a
source able to efficiently accelerate both protons and electrons, or by 2
sources lying in the same sky region.

In the following, we use the Fermi-LAT and XMM-Newton results derived as
described above and the TeV spectra obtained with VERITAS (Aliu et al., 2014)
and MILAGRO (Abdo et al., 2007). We have not used the HAWC data reported in
(Abeysekara et al., 2020) for our fit, because they are inconsistent with both
the VERITAS and HESS data in an energy range (1-10 TeV) where these
instruments proved to show reliable results. Also HESS data were not
considered in the fitting procedure as they are fully compatible with the
VERITAS ones. This discrepancy is explained in (Abeysekara et al., 2020) as
consequence a the larger source extent observed by HAWC. However, we note that
our model reproduces quite well the slope of the HAWC spectrum at energies of
around 100 TeV, hence in the range where there is consistency between the
results reported in (Abeysekara et al., 2020) and (Abeysekara et al., 2017).

The X-ray upper limit, derived for a region of about 75 arcmin2, was rescaled
for the larger area enclosed in the 3$\sigma$ contours of the TeV emission
($\sim$ 1400 arcmin2). This is a conservative assumption, because any diffuse
X-ray emission produced by particles accelerated by the pulsar would likely
decrease in intensity at larger distances.

###  4.1 1-zone leptonic Hypothesis

In the leptonic scenario, we assume that the whole emission from MGRO J1908+06
is due to a population of relativistic electrons interacting, through IC, with
the ISRF photons. These electrons can be supplied by the energetic pulsar PSR
J1907+0602 or by the SNR G40.5-0.5.

In order to infer the properties of the parent particle distribution, we fit
the multiwavelength SED through a Markov Chain Monte Carlo (MCMC) procedure
using the naima python package (see Zabalza 2015 for a detailed description of
the fitting procedure and of the IC radiative model). We modeled the three
dominant photon fields with energy densities fixed at $\epsilon_{\rm CMB}$ =
0.261 erg cm-3 for the Cosmic Microwave Background, $\epsilon_{\rm FIR}$ = 0.5
erg cm-3 for the far-infrared dust emission, and $\epsilon_{\rm NIR}$ = 1.0
erg cm-3 for the near-infrared stellar emission. The electron distribution was
modeled as a broken power law with an exponential cut-off with all the
parameters free to vary during the fitting procedure.

We obtained a best fit (see Fig. 4) with the parameters for the electron
distribution shown in Table 1.

The high-energy electrons responsible for the TeV emission interact also with
the ambient magnetic field producing synchrotron radiation. We found that the
electron population obtained in our best fit can be reconciled with the XMM-
Newton upper limit in the X-ray band only for ambient magnetic fields (B)
smaller than 1.2 $\mu$G.

This limit is rather small compared to the typical value of the Galactic
magnetic field of 5 $\mu$G (Haverkorn, 2015). If we assume a magnetic field of
5 $\mu$G and the same spectral slope of the best fit, the normalization of the
electrons population must be reduced by a factor $\sim$15 to be consistent
with the X-ray upper limit (see Fig. 5). This results shows that a 1-zone
leptonic model alone cannot explain the whole multiwavelength set of data.

Figure 4: Multiwavelength SED of MGRO J1908+06 fitted with the leptonic
scenario (Inverse Compton + synchrotron) with B=1.2 $\mu$G. The black line
shows the total IC emission, whose components are indicated by the thin solid,
dashed and dotted black line for the CMB, FIR and NIR, respectively. The red
solid line is for the synchrotron emission. The XMM-Newton upper limit (black
arrow) as well as the Fermi-LAT data (red squares) are from this work. The
other data are taken from Aliu et al. (2014) (VERITAS - green diamonds) and
Abdo et al. (2007) (MILAGRO - pink hexagons). The blue and gray butterflies
are the HAWC models from Abeysekara et al. (2020) and from Abeysekara et al.
(2017), plotted for comparison, but not used in the fit.

Figure 5:  Same as Fig. 4, but with B=5 $\mu$G and the electron population
normalization reduced by a factor 15.

.

###  4.2 1-zone hadronic hypothesis

In the hadronic scenario we assume that a population of relativistic protons
interacts with dense interstellar material and produces TeV photons via pion
decay. A good candidate for the acceleration of these protons is the SNR
G40.5-0.5. In fact, our analysis of the molecular gas around the SNR
demonstrated the presence of molecular clouds in good spatial correlation with
the SNR shell. We used naima to fit the $\gamma$-ray SED, assuming a proton
distribution described by a broken power law with exponential cut-off and pion
decay as radiative model. The parameters of the best fit, shown in Fig. 6, are
given in Table 1.

With the average densities of the clouds derived in section 3.1, $180$ or $60$
particles cm-3 depending on the considered distance, the total proton energy
required by the fit is $7\times 10^{47}$ erg or $2\times 10^{49}$ erg,
respectively.

Figure 6: MGRO J1908+06 hadronic emission model. The black line shows the Pion
Decay best fit model. Data are the same as in Fig. 4

###  4.3 2-zone model

Considering the complex spatial distribution of the VHE emission, it can not
be excluded that both a TeV PWN powered by PSR J1907+0602 and hadronic
processes associated to the SNR G40.5-0.5 contribute to the $\gamma$-ray
emission observed from this sky region. Therefore, we have also explored a
hybrid emission model in which the TeV emission is due to the superposition of
leptonic and hadronic components from these two sources. Of course, in this
scenario, MGRO J1908+06 might consist of two physically separated sources, not
necessarily at the same distance.

Figure 7: MGRO J1908+06 2-zone emission model. The black line shows the total
emission while the gray lines and the orange lines are for the IC and pion
decay emission respectively. Other colors are the same as in Fig. 4. The
magnetic field is B=4 $\mu$G.

To fit the spectral energy distribution we assume that the steep spectrum at
GeV energies has a leptonic origin, while the hadronic emission is responsible
for the softer part at TeV energies. We used naima to recover the radiative
models from a particle distribution. We used an exponential cut-off broken
power law both for electrons and protons. The resulting emission model is
plotted in Fig. 7, while the assumed parameters are reported in Table 1. The
recovered total protons energy for distances of 3 and 9 kpc is $4\times
10^{47}$ erg and $10^{49}$ erg, respectively, while the electrons energy is
$9\times 10^{46}$ erg.

Model | Component | d | $\Gamma_{1}$ | $\Gamma_{2}$ | $W$ | $E_{0}$ | $E_{b}$ | $E_{c}$
---|---|---|---|---|---|---|---|---
|  | (kpc) |  |  | (erg) | (TeV) | (TeV) | (PeV)
1-zone | Leptonic | 3 | 1.0 $\pm$ 0.4 | 2.6 $\pm$ 0.1 | 2$\times 10^{47}$ | 10 | 2.7 $\pm$ 0.7 | 7.1 $\pm$ 6.0
1-zone | Hadronic | 3 | 1.0 $\pm$ 0.1 | 2.1 $\pm$ 0.1 | 7$\times 10^{47}$ | 30 | 2.8 $\pm$ 0.8 | 3.0 $\pm$ 0.9
1-zone | Hadronic | 9 | 1.1 $\pm$ 0.2 | 2.1 $\pm$ 0.1 | 2$\times 10^{49}$ | 30 | 3.4 $\pm$ 1.2 | 1.9 $\pm$ 0.5
2-component | Leptonic | 3 | 1.2 | 1.2 | 9$\times 10^{46}$ | 10 | 0.2 | 0.011
| Hadronic | 3 | 1.6 | 2.0 | 4$\times 10^{47}$ | 30 | 200 | >1
| Hadronic | 9 | 1.6 | 2.0 | 1$\times 10^{49}$ | 30 | 200 | >1
Table 1: Parameters for all emission models considered. $\Gamma_{1}$ and
$\Gamma_{2}$ are the indices before and after the break ($E_{b}$), W is the
particles total energy, $E_{0}$ is the reference energy while $E_{c}$ is the
cut-off energy.

##  5 Discussion

The association of MGRO J1908+06 with an offset relic PWN driven by PSR
J1907+0602 was initially considered the most likely origin of the VHE (Abdo et
al., 2010). Our results show that, for reasonable values of the ambient
magnetic field, a leptonic emission model fitting the $\gamma$-ray spectrum
from the whole source would produce a synchrotron X-ray flux incompatible with
the upper limit in the few keV region. The leptonic interpretation is also
disfavored by the spatial shape of the TeV emission, extending far from the
pulsar position and without evident sign of spectral softening with distance
from the pulsar, as it would be expected from electron cooling (Aliu et al.,
2014). We further note that, due to the Klein-Nishina suppression of the IC
cross section at high energies, a rather large value of the electron maximum
energy is required to fit the $\gamma$-ray spectrum.

Our analysis of the molecular gas around SNR G40.5-0.5 demonstrates the
presence of molecular clouds in good spatial correlation with the SNR (see
also Duvidovich et al. 2020) and motivates the exploration of a hadronic
scenario. The cloud densities required by the best fit proton distribution
used to reproduce the observed $\gamma$-ray spectrum are consistent with the
ones that we derived from an analysis of the CO data, independent of the
source association with the near (3 kpc) or far (9 kpc) source distance.
However, the VHE emission of MGRO J1908+06 extends beyond the spatial
distribution of the target material, with no obvious molecular cloud
counterparts in the southern region. The very hard photon index $\Gamma_{1}$
required to fit the Fermi-LAT data (not seen in other TeV sources associated
with SNRs) also disfavors a fully hadronic model for MGRO J1908+06.

The difficulties discussed above for single zone models are easily solved
assuming that the MGRO J1908+06 spectrum is a sum of two components (hadronic
and leptonic). The available data do not allow us to perform a spatially-
resolved spectral analysis. Thus we cannot exclude that these two components
arise from different zones of the source (for example, a TeV PWN powered by
PSR J1907+0602 could be responsible for the southern lobe, and the interaction
between the SNR G40.5-0.5 and the molecular clouds for the northern part).

An important result of this analysis is that for both the 1-zone hadronic
model and the 2-component model the maximum energy of the emitting protons
required to fit the gamma-ray spectrum is greater than 1 PeV.

##  6 Conclusions

Our multiwavelength modeling of MGRO J1908+06 confirms that this source is one
of the best Galactic PeVatron candidates. We found that single-zone models,
although in principle justified by the presence of plausible counterparts in
both the leptonic and hadronic scenarios (a pulsar wind nebula powered by PSR
J1907+0602 or SNR G40.05-0.5 interacting with molecular clouds, respectively),
run into problems to explain the multiwavelength and spatial morphology
properties of MGRO J1908+06. Therefore a 2-zone model is preferred to describe
the emission from this source. Spatially resolved data, as those that will be
provided by the next generation of Cherenkov telescopes such as the upcoming
ASTRI Mini-array and the Cherenkov Telescope Array (CTA), are needed to
separate the emission components of this source.

## Acknowledgements

This work made use of data from FUGIN, FOREST Unbiased Galactic Plane Imaging
survey with the Nobeyama 45-m telescope, a legacy project in the Nobeyama 45-m
radio telescope.

## DATA Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abdo et al. (2007) Abdo A. A., et al., 2007, ApJ, 664, L91
  * Abdo et al. (2010) Abdo A. A., et al., 2010, ApJ, 711, 64
  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Abeysekara et al. (2017) Abeysekara A. U., et al., 2017, The Astrophysical Journal, 843, 40
  * Abeysekara et al. (2020) Abeysekara A. U., et al., 2020, Phys. Rev. Lett., 124, 021102
  * Aharonian et al. (2009) Aharonian F., et al., 2009, A&A, 499, 723
  * Aliu et al. (2014) Aliu E., et al., 2014, ApJ, 787, 166
  * Bolatto et al. (2013) Bolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, 51, 207
  * Clemens (1985) Clemens D. P., 1985, ApJ, 295, 422
  * Cordes & Lazio (2002) Cordes J. M., Lazio T. J. W., 2002, arXiv e-prints, pp astro–ph/0207156
  * Downes et al. (1980) Downes A. J. B., Pauls T., Salter C. J., 1980, A&A, 92, 47
  * Duvidovich et al. (2020) Duvidovich L., Petriella A., Giacani E., 2020, MNRAS, 491, 5732
  * Gallant (2007) Gallant Y. A., 2007, Astrophysics and Space Science, 309, 197–202
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration et al., 2018, A&A, 612, A9
  * HESS Collaboration et al. (2016) HESS Collaboration et al., 2016, Nature, 531, 476
  * Haverkorn (2015) Haverkorn M., 2015, in Lazarian A., de Gouveia Dal Pino E. M., Melioli C., eds, Vol. 407, Magnetic Fields in Diffuse Media. p. 483 (arXiv:1406.0283), doi:10.1007/978-3-662-44625-6_17
  * Hobbs et al. (2004) Hobbs G., et al., 2004, MNRAS, 352, 1439
  * Lyne et al. (2017) Lyne A. G., et al., 2017, ApJ, 834, 137
  * MAGIC Collaboration et al. (2020) MAGIC Collaboration et al., 2020, A&A, 642, A190
  * Ranasinghe & Leahy (2017) Ranasinghe S., Leahy D. A., 2017, ApJ, 843, 119
  * Rosolowsky & Leroy (2006) Rosolowsky E., Leroy A., 2006, PASP, 118, 590
  * Rosolowsky et al. (2008) Rosolowsky E. W., Pineda J. E., Kauffmann J., Goodman A. A., 2008, ApJ, 679, 1338
  * Stil et al. (2006) Stil J. M., et al., 2006, The Astronomical Journal, 132, 1158–1176
  * Zabalza (2015) Zabalza V., 2015, in 34th International Cosmic Ray Conference (ICRC2015). p. 922 (arXiv:1509.03319)