# Machine learning on knowledge graphs for context-aware security monitoring
††thanks:  This work was partially funded by the Federal Ministry for Economic
Affairs and Energy of Germany (BMWi) within the IIP-Ecosphere Project.
††thanks:  †These authors contributed equally to this work.

Josep Soler Garrido†, Dominik Dold† and Johannes Frank  Siemens Technology
Munich, Germany
{josep.soler_garrido, dominik.dold<EMAIL_ADDRESS>
###### Abstract

Machine learning techniques are gaining attention in the context of intrusion
detection due to the increasing amounts of data generated by monitoring tools,
as well as the sophistication displayed by attackers in hiding their activity.
However, existing methods often exhibit important limitations in terms of the
quantity and relevance of the generated alerts. Recently, knowledge graphs are
finding application in the cybersecurity domain, showing the potential to
alleviate some of these drawbacks thanks to their ability to seamlessly
integrate data from multiple domains using human-understandable vocabularies.
We discuss the application of machine learning on knowledge graphs for
intrusion detection and experimentally evaluate a link-prediction method for
scoring anomalous activity in industrial systems. After initial unsupervised
training, the proposed method is shown to produce intuitively well-calibrated
and interpretable alerts in a diverse range of scenarios, hinting at the
potential benefits of relational machine learning on knowledge graphs for
intrusion detection purposes.

###### Index Terms:

artificial intelligence, cybersecurity, industrial control systems, knowledge
graphs, machine learning

††publicationid: pubid: 978-1-7281-5684-2/20/$31.00 ©2021 IEEE

85(17.5,256.75) ©2021 IEEE. Personal use of this material is permitted.
Permission from IEEE must be obtained for all other uses, in any current or
future media, including reprinting/republishing this material for advertising
or promotional purposes, creating new collective works, for resale or
redistribution to servers or lists, or reuse of any copyrighted component of
this work in other works.

##  I Introduction

Artificial Intelligence (AI) methods and techniques such as Machine Learning
(ML) are bound to play an increasingly prominent role in the cybersecurity
domain [1]. Potential applications are found both in defensive and offensive
scenarios across most cybersecurity functions [2, 3]. Among these, detection
tasks are particularly poised to benefit from the ability to automatically
analyze and learn from vast quantities of data. Indeed, there are many
relevant examples of the application of deep learning and similar techniques
for intrusion detection systems (IDS) [4], based either on supervised learning
techniques targeting detection of known malicious signatures, or on anomaly
detection methods able to detect unknown attack patterns and find deviations
from a previously learned baseline. However, it has been argued that, while
advantageous in terms of detection capabilities and data requirements, the
latter class of machine learning-based IDS systems must address some drawbacks
in order to improve their usefulness [4], including, among others, their
tendency to produce an unmanageable amount of alerts [2, 5], their limited
ability to learn from and interact with human operators [5, 6], and the lack
of interpretability of the alerts generated [7].

Addressing these drawbacks would be of great benefit not only to defend
conventional IT systems, but also in the context of modern operational
technology (OT) systems, such as those used in factories and other industrial
automation settings. Indeed, while industrial control systems have
traditionally been extremely deterministic in their operation—and their
security has been at least partially implemented by means of isolation
[8]—modern Industry 4.0 automation settings are characterized by a convergence
of IT and OT infrastructure [9], displaying increasingly complex activity
patterns and network topologies and making extensive use of autonomous systems
and components (e.g. robots, or AI-enabled software applications) that
constantly interact with one another. This new paradigm has the potential to
substantially improve the flexibility, reliability and efficiency of
industrial systems, but also poses new cybersecurity challenges [10] and
demands a high degree of domain-specific knowledge from analysts assessing
potential integrity issues or indications of security compromises.

It is in this context that we propose the application of machine learning on
knowledge graphs in order to improve the quality and relevance of IDS-
generated alerts in modern industrial systems, increasing their usefulness for
human operators. Intuitively, this can be achieved by leveraging the semantic
integration of domain-specific data on a single knowledge base, which enables
us to better contextualize and enrich cybersecurity-relevant observations, and
allows machine learning methods to leverage this additional context, i.e., to
learn from these observations in a way that makes use of the rich set of
interconnections and relations between different entities in the graph. In our
case, this is shown to result in a model of the behavior of the different
actors in an industrial automation system that exhibits particularly useful
generalization properties and can be employed to reliably detect and score
unexpected activity patterns.

The rest of the paper is organized as follows. Section II introduces knowledge
graphs and presents a brief overview of existing applications in the
cybersecurity domain. In section III we describe the graph learning algorithm
employed for our specific security monitoring application. Section IV
describes an empirical evaluation setup based on an industrial automation
hardware prototype. The experiments performed and results obtained are
presented and discussed in section V, before finally drawing some conclusions
in section VI.

##  II Background and related work

###  II-A Knowledge graph model

A knowledge graph is a specific type of knowledge base where information is
encoded in the form of a directed labeled graph, with nodes representing
entities and edges representing different types of possible relationships
between entities. The information contained in a knowledge graph can be
alternatively represented as a list of statements with the form subject-
predicate-object, or $\\{s,p,o\\}$ for short, where $s$ corresponds to an
entity in the graph, $p$ to a specific relationship, and $o$ can either be an
entity or a literal, e.g., a specific data value. A simple knowledge graph
illustrating these two alternative representations is depicted in Figure 1.

Knowledge graphs are particularly useful structures to integrate data from
multiple areas of knowledge, typically making use of domain-specific
vocabularies and ontologies that model the different categories,
relationships, rules and constraints in a specific area of knowledge.
Sophisticated data discovery and analytics tasks can be performed on the
resulting graphs, e.g., by means of semantic queries using languages such as
SPARQL in RDF graphs [11]. Use of machine learning methods is also possible on
knowledge graphs, typically by means of so-called node embeddings: vector
representations of graph entities which are more suitable for processing via
neural networks and similar methods than their original symbolic
representations.

Many different methods have been proposed to learn embeddings [12]. A specific
method used in the context of our work is RESCAL [13]. Assuming a knowledge
graph with $n$ entities and $m$ different types of relations, a suitable
representation is given by $\bf{X}$, a binary tensor of size $n\times m\times
n$ where entry $X_{s,p,o}$ contains a $1$ if the corresponding triple
$\\{s,p,o\\}$ is present in the graph and $0$ otherwise. Given such a
representation, RESCAL formulates the graph embedding problem as a tensor
factorization operation

$X_{s,p,o}\,{\approx}\,\bm{e}^{T}_{s}\bm{R}_{p}\bm{e}_{o},$ (1)

where $\bm{e}_{j}$ is the $R\times 1$ embedding vector of entity $j$ and
$\bm{R}_{k}$ is a $R\times R$ matrix modelling the interactions between entity
embeddings for the k-th relationship type.

Embeddings are obtained through minimization of the reconstruction loss

$L_{\mathrm{MSE}}=\sum_{s,p,o}\|X_{s,p,o}-\bm{e}^{T}_{s}\bm{R}_{p}\bm{e}_{o}\|^{2}\,,$
(2)

which can be performed, for instance, by stochastic gradient descent (SGD),
sequentially sampling batches of positive and negative (i.e., not present in
the graph) triple statements in order to iteratively update the learned
representations. Once the learning step is completed, the resulting embeddings
can be used for a wide range of downstream applications. For instance, in a
knowledge graph completion task, non-observed relationships can be inferred.
In the RESCAL case this is done by means of the score
$\theta_{s,p,o}=\bm{e}^{T}_{s}\bm{R}_{p}\bm{e}_{o}$, which captures the
likelihood that, given the set of observed triples, $\\{s,p,o\\}$ also
represents a true fact.

Figure 1: Simple graph describing some elements in an industrial network, and
list of statements represented in the graph.

Figure 2: Hardware demonstrator used for evaluation (left) and simplified
network diagram including main traffic flows (right).

###  II-B Cybersecurity applications

Analysis of security incidents typically requires consideration of multiple
data sources, some of which are often exchanged between organizations. In
order to facilitate this, common schemas and data representation formats, such
as STIX [14], have been introduced that enable organizations to exchange
threat intelligence in a consistent way. More recently, these have evolved
into fully-fledged ontologies enabling inference and reasoning [15]. These
ontologies model a wide range of cybersecurity-relevant knowledge such as
product information, known vulnerabilities and attack patterns, and can
additionally be linked to domain-specific knowledge, e.g. coming from
industrial automation systems [16].

Construction of high-quality knowledge graphs is a challenging task,
especially when it requires extraction of information from unstructured
textual or heterogeneous data. Therefore, the use of machine learning to
perform automatic extraction of cybersecurity knowledge from various sources
has been an area of strong academic interest [17, 18].

Once constructed, these knowledge graphs find a wide range of applications,
e.g., vulnerability assessment [19], intrusion and threat detection [19, 17,
20], malware analysis [21], and other forms of AI-enabled assistance to
security analysts [22]. These applications make use of various methods in
order to extract insights from cybersecurity knowledge graphs, ranging from
simple SPARQL queries to complex rule-based systems, ranking algorithms or
logical reasoning methods.

##  III Learning method

We apply machine learning on knowledge graphs to detect unexpected activity in
industrial automation systems. Multiple studies have considered detection of
anomalies using graph data [23]. More recently, several of these methods have
focused on multi-relational data, as found in knowledge graphs [24, 25, 26].
In fact, knowledge graph completion methods, able by definition to infer
missing facts, can be readily applied to anomaly detection problems [27]. That
is the approach we follow, learning a baseline of normal behavior by training
a generative graph embedding algorithm on a graph built from a training
dataset, and using it thereafter in a link-prediction setting to rank the
likelihood of triple statements resulting from events observed at test time
and determine whether they represent a significant deviation from the
baseline.

We employ a tensor factorization method inspired by RESCAL. In our case, a
novel probabilistic formulation is employed [28], where the joint probability
of a graph represented by $\bm{X}$ is modelled as

$\displaystyle p(\boldsymbol{X})=\frac{1}{Z}\,e^{-E(\boldsymbol{X})}\,,$ (3)

where $E(\boldsymbol{X})$ is an energy function given by

$E(\boldsymbol{X})=-\sum_{s,p,o}X_{s,p,o}\,\theta_{s,p,o}\,,$ (4)

and $Z$ is a normalization constant given by a summation over all possible
graph realizations

$Z=\sum_{\boldsymbol{X}^{\prime}}e^{-E(\boldsymbol{X}^{\prime})}\,.$ (5)

Such a stochastic description of graph data allows us to handle data streams
where identical events appear multiple times or underlie statistical
variations—something that is not covered by traditional graph embedding
algorithms like RESCAL. By employing wake-sleep learning [29], where the model
learns in an unsupervised mode by only observing positive triples, embeddings
are found that assign low energies to triple statements consistent with the
training graph, and high energies otherwise. This energy-based formulation is
experimentally shown to produce a broad range of triple scores
$\theta_{s,p,o}$ compared to standard RESCAL factorization and other similar
graph embedding methods, which in our case have been observed to result in a
more bimodal distribution of scores, i.e., they tend to assign either very low
or very high scores to potentially anomalous triples. Having a system that
produces a wide range of scores is beneficial in the proposed IDS context, as
it allows human analysts to better prioritize alerts. These benefits are
further illustrated in the next sections which describe a thorough
experimental evaluation campaign.

##  IV Experimental setup

The hardware prototype is described in Figure 2, following the design of
modern industrial systems integrating IT and OT elements. The automation side
is equipped with a Siemens S7-1500 PLC connected to peripherals via an
industrial network. These include a drive subsystem controlling the motion of
a conveyor belt, an industrial camera, a human-machine-interface (HMI), and a
distributed I/O subsystem with modules interfacing with various sensors for
object positioning and other measurements. The PLC exposes values reported by
these sensors as well as information about the state of the system by means of
an OPC-UA server [30].

The variables exposed by the server are consumed on the IT part of the
demonstrator by applications hosted on edge computing servers, i.e., computing
infrastructure directly located at the shop floor which is typically devoted
to data-driven tasks that require close integration with the automation system
and short response times, such as real-time system monitoring, fault
prediction or production optimization.

Industrial edge applications have dynamic lifecycles, and this is captured in
the prototype by recreating a development environment. This cycle starts with
development hosts consuming potentially high volumes of data from a historian,
a database that constantly stores process data from the automation system.
Thereafter, containerized applications—ostensibly developed on the basis of
the consumed data—are pushed into an application repository. Finally, edge
computing hosts fetch application updates periodically. To make the behavior
more realistic, development hosts occasionally access the internet with low
traffic volumes. The environment is fully virtualized and performs these
activities in an autonomous manner, with an option to manually induce
different types of anomalous behaviors in order to test the response of our
IDS system.

A knowledge graph is built out of the running prototype by integrating three
main sources of knowledge: information about the automation system,
observations at the network level (e.g., connections between hosts), and
observations at the application level (e.g., data access events). A sizeable
portion of the information is related to the automation system, which is
extracted from engineering tools in the AutomationML format and ingested into
the graph using a readily available ontology [31]. Information about
application activity is obtained from the OPC-UA server logs, including
session information, which variables are accessed and in which way. Finally,
all network traffic is passed through the security monitoring tool Zeek [32],
which produces a stream of observed connections that are ingested using a
simple custom data model.

(a) sub-graph containing variable access events

(b) sub-graph containing network connections

Figure 3: Illustration of two evaluation scenarios, including baseline activity observed during training (black edges) and deviations of different severity (dashed colored edges: Expected, Unexpected, Suspicious, Highly Suspicious). (a) Expected scores for data-access events, depending on the context provided by observed ones, e.g., which modules or subsystems previously accessed variables are related to. (b) Expected scores for new HTTPS connections (4–5), depending on the context provided by observed ones (1–3), e.g. whether hosts of the same type are involved in similar traffic patterns. Test case | Scenario description | Expected severity
---|---|---
Variable access | App accesses the same variables as during training. | Observed
App accesses variables from the same (or a closely related) device or module as those observed during training. | Expected
App changes the way it accesses some variables (e.g. writes instead of reads). | Unexpected
App accesses variables completely unrelated to those accessed during training. | Highly Suspicious
HTTPS access | A dev. host makes an HTTPS access (internet or historian) as observed during training for the same host. | Observed
A dev. host makes an HTTPS access for the first time, but in a way consistent with the actions of other dev. hosts. | Expected
A dev. host observed to access the internet during training accesses a previously unseen public IP address. | Expected
A local address not corresponding to a dev. host (e.g. an edge server) accesses the historian. | Suspicious
A local address not corresponding to a dev. host (e.g. an edge server) accesses a public IP address. | Highly Suspicious
A high-volume HTTP access is made to a public IP address (high volumes only from historian in baseline). | Highly Suspicious
SSH access | A dev. host makes an SSH access to the app repository as seen during training. | Observed
A dev. host makes an SSH access to the app repository for the first time, but exactly like other dev. hosts. | Expected
A dev. host accesses the app repository via SSH as during training but with a slightly higher data volume. | Expected
The historian host (not a dev. host but on the same network) accesses the app repository via SSH. | Unexpected
A dev. host accesses an edge server via SSH, but during training no edge servers received SSH connections. | Suspicious
SSH connection between two edge servers. During training no edge servers started or received SSH connections. | Highly Suspicious
Credential use | Access to OPC-UA server from the same IP address as observed for the corresponding app during training. | Observed
Access to OPC-UA server from a different IP address but which also corresponds to an edge server. | Expected
Access to OPC-UA server from an IP address that corresponds to a development host. | Suspicious
Network scan | Connections (source, dest., port) matching those observed during training. | Observed
Connection matching source-destination pairs observed during training, but on a different port. | Expected
Connection which does not match any source-destination pair observed during training. | Suspicious
Attempt to connect to an IP which is not assigned to any host. | Highly Suspicious
TABLE I: List of scenarios and baseline deviations generated for evaluation,
and their a-priori severity assigned. Figure 4: Probabilities—derived from
scores—for test-set triples corresponding to (a) the different scenarios
described in Table I and (b) specific graph entities, from left to right:
variables accessed by application #1, HTTPS events performed by development
host #2, HTTPS events with the historian as destination, SSH events involving
edge computing host #1 and SSH events involving the app repository. Markers
and colors denote the expected severity: Highly Suspicious (Y), Suspicious
($\times$), Unexpected ($\bullet$), Expected ($+$) and Observed ($\star$). To
increase readability, only the average probability of Observed triples is
shown.

##  V Evaluation

###  V-A Methodology

Initially, a baseline is captured with the system under normal operating
conditions, and the collected data is used to train the link prediction
algorithm of section III in an unsupervised manner. Thereafter, in order to
qualitatively evaluate its predictions, we trigger a set of actions which
result in events not observed during normal operation, but which would be
assigned a wide range of severity levels by a human expert upon detailed
analysis of the available contextual information.

Table I contains the full list of scenarios evaluated, ranging from subtle
baseline deviations which are likely benign in nature to more severe ones
which may be indicative of malicious activity such as reconnaissance, app and
host compromise, credential stealing, lateral movement or data exfiltration.

Some of these scenarios, related to app and network activity, are illustrated
in Figure 3. The examples highlight different expected scores for specific
statements based on observed activity, e.g. whether similar apps or hosts
behave in similar ways. It should be noted that any single event, e.g. a
connection between two hosts, can result in multiple statements on the graph,
with potentially different associated scores. For example, source, destination
and port in a connection may be ranked as expected while the edges
corresponding to traffic volume or to the service name may be flagged as
suspicious, giving the resulting alerts a high degree of interpretability when
examined by a human security analyst.

###  V-B Results

Results for the different test cases are shown in Figure 4a in probabilistic
form by applying the logistic function to obtained scores
$p_{s,p,o}=\sigma\left(\theta_{s,p,o}\right)$. In general, it is observed that
they follow the expected trend, with empirical score distributions showing
lower expected values for scenarios with higher severity levels. Perhaps
unsurprisingly, there are some exceptions, e.g., the statements related to
variable write access which are ranked higher than initially expected.
However, this can be due to the fact that write operations are rather
infrequent in our dataset, making it more difficult for the algorithm to learn
the general expected behaviour.

It should be noted that data modelling choices, i.e., how the raw captured
data is represented in the form of entities and relationships, can have an
impact on the resulting scores. This means that some representation choices
can be beneficial for certain scenarios and detrimental to others. Despite
this, for practical reasons, our evaluation campaign was carried out using a
single model trained on a single graph and tested on all the different
scenarios.

Results for statements involving specific graph entities are shown in Figure
4b. Compared to the previous aggregated plots, these typically show a clearer
trend, i.e., a lower degree of score overlap for different severity levels.
This indicates that it may be beneficial for analysts using our IDS system to
compare the observed scores with individual expected baselines for the
specific entities involved, rather than defining general thresholds for all
events to assess their severity.

Finally, it should be noted that all the events generated by the system
prototype at test-time are represented in the figures for the corresponding
test cases. This includes any statements resulting from spurious events unseen
during training but unrelated to the test-cases themselves. These can result
from non-deterministic behaviors built into the different software components
or unexpected network issues, and typically represent small deviations from
the baseline which are ranked with high scores. This indicates that the
proposed method is reasonably robust in realistic operation conditions, i.e.,
does not produce an unreasonable amount of false positives.

##  VI Summary and Conclusions

We propose the application of relational learning on knowledge graphs to
security monitoring and intrusion detection. The collective learning
properties of graph embedding methods allow the resulting models to generalize
beyond individual observations, benefiting from the context provided by a rich
set of entity and relationship types. This in turn translates into an
efficient use of training data, inherent robustness to false alarms in the
presence of previously unobserved events, and potentially shorter baselining
periods.

The proposed method is experimentally evaluated in an industrial automation
system prototype across a wide range of scenarios, and it is shown to leverage
context information to produce a meaningful range of severity scores. This is
particularly useful in the IDS setting, as observations are often not easily
categorized as completely benign or malicious a priori, rather evoking varying
degrees of interest or suspicion.

While these results are produced on a reduced-scale prototype, evaluation on a
real-size multi-station manufacturing environment is underway. In addition,
ongoing work aims to systematically assess the method’s robustness to changing
data modelling choices, and initial empirical results using various intuitive
alternatives (e.g., representing network connections as either nodes or as
edges connecting IP addresses) are encouraging, showing limited overall
variations in detection performance, as improvements in some scenarios tend to
compensate a slight degradation in others.

Finally, future work includes creating synergistic combinations of the
proposed graph-based method with other conventional solutions. Security
analysts stand to benefit from a well-calibrated range of scores attached to
IDS events, helping them focus their efforts on the most promising leads.

## Acknowledgment

The authors would like to thank their colleagues at the Semantics and
Reasoning group and the Siemens AI Lab for the support and the many useful
discussions in the course of this project. They would also like to thank Lukas
Vogtmann for his assistance in building the evaluation prototype.

## References

  * [1] National Science and Technology Council  “AI and Cybersecurity: Opportunities and Challenges”  In _Network & Information Technology R&D and Machine Learning & Artificial Intelligence subcommittees_, 2020  URL: https://www.nitrd.gov/pubs/AI-CS-Tech-Summary-2020.pdf
  * [2] Cong Truong, Diep Quoc Bao and Ivan Zelinka  “Artificial Intelligence in the Cyber Domain: Offense and Defense”  In _Symmetry_ 12, 2020, pp. 410  DOI: 10.3390/sym12030410
  * [3] “Framework for Improving Critical Infrastructure Cybersecurity”  National Institute of StandardsTechnology, 2018  URL: https://books.google.de/books?id=uojXvQEACAAJ
  * [4] Arwa Aldweesh, Abdelouahid Derhab and Ahmed Z. Emam  “Deep learning approaches for anomaly-based intrusion detection systems: A survey, taxonomy, and open issues”  In _Knowledge-Based Systems_ 189, 2020, pp. 105124  DOI: https://doi.org/10.1016/j.knosys.2019.105124
  * [5] K. Bresniker et al.  “Grand Challenge: Applying Artificial Intelligence and Machine Learning to Cybersecurity”  In _Computer_ 52.12, 2019, pp. 45–52  DOI: 10.1109/MC.2019.2942584
  * [6] K. Veeramachaneni et al.  “AI2: Training a Big Data Machine to Defend”  In _2016 IEEE 2nd International Conference on Big Data Security on Cloud (BigDataSecurity), IEEE International Conference on High Performance and Smart Computing (HPSC), and IEEE International Conference on Intelligent Data and Security (IDS)_ , 2016, pp. 49–54  DOI: 10.1109/BigDataSecurity-HPSC-IDS.2016.79
  * [7] M. Wang, K. Zheng, Y. Yang and X. Wang  “An Explainable Machine Learning Framework for Intrusion Detection Systems”  In _IEEE Access_ 8, 2020, pp. 73127–73141  DOI: 10.1109/ACCESS.2020.2988359
  * [8] Deval Bhamare et al.  “Cybersecurity for industrial control systems: A survey”  In _Computers & Security_ 89, 2020, pp. 101677  DOI: https://doi.org/10.1016/j.cose.2019.101677
  * [9] A. Chemudupati et al.  “The convergence of IT and Operational Technology”  In _Atos White Paper_ , 2012
  * [10] R. Paes, D.. Mazur, B.. Venne and J. Ostrzenski  “A Guide to Securing Industrial Control Networks: Integrating IT and OT Systems”  In _IEEE Industry Applications Magazine_ 26.2, 2020, pp. 47–53  DOI: 10.1109/MIAS.2019.2943630
  * [11] Eric Prud’hommeaux and Andy Seaborne  “SPARQL Query Language for RDF” http://www.w3.org/TR/rdf-sparql-query/, W3C Recommendation, 2008  URL: http://www.w3.org/TR/rdf-sparql-query/
  * [12] H. Cai, V.. Zheng and K.. Chang  “A Comprehensive Survey of Graph Embedding: Problems, Techniques, and Applications”  In _IEEE Transactions on Knowledge and Data Engineering_ 30.9, 2018, pp. 1616–1637  DOI: 10.1109/TKDE.2018.2807452
  * [13] Maximilian Nickel, Volker Tresp and Hans-Peter Kriegel  “A three-way model for collective learning on multi-relational data.”  In _Icml_ 11, 2011, pp. 809–816
  * [14] Sean Barnum  “Standardizing cyber threat intelligence information with the structured threat information expression (STIX)”  In _Mitre Corporation_ 11, 2012, pp. 1–22
  * [15] Zareen Syed et al.  “UCO: A Unified Cybersecurity Ontology”  In _AAAI Workshop: Artificial Intelligence for Cyber Security_ , 2016
  * [16] M. Eckhart, A. Ekelhart and E.. Weippl  “Automated Security Risk Identification Using AutomationML-based Engineering Data”  In _IEEE Transactions on Dependable and Secure Computing_ , 2020, pp. 1–1  DOI: 10.1109/TDSC.2020.3033150
  * [17] S.. Narayanan et al.  “Early Detection of Cybersecurity Threats Using Collaborative Cognition”  In _2018 IEEE 4th International Conference on Collaboration and Internet Computing (CIC)_ , 2018, pp. 354–363  DOI: 10.1109/CIC.2018.00054
  * [18] A. Pingle et al.  “RelExt: Relation Extraction using Deep Learning approaches for Cybersecurity Knowledge Graph Improvement”  In _2019 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM)_ , 2019, pp. 879–886  DOI: 10.1145/3341161.3343519
  * [19] Elmar Kiesling, Andreas Ekelhart, Kabul Kurniawan and Fajar Ekaputra  “The SEPSES Knowledge Graph: An Integrated Resource for Cybersecurity”  In _The Semantic Web – ISWC 2019_ Cham: Springer International Publishing, 2019, pp. 198–214
  * [20] Pejman Najafi et al.  “MalRank: A Measure of Maliciousness in SIEM-Based Knowledge Graphs”  In _Proceedings of the 35th Annual Computer Security Applications Conference_ , ACSAC ’19  San Juan, Puerto Rico, USA: Association for Computing Machinery, 2019, pp. 417–429  DOI: 10.1145/3359789.3359791
  * [21] A. Piplai et al.  “Creating Cybersecurity Knowledge Graphs From Malware After Action Reports”  In _IEEE Access_ 8, 2020, pp. 211691–211703  DOI: 10.1109/ACCESS.2020.3039234
  * [22] Sudip Mittal, Anupam Joshi and Tim Finin  “Cyber-All-Intel: An AI for Security related Threat Intelligence”  In _CoRR_ abs/1905.02895, 2019  arXiv: http://arxiv.org/abs/1905.02895
  * [23] Leman Akoglu, Hanghang Tong and Danai Koutra  “Graph based anomaly detection and description: A survey”  In _Data Mining and Knowledge Discovery_ 29.3 Springer Netherlands, 2015, pp. 626–688  DOI: 10.1007/s10618-014-0365-y
  * [24] Li Zheng et al.  “AddGraph: Anomaly Detection in Dynamic Graph Using Attention-based Temporal GCN”  In _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19_ International Joint Conferences on Artificial Intelligence Organization, 2019, pp. 4419–4425  DOI: 10.24963/ijcai.2019/614
  * [25] J. Jiang et al.  “Anomaly Detection with Graph Convolutional Networks for Insider Threat and Fraud Detection”  In _MILCOM 2019 - 2019 IEEE Military Communications Conference (MILCOM)_ , 2019, pp. 109–114  DOI: 10.1109/MILCOM47813.2019.9020760
  * [26] B. Jia et al.  “Pattern Discovery and Anomaly Detection via Knowledge Graph”  In _2018 21st International Conference on Information Fusion (FUSION)_ , 2018, pp. 2392–2399  DOI: 10.23919/ICIF.2018.8455737
  * [27] Heiko Paulheim  “Knowledge graph refinement: A survey of approaches and evaluation methods”  In _Semantic Web_ 8, 2016, pp. 489–508  DOI: 10.3233/SW-160218
  * [28] Dominik Dold and Josep Soler Garrido  “An energy-based model for neuro-symbolic reasoning on knowledge graphs”  In _unpubl._ , in prep.
  * [29] Geoffrey E Hinton, Peter Dayan, Brendan J Frey and Radford M Neal  “The” wake-sleep” algorithm for unsupervised neural networks”  In _Science_ 268.5214 American Association for the Advancement of Science, 1995, pp. 1158–1161
  * [30] Wolfgang Mahnke, Stefan-Helmut Leitner and Matthias Damm  “OPC unified architecture”  Springer Science & Business Media, 2009
  * [31] O. Kovalenko et al.  “Modeling AutomationML: Semantic Web technologies vs. Model-Driven Engineering”  In _2015 IEEE 20th Conference on Emerging Technologies Factory Automation (ETFA)_ , 2015, pp. 1–4  DOI: 10.1109/ETFA.2015.7301643
  * [32] “Zeek documentation” Accessed: 2021-02-04, https://docs.zeek.org/en/stable/intro/


# A new method of reconstructing Galactic three-dimensional structures using
ultralong-wavelength radio observations

Yanping Cong  National Astronomical Observatories, Chinese Academy of
Sciences, 20A Datun Road, Chaoyang District, Beijing 100101, China  School of
Astronomy and Space Science, University of Chinese Academy of Sciences,
Beijing 100049, China  Bin Yue  National Astronomical Observatories, Chinese
Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100101, China
Yidong Xu  National Astronomical Observatories, Chinese Academy of Sciences,
20A Datun Road, Chaoyang District, Beijing 100101, China  Yuan Shi  National
Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road,
Chaoyang District, Beijing 100101, China  School of Astronomy and Space
Science, University of Chinese Academy of Sciences, Beijing 100049, China
Xuelei Chen  National Astronomical Observatories, Chinese Academy of Sciences,
20A Datun Road, Chaoyang District, Beijing 100101, China  School of Astronomy
and Space Science, University of Chinese Academy of Sciences, Beijing 100049,
China  Department of Physics, College of Sciences, Northeastern University,
Shenyang, 110819, China  Center of High Energy Physics, Peking University,
Beijing 100871, China

###### Abstract

The free-free absorption of low frequency radio waves by thermal electrons in
the warm ionized medium of our Galaxy becomes very significant at $\lesssim
10$ MHz (ultralong-wavelength), and the absorption strength depends on the
radio frequency. Upcoming space experiments such as the Discovering Sky at the
Longest wavelength (DSL) and Farside Array for Radio Science Investigations of
the Dark ages and Exoplanets (FARSIDE) will produce high-resolution multi-
frequency sky maps at the ultralong-wavelength, providing a new window to
observe the Universe. In this paper we propose that from these ultralong-
wavelength multi-frequency maps, the three-dimensional distribution of the
Galactic electrons can be reconstructed. This novel and robust reconstruction
of the Galactic electron distribution will be a key science case of those
space missions. Ultralong-wavelength observations will be a powerful tool for
studying the astrophysics relevant to the Galactic electron distribution, for
example, the impacts of supernova explosions on electron distribution, and the
interaction between interstellar atoms and ionizing photons escaped from the
HII regions around massive stars. An animation shows the reconstructed results
using NE2001 model as input test. On ArXiv, it is given in the directory:
Ancillary files. In the paper the animation is linked to Fig. 5.

Interstellar medium (847) — Interstellar plasma (851)— Interstellar absorption
(831)— Milky Way Galaxy (1054) — Radio interferometers (1345)

##  1 Introduction

Understanding our neighborhood and environment has long been a human pursuit,
and the vicinity of the Solar system is of great interest, as it has direct
impacts on the Solar system. Besides the well-observed nearby stars (Gaia
Collaboration et al., 2022) the space is permeated with the interstellar
medium (ISM), which serves as the reservoir of the material from which stars
formed. Our Sun resides in a low-density cavity named “Local Bubble” (LB) or
“Local Hot Bubble” (LHB, Cox & Reynolds 1987), which is filled with X-ray-
emitting hot gas (Snowden et al., 1998). The LB has a typical electron number
density of $\sim 5\times 10^{-3}$ cm-3, which is only $\sim 20\%-50\%$ of the
average value at that radius on the Galactic disk (Ocker et al., 2020), and it
extends to $\sim$100 pc (Frisch et al., 2011). This low density cavity in the
ISM is thought to be created by a series of supernova explosions in
prehistoric time (Zucker et al., 2022). The spatial distribution and
morphology of the ISM reflects the complex interactions between the gas and
stars in the Galactic “ecosystem", but our knowledge of this neighborhood is
far from complete. For example, the Loop I bubble is another low density
cavity adjacent to the LB, and the well-known giant North Polar Spur (NPS)
feature in low-frequency radio sky maps is believed by many to be the
brightest part of its bubble shell (Wolleben, 2007). However, there are also
arguments that the NPS is located at the much more distant Galactic Center and
is related to the Fermi Bubble (Sarkar, 2019).

The ISM has a number of different components or phases, which are observed or
probed by various methods (Draine, 2011). The atomic and molecular hydrogen,
accounting for $\sim 80\%$ of the ISM hydrogen in mass, are traced by the 21cm
and CO lines respectively. The dust ($\sim 1$% of ISM in mass) is measured
using interstellar extinction (Lallement et al., 2014; Capitanio et al.,
2017). The hot gaseous halo ($\lesssim$ 5% of ISM mass within Galactocentric
distance of $\sim 20$ kpc) is measured from the X-ray absorption or emission
lines of highly ionized metals (Miller & Bregman, 2013). The diffuse warm
ionized medium (WIM), which occupies $\sim 20-30\%$ of the volume near the
Galactic plane and accounts for $\sim 30\%$ of the total ISM mass (Reynolds,
1991a; Kulkarni & Heiles, 1988), has been observed via the diffuse H$\alpha$
or other faint nebular emission lines (Reynolds, 1977; Reynolds & Haffner,
2000; Haffner et al., 2003; Finkbeiner, 2003; Dickinson et al., 2003), or by
the pulsar dispersion measures (DM) which could be used to reconstruct the 3D
distribution of the electrons (Cordes et al., 1991; Taylor & Cordes, 1993;
Ferrière, 2001; Gómez et al., 2001; Schnitzeler, 2012; Greiner et al., 2016;
Ocker et al., 2020), though the spatial resolution is limited by the available
pulsars. Moreover, the scattering measure (SM) can provide information on the
intrinsic fluctuations of the WIM (Williamson, 1972).

The ISM can be probed using low-frequency radio observations, based on the
fact that the ionized gas has a stronger absorption of radio waves at lower
frequencies. The radiative transfer equation of the radiation is

$\frac{dI_{\nu}}{ds}=-\alpha_{\nu}I_{\nu}+j_{\nu}$ (1)

where $I_{\nu}$ is the intensity at frequency $\nu$, and the absorption
coefficient $\alpha_{\nu}$ of the ISM at low radio frequency is dominated by
the electron free-free process111The synchrotron self-absorption is
unimportant for the typical Galactic magnetic field ($\sim 1~{}\mu$G) (Sun et
al., 2008) and cosmic ray density ($\sim 10^{-13}$ GeV-1 cm-3) (Peron et al.,
2021) above $\sim 1$ kHz (Ghisellini, 2013). The radio wave cannot propagate
below the ISM plasma frequency $\sim 2$ kHz (Jester & Falcke, 2009), so in the
frequency range considered in this paper, free-free absorption is the dominant
absorption mechanism., given by (Condon & Ransom, 2016)

$\alpha_{\nu}\approx 3.28\times 10^{-7}\left(\frac{T_{e}}{10^{4}{\rm
K}}\right)^{-1.35}\left(\frac{\nu}{\rm
GHz}\right)^{-2.1}\left(\frac{n_{e}}{\rm cm^{-3}}\right)^{2}{\rm pc}^{-1},$
(2)

where $T_{e}$ and $n_{e}$ are the free electron temperature and number density
respectively. For diffuse Galactic electrons the absorption effect becomes
significant at $\nu\lesssim 10$ MHz, while for dense HII regions the
absorption becomes significant at even higher frequencies. The absorption
strength contains information of electron densities at different distances. By
using the multi-frequency data, the density distribution of the electron can
be reconstructed, if the radiation emissivity $j_{\nu}$ (or in terms of
brightness temperature $T=\frac{c^{2}I_{\nu}}{2k\nu^{2}}$,
$\epsilon=\frac{c^{2}j_{\nu}}{2k\nu^{2}}$) is known. As this method is
sensitive to free electron density, it probes primarily the WIM component of
the ISM. Given the large volume and mass fraction of the WIM, it is obviously
highly interesting and useful to reconstruct its 3D structures.

The data at frequencies below 30 MHz, which we shall refer to as the
ultralong-wavelength band, and especially below 10 MHz, is very scarce, as
ground-based observation is hindered by ionosphere and radio frequency
interference (RFI)(Jester & Falcke, 2009). There were some ground-based
observations at Tasmania of Australia and Canada (e.g. Reber & Ellis 1956;
Ellis et al. 1962; Ellis & Hamilton 1966a; Cane & Whitham 1977; Cane 1979;
Reber 1994, for its history see also Orchiston et al. 2015a; George et al.
2015a, b; Orchiston et al. 2015b; George et al. 2015c, 2016); and some early
space observation (e.g. Hartz 1964; Alexander & Stone 1965; Smith 1965;
Alexander et al. 1969; Brown 1973; Alexander & Novaco 1974; Novaco & Brown
1978; Manning & Dulk 2001). From these observations, it is noted that there is
a downturn in the global brightness below $\sim 3$ MHz, which is believed to
be caused by the free-free absorption mechanism of the ISM (Ellis et al.,
1962; Alexander & Stone, 1965; Ellis & Hamilton, 1964; Novaco & Brown, 1978).
It has also been noted that at frequencies below a few MHz, in contrast to
higher frequencies, the Galactic poles appear to be brighter than the Galactic
plane, which is attributed to the stronger absorption on the Galactic plane
(Alexander et al., 1970; Brown, 1973; Novaco & Brown, 1978; Ellis, 1982;
Manning & Dulk, 2001).

Information on the Galactic electron distribution had been inferred from the
ultralong-wavelength observations, even with the crude data obtained in the
early observations. From the measured ultralong-wavelength global spectrum, by
assuming an absorption length, electron density were derived (Hoyle & Ellis,
1963; Smith, 1965; Ellis & Hamilton, 1966b; Alexander et al., 1970; Brown,
1973; Cane, 1979; Kassim, 1989; Fleishman & Tokarev, 1995), which is found to
be $\sim 0.03-0.1$ cm-3 near the Galactic plane. These early works assumed a
simple constant electron density. A slightly more sophisticated model assumes
that the electron density depends on the vertical distance to the Galactic
plane, which is constrained by pulsar DM (e.g. Reynolds 1991b; Nordgren et al.
1992; Taylor & Cordes 1993; Gómez et al. 2001). Peterson & Webber (2002) found
however that in this model the electron density derived from the pulsar
observations is too low to give raise the observed downturn of the ultralong-
wavelength global spectrum. To fit both observations, they suggested that the
electrons in the WIM is clumpy, which would induce stronger free-free
absorption.

Jones et al. (2016) proposed that the free-free absorption can change the
spectral index of the extragalactic background sources, and this can be used
to construct the 2D distribution of electrons (the column density along the
line of sight). However, by using discrete point sources as the background,
this approach would be limited by the number of available point sources, just
as in the pulsar observations.

The existing ultralong-wavelength observations have relatively low resolution.
However, reconstruction of the full 3D Galactic electron structures will be
feasible when high resolution ultralong-wavelength sky maps are available.
Recently, a number of ultralong-wavelength space missions have been proposed
(Chen et al., 2019), such as the Discovering the Sky at the Longest wavelength
(DSL) lunar orbit array (Chen et al., 2020; Shi et al., 2022b, a), and the
Farside Array for Radio Science Investigations of the Dark ages and Exoplanets
(FARSIDE) (Burns et al., 2019), which have the capability of producing high
resolution maps, and will enable the full 3D reconstruction. Anticipating high
resolution multi-frequency sky maps in the near future, we propose here to
reconstruct the full 3D distribution of electrons using the ultralong-
wavelength observation. The spatial distribution of the electrons is encoded
in the spectrum: at the lower frequency, there is a higher contribution of
absorption from more nearby absorption. If the emissivity along a line of
sight (LoS) is known, then the electron density profile along it can be
uniquely determined from a high resolution spectrum, without requiring any
prior information about the density profile. Combining the observations of
many LoS toward different directions one can obtain the full 3D distribution
of the electrons. To our knowledge, such multi-frequency tomographic method
has not been applied to the full reconstruction of the ISM before.

The sky radio radiation at different frequencies is produced by various
physical processes. Besides the Cosmic Microwave Background (CMB), the thermal
emission of the Galactic dust dominates the radiation above $\sim 60$ GHz. The
sky between $\sim 10-60$ GHz is dominated by the bremsstrahlung radiation from
Galactic thermal electrons. In this frequency range, there is also the
spinning dust emission that peaks at $\sim 20$ GHz. Below $\sim 10$ GHz, the
synchrotron radiation produced by the Galactic cosmic ray electrons increases
rapidly and becomes dominant (de Oliveira-Costa et al., 2008). In the
synchrotron background, roughly one-third (Seiffert et al., 2011) is from
extragalactic sources such as star-forming galaxies, AGNs, galaxy clusters,
and so on (Singal et al., 2018). Our proposed reconstruction depends on the
observation of the ultralong-wavelength radiation, and the understanding of
the unabsorbed synchrotron emissivity. The unabsorbed synchrotron radiation
has a nearly power-law spectrum, and the variation of the spectral index is
relatively small, which simplifies the reconstruction.

The outline of this paper is: In Sec. 2 we introduce our methods. In Sec. 3 we
present the statistics of our results, and show the reconstructed 1D, 2D
(Galactic plane) and 3D structures in simulations and compare them with the
input electron model. In Sec. 4 we summarize the results and discuss some
potential uncertainties.

##  2 Methods

At the frequency of interest, the emissivity is dominated by the synchrotron
radiation from cosmic ray electrons, which are believed to be produced
primarily by supernova remnants and other accelerators and propagate over the
Galaxy (Orlando & Strong, 2013a). Due to the diffusive propagation, its
intensity is relatively smooth, though it can be enhanced near the injecting
source or at regions with stronger magnetic field which can trap the
particles. Here, we model the emissivity as (Cong et al., 2021)

$\epsilon(\nu,R,Z)=A\left(\frac{R+r_{1}}{R_{0}}\right)^{\alpha}e^{-R/R_{0}}e^{-|Z/Z_{0}|^{\gamma}}\left(\frac{\nu}{\nu_{*}}\right)^{\beta_{\rm
G}}.$ (3)

where $R,Z$ are the Galactic cylindrical coordinates, $r_{1}=0.1$ kpc is a
small cut-off radius, $\nu_{*}=408$ MHz, $A$, $R_{0}$, $\alpha$, $Z_{0}$ and
$\gamma$ are model parameters, they are obtained by fitting the Haslam 408 MHz
map (Haslam et al., 1982; Remazeilles et al., 2015). That is, we integrate
this emissivity along each LoS to obtain a sky map at 408 MHz, then find the
best-fit emissivity parameters for which the sky map is closest (measured by
the Euclidean distance) to the observed map. Details can be found in Cong et
al. (2021). We use the constant $\beta_{\rm G}=-2.51$ that is fitted from
multi-frequency observations between 45 - 408 MHz.

The purpose of this paper is to demonstrate the feasibility of the algorithm,
therefore the simple emissivity model of Eq. (3) is used. One may also use
more sophisticated models, for example the emissivity constructed by cosmic
ray propagation and magnetic distribution model, and constrained by
observations like gamma-ray sky maps, e.g. GALPROP (Porter et al., 2021).
Moreover, there have been some attempts to construct the synchrotron
emissivity from observations. The opaque HII regions can fully absorb the
synchrotron radiation behind them. Therefore they can be used to separate the
synchrotron emission from regions in front of and behind the HII regions (Nord
et al., 2006; Hindson et al., 2016; Su et al., 2017, 2018; Polderman et al.,
2019).

We make the mock multi-frequency sky maps $T_{\rm sky}^{\rm obs}(\nu_{j},l,b)$
using the ULSA sky model (Cong et al., 2021). The ULSA generates the sky maps
at ultralong-wavelengths including the free-free absorption effect of the
Galactic electrons. There are a number of Galactic electron distribution
models by synthesizing different observations (Cordes et al., 1991; Gaensler
et al., 2008; Schnitzeler, 2012; Yao et al., 2017), here the NE2001 model
(Cordes & Lazio, 2002, 2003) is adopted. The basic distribution is described
by a few components with analytical formulae, for example the thin and thick
Galactic disks, and the spiral arms. Known structures such as voids and dense
HII regions are also added. Fluctuations caused by unresolved small-scale
($<1$ pc) electron structures are described by fluctuation parameters. When
simulating the mock sky maps, we adopt the same fluctuation parameter for all
the electron components in NE2001 for simplicity. Cong et al. (2021) found
that using the default fluctuation parameters in NE2001 would over-predict the
sky brightness below $\sim 3$ MHz. Instead, if increasing the fluctuation
parameter for the thick disk from 0.2 pc-2/3 to 3.0 pc-2/3, the predicted
global spectrum is consistent with the current observations. Since the thick
disk dominates the all-sky absorption, we adopt the fluctuation parameter 3.0
pc-2/3 for all components. The produced sky maps are similar to those in Fig.
8 of Cong et al. (2021). Meanwhile, as our reconstruction relies on the free-
free absorption, adopting a larger fluctuation parameter, and hence a more
efficient absorption effect, is helpful for reconstructing the electron
content. Moreover, we adopt a constant electron temperature of 8000 K. In Fig.
1 we plot the mock sky maps $T_{\rm sky}^{\rm obs}(\nu_{j},l,b)$ at $\nu=0.1$
and 0.5 MHz respectively as examples. For sky maps at $\geq 1$ MHz, we refer
interested readers to the Fig. 8 of Cong et al. (2021).

Figure 1:  The mock sky maps at 0.1 MHz (top) and 0.5 MHz (bottom) used in
this paper.

We demonstrate the reconstruction of electron density distribution from the
mock maps. The sky is pixelized with the HealPix scheme, with NSIDE$=$32, so
there are a total of 12,288 pixels for the full sky. Better performance is
expected if observations have higher resolutions. Along each pixel direction,
the LoS is divided into $N$ bins, each with an emissivity and electron
density. The emissivity which is dominated by the synchrotron is assumed to be
known. The sky brightness temperature of this LoS is modeled as

$\displaystyle T_{\rm sky}(\nu,l,b)\approx$
$\displaystyle\sum_{i=1}^{N}\bar{\epsilon}_{i}(\nu,R,Z)\Delta
x_{i}\frac{1-\exp(-\Delta\tau_{i})}{\Delta\tau_{i}}\exp(-\tau_{i-1})$
$\displaystyle+T^{\rm iso}_{E}(\nu)\exp(-\tau_{N}),$ (4)

where $\bar{\epsilon}_{i}$, $\Delta x_{i}$ and $\Delta\tau_{i}$ are the mean
emissivity, length and optical depth of the $i$-th bin, and

$\tau_{i-1}=\sum_{i^{\prime}=1}^{i-1}\Delta\tau_{i^{\prime}}.$ (5)

$T^{\rm iso}_{E}(\nu)$ is the isotropic extragalactic background used in ULSA
model, given by $T^{\rm iso}_{E}=1.2(\nu/{\rm GHz})^{-2.58}$ K (Seiffert et
al., 2011).

We make log-spaced bins along the LoS, starting from $x_{\rm min}=0.005$ kpc
with a step-interval of $\Delta\log x=0.2$, as the brightness temperature is
related to the exponential of optical depth. We have limited the LoS to a
distance at most 20 kpc from the Galactic center. The free-free optical depth
of each bin is

$\displaystyle\Delta\tau_{i}$ $\displaystyle=\alpha_{\nu,i}\Delta x_{i}$
$\displaystyle=3.28\times 10^{-7}\left(\frac{T_{e}}{10^{4}{\rm
K}}\right)^{-1.35}\left(\frac{\nu}{\rm
GHz}\right)^{-2.1}\left(\frac{\left<n^{2}_{e,i}\right>}{\rm
cm^{-6}}\right)\frac{\Delta x_{i}}{\rm pc}.$ (6)

If $T_{e}$ is known, then the reconstructed quantity is
$\left<n^{2}_{e,i}\right>$. Here the symbol “$\left<\right>$” denotes the mean
inside each bin. From $\left<n^{2}_{e,i}\right>=F_{\rm
fluc}\left<n_{e,i}\right>^{2}$, and if the fluctuation parameter222 When
converting the $\left<n_{e}^{2}\right>$ into $\left<n_{e}\right>$, we define
fluctuation parameter $F_{\rm
fluc}=\left<n_{e}^{2}\right>/\left<n_{e}\right>^{2}$. It is a dimensionless
and phenomenological parameter that describes the connection between the
electron density and the free-free absorption strength, without specifying a
physical picture of the electron distribution. In NE2001, the physically
motivated fluctuation parameter is defined as $F_{\rm
NE2001}=\zeta\omega^{2}\eta^{-1}l_{0}^{-2/3}$, where $\omega^{2}$ (to avoid
confusing the variable emissivity, in this paper we replace the
“$\epsilon^{2}$” in NE2001 paper with “$\omega^{2}$”) describes the electron
fluctuations inside electron clouds, $\zeta$ describes the cloud-to-cloud
fluctuations, $\eta$ is the filling factor of electron clouds, and $l_{0}$ is
their outer scale. Therefore $F_{\rm fluc}$ performs as $l_{0}^{2/3}F_{\rm
NE2001}\omega^{-2}(1+\omega^{2})$ in formula. We have checked that when
$l_{0}=1$ pc is adopted as in NE2001, the same $F_{\rm fluc}$ and $F_{\rm
NE2001}$ values give the same conversion factors. For this reason, we keep the
dimensionless definition of the fluctuation parameter as it is only used to
display the reconstruction results more intuitively.  $F_{\rm fluc}$ is
further known, it is translated into $\left<n_{e,i}\right>$. For demonstration
purpose we assume $T_{e}=8000$ K (Reynolds, 1990) and $F_{\rm fluc}\equiv 3.0$
is known, to show the reconstructed electron density as result. Otherwise, our
method only reconstructs a combined quantity $T_{e}^{-1.35}F_{\rm
fluc}\left<n_{e,i}\right>^{2}$. In observations $T_{e}$ varies between 5500 K
and 20000 K (Reynolds, 1990) or 6000 K - 10000 K (Haffner et al. 2009 and
references therein), and $F_{\rm fluc}$ can be $\gtrsim 3$ (Gaensler et al.,
2008). We assume that the emissivity $\epsilon(\nu,R,Z)$ is known and take the
form as given by Eq. (3).

For the multi-frequency maps, we assume the frequency interval is 0.1 MHz for
the frequency range of 0.1-1.0 MHz, and 0.2 MHz for the frequency range of 1.0
to 10.0 MHz, so there are totally 55 mock maps at 55 frequencies. We treat
$\left<n_{e,i}\right>$ of each bin as an unknown parameter. By comparing the
multi-frequency sky brightness temperature Eq. (4) that is a function of
$\left<n_{e,i}\right>$, with observations (mock maps), we can simultaneously
find $\left<n_{e,i}\right>$ of all bins that minimize

$\chi^{2}(l,b)=\sum_{j=1}^{N_{\rm freq}}\frac{[T_{\rm sky}(\nu_{j},l,b)-T_{\rm
sky}^{\rm obs}(\nu_{j},l,b)]^{2}}{\sigma^{2}_{\rm noise}},$ (7)

using the MCMC procedure (Foreman-Mackey et al., 2013), where $N_{\rm freq}$
is the number of frequency bins.

In this work we use fixed grids for solving the electron densities, which may
not accurately capture some small and dense clumps. The absorption by a single
dense clump may obscure the whole LoS, but if we already know some information
about these clumps, i.e. their location, size, and density from other
observations, we can use such information to improve the fitting. If a bin
contains a dense clump, then the density contribution from this clump to the
bin is known, leaving only the diffuse electron density as an unknown
quantity. For those dense clumps with extremely strong absorption, data at
higher frequencies can also be employed. We assume that clumps with density
$>0.1$ cm-3 have been known. For these LoS we add 20 frequency points between
10 and 30 MHz, with an interval of 1 MHz. The solution may also be improved
with adaptive grids, but we leave such technical improvements to future works.

Regarding the noise level of the observations, at low-frequency the system
temperature is dominated by the sky temperature (Shi et al., 2022b),

$\displaystyle\sigma_{T_{\rm noise}}$ $\displaystyle\sim\frac{D^{2}T_{\rm
sky}}{A_{\rm eff}\sqrt{N_{a}(N_{a}-1)t_{\rm obs}\Delta\nu}}$
$\displaystyle=\frac{4\pi T_{\rm sky}}{\theta^{2}_{\rm
res}\sqrt{N_{a}(N_{a}-1)t_{\rm obs}\Delta\nu}},$ (8)

where $A_{\rm eff}=\frac{\lambda^{2}}{4\pi}$ is the effective area of the
antenna with physical size $\ll\lambda$; $N_{a}$ is the number of antennas;
$D$ is the baseline length, and $\theta_{\rm res}\sim\frac{\lambda}{D}$ is the
angular resolution; $t_{\rm obs}$ is the integration time and $\Delta\nu$ is
the width of the frequency point. For $\theta_{\rm res}=1^{\circ}$, $N_{a}=8$
and $\Delta\nu=0.1$ MHz, we have $\sigma_{T_{\rm noise}}/T_{\rm sky}\sim 1\%$
if $t_{\rm obs}=0.1$ year; and $\sigma_{T_{\rm noise}}/T_{\rm sky}\sim 10\%$
if $t_{\rm obs}=0.5$ day. We adopt 1% noise level for all simulations in this
paper.

##  3 Results

In Fig. 2 we show the minimum $\chi^{2}$ normalized by the number of frequency
bins for each LoS as HealPix sky map (top), and the probability density (PD)
of their distribution (bottom). About 90% of the LoS have $\chi^{2}_{\rm
min}/N_{\rm freq}<1.5$, so the MCMC fitting does indeed have a good
performance. About 2% of the LoS have $\chi^{2}_{\rm min}/N_{\rm freq}>3.0$,
for displaying purpose they are not shown in Fig. 2. The distribution of
$\chi^{2}_{\rm min}/N_{\rm freq}$ is well-fitted by a Gaussian function with
the peak 1.72, the center 1.12, and the radius 0.23. We check that most LoS
with $\chi^{2}_{\rm min}/N_{\rm freq}>2.0$ have dense clumps. As described in
Sec. 2, if there is a dense clump in a LoS, we assume that its contribution to
the density in the relevant bin is known, and only fit the density of diffuse
electrons in that bin. This can improve the reconstructed density profile, see
Fig. 4 later on. But it cannot solve the problem completely. Because usually
the size of the dense clump is much smaller than the bin width, its absorption
to the emissivity in the bin also depends on its exact location inside the
bin, but this information is not used in the fix grid computation. In the
northern hemisphere between $-30^{\circ}\lesssim l\lesssim 30^{\circ}$ there
is the Loop I bubble, it is a low density cavity with huge angular size, the
absorption within the bubble is weak, see the sky maps in Fig. 1. Moreover,
such a bubble has an edge that is modeled as a sharp break in the NE2001. The
sharp break feature is hard to capture because of the bin resolution limit.
Because of the above reasons, the $\chi^{2}_{\rm min}/N_{\rm freq}$ in this
region is biased to higher values.

Figure 2:  Top: The sky map of $\chi^{2}_{\rm min}/N_{\rm freq}$ for all LoS.
Bottom: The probability density of the $\chi^{2}_{\rm min}/N_{\rm freq}$
distribution. The solid curve is a Gaussian fit, with the peak 1.72, the
center 1.12 and the radius 0.23. There are about $2\%$ of LoS have
$\chi^{2}_{\rm min}/N_{\rm freq}>3.0$. For displaying purpose, in the top
panel they are shown by color same to 3.0, while in the bottom panel these are
not counted.

We then compare our reconstructed electron density with NE2001 333 By
definition, the fluctuation parameter in NE2001 only accounts for the
fluctuations on scales smaller than 1 pc. However, our bin width is much
larger than 1 pc. So the fluctuations on scales larger than 1 pc and smaller
than bin width are absorbed in the reconstructed quantity
$\left<n_{e,i}\right>$ for the $i$-th bin. For this bin, if we divide it into
$N_{\rm sub}$ sub-bins with width of 1 pc, then in principle the reconstructed
$\left<n_{e,i}\right>$ should be more close to
$\sqrt{\sum_{j}n^{2}_{e,j}/N_{\rm sub}}$ rather than $\sum_{j}n_{e,j}/N_{\rm
sub}$, where $n_{e,j}$ is the NE2001 electron density of the $j$-th sub-bin in
the $i$-th bin. However, in most cases these two quantities are close to each
other, unless the bin width $\gtrsim 1$ kpc, or the bin contains dense clumps.
Throughout this paper, when comparing with NE2001 model we take the quantity
$\sqrt{\sum_{j}n^{2}_{e,j}/N_{\rm sub}}$ of NE2001. . The same fluctuation
parameter of 3.0 is assumed for all gas components in order to translate the
reconstructed $\left<n_{e}^{2}\right>$ into the mean electron density that is
more intuitive for displaying. We plot the mean absolute value of relative
errors between the reconstructed electron density and that of NE2001
$\overline{|f_{\rm err}|}$ for all LoS as a map in the top panel of Fig. 3. To
reduce the discrete errors, both the reconstructed density and the NE2001
density are smoothed by a Gaussian kernel with a radius equal to one bin
width. For each LoS we only count the bins with $R<12$ kpc and $|Z|<2$ kpc,
beyond that the electron density is too small and the absorption too weak for
relative error to be a good indicator of the reconstruction performance. In
the bottom panel of Fig. 3, we show the probability density distribution of
$\overline{|f_{\rm err}|}$. About 70% LoS have $\overline{|f_{\rm
err}|}<20\%$. Most of them are near the Galactic plane, see top panel of Fig.
3. This is natural because in the Galactic plane the absorption is stronger,
making it possible to reconstruct the electron density profile more
accurately. High Galactic latitude regions have larger errors, not only
because the absorption is weaker, but also because the extragalactic
background has a larger fraction in the total received flux. Around $l\sim
330^{\circ}$, there is a patch where the relative error is much higher than
the surrounding regions. We check that in NE2001 there is a low density void
with density as low as $0.001$ cm-3. In our reconstruction usually the low
density regions have larger relative errors. Nevertheless, our reconstruction
has good performance, considering the fact that in our method not only each
LoS, but also the electron density in each bin of a LoS, are all independent
of the others.

Figure 3:  Top: The sky map for the mean absolute value of the relative error
of all LoS. Bottom: The probability density of the mean absolute value of the
relative error distribution.

We show the reconstructed 1D density profiles for some LoS in Fig. 4 as
examples. In Galactic plane our algorithm correctly reproduces the electron
density profiles up to $\sim 10$ kpc, in high Galactic latitude regions up to
$\sim 3$ kpc, although there are deviations on small scales. For all these
LoS, within $\sim 0.03-0.2$ kpc, the electron density is relatively low,
indicating that our Sun is located in a bubble. For the LoS pointing to $l\sim
90^{\circ}$ however, the low density environment extends to $\sim 2$ kpc, this
is the low density region simply called “LDR” in the NE2001 model (Toscano et
al., 1999; Cordes & Lazio, 2003). A LoS penetrating the Gum Nebula is shown in
(g) and (h) panels. Panel (g) is the result directly reconstructed, and panel
(h) is the result when the location, size and density of Gum Nebula is known
and its density is added to the relevant bin. Indeed, the density profile in
panel (h) is much improved than in panel (g). Nevertheless, except for the Gum
Nebula, other dense clumps have much smaller angular sizes on the sky, the
large-scale morphology of our reconstructed maps almost does not change, if we
simply remove those LoS penetrating dense clumps.

Figure 4:  The reconstructed electron density profiles along 7 LoS in the
Galactic plane and in high Galactic regions, pointing to different directions:
(a) $l=2^{\circ}$, $b=0^{\circ}$; (b) $l=90^{\circ}$, $b=0^{\circ}$; (c)
$l=180^{\circ}$, $b=0^{\circ}$; (d) $l=101.67^{\circ}$, $b=49.70^{\circ}$; (e)
$l=5.62^{\circ}$, $b=34.23^{\circ}$; (f) $l=260.16^{\circ}$,
$b=-41.81^{\circ}$; (g) $l=270^{\circ}$, $b=0^{\circ}$, the bump at $\sim 0.5$
kpc is the Gum Nebula; (h) same as (g) but assume information on Gum Nebula on
this LoS is already known, and $0.1\leq\nu\leq 30$ MHz data is used for
fitting. As comparison in each panel we also plot the NE2001 density of same
$r$ bins. In all panels the reconstructed results and the NE2001 model are
smoothed by Gaussian kernel with radius equal to one bin.  Figure 5:  Top: The
reconstructed electron distribution. Bottom : The NE2001 model. Both the top
and bottom panel are smoothed by a Gaussian kernel with a radius of 0.05 kpc.
An associated animation for the reconstructed electron density is available
online. In the animation, the 3D box rotates about the $Z$-axis in the first
30 seconds duration, so one can see the structures from different
perspectives. In the last 3 seconds, we zoom in the 3D box until the innermost
LB is seen. (On ArXiv, please find the animation in the directory: Ancillary
files)

With the mock full-sky map, we can reconstruct the 3D distribution of Galactic
electrons by synthesizing all LoS. The tomographic reconstruction result of
our Sun’s vicinity is shown in the top panel of Fig. 5, and for comparison the
original NE2001 model distribution is shown in the bottom panel. This figure
uses 12288 LoS, and the density volume is a heliocentric cuboid of $6\times
6\times 4$ kpc3, with pixel size $0.05$ kpc. However, the poorest resolution
at the edge is $\sim 0.5$ kpc, limited by the bin width. The reconstructed
volume map is also shown in a video linked to Fig. 5. Both the reconstructed
density and the NE2001 density are smoothed by a Gaussian kernel with radius
0.05 kpc. Our method reconstructs major features from NE2001, including the
LB, the LSB, the Gum Nebula, and so on.

In our reconstructed electron density field, there are some streaks as shown
in the top panel of Fig. 5. This is a consequence of our solving the density
along the radial direction from our position as the center. Since we adopt a
grid along the radial direction, the grid cells will have shapes that either
appear elongated if their size along the radial direction is larger than the
tangential direction, or oblate if the reverse is true. In our case we adopted
a log-spaced grid which is good for the numerical problem, but at the far end
the streak is more apparent. Increasing the number of LoS and reducing the bin
width can relieve the steak problem. However, if the bin width is too small
then the noise of each bin is large since the contribution of a single bin to
the total absorption is small.

In Fig. 6 we show the probability density as a function of $|f_{\rm err}|$ and
$n_{e}$, where $|f_{\rm err}|$ is the absolute value of the relative error
between the top and bottom panels of Fig. 5 for each pixel, and $n_{e}$ is the
density of relevant pixel in the bottom panel of Fig. 5. About 70% pixels have
$|f_{\rm err}|<30\%$, and the median of $|f_{\rm err}|$ is 17%. About 8%
pixels have $|f_{\rm err}|>100\%$, most of them are low-density regions, see
the top left part of Fig. 6. This is not surprising because the absorption in
low-density regions is weaker, therefore hard to constrain the density
accurately. Considering the above, any structures with size comparable to the
resolution limit and/or with density contrast comparable to the uncertainty
level are hard to see from the figure.

Figure 6:  The $|f_{\rm err}|$ vs. $n_{e}$ in the reconstructed 3D Cartesian
coordinate density field. Colorbar indicates the probability density with
normalization $\int{\rm PD}d|f_{\rm err}|dn_{e}=1$.

Moving further from our vicinity on to larger scales, we meet two spiral arms
that bracket our Sun, as shown in the reconstructed electron density on the
Galactic plane in Fig. 7. The reconstructed electron map clearly shows the
Carina-Sagittarius arm in the direction towards the Galactic center, and the
Perseus arm in the anti-center direction. Previously the spiral arms have been
mapped via neutral hydrogen, CO (trace the molecular hydrogen, Englmaier et
al. 2011; Nakanishi & Sofue 2016), dust extinction survey (Rezaei Kh. et al.,
2018), and other tracers for high-mass star formation (Hou & Han, 2014). The
Ultralong-wavelength observations will allow reconstruction of the diffuse
free electrons in the spiral arms.

Figure 7:  Top: The face on view of the electron distribution on Galactic
plane reconstructed from 900 equally separated LoS. To suppress the noise and
highlight the large-scale structures, the density fields are smoothed by
Gaussian kernal with a radius 0.25 kpc. Bottom: The electron density in the
NE2001 model. For comparison purpose, the NE2001 data is binned and smoothed
same to reconstructed density. To guide the eye we plot the locations of
spiral arms modeled in NE2001.

##  4 Summary and Discussions

The ultralong-wavelength radio spectrum below $\sim$10 MHz has largely been
constrained from the ground and space at poor spatial resolution. Some space
missions have been proposed to obtain high-resolution sky maps in this band.
In this paper we proposed that, thanks to the frequency-dependent absorption
to electromagnetic wave, from the upcoming ultralong-wavelength observations
the three-dimensional structures of the Galactic electrons can be
reconstructed. Using mock sky maps we performed simulations to prove the
feasibility of this method.

There are dense clumps in the sky, and they often have active star formation
and also emit synchrotron radiation. If such emission is not accurately known,
then it may cause errors in reconstruction. To assess this uncertainty, we
consider the nearby Gum Nebula. According to Woermann (1998), at 408 MHz the
Gum Nebula has typical surface brightness temperature of roughly 10 K, and the
typical spectrum index derived from 408 and 2326 MHz is roughly $-2.5$. Taking
the Gum Nebula as a sphere with radius 0.14 kpc and electron density 0.43
cm-3, at a distance 0.5 kpc from us (Cordes & Lazio, 2002), and the emissivity
is uniformly distributed in the sphere, we expect $\sim 5\times 10^{4}$ K
surface brightness at 1 MHz if $F_{\rm fluc}=3.0$. This is actually smaller
than the mean brightness of the Galactic plane in our sky model at the same
frequency. According to this simple estimate, a large fraction of the clump
emission might be absorbed by itself, since it has density much higher than
the ambient medium. So at least in this case, the influence of the clump
emission is modest.

Our synchrotron emissivity model is cylindrical. However the real emissivity
could have non-cylindrical large-scale structures, induced by the spiral
structures in the regular magnetic field (Orlando & Strong, 2013b). We have
tested that for such a spiral emissivity model, the electron density can be
still reconstructed and all conclusions do not change. The reconstruction
errors could come from the unknown random small-scale fluctuations in the
synchrotron emissivity. To estimate its impact, we add fluctuations to the
emissivity model of Eq. (3), and then reconstruct the electron density profile
assuming the small-scale fluctuations are not present. Motivated by the
measured angular power spectrum of Galactic synchrotron radiation
$C_{l}\propto l^{-2}$ (Iacobelli et al., 2013), we model the small-scale
fluctuations of emissivity to be Gaussian with the power spectrum

$P_{\epsilon}(k)\propto k^{-3},~{}~{}k>k_{\rm min},$ (9)

where $k$ is the spatial wavenumber. We assume that fluctuations only exist on
scales smaller than $2\pi/k_{\rm min}$ ($k>k_{\rm min}$). On larger scales, we
assume the Eq. (3) correctly describes the cosmic ray emissivity. We normalize
$P_{\epsilon}$ so that the emissivity has a mean relative fluctuations level
$\sim 50\%$ at the smoothing scale 0.1 kpc. That is, when the fluctuations
field following the power spectrum of Eq. (9) is added to the cylindrical
emissivity field, and the new emissivity field is smoothed with Gaussian
kernel with radius $=0.1$ kpc, then
$\sqrt{\overline{(\delta\epsilon)^{2}}}/\overline{\epsilon}\sim 50\%$. We fond
that, despite this perturbation, which induces errors in the reconstruction at
small scales, the reconstruction result remain largely unchanged at large
scales. The results of this reconstruction are shown in Fig. 8, for $k_{\rm
min}=5.0~{}\mathrm{kpc}^{-1}$ and $k_{\rm min}=0.5~{}\mathrm{kpc}^{-1}$
respectively.

Figure 8:  The emissivity with small-scale fluctuations and its known smooth
component (top panels), and the reconstructed electron density profiles for
the LoS with $l=2.0^{\circ}$ on Galactic plane (bottom panels). For left
panels we adopt $k_{\rm min}=5.0$ kpc-1 and for right panels we adopt $k_{\rm
min}=0.5$ kpc-1. All density profiles, including the NE2001, have been
smoothed by a Gaussian kernel with radius equal to one bin.

Our method directly reconstructs $\left<n_{e}^{2}\right>$ of each bin, a
quantity that produces absorption equivalent to real electron distribution
inside this bin. Since $\left<n_{e}^{2}\right>=F_{\rm
fluc}\left<n_{e}\right>^{2}$, we actually cannot distinguish the cases if a
region has smaller electron density but larger fluctuation parameter, or the
reverse, unless the $F_{\rm fluc}$ is derived from other methods. Just for
convenience of translating the results into a more intuitive quantity
$\left<n_{e}\right>$ we assume a constant $F_{\rm fluc}$. The $F_{\rm fluc}$
may vary bin-by-bin since the fluctuations of the WIM can depend on locations
(Reynolds, 1991a; Gaensler et al., 2008), however we have checked that it will
not change our reconstruction of $\left<n^{2}_{e}\right>$. Just, if the
correct $F_{\rm fluc}$ is not known, the translated $\left<n_{e}\right>$ would
deviate from the real mean electron density. For example, in Gaensler et al.
(2008) they investigated the volume filling factor of electron clouds. If the
electron fluctuations inside clouds are negligible, then the fluctuation
parameter can be well approximated as the reciprocal of the filling factor. In
this case the fluctuation parameter is $\sim 25$ in the Galactic plane and
$\sim 3$ at $|Z|=1.4$ kpc. For such a fluctuation parameter then our
$\left<n_{e}\right>$ at Galactic plane is overestimated by a factor $\sim 3$.
In Ocker et al. (2021) they considered a complicated electron cloud model
rather than uniform density. According to the observations of scattering and
dispersion measures of the Galactic pulsars, they found that the fluctuation
parameters at the low Galactic latitudes and inner Galaxy are $\gtrsim
10-1000$ times larger than the typical thick disk at high latitudes.
Considering this, the error on the reconstructed electron density by assuming
a constant fluctuation parameter can be even much larger, particularly at low
Galactic latitudes and for the LoS towards the inner Galaxy. Future studies
trying to disentangle the electron density would merit a more sophisticated
model for the spatial variation in the fluctuation parameter. Perhaps values
independently inferred from pulsar observations along similar LoS could be
employed, although this would be limited by the sparseness of pulsar spatial
distribution. In some biased regions, the WIM has larger fluctuations. The
synchrotron emissivity from such regions will be also biased, because usually
such regions are more close to the sources of the cosmic ray– supernova
remnants. We tested that, if the biased emissivity is correctly involved in
the emissivity model, then no matter if the fluctuation parameter is biased or
not in the same regions, the $\left<n^{2}_{e}\right>$ can still be correctly
reconstructed. However if the biased emissivity is missed in the model, then
like what we have discussed in last paragraph, it would induce small-scale
uncertainties on the reconstructed $\left<n^{2}_{e}\right>$ and the final
$\left<n_{e}\right>$.

Despite some modeling uncertainties, the reconstruction from ultralong-
wavelength observation is fairly robust. Once the ultralong-wavelength sky is
systematically surveyed, we will gain new knowledge about our residence in the
Milky Way, and pin down the long-lasting debate on the distance of the NPS. It
will be an important step in obtaining a synthetic picture of our Galaxy.

## Acknowledgments

We thank the anonymous referee for the very helpful comments. This work is
supported by National SKA Program of China Nos. 2020SKA0110402 and
2020SKA0110401, Inter-government cooperation Flagship program grant No.
2018YFE0120800, the National Natural Science Foundation of China grant Nos.
11973047, 11633004, the Chinese Academy of Sciences (CAS) Key Instrument grant
ZDKYYQ20200008, the Strategic Priority Research Program XDA15020200, and the
CAS Frontier Science Key Project QYZDJ-SSW-SLH017.

## References

  * Alexander et al. (1970) Alexander, J. K., Brown, L. W., Clark, T. A., & Stone, R. G. 1970, A&A, 6, 476
  * Alexander et al. (1969) Alexander, J. K., Brown, L. W., Clark, T. A., Stone, R. G., & Weber, R. R. 1969, ApJ, 157, L163, doi: 10.1086/180411
  * Alexander & Novaco (1974) Alexander, J. K., & Novaco, J. C. 1974, AJ, 79, 777, doi: 10.1086/111608
  * Alexander & Stone (1965) Alexander, J. K., & Stone, R. G. 1965, ApJ, 142, 1327, doi: 10.1086/148418
  * Brown (1973) Brown, L. W. 1973, ApJ, 180, 359, doi: 10.1086/151968
  * Burns et al. (2019) Burns, J. O., Hallinan, G., Lux, J., et al. 2019, arXiv e-prints, arXiv:1911.08649.  https://arxiv.org/abs/1911.08649
  * Cane (1979) Cane, H. V. 1979, MNRAS, 189, 465, doi: 10.1093/mnras/189.3.465
  * Cane & Whitham (1977) Cane, H. V., & Whitham, P. S. 1977, MNRAS, 179, 21, doi: 10.1093/mnras/179.1.21
  * Capitanio et al. (2017) Capitanio, L., Lallement, R., Vergely, J. L., Elyajouri, M., & Monreal-Ibero, A. 2017, A&A, 606, A65, doi: 10.1051/0004-6361/201730831
  * Chen et al. (2020) Chen, X., Yan, J., Deng, L., et al. 2020, Phil. Trans. Roy. Soc. Lond. A, 379, 20190566, doi: 10.1098/rsta.2019.0566
  * Chen et al. (2019) Chen, X., Burns, J., Koopmans, L., et al. 2019, arXiv e-prints, arXiv:1907.10853.  https://arxiv.org/abs/1907.10853
  * Condon & Ransom (2016) Condon, J. J., & Ransom, S. M. 2016, Essential Radio Astronomy
  * Cong et al. (2021) Cong, Y., Yue, B., Xu, Y., et al. 2021, ApJ, 914, 128, doi: 10.3847/1538-4357/abf55c
  * Cordes & Lazio (2002) Cordes, J. M., & Lazio, T. J. W. 2002, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0207156
  * Cordes & Lazio (2003) —. 2003, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0301598
  * Cordes et al. (1991) Cordes, J. M., Weisberg, J. M., Frail, D. A., Spangler, S. R., & Ryan, M. 1991, Nature, 354, 121, doi: 10.1038/354121a0
  * Cox & Reynolds (1987) Cox, D. P., & Reynolds, R. J. 1987, ARA&A, 25, 303, doi: 10.1146/annurev.aa.25.090187.001511
  * de Oliveira-Costa et al. (2008) de Oliveira-Costa, A., Tegmark, M., Gaensler, B. M., et al. 2008, MNRAS, 388, 247, doi: 10.1111/j.1365-2966.2008.13376.x
  * Dickinson et al. (2003) Dickinson, C., Davies, R. D., & Davis, R. J. 2003, MNRAS, 341, 369, doi: 10.1046/j.1365-8711.2003.06439.x
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium
  * Ellis (1982) Ellis, G. R. A. 1982, Australian Journal of Physics, 35, 91, doi: 10.1071/PH820091
  * Ellis & Hamilton (1964) Ellis, G. R. A., & Hamilton, P. A. 1964, Nature, 204, 272, doi: 10.1038/204272a0
  * Ellis & Hamilton (1966a) —. 1966a, ApJ, 143, 227, doi: 10.1086/148493
  * Ellis & Hamilton (1966b) —. 1966b, ApJ, 146, 78, doi: 10.1086/148860
  * Ellis et al. (1962) Ellis, G. R. A., Waterworth, M. D., & Bessell, M. 1962, Nature, 196, 1079, doi: 10.1038/1961079a0
  * Englmaier et al. (2011) Englmaier, P., Pohl, M., & Bissantz, N. 2011, Memorie della Societa Astronomica Italiana Supplementi, 18, 199.  https://arxiv.org/abs/0812.3491
  * Ferrière (2001) Ferrière, K. M. 2001, Reviews of Modern Physics, 73, 1031, doi: 10.1103/RevModPhys.73.1031
  * Finkbeiner (2003) Finkbeiner, D. P. 2003, ApJS, 146, 407, doi: 10.1086/374411
  * Fleishman & Tokarev (1995) Fleishman, G. D., & Tokarev, Y. V. 1995, A&A, 293, 565
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Frisch et al. (2011) Frisch, P. C., Redfield, S., & Slavin, J. D. 2011, ARA&A, 49, 237, doi: 10.1146/annurev-astro-081710-102613
  * Gaensler et al. (2008) Gaensler, B. M., Madsen, G. J., Chatterjee, S., & Mao, S. A. 2008, PASA, 25, 184, doi: 10.1071/AS08004
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Drimmel, R., Romero-Gomez, M., et al. 2022, arXiv e-prints, arXiv:2206.06207.  https://arxiv.org/abs/2206.06207
  * George et al. (2015a) George, M., Orchiston, W., Slee, B., & Wielebinski, R. 2015a, Journal of Astronomical History and Heritage, 18, 14
  * George et al. (2015b) —. 2015b, Journal of Astronomical History and Heritage, 18, 177
  * George et al. (2016) —. 2016, Journal of Astronomical History and Heritage, 19, 185
  * George et al. (2015c) George, M., Orchiston, W., Wielebinski, R., & Slee, B. 2015c, Journal of Astronomical History and Heritage, 18, 312
  * Ghisellini (2013) Ghisellini, G. 2013, Radiative Processes in High Energy Astrophysics, Vol. 873, doi: 10.1007/978-3-319-00612-3
  * Gómez et al. (2001) Gómez, G. C., Benjamin, R. A., & Cox, D. P. 2001, AJ, 122, 908, doi: 10.1086/321180
  * Greiner et al. (2016) Greiner, M., Schnitzeler, D. H. F. M., & Enßlin, T. A. 2016, A&A, 590, A59, doi: 10.1051/0004-6361/201526717
  * Haffner et al. (2003) Haffner, L. M., Reynolds, R. J., Tufte, S. L., et al. 2003, ApJS, 149, 405, doi: 10.1086/378850
  * Haffner et al. (2009) Haffner, L. M., Dettmar, R. J., Beckman, J. E., et al. 2009, Reviews of Modern Physics, 81, 969, doi: 10.1103/RevModPhys.81.969
  * Hartz (1964) Hartz, T. R. 1964, Nature, 203, 173, doi: 10.1038/203173a0
  * Haslam et al. (1982) Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1
  * Hindson et al. (2016) Hindson, L., Johnston-Hollitt, M., Hurley-Walker, N., et al. 2016, PASA, 33, e020, doi: 10.1017/pasa.2016.19
  * Hou & Han (2014) Hou, L. G., & Han, J. L. 2014, A&A, 569, A125, doi: 10.1051/0004-6361/201424039
  * Hoyle & Ellis (1963) Hoyle, F., & Ellis, G. R. A. 1963, Australian Journal of Physics, 16, 1, doi: 10.1071/PH630001
  * Iacobelli et al. (2013) Iacobelli, M., Haverkorn, M., Orrú, E., et al. 2013, A&A, 558, A72, doi: 10.1051/0004-6361/201322013
  * Jester & Falcke (2009) Jester, S., & Falcke, H. 2009, New A Rev., 53, 1, doi: 10.1016/j.newar.2009.02.001
  * Jones et al. (2016) Jones, D. I., Igoshev, A. P., & Haverkorn, M. 2016, MNRAS, 460, 3298, doi: 10.1093/mnras/stw1192
  * Kassim (1989) Kassim, N. E. 1989, ApJ, 347, 915, doi: 10.1086/168183
  * Kulkarni & Heiles (1988) Kulkarni, S. R., & Heiles, C. 1988, in Galactic and Extragalactic Radio Astronomy, ed. K. I. Kellermann & G. L. Verschuur, 95–153
  * Lallement et al. (2014) Lallement, R., Vergely, J. L., Valette, B., et al. 2014, A&A, 561, A91, doi: 10.1051/0004-6361/201322032
  * Manning & Dulk (2001) Manning, R., & Dulk, G. A. 2001, A&A, 372, 663, doi: 10.1051/0004-6361:20010516
  * Miller & Bregman (2013) Miller, M. J., & Bregman, J. N. 2013, ApJ, 770, 118, doi: 10.1088/0004-637X/770/2/118
  * Nakanishi & Sofue (2016) Nakanishi, H., & Sofue, Y. 2016, PASJ, 68, 5, doi: 10.1093/pasj/psv108
  * Nord et al. (2006) Nord, M. E., Henning, P. A., Rand, R. J., Lazio, T. J. W., & Kassim, N. E. 2006, AJ, 132, 242, doi: 10.1086/504407
  * Nordgren et al. (1992) Nordgren, T. E., Cordes, J. M., & Terzian, Y. 1992, AJ, 104, 1465, doi: 10.1086/116331
  * Novaco & Brown (1978) Novaco, J. C., & Brown, L. W. 1978, ApJ, 221, 114, doi: 10.1086/156009
  * Ocker et al. (2020) Ocker, S. K., Cordes, J. M., & Chatterjee, S. 2020, ApJ, 897, 124, doi: 10.3847/1538-4357/ab98f9
  * Ocker et al. (2021) —. 2021, ApJ, 911, 102, doi: 10.3847/1538-4357/abeb6e
  * Orchiston et al. (2015a) Orchiston, W., George, M., Slee, B., & Wielebinski, R. 2015a, Journal of Astronomical History and Heritage, 18, 3
  * Orchiston et al. (2015b) Orchiston, W., Slee, B., George, M., & Wielebinski, R. 2015b, Journal of Astronomical History and Heritage, 18, 285
  * Orlando & Strong (2013a) Orlando, E., & Strong, A. 2013a, MNRAS, 436, 2127, doi: 10.1093/mnras/stt1718
  * Orlando & Strong (2013b) Orlando, E., & Strong, A. W. 2013b, in International Cosmic Ray Conference, Vol. 33, International Cosmic Ray Conference, 841.  https://arxiv.org/abs/1307.2264
  * Peron et al. (2021) Peron, G., Aharonian, F., Casanova, S., Yang, R., & Zanin, R. 2021, ApJ, 907, L11, doi: 10.3847/2041-8213/abcaa9
  * Peterson & Webber (2002) Peterson, J. D., & Webber, W. R. 2002, ApJ, 575, 217, doi: 10.1086/341258
  * Polderman et al. (2019) Polderman, I. M., Haverkorn, M., Jaffe, T. R., & Alves, M. I. R. 2019, A&A, 621, A127, doi: 10.1051/0004-6361/201834405
  * Porter et al. (2021) Porter, T. A., Johannesson, G., & Moskalenko, I. V. 2021, arXiv e-prints, arXiv:2112.12745.  https://arxiv.org/abs/2112.12745
  * Reber (1994) Reber, G. 1994, JRASC, 88, 297
  * Reber & Ellis (1956) Reber, G., & Ellis, G. R. 1956, J. Geophys. Res., 61, 1, doi: 10.1029/JZ061i001p00001
  * Remazeilles et al. (2015) Remazeilles, M., Dickinson, C., Banday, A. J., Bigot-Sazy, M. A., & Ghosh, T. 2015, MNRAS, 451, 4311, doi: 10.1093/mnras/stv1274
  * Reynolds (1977) Reynolds, R. J. 1977, ApJ, 216, 433, doi: 10.1086/155484
  * Reynolds (1990) —. 1990, in Low Frequency Astrophysics from Space, ed. N. E. Kassim & K. W. Weiler, Vol. 362, 121, doi: 10.1007/3-540-52891-1
  * Reynolds (1991a) —. 1991a, ApJ, 372, L17, doi: 10.1086/186013
  * Reynolds (1991b) Reynolds, R. J. 1991b, in The Interstellar Disk-Halo Connection in Galaxies, ed. H. Bloemen, Vol. 144, 67
  * Reynolds & Haffner (2000) Reynolds, R. J., & Haffner, L. M. 2000, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0010618
  * Rezaei Kh. et al. (2018) Rezaei Kh., S., Bailer-Jones, C. A. L., Hogg, D. W., & Schultheis, M. 2018, A&A, 618, A168, doi: 10.1051/0004-6361/201833284
  * Sarkar (2019) Sarkar, K. C. 2019, MNRAS, 482, 4813, doi: 10.1093/mnras/sty2944
  * Schnitzeler (2012) Schnitzeler, D. H. F. M. 2012, MNRAS, 427, 664, doi: 10.1111/j.1365-2966.2012.21869.x
  * Seiffert et al. (2011) Seiffert, M., Fixsen, D. J., Kogut, A., et al. 2011, ApJ, 734, 6, doi: 10.1088/0004-637X/734/1/6
  * Shi et al. (2022a) Shi, Y., Deng, F., Xu, Y., et al. 2022a, ApJ, 929, 32, doi: 10.3847/1538-4357/ac5965
  * Shi et al. (2022b) Shi, Y., Xu, Y., Deng, L., et al. 2022b, MNRAS, 510, 3046, doi: 10.1093/mnras/stab3623
  * Singal et al. (2018) Singal, J., Haider, J., Ajello, M., et al. 2018, PASP, 130, 036001, doi: 10.1088/1538-3873/aaa6b0
  * Smith (1965) Smith, F. G. 1965, MNRAS, 131, 145, doi: 10.1093/mnras/131.1.145
  * Snowden et al. (1998) Snowden, S. L., Egger, R., Finkbeiner, D. P., Freyberg, M. J., & Plucinsky, P. P. 1998, ApJ, 493, 715, doi: 10.1086/305135
  * Su et al. (2017) Su, H., Hurley-Walker, N., Jackson, C. A., et al. 2017, MNRAS, 465, 3163, doi: 10.1093/mnras/stw2959
  * Su et al. (2018) Su, H., Macquart, J. P., Hurley-Walker, N., et al. 2018, MNRAS, 479, 4041, doi: 10.1093/mnras/sty1732
  * Sun et al. (2008) Sun, X. H., Reich, W., Waelkens, A., & Enßlin, T. A. 2008, A&A, 477, 573, doi: 10.1051/0004-6361:20078671
  * Taylor & Cordes (1993) Taylor, J. H., & Cordes, J. M. 1993, ApJ, 411, 674, doi: 10.1086/172870
  * Toscano et al. (1999) Toscano, M., Britton, M. C., Manchester, R. N., et al. 1999, ApJ, 523, L171, doi: 10.1086/312276
  * Williamson (1972) Williamson, I. P. 1972, MNRAS, 157, 55, doi: 10.1093/mnras/157.1.55
  * Woermann (1998) Woermann, B. 1998, PhD thesis, Rhodes University, South Africa
  * Wolleben (2007) Wolleben, M. 2007, ApJ, 664, 349, doi: 10.1086/518711
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29, doi: 10.3847/1538-4357/835/1/29
  * Zucker et al. (2022) Zucker, C., Goodman, A. A., Alves, J., et al. 2022, Nature, 601, 334, doi: 10.1038/s41586-021-04286-5