# STAL: Spike Threshold Adaptive Learning Encoder for Classification of Pain-
Related Biosignal Data

Freek Hens  0009-0005-0405-0560 Mohammad Mahdi Dehshibi Senior Member, IEEE
0000-0001-8112-5419 Leila Bagheriye Member, IEEE  0000-0003-1605-5850 Mahyar
Shahsavari  0000-0001-7703-6835 Member, IEEE  Ana Tajadura-Jiménez The paper
is submitted for review on . “This work was supported by the European Union’s
Horizon 2020 program under grant No 101002711 for the BODYinTRANSIT project.”
Freek Hens (e-mail: freek.hens@outlook.com), Mahyar Shahsavari, and Leila
Bagheriye are with the Donders Institute for Brain, Cognition and Behavior at
Radboud University in Nijmegen, The Netherlands. Mohammad Mahdi Dehshibi and
Ana Tajadura-Jiménez are with the Department of Computer Science and
Engineering at Universidad Carlos III de Madrid, Léganes, Spain.
0000-0003-3166-3512

###### Abstract

This paper presents the first application of spiking neural networks (SNNs)
for the classification of chronic lower back pain (CLBP) using the EmoPain
dataset. Our work has two main contributions. We introduce Spike Threshold
Adaptive Learning (STAL), a trainable encoder that effectively converts
continuous biosignals into spike trains. Additionally, we propose an ensemble
of Spiking Recurrent Neural Network (SRNN) classifiers for the multi-stream
processing of sEMG and IMU data. To tackle the challenges of small sample size
and class imbalance, we implement minority over-sampling with weighted sample
replacement during batch creation. Our method achieves outstanding performance
with an accuracy of 80.43%, AUC of 67.90%, F1 score of 52.60%, and Matthews
Correlation Coefficient (MCC) of 0.437, surpassing traditional rate-based and
latency-based encoding methods. The STAL encoder shows superior performance in
preserving temporal dynamics and adapting to signal characteristics.
Importantly, our approach (STAL-SRNN) outperforms the best deep learning
method in terms of MCC, indicating better balanced class prediction. This
research contributes to the development of neuromorphic computing for
biosignal analysis. It holds promise for energy-efficient, wearable solutions
in chronic pain management.

Spike Train Encoder, Spiking Recurrent Neural Network, EmoPain, Chronic Pain

##  1 Introduction

Biosignals are pivotal in non-invasively analysing and understanding
physiological processes and human body function [1]. Among these, surface
electromyography (sEMG) and data from inertial measurement units (IMUs) offer
a comprehensive representation of an individual’s motor behaviour [2]. sEMG
measures the electrical activity of skeletal muscles, while IMUs directly
quantify body segment movements. Analysing human motor behaviour has found
widespread applications in several domains, particularly in chronic pain
management [3].

Chronic lower back pain (CLBP) is a challenging condition that impacts
millions of individuals worldwide, often resulting in reduced quality of life
and significant economic strain. The precise detection and monitoring of CLBP
are essential for developing effective pain management strategies [4].
However, the analysis of biosignals for this purpose presents several
challenges. One primary challenge is signal contamination from various
sources, including motion artefacts and crosstalk from nearby physiological
processes [5]. Another lies in inter- and intra-individual variability due to
physiological differences or factors like fatigue and stress [6]. Beyond these
challenges, sEMG data are high-dimensional and multi-channel, with a non-
stationary nature due to variations in muscle anatomy between individuals.
Similarly, IMU signals present challenges such as drift and integration
errors, sensitivity to sensor placement, and the complexity of interpreting 3D
motion data [7].

Conventionally, researchers have addressed these challenges using handcrafted
features and filtering techniques. While promising, these approaches often
require domain expertise and may fail to capture the intricate dynamics of
complex biosignals [8]. More recently, deep learning models such as
convolutional and recurrent neural networks [9, 7] have demonstrated the
potential to derive intricate features from raw data for CLBP classification.
However, these models typically demand large, high-quality datasets for
training and are computationally intensive.

Spiking Neural Networks (SNNs) offer a promising alternative for biosignal
analysis. Inspired by biological neural systems, SNNs process information
through discrete spike events, mimicking the communication mechanism of
biological neurons. This biological plausibility is particularly relevant for
sEMG and IMU analysis, as SNNs can capture the complex temporal dynamics of
these signals more effectively [10]. Moreover, when deployed on neuromorphic
hardware, SNNs can operate with significantly lower power consumption compared
to other neural network architectures. This feature makes them ideal for
resource-constrained environments like wearable devices [11].

Despite these advantages, a significant challenge in using SNNs for biosignal
analysis lies in converting continuous analogue signals into binary spike
trains – the native language of SNNs. Various encoding methods have been
proposed, including rate-based, time-to-first-spike, and delta-modulated
approaches. However, these methods often struggle to preserve the temporal
dynamics of biosignals, resist inherent noise, and adapt to signal
characteristics without violating biological constraints [12, 13].

In this paper, we present the first application of neuromorphic computing
techniques to CLBP classification using sEMG and IMU data provided in the
EmoPain dataset [3]. This dataset includes recordings from 25 healthy
individuals and 21 subjects with CLBP performing a sequence of physical
exercises. We introduce two key innovations:

  1. 1.

Spike Threshold Adaptive Learning (STAL): A learnable encoder designed to
convert continuous sEMG and IMU signals into spike trains that efficiently
preserve temporal information.

  2. 2.

An ensemble of Spiking Recurrent Neural Network (SRNN) classifiers: This
multi-stream SRNN uses recurrent leaky integrate-and-fire neurons to process
the spike trains and improves accuracy and robustness in classifying
individuals as healthy or having CLBP with an ensemble of Random Forest
classifiers.

By combining the biological plausibility and efficiency of SNNs with the
proven performance of ensemble learning methods, our two-fold contribution
provides opportunities for developing more effective, personalised
rehabilitation strategies for individuals with CLBP. Furthermore, the proposed
approach has potential broader impacts beyond CLBP classification. It could be
adapted for other biosignal analysis tasks in healthcare, such as detecting
neurodegenerative disorders or monitoring cardiovascular health. The low-power
consumption of SNNs also makes this approach promising for long-term,
continuous health monitoring using wearable devices.

The remainder of this paper is organised as follows: Section 2 reviews related
work in biosignal encoding and studies that addressed challenges defined in
the EmoPain database. Section 3 details our proposed method, including the
STAL encoder and the ensemble of SRNN classifiers. Section 4 presents our
experimental setup and results, comparing our method’s performance to existing
approaches. Finally, Section 5 concludes the paper and discusses future
research directions.

##  2 Literature Review

The use of spiking neural networks in biosignal processing has shown great
promise, thanks to their biological plausibility, lower power consumption, and
reduced training requirements compared to conventional deep learning methods
[14]. However, their potential application in chronic pain detection and
identification of pain-related protective behaviour is an unexplored area,
presenting a unique opportunity in this domain.

Figure 1: A schematic of the proposed architecture, including the Spike
Threshold Adaptive Learning (STAL) encoder and Spiking Recurrent Neural
Network (SRNN), which are responsible for classifying data from EmoPain [3]
into Healthy and Chronic Lower Back Pain (CLBP) classes. The multimodal
biosignals consist of sEMG sensors and derived features (i.e., Joint Energy
and Joint Angle) from IMU sensors. These continuous signals convert into spike
trains using the STAL encoder and are subsequently classified using an
ensemble of SRNNs. The ensemble prediction is generated using a Random Forest
meta classifier.

A critical component of using SNNs in biosignal processing is the conversion
of continuous signals into discrete spike trains, known as spike encoding.
This process is particularly challenging for EMG and IMU data due to their
rapidly changing patterns and susceptibility to noise [7]. Several encoding
methods have been proposed, each with distinct trade-offs in temporal
resolution, noise robustness, and computational efficiency.

The use of rate-based encoding, though energy-efficient, often faces
challenges in capturing the rapid changes inherent in EMG and IMU signals
[15]. Conversely, phase coding provides noise resilience at the cost of
complex synchronisation mechanisms [16]. Other techniques, including Send-on-
Delta, Ben’s spiker algorithm [17], and Variational mode decomposition (VMD)
[18], have their unique advantages in precise spike generation or direct
encoding of analogue signals but frequently encounter difficulties in handling
the intricate dynamics of _in vivo_ EMG and IMU data and fail to fully
leverage the spatial information in high-density EMG arrays.

Recent research has explored learnable encoders for multimodal biosignal
analysis in SNNs [19], offering adaptability to specific signal
characteristics. However, these approaches often need more focus on high
temporal resolution, noise robustness, and computational efficiency, which are
crucial for real-world CLBP applications using wearable sensors [20].

The limitations of existing encoding methods in handling the unique challenges
of EMG and IMU data for CLBP classification motivate the development of a new
approach. An ideal encoder for this purpose should retain temporal dynamics,
resist inherent noise, and adjust to signal characteristics without violating
biological constraints. Hence, we introduce Spike Threshold Adaptive Learning
in this study, which is designed to address the unique requirements of CLBP
classification in the EmoPain dataset.

##  3 Proposed Method

This section details the proposed two-stage approach for classifying chronic
lower back pain using biosignals. In the first stage, a novel learnable
encoder – Spike Threshold Adaptive Learning (STAL) – converts biosignals
(sEMG, Joint Angle, and Joint Energy) into spike trains. In the second stage,
independent spiking recurrent neural networks are used within a multi-stream
ensemble classification framework to analyse the encoded spike trains. This
ensemble approach combines the predictions from each modality-specific SRNN
using a Random Forest classifier to perform classification. Figure 1
illustrates a schematic of the proposed architecture.

Let $\mathbf{X}\in\mathbb{R}^{T\times c}$ denote the input data, where $T$ is
the number of time steps and $c$ is the number of channels. To maintain
consistency in notation, we refer to the biosignal input dimensions as
“channels,” irrespective of whether they represent physical sensor channels or
derived features. The data is segmented into windows of size $\omega$ time
steps, resulting in a collection of windowed data points denoted by
$\mathbf{X}^{(\omega)}\in\mathbb{R}^{\omega\times c}$ for each window
$\omega$.

###  3.1 STAL: Spike Threshold Adaptive Learning

The STAL architecture, depicted in Fig. 2, consists of feature extraction and
feature-to-spike conversion modules. The feature extraction module comprises
two consecutive blocks, each containing a dense layer, dropout, ReLU
activation, and batch normalisation. The ReLU activation function
($\mathrm{ReLU}(\cdot)=\max(0,\cdot)$) introduces non-linearity, allowing the
network to learn complex patterns. Dropout helps prevent overfitting by
randomly setting a fraction of input units to 0 during training. Batch
normalisation improves the stability and performance of the neural network by
normalising the inputs to each layer.

Figure 2: The STAL architecture consists of feature extraction and feature-to-
spike conversion modules. Here, $\mathbf{\mathcal{L}}$ is the encoder loss
function, formulated in Eq. (3.1).

To prepare for feature extraction, we reshape $\mathbf{X}^{(\omega)}$ into a
feature vector $\mathbf{X}_{f}\in\mathbb{R}^{1\times\omega c}$. The first
block processes $\mathbf{X}_{f}$ through a dense layer $\mathbf{Z}^{(1)}$,
resulting in $\mathbf{Z}^{(1)}\in\mathbb{R}^{1\times\ell}$, where $\ell$
denotes the number of neurons in the layer. The second block applies similar
operations to the output of the first block, producing
$\mathbf{Z}^{(2)}\in\mathbb{R}^{1\times\ell}$.

The feature-to-spike conversion module of the STAL converts the extracted
features into surrogate spiking activity. This process begins with an element-
wise repetition operation on the output of the second block using Eq. (1),
expanding it to $\mathbf{h}\in\mathbb{R}^{1\times\omega\psi c}$, where $\psi$
represents the number of spikes allowed at each time step.

$\mathbf{h}=\mathrm{Repeat}(\mathbf{Z}^{(2)},\psi)$ (1)

A learnable threshold vector, denoted as
$\mathbf{\Phi}\triangleq\left[\phi_{1},\phi_{2},\cdots,\phi_{\omega\psi
c}\right]^{\top}$, is incorporated into the encoding function
$f:\mathbb{R}\rightarrow\\{0,1\\}$. This function is designed as a
differentiable surrogate in Eq. (2) to generate binary spike trains
$\mathbf{B}\in\\{0,1\\}^{\omega\times\psi\times c}$. Each threshold, denoted
by $\phi_{i}\sim\mathcal{U}(a,b),~{}\forall i\in\\{1,2,\cdots,\omega\psi
c\\}$, is initialised from a uniform distribution with hyperparameters
$a,b\in\mathbb{R}^{+}$.

$f_{\mathrm{surr}}(\mathbf{h};\mathbf{\Phi})=\sigma\left(\alpha\left(\mathbf{h}-\mathbf{\Phi}\right)\right)$
(2)

where $\sigma(\cdot)$ is the sigmoid function and $\alpha$ controls the slope
of the surrogate function. When $\alpha\rightarrow\infty$, the function
approaches a true binary step function, mimicking the behaviour of real
spiking neurons.

To align the dimensions of the spike train with the input for subsequent
computations, we reshape $\mathbf{B}$ into
$\mathbf{\hat{B}}\in\\{1,2,\cdots,\mathcal{K}(\psi)\\}^{1\times\omega c}$
using Eq. (3). Here, $\mathcal{K}(\psi)$ denotes the number of possible
combinations arising from summing the activations across $\psi$ spike
positions at each time step and channel, capturing the temporal information
about the order of spikes within $\mathbf{\hat{B}}$.

$\mathbf{\hat{B}}_{i}=\sum_{j=1}^{\psi}p_{j}\mathbf{B}_{i,j,k},\quad\forall
i\in\\{1,\cdots,\omega\\},k\in\\{1,\cdots,c\\}$ (3)

where $p_{j}$ represents the weight associated with the $j$-th spike position.
This process incorporates the temporal information about the order of spikes
into the final representation, as different spike orders will lead to
different weighted sums and, consequently, different entries in
$\mathbf{\hat{B}}$.

We introduce a customised loss function to update the parameters of the STAL
encoder. This function consists of the Mutual Information term, inspired by
[21], and the Sparsity term. The Mutual Information term
($\mathcal{L}_{\mathcal{MI}}$) measures the extent to which the encoder
captures relevant information from the input data ($\mathbf{X}_{f}$) at
multiple levels of representation, i.e., the encoded information in the hidden
layers ($\mathbf{Z}^{(1)}$ and $\mathbf{Z}^{(2)}$) and encoded spike train
($\mathbf{\hat{B}}$). Meanwhile, the Sparsity term
($\mathcal{L}_{\mathcal{S}}$) serves as a regularisation term, encouraging the
generation of sparse spike trains that efficiently represent the input data,
promoting resource efficiency and robustness against overfitting. Formally,
the loss function is formulated as in Eq. (3.1).

$\displaystyle\mathcal{L}_{\mathcal{MI}}$
$\displaystyle=-\dfrac{1}{3}\left(\mathcal{I}\left[\mathbf{X}_{f};\mathbf{Z}^{(1)}\right]+\mathcal{I}\left[\mathbf{X}_{f};\mathbf{Z}^{(2)}\right]+\mathcal{I}\left[\mathbf{X}_{f};\mathbf{\hat{B}}\right]\right),$
$\displaystyle\mathcal{L}_{\mathcal{S}}$
$\displaystyle=\lambda\cdot||\mathbf{X}_{f}-\mathbf{\hat{B}}||_{1},$
$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{L}_{\mathcal{MI}}+\mathcal{L}_{\mathcal{S}}$ (4)

where $\mathcal{I}\left[\cdot;\cdot\right]$ computes the Mutual Information
between tensors, $||\cdot||_{1}$ is the L1 norm, and $\lambda$ is the strength
coefficient controlling the level of sparseness imposed on the internal
representations.

To calculate Mutual Information, the encoded information in the hidden layers
and spike train need to have the same dimensions as the input data. We reshape
the hidden layer outputs when their dimensions differ. If $\ell<\omega c$, the
features will be repeated $\left\lfloor\frac{\ell}{\omega c}\right\rfloor$
times. Conversely, if $\ell>\omega c$, average pooling will be used with a
pool size and stride of $\left\lfloor\frac{\omega c}{\ell}\right\rfloor$.

The role of $\mathcal{L}_{\mathcal{S}}$ in regulating the sparsity of the
spike train is essential. It uses L1 norm regularisation to calculate the
difference between the input data ($\mathbf{X}_{f}$) and the surrogate spike
train ($\mathbf{\hat{B}}$). This difference serves as an approximation of the
number of spikes in the spike train. To limit the number of spikes to align
approximately with the _area_ of the input signal, Eq. (5) defines the impact
of the sparsity term on the training process. This control over sparsity is
achieved by taking the derivative of the sparsity term with respect to
$\mathbf{\hat{B}}$.

$\dfrac{\partial\mathcal{L}_{\mathcal{S}}}{\partial\mathbf{\hat{B}}}=-\lambda\cdot\mathrm{sign}(\mathbf{X}_{f}-\mathbf{\hat{B}}).$
(5)

where $\mathrm{sign}(\cdot)$ is the sign function. This derivative encourages
the output values to be close to 0 or 1, leading to a sparse output:

  1. 1.

When $\mathbf{\hat{B}}>\mathbf{X}_{f}$, the derivative is negative, pushing
output values towards 0.

  2. 2.

When $\mathbf{\hat{B}}<\mathbf{X}_{f}$, the derivative is positive, pushing
output values towards 1.

  3. 3.

When $\mathbf{\hat{B}}=\mathbf{X}_{f}$, the derivative is 0, leaving output
values unchanged.

###  3.2 Spiking Recurrent Neural Network

The proposed spiking recurrent neural network, depicted in Fig. 3, leverages
the temporal dynamics of spike trains for chronic lower back pain
classification. The SRNN processes the output of the STAL encoder
($\mathbf{\hat{B}}$) through a network of Recurrent Leaky Integrate-and-Fire
(R-LIF) neurons.

The network comprises $L$ layers of R-LIF neurons, where $L=2$ in our current
implementation. The first layer (i.e., hidden layer) contains 500 neurons, and
the second layer (i.e., readout layer) contains two neurons corresponding to
the healthy and CLBP classes. All layers implement recurrent connections,
allowing the network to capture complex temporal dependencies in the input
spike trains.

Each R-LIF neuron is implemented using an exponentially decaying membrane
potential, represented in discrete time steps. For each layer
$l\in\\{1,\cdots,L\\}$, the membrane potential $\mathbf{U}^{l}[t]$ and output
spikes $\mathbf{S}^{l}[t]$ at time step $t$ are computed by Eq. (3.2).

$\displaystyle\mathbf{U}^{l}[t+1]$
$\displaystyle=\beta\mathbf{U}^{l}[t]+\mathbf{I}^{l}[t+1]+$
$\displaystyle\mathbf{V}^{l}(\mathbf{S}^{l}[t])-\mathbf{R}^{l}[t]\mathbf{U}_{\mathrm{thr}},$
$\displaystyle\mathbf{S}^{l}[t+1]$
$\displaystyle=\Theta(\mathbf{U}^{l}[t+1]-\mathbf{U}_{\mathrm{thr}}),$ (6)

where $\beta=0.99$ is the membrane potential decay rate chosen to enhance the
LIF neurons’ capacity to capture temporal dependencies in input spike trains
by promoting longer memory. $\mathbf{I}^{l}[t]$ is the input current for layer
$l$. The recurrent connections, $\mathbf{V}^{l}(\mathbf{S}^{l}[t])$, capture
the impact of spikes from all interconnected neurons within layer $l$. The
binary vector $\mathbf{S}^{l}[t]$ represents the output spikes of all neurons
in layer $l$ at time $t$. We set the membrane threshold
$\mathbf{U}_{\mathrm{thr}}$ to 1 in this study. The binary reset-by-
subtraction mechanism $\mathbf{R}^{l}[t]$ is set to 1 when the threshold is
reached and 0 otherwise, effectively subtracting $\mathbf{U}_{\mathrm{thr}}$
(which is 1 in this case) from the membrane potential. Lastly, $\Theta(\cdot)$
refers to the Heaviside step function.

Figure 3:  The SRNN architecture. The spike trains $\mathbf{\hat{B}}$ of each
data type (sEMG, Energy, Angle) are classified separately. $T$ represent the
total time steps, $c$ are the channels, $\tau$ is the time dimension of the
Recurrent Leaky Integrate-and-Fire (R-LIF) neurons.

All layers implement all-to-all connectivity. The input currents and recurrent
connections for each layer are computed using Eq. (3.2).

$\displaystyle\mathbf{I}^{1}[t]$
$\displaystyle=\mathbf{W}^{\\{0,1\\}}\mathbf{\hat{B}}[t],$
$\displaystyle\mathbf{I}^{l}[t]$
$\displaystyle=\mathbf{W}^{\\{l-1,l\\}}\mathbf{S}^{l-1}[t],\quad\text{for
}l>1,$ $\displaystyle\mathbf{V}^{l}(\mathbf{S}^{l}[t])$
$\displaystyle=\mathbf{W}_{\mathrm{rec}}^{l}\mathbf{S}^{l}[t],\quad\text{for
}l\in{1,\cdots,L},$ (7)

where $\mathbf{W}^{\\{0,1\\}}\in\mathbb{R}^{n_{1}\times c}$,
$\mathbf{W}^{\\{l-1,l\\}}\in\mathbb{R}^{n_{l}\times n_{l-1}}$, and
$\mathbf{W}_{\mathrm{rec}}^{l}\in\mathbb{R}^{n_{l}\times n_{l}}$ are learnable
weight matrices. Here, $n_{l}$ represents the number of neurons in layer $l$,
with $n_{1}=500$ and $n_{2}=2$ in our architecture, and $\mathbf{\hat{B}}[t]$
represents the input at time step $t$.

The SRNN processes the input over $\tau$ time steps, where $\tau=5$ is chosen
to capture the temporal dynamics of the R-LIF neurons. The network’s output is
computed using Eq. (3.2), which accumulates the spikes and membrane potentials
of the output layer over all time steps.

$\displaystyle\mathbf{S}_{\mathrm{out}}$
$\displaystyle=\left[\mathbf{S}^{L}[1],\mathbf{S}^{L}[2],\cdots,\mathbf{S}^{L}[\tau]\right],$
$\displaystyle\mathbf{U}_{\mathrm{out}}$
$\displaystyle=\left[\mathbf{U}^{L}[1],\mathbf{U}^{L}[2],\cdots,\mathbf{U}^{L}[\tau]\right],$
(8)

where the final classification is determined by the neuron with the highest
cumulative activity across these $\tau$ time steps in the output layer.

To train the SRNN, we use the cross-entropy loss function, which is
particularly suitable for our binary classification task between healthy and
CLBP classes. The cross-entropy loss $\mathcal{L}$ is defined as in Eq. (9).

$\mathcal{L}_{\mathrm{SRNN}}=-\sum_{i=1}^{C}y_{i}\log(\hat{y}_{i})$ (9)

where $C$ is the number of classes (in our case, $C=2$), $y_{i}$ is the true
label (0 or 1) for class $i$, and $\hat{y}_{i}$ is the predicted probability
for class $i$. The predicted probabilities are obtained by applying a softmax
function to the accumulated output $\mathbf{S}_{\mathrm{out}}$ and
$\mathbf{U}_{\mathrm{out}}$. To handle the non-differentiable nature of spikes
during backpropagation, we employ a surrogate gradient function. Specifically,
we use the gradient of the shifted arc-tan function in Eq. (10).

$\frac{\partial S}{\partial U}=\frac{1}{\pi}\frac{1}{\left(1+\left(\pi
U\frac{\alpha}{2}\right)^{2}\right)}.$ (10)

This surrogate gradient allows us to backpropagate the loss through the
network and update the weights of the SRNN, enabling effective training
despite the discrete nature of spike events.

###  3.3 Multi-Stream Ensemble Classification

We have implemented a multi-stream ensemble classification approach using
separate STAL-SRNN pipelines to leverage the diverse information captured by
various biosignal modalities. This tailored approach provides customised
encoding and classification strategies for each type of data to compensate for
the information loss during continuous-to-spike conversion and leverages
cross-modal interactions. For each data modality
$m\in\\{\mathrm{sEMG},\mathrm{Angle},\mathrm{Energy}\\}$, we train an
independent STAL-SRNN pipeline formulated in Eq. (11).

$\mathbf{Y}_{m}=\mathcal{P}_{m}(\mathbf{\hat{B}}_{m};\boldsymbol{\theta}_{m}),$
(11)

where $\mathbf{\hat{B}}_{m}$ is the output of the STAL encoder for modality
$m$, as defined in Eq. (3). Here, $\mathcal{P}_{m}(\cdot)$ denotes the STAL-
SRNN pipeline for that modality, $\boldsymbol{\theta}_{m}$ represents the
learnable parameters of the pipeline, and $\mathbf{Y}_{m}\in\mathbb{R}^{2}$ is
the output prediction (i.e., probabilities for Healthy and CLBP classes) for
modality $m$. Each STAL-SRNN pipeline is trained independently using the
cross-entropy loss function, as described in Section 3.2.

The outputs from these individual pipelines serve as features for the meta-
classifier, which is implemented using a Random Forest classifier [22]. This
classifier is comprised of multiple decision trees, each trained on a
bootstrap sample of the training data using the Gini impurity function to
measure the quality of a split. The final prediction is determined by majority
voting among the trees using Eq. (12).

$\hat{y}=\underset{C}{\mathrm{mode}}\left\\{\mathrm{RF}_{i}\left(\mathbf{Y}_{\mathrm{sEMG}},\mathbf{Y}_{\mathrm{Angle}},\mathbf{Y}_{\mathrm{Energy}}\right)\right\\}_{i=1}^{N_{t}},$
(12)

where $\mathrm{RF}_{i}(\cdot)$ represents the prediction of the $i$-th
decision tree in the Random Forest. The final prediction $\hat{y}$ is
determined by majority voting among the $N_{t}$ trees.

##  4 Experiments and Results

###  4.1 Dataset and Preparation

The EmoPain dataset [3] is a multimodal collection of data from 46
participants, encompassing healthy individuals and those with CLBP.
Participants performed various exercises (one-leg-stand, reach-forward, stand-
to-sit, sit-to-stand, and bend-down) at normal and difficult levels, with
transitions included to simulate real-life movements.

Data were collected using 18 wearable IMU sensors and four sEMG sensors. The
IMU sensors recorded full-body 3D motion, from which 13 joint angle data
(Joint Angle) and 13 angular energy data (Joint Energy) were derived. The sEMG
sensors captured muscle activity data from the right and left upper and lower
back muscle groups. The dataset includes self-reported pain intensity ratings
(0-10) from CLBP individuals and labels for six pain-related behaviours
provided by physiotherapists. We assigned positive labels to individuals
reporting pain intensity greater than 5, resulting in 12 positive and 34
negative labels.

To prepare the data for analysis, we addressed the variation in recording
lengths by trimming recordings that exceeded 18,000 frames by 6,000 frames
from both ends. This method allowed us to balance standardisation with data
preservation, reducing information loss and the impact of stress and fatigue
[5]. We used normalisation techniques for all data modalities (sEMG, Angles,
and Energies) to ensure consistent scaling and facilitate model training.
Specifically, we normalised values to a 0-1 range by subtracting the minimum
value and dividing by the maximum value within each sample. Additionally, we
zero-padded the data to match the longest recording length (17,995 time
steps). Lastly, we employed a sliding window approach with a 3,000 time-step
window (50 seconds at 60 Hz), allowing the STAL to capture longer-term
temporal dependencies and enhance robustness to noise.

###  4.2 Implementation Details

The experiments were conducted on a high-performance computing system equipped
with an AMD EPYC 7302 16-core processor, 504GB of RAM, and an NVIDIA RTX A6000
GPU. The deep learning framework used was PyTorch version 2.3.1 [23]. The
AdamW optimiser [24] was employed for both the STAL encoder and the SRNN
classifier, with fixed learning rates of $5.0\times 10^{-3}$ and $7.5\times
10^{-4}$, respectively, due to its effectiveness in handling sparse gradients,
a common characteristic in SNN training.

The number of neurons in the hidden layers was chosen based on the size of the
sliding window. Therefore, both $\mathbf{Z^{(1)}}$ and $\mathbf{Z^{(2)}}$ were
designed with 3000 neurons each to effectively process the temporal
information within the data. The SRNN classifier’s hidden layer was composed
of 500 recurrent LIF neurons. This decision was deliberate to find a balance
between computational efficiency, model performance, and generalisability.
Additionally, to ensure that the model learns from all classes in the context
of this imbalanced dataset with a small sample size and to mitigate any
biasing effects, we applied minority over-sampling with weighted sample
replacement during batch creation.

Table 1: performance of four encoders on the EmoPain Database. |  |  |  |  |  | Spike density $\downarrow$
---|---|---|---|---|---|---
Encoder | Classifier | Accuracy $\uparrow$ | AUC $\uparrow$ | F1 $\uparrow$ | MCC $\uparrow$ | Ensemble | sEMG | Energy | Angle
Baseline methods
Rate-coding [17] | SRNN | 73.91% | 0.500 | 0.000 | 0.000 | 0.167 | 0.068 | 0.584 | 0.703
Latency-coding [17] | SRNN | 71.74% | 0.512 | 0.133 | 0.044 | 0.200 | 0.200 | 0.200 | 0.200
Proposed method
STAL–Stacked | SRNN | 80.43% | 0.679 | 0.526 | 0.437 | 0.549 | 0.360 | 0.703 | 0.793
STAL–Vanilla | SRNN | 76.09% | 0.596 | 0.353 | 0.270 | 0.021 | 0.007 | 0.439 | 0.528

An empirical study identified optimal hyperparameters for our STAL-SRNN
architecture. The best configuration (see Fig. 4) included a batch size of 32
for sEMG, 8 for Energy, and 16 for Angle, along with a dropout rate of 0.5 and
$\psi=5$.

Figure 4: (a) The influence of batch size on the F1-score of the proposed
STAL-SRNN for each modality. (b) The Spike density for the proposed
architecture as a function of $(\psi,p)$. While the difference in spike
density between $(\psi,p)=(5,0.5)$ and $(\psi,p)=(10,0.5)$ is less than 1%
(0.549 vs 0.540), the former achieves better performance across other
classification metrics.

The training process consisted of 30 epochs for the STAL encoder and 25 epochs
for the SRNN classifier, with early stopping implemented to prevent
overfitting. This training regime was carefully chosen to balance model
convergence and computational efficiency. All the hyperparameters and settings
determined through this empirical search were subsequently used consistently
across all experiments in this study. This approach ensures the
reproducibility of our results and allows for a fair comparison across
different experimental conditions111The implementation is accessible for
reproducibility purposes on https://github.com/freek1/emopain-stl..

###  4.3 Evaluation Metrics

We used both standard classification metrics and a sparsity measure tailored
to SNNs to evaluate the effectiveness of the STAL-SRNN framework for CLBP
classification on the EmoPain dataset. All metrics were assessed using leave-
one-subject-out cross-validation to ensure a rigorous analysis, especially
given the relatively small size of the EmoPain dataset.

The standard classification metrics employed in this study encompassed
accuracy, Macro F1 score, Area Under the Curve (AUC), and Matthews Correlation
Coefficient (MCC). These metrics were suggested in the EmoPain challenge [25]
for comprehensive insights into the model’s performance across various
aspects.

In addition to these standard metrics, we introduced a spike density measure
to assess the efficiency of our encoded spike trains. This measure provides
insight into the sparsity of our neural representations across our multimodal
approach. The overall spike density is computed as the harmonic mean of the
individual densities for each data modality, as shown in Eq. (4.3).

$\displaystyle\mathrm{SpikeDensity}$
$\displaystyle=\frac{3}{\sum_{m}\frac{1}{D_{m}}},$ $\displaystyle D_{m}$
$\displaystyle=\frac{\sum_{i=1}^{T}\sum_{j=1}^{\psi}\sum_{k=1}^{c_{m}}\mathbf{B}_{m,i,j,k}}{T\psi
c_{m}},$ (13)

where, $\mathbf{B}_{m}\in\\{0,1\\}^{T\times\psi\times c_{m}}$ represents the
spike train for modality $m$, summed over $T$ time steps, $c_{m}$ channels,
and $\psi$ spikes per time step. The denominator represents the maximum
possible number of spikes for that modality.

###  4.4 Experimental Results

Our analysis focuses on three main aspects: (1) the performance of different
spike train encoding strategies with the SRNN classifier, (2) an ablation
study to understand the impact of architectural choices and hyperparameters,
and (3) a comparison with state-of-the-art deep learning models for chronic
pain assessment.

####  4.4.1 Spike Train Encoders

We conducted an assessment of four spike train encoders using the SRNN
classifier with leave-one-subject-out cross-validation on the EmoPain dataset.
The term STAL–Stacked refers to the proposed encoder outlined in Section 3.1.
STAL–Vanilla represents a simplified version of STAL–Stacked, wherein the
input is directly connected to the Repeat layer, excluding the two feature
extraction blocks. Rate-coding and Latency-coding are considered baseline
methods that exclusively encode information in the spatial and temporal
domains, respectively. The experimental results are presented in Table 1, and
Fig. 5 shows the confusion matrices for each encoder.

Figure 5: Confusion matrices of the encoders with the SRNN classifier.

The results indicate that STAL–Stacked outperformed all other models,
achieving the highest accuracy of 80.43%, AUC of 0.679, F1 score of 0.526, and
Matthews Correlation Coefficient (MCC) of 0.437. This suggests that
STAL–Stacked offers the most robust and balanced performance among the
encoders tested. STAL–Vanilla emerged as the second-best performer, achieving
an accuracy of 76.09%, AUC of 0.596, F1 score of 0.353, and MCC of 0.270.
Notably, STAL–Vanilla achieved this with the least dense spike trains (spike
density of 0.021), indicating its potential for efficient computation.

In comparison, rate-coding and latency-coding methods showed lower performance
across all metrics. Rate-coding achieved 73.91% accuracy but did not
effectively differentiate between classes (F1 score and MCC of 0). Latency-
coding performed slightly better in terms of F1 score and MCC but had the
lowest overall accuracy at 71.74

####  4.4.2 Ablation Study

To assess the effectiveness of the STAL architecture and the impact of
different design choices, we conducted an ablation study comparing
STAL–Stacked and STAL–Vanilla and examining the effects of dropout
probability.

The primary distinction between STAL–Stacked and STAL–Vanilla lies in the
inclusion of feature extraction blocks. STAL–Stacked integrates two feature
extractor blocks before the Repeat layer, while STAL–Vanilla directly connects
the input to the Repeat layer. The superior performance of STAL–Stacked
(80.43% accuracy compared to 76.09% for STAL–Vanilla) suggests that these
additional layers contribute to a more effective feature representation.
Additionally, STAL–Stacked employed a higher spike density (0.549) compared to
STAL–Vanilla (0.021). This outcome indicates that the improved performance of
STAL–Stacked is attributed to a more information-rich encoding, albeit at the
cost of increased computational resources.

Both variants of the STAL method performed better than the baseline rate- and
latency-based encoding approaches. This improvement can be attributed to
STAL’s unique ability to encode spike trains both spatially and temporally at
the same time. Specifically, STAL encodes spatially by learning a distinct
threshold value $\phi_{i}\in\mathbf{\Phi}$ for each feature
$h_{i}\in\mathbf{h}$, as shown in Eq. (2). At the same time, the loss function
encourages temporal encoding of information by promoting spikes at the end of
the spike train to encode a high signal, as demonstrated by the positional
encoding in Eq. (3). This dual encoding capability enables STAL to capture
more complex patterns in the data, resulting in improved performance compared
to the baseline methods.

In contrast, rate-coding encodes information in the spatial domain by
proportionally correlating the probability sum of generated spikes for each
input value while randomly positioning the spikes in the temporal domain.
Similarly, latency-coding only encodes information temporally. This is evident
from the uniform distribution in the spatial domain, where 1 spike is
generated for each input, and the spike is strategically positioned in the
temporal domain.

In our ablation study, we focused on examining the impact of dropout
probability on the performance of our models. We found that in the
STAL–Stacked architecture, a higher dropout probability $(p=0.5)$ had a
positive effect on performance. This can be attributed to the presence of
fully connected layers in the feature extraction blocks, which facilitated
information distribution during dropout. Conversely, in STAL–Vanilla, a lower
dropout probability $(p=0.0)$ was found to be optimal. This is because the
absence of feature extraction layers in this architecture made it more
susceptible to information loss during dropout.

####  4.4.3 Comparison with Deep Learning Models

To contextualise our results within the broader landscape of pain assessment
methods, we compared our STAL-SRNN approach with state-of-the-art deep
learning models on the EmoPain dataset (see Table 2). This comparison serves
multiple purposes:

  1. 1.

It establishes a performance baseline for spiking neural networks in chronic
pain assessment, an area predominantly explored using traditional deep
learning methods.

  2. 2.

It highlights the trade-offs between deep learning models’ high performance
and the potential efficiency advantages of spiking neural networks.

  3. 3.

It shows that despite the inherent differences in input processing (continuous
vs. spike-based), our STAL-SRNN can achieve competitive results in certain
aspects.

Table 2: Comparison of the proposed spiking architecture with state-of-the-art deep learning models on the EmoPain dataset. Methods | AUC | F1 | MCC
---|---|---|---
Deep Learning
MiMT [26] | 0.648 | 0.722 | 0.327
LSTM+GCN [27] | 0.690 | 0.731 | 0.365
GRU-RNN [7] | 0.855 | 0.765 | 0.411
L-SFAN [9] | 0.849 | 0.772 | 0.510
Proposed Spiking Method
STAL-SRNN | 0.679 | 0.526 | 0.437

Our STAL-SRNN demonstrates competitive performance on the EmoPain dataset,
achieving an AUC of 0.679 for CLBP classification. While this AUC is lower
than the deep learning approaches (GRU-RNN: 0.855, L-SFAN: 0.849), it is
important to highlight that our method outperforms the GRU-RNN model in MCC
(0.437 vs 0.411). The competitive MCC suggests that STAL-SRNN has the
potential for balanced classification in this domain.

##  5 Conclusion and Future Work

In this study, we presented the first application of spiking neural networks
for chronic lower back pain classification using the EmoPain dataset. Our work
introduces two significant contributions to the field of neuromorphic
computing and biosignal analysis. We introduced Spike Threshold Adaptive
Learning (STAL), a novel, trainable encoder that effectively converts
continuous biosignals into spike trains. Experimental results demonstrate that
STAL could preserve crucial temporal dynamics and adapt to signal
characteristics. We proposed an ensemble of Spiking Recurrent Neural Network
(SRNN) classifiers designed to process surface Electromyography (sEMG) and
Inertial Measurement Unit (IMU) data within a multi-stream architecture. This
ensemble leverages a random forest classifier to efficiently harness the
temporal processing capabilities inherent in SRNNs.

Our STAL-SRNN framework demonstrates outstanding performance, achieving an
accuracy of 80.43%, AUC of 67.90%, F1 score of 52.60%, and Matthews
Correlation Coefficient (MCC) of 0.437. Notably, it outperforms traditional
rate-based and latency-based encoding methods and surpasses the best-
performing deep learning approach in terms of MCC, indicating better-balanced
class prediction. This achievement was particularly significant given the
challenges posed by the small sample size and class imbalance in the EmoPain
dataset.

The success of our approach could represent a significant step forward in
applying neuromorphic computing techniques to chronic pain management
challenges. By leveraging the inherent efficiency of SNNs, our method bridged
the gap between efficient spike-based computation and complex biosignal
analysis. The STAL-SRNN provides new possibilities for personalised,
continuous health monitoring and intervention in chronic pain management.

Future work could explore extending the framework to simultaneously classify
pain levels and identify specific pain-related behaviours, providing more
comprehensive insights for pain management. Furthermore, exploring the
deployment of the STAL-SRNN framework on neuromorphic hardware could fully
leverage its energy efficiency for long-term, wearable applications.

## References

  * [1] A. Supratak, C. Wu, H. Dong, K. Sun, and Y. Guo, _Survey on Feature Extraction and Applications of Biosignals_. Springer International Publishing, 2016, pp. 161–182.
  * [2] N. Li, R. Zhou, B. Krishna, A. Pradhan, H. Lee, J. He, and N. Jiang, “Non-invasive Techniques for Muscle Fatigue Monitoring: A Comprehensive Survey,” _ACM Computing Survey_ , vol. 56, no. 9, 2024.
  * [3] M. S. H. Aung, S. Kaltwang, B. Romera-Paredes, B. Martinez, A. Singh, M. Cella, M. Valstar, H. Meng, A. Kemp, M. Shafizadeh, A. C. Elkins, N. Kanakam, A. De Rothschild, N. Tyler, P. J. Watson, A. C. D. C. Williams, M. Pantic, and N. Bianchi-Berthouze, “The Automatic Detection of Chronic Pain-Related Expression: Requirements, Challenges and the Multimodal EmoPain Dataset,” _IEEE Transactions on Affective Computing_ , vol. 7, no. 4, pp. 435–451, 2016.
  * [4] R. J. Yong, P. M. Mullins, and N. Bhattacharyya, “Prevalence of chronic pain among adults in the united states,” _Pain_ , vol. 163, no. 2, pp. e328–e332, 2022.
  * [5] R. Fernandez Rojas, N. Brown, G. Waddington, and R. Goecke, “A systematic review of neurophysiological sensing for the assessment of acute pain,” _npj Digital Medicine_ , vol. 6, no. 1, p. 76, 2023.
  * [6] N. Makaram, P. A. Karthick, and R. Swaminathan, “Analysis of Dynamics of EMG Signal Variations in Fatiguing Contractions of Muscles Using Transition Network Approach,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 70, pp. 1–8, 2021.
  * [7] M. M. Dehshibi, T. Olugbade, F. Diaz-de Maria, N. Bianchi-Berthouze, and A. Tajadura-Jiménez, “Pain Level and Pain-Related Behaviour Classification Using GRU-Based Sparsely-Connected RNNs,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 17, no. 3, pp. 677–688, 2023.
  * [8] E. Farago, D. MacIsaac, M. Suk, and A. D. C. Chan, “A Review of Techniques for Surface Electromyography Signal Quality Analysis,” _IEEE Reviews in Biomedical Engineering_ , vol. 16, pp. 472–486, 2023.
  * [9] J. Ortigoso-Narro, F. D. de Maria, M. M. Dehshibi, and A. Tajadura-Jiméne, “L-SFAN: Lightweight Spatially-focused Attention Network for Pain Behavior Detection,” pp. 1–9, 2024.
  * [10] F. Wang, W. M. Severa, and F. Rothganger, “Acquisition and Representation of Spatio-Temporal Signals in Polychronizing Spiking Neural Networks,” in _Proceedings of the 7th Annual Neuro-Inspired Computational Elements Workshop_. Association for Computing Machinery, 2019.
  * [11] A. Sun, X. Chen, M. Xu, X. Zhang, and X. Chen, “Feasibility study on the application of a spiking neural network in myoelectric control systems,” _Frontiers in Neuroscience_ , vol. 17, p. 1174760, 2023.
  * [12] D. Auge, J. Hille, E. Mueller, and A. Knoll, “A survey of encoding techniques for signal processing in Spiking Neural Networks,” _Neural Processing Letters_ , vol. 53, no. 6, pp. 4693–4710, 2021.
  * [13] N. Perez-Nieves and D. F. M. Goodman, “Sparse Spiking Gradient Descent,” in _Advances in Neural Information Processing Systems_ , 2021, pp. 11 795–11 808.
  * [14] M. Shahsavari, D. Thomas, A. Brown, and W. Luk, “Neuromorphic Design Using Reward-based STDP Learning on Event-Based Reconfigurable Cluster Architecture,” in _International Conference on Neuromorphic Systems_. Association for Computing Machinery, 2021.
  * [15] J. M. Fabian, D. C. O’Carrol, and S. D. Wiederman, “Sparse spike trains and the limitation of rate codes underlying rapid behaviours,” _Biology Letters_ , vol. 19, no. 5, p. 20230099, 2023.
  * [16] M. Shahsavari, D. Thomas, M. van Gerven, A. Brown, and W. Luk, “Advancements in spiking neural network communication and synchronization techniques for event-driven neuromorphic systems,” _Array_ , vol. 20, p. 100323, 2023.
  * [17] J. K. Eshraghian, M. Ward, E. O. Neftci, X. Wang, G. Lenz, G. Dwivedi, M. Bennamoun, D. S. Jeong, and W. D. Lu, “Training Spiking Neural Networks Using Lessons From Deep Learning,” _Proceedings of the IEEE_ , vol. 111, no. 9, pp. 1016–1054, 2023.
  * [18] E. Donati, M. Payvand, N. Risi, R. Krause, and G. Indiveri, “Discrimination of EMG Signals Using a Neuromorphic Implementation of a Spiking Neural Network,” _IEEE Transactions on Biomedical Circuits and Systems_ , vol. 13, no. 5, pp. 795–803, 2019.
  * [19] P. Gong, P. Wang, Y. Zhou, and D. Zhang, “A Spiking Neural Network With Adaptive Graph Convolution and LSTM for EEG-Based Brain-Computer Interfaces,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 31, pp. 1440–1450, 2023.
  * [20] R. D. Gurchiek, N. Donahue, N. M. Fiorentino, and R. S. McGinnis, “Wearables-Only Analysis of Muscle and Joint Mechanics: An EMG-Driven Approach,” _IEEE Transactions on Biomedical Engineering_ , vol. 69, no. 2, pp. 580–589, 2021.
  * [21] A. El Ferdaoussi, J. Rouat, and E. Plourde, “Efficiency metrics for auditory neuromorphic spike encoding techniques using information theory,” _Neuromorphic Computing and Engineering_ , vol. 3, no. 2, p. 024007, 2023.
  * [22] T. K. Ho, “Random decision forests,” in _Proceedings of 3rd International Conference on Document Analysis and Recognition_ , vol. 1. IEEE, 1995, pp. 278–282.
  * [23] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer, “Automatic differentiation in PyTorch,” in _NIPS Autodiff Workshop_ , 2017.
  * [24] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” in _International Conference on Learning Representations_ , 2019, pp. 1–18.
  * [25] J. O. Egede, S. Song, T. A. Olugbade, C. Wang, C. D. C. Amanda, H. Meng, M. Aung, N. D. Lane, M. Valstar, and N. Bianchi-Berthouze, “Emopain challenge 2020: Multimodal pain evaluation from facial and bodily expressions,” in _2020 15th IEEE International Conference on Automatic Face and Gesture Recognition (FG 2020)_. IEEE, 2020, pp. 849–856.
  * [26] T. Olugbade, N. Gold, A. C. d. C. Williams, and N. Bianchi-Berthouze, “A Movement in Multiple Time Neural Network for Automatic Detection of Pain Behaviour,” in _International Conference on Multimodal Interaction_. Association for Computing Machinery, 2021, pp. 442––445.
  * [27] C. Wang, Y. Gao, A. Mathur, A. C. De C. Williams, N. D. Lane, and N. Bianchi-Berthouze, “Leveraging Activity Recognition to Enable Protective Behavior Detection in Continuous Data,” _Proc. ACM Interact. Mob. Wearable Ubiquitous Technol._ , vol. 5, no. 2, 2021.


# Steeper Scattered Disks Buckle Faster

Alexander Zderic JILA and Department of Astrophysical and Planetary Sciences,
CU Boulder, Boulder, CO 80309, USA<EMAIL_ADDRESS>Ann-Marie
Madigan JILA and Department of Astrophysical and Planetary Sciences, CU
Boulder, Boulder, CO 80309, USA

###### Abstract

Disks of low-mass bodies scattered by giant planets to large semi-major axis
and constant periapsis orbits are vulnerable to a buckling instability. This
instability exponentially grows orbital inclinations, raises periapsis
distances, and coherently tilts orbits resulting in clustering of arguments of
periapsis. The dynamically hot system is then susceptible to the formation of
a lopsided mode. Here we show that the timescale of the buckling instability
decreases as the radial surface density of the population becomes more
centrally dense, i.e., steeper scattered disks buckle faster. Accounting for
differential apsidal precession driven by giant planets, we find that
$\sim\\!10\,M_{\oplus}$ is sufficient for a primordial scattered disk in the
trans-Neptunian region to have been unstable if $dN\propto a^{-2.5}da$.

††software: REBOUND (Rein & Liu, 2012)

##  1 Introduction

Axisymmetric disks of bodies on eccentric orbits in a near-Keplerian potential
are collectively unstable to an out-of-plane buckling (Madigan & McCourt,
2016). The dynamics are driven by orbit-averaged torques between the massive
disk bodies and result in exponential growth of orbit inclinations, a decrease
in orbital eccentricities, and clustering in arguments of periapsis (Madigan
et al., 2018). This ‘inclination instability’ results in a hotter system which
is susceptible to the formation of a lopsided mode or clustering in longitude
of periapsis (Zderic et al., 2020, 2021).

We have suggested that this collective instability could explain the orbital
anomalies observed in the trans-Neptunian outer solar system which include
detached objects, extreme inclination orbits, and possible clustering in
arguments of periapsis and longitude of periapsis (Trujillo & Sheppard, 2014;
Batygin & Brown, 2016; Shankman et al., 2017; Lawler et al., 2017; Becker et
al., 2018; Brown & Batygin, 2019; Bernardinelli et al., 2020; Kavelaars et
al., 2020; Napier et al., 2021). Migrating giant planets interacting with
planetesimals produce massive scattered disks early in the solar system’s
evolution (Duncan & Levison, 1997; Vokrouhlický & Nesvorný, 2019). These
eccentric primordial scattered disks are susceptible to the inclination
instability through mutual orbit-averaged torques.

The gravitational torques that drive the instability are sensitive to external
perturbations. The differential precession induced in the trans-Neptunian
region by the outer giant planets reduces the strength of inter-orbit torques
and lengthens the instability timescale. Zderic & Madigan (2020) found that
the total disk mass between $100-1000$ au must have been $\gtrsim
20\,M_{\oplus}$ at some time in the past for the instability to occur. This
was based on $N$-body simulations of a primordial scattered disk with a
relatively shallow surface density profile, $dN\sim a^{-1}da$. Here we show
that the steepness of the radial mass distribution strongly affects the
timescale of the inclination instability, with more centrally dense mass
distributions having significantly faster timescales. We run a suite of
$N$-body simulations of disks undergoing the inclination instability with a
range of radial mass distributions, and find that the timescale of the
inclination instability varies with the median period within the disk. There
is a factor of $\gtrsim\\!6$ difference in timescales across the range of
power law indices explored. In Zderic & Madigan (2020), we found that
timescale of the inclination instability plays an important role in
determining the instability’s robustness to external sources of differential
orbital precession. Here we find that disks with steeper radial distributions
require less mass to be unstable to the inclination instability. For example,
Huang et al. (2022) generate a primordial scattered disk using simple Neptune-
scattering simulations with test particles finding $dN\sim a^{-2.5}da$ after
$50\,{\rm Myr}$. With this steeper distribution, $\sim\\!10\,M_{\oplus}$ is
sufficient for the inclination instability to have occurred in the primordial
trans-Neptunian region.

##  2 Simulations and Timescales

Here we run $N$-body simulations of unstable scattered disks to quantify the
effect of varying the radial mass profile. Our simulations are run using the
REBOUND $N$-body integrator framework (Rein & Liu, 2012; Tamayo et al., 2020).
We use the IAS15 integrator with interactive, massive particles. We run simple
simulations to directly compare with our previous results. The $N$-body disks
have $M_{\rm disk}=10^{-3}\,M_{\odot}$ and $N=400$ with each particle having
identical mass, $m=M_{\rm disk}/N$. The unrealistically large disk mass
accelerates the secular dynamics which reduces the computing time required to
capture the instability, and the low particle number is required for these
relatively inefficient but highly accurate simulations.111Fixed-timestep,
hybrid-symplectic integrators fail to resolve the highly eccentric orbits in
our simulations. However, walltimes for $N\gtrsim 10^{3}$ with the IAS15
integrator for our disks are prohibitively long. The artificially large
$M_{\rm disk}$ and low $N$ means we cannot add the giant planets to the
simulation as the back-reaction of the disk on the giant planets would be
unphysically large. We discuss how these results apply to the actual
primordial trans-Neptunian population, including the influences of the giant
planets, in Section 4.

The initial disk orbital parameters are as follows: the periapsis distances
are $30\,{\rm au}$, inclinations ($i$) are drawn from a Rayleigh distribution
with a mean inclination of $5^{\circ}$, and the arguments of periapsis
($\omega$), longitudes of the ascending node ($\Omega$), and mean anomalies
are drawn from a uniform distribution in $[0,360^{\circ})$. The semi-major
axes of the orbits are drawn from a 1D power law, $dN\propto a^{-\alpha}da$,
in the range $[a_{i},10a_{i}]$ where $a_{i}=100\,{\rm au}$ (with constant
periapsis, $e\in[0.7,0.97]$). As the disks are composed of eccentric orbits,
the semi-major axis distribution is different from the radial distribution.
However, $a^{-\alpha}\approx r^{-\alpha}$ outside of the apoapsis of the
innermost orbit. The simplicity of these simulations allow us to directly
compare with our previous work and to apply these results at a later date to
exoplanet systems of varying orbital architecture. Future work applying the
instability to the outer solar system will need to incorporate scattered disks
created by the four (or more) migrating gas and ice giants.

Figure 1: Absolute value of the mean $i_{\text{b}}$ for three disks with
different semi-major axis distribution indices, $\alpha$, and identical $N$
and $M_{\rm disk}$. We’ve truncated points after the linear stage of the
instability to emphasize the growth rate differences. The timescale for the
instability clearly decreases with increasing $\alpha$ (steeper density
profiles).

The collective dynamics behind the instability are described in Madigan et al.
(2018). In brief, we can understand the mechanism by considering just two
eccentric orbits in a disk. A small vertical perturbing force will induce a
torque on an eccentric orbit which causes it to roll over its major axis. This
rolling in turn induces a torque on a nearby orbit causing it to pitch over
its semi-latus rectum. For sufficiently eccentric orbits ($e\gtrsim 0.5$), the
pitching acts to raise or lower the apocenter of the second orbit such that it
reinforces the initial vertical force. One can solve linearized equations of
motion of the system and demonstrate exponential growth of orbital
inclinations in the initially thin disk. The coherent pitching and rolling of
the orbits in the same direction means that the apocenters of the orbits all
rise above or drop below the mid-plane together. The disk buckles into a cone
or bowl shape rather than just expanding vertically. It is instructive to look
at roll and pitch angles of the orbits, $i_{\text{a}}$ and $i_{\text{b}}$,
which are related to the Kepler angles as follows (valid for prograde orbits,
see Madigan et al. (2018)):

$\displaystyle i_{\text{a}}$
$\displaystyle=\arcsin{\left(\sin{i}\cos{\omega}\right)},$ (1) $\displaystyle
i_{\text{b}}$ $\displaystyle=\arcsin{\left(-\sin{i}\sin{\omega}\right)},$ (2)
$\displaystyle i$ $\displaystyle=\arccos\left[\cos i_{\text{a}}\cos
i_{\text{b}}\left[1-(\tan i_{\text{a}}\tan
i_{\text{b}})^{2}\right]^{1/2}\right].$ (3)

The instability is characterized by the exponential growth of the mean
$i_{\text{a}}$ and $i_{\text{b}}$ in opposite directions. We define the
timescale of the instability, $t_{\rm insta}$, in our $N$-body simulations as
the inverse of the exponential growth rate of the mean $i_{\text{b}}$ of the
disk given in units of the period of the innermost orbit,
$P_{i}=2\pi\mu^{-\nicefrac{{1}}{{2}}}a_{i}^{\nicefrac{{3}}{{2}}}$, where
$\mu=GM_{\odot}$. The growth rate is given by the slope of the mean
$i_{\text{b}}$ in the linear regime, and it is obtained from with a simple
linear, least-squares fit to the log of the mean $i_{\text{b}}$. An example is
shown in Figure 1 where we show the mean $i_{\text{b}}$ for three disks with
different $\alpha$. The linear stage lasts from $\sim\\!0.2^{\circ}$ to
$\sim\\!10^{\circ}$, and instability timescale clearly varies with $\alpha$.
When the instability saturates, mean $i_{b}$ in the disks is $\sim 20^{\circ}$
while the mean inclination has reached $35^{\circ}-50^{\circ}$. Coinciding
with the buckling, the mean eccentricity of the disk drops to conserve vector
angular momentum of the disk. This drop is most pronounced at inner edge of
the disk with the mean eccentricity dropping to $0.4-0.6$ in this region. This
drop at approximately fixed semi-major axes results in an increase in mean
perihelion, $q=a(1-e)$.

##  3 Results

We run 10 $N$-body scattered disk simulations for each
$\alpha\in[0,0.5,1,1.5,2,2.5]$ and calculate the timescale of the instability
as described in Section 2. We show the measured timescales as a function of
$\alpha$ in Figure 2, where we’ve added jitter to the data to reveal
overlapping points. The instability timescale monotonically decreases with
increasing $\alpha$ (steeper density profiles). There is a factor of
$\gtrsim\\!6$ difference in timescales across the range in $\alpha$.

Figure 2: Timescale of the inclination instability, $t_{\rm insta}$, as a
function of semi-major axis distribution index, $\alpha$. $N=400$ and $M_{\rm
disk}=10^{-3}\,M_{\oplus}$ for all simulations. The instability timescale
monotonically decreases with increasing $\alpha$. Jitter is added to reveal
overlapping points.

Being a secular (orbit-averaged) phenomenon, the inclination instability
scales with the secular timescale,

$t_{\rm sec}\sim\frac{M_{\odot}}{M_{\rm disk}}\,\frac{P}{2\pi}.$ (4)

Here we define $t_{\rm sec}$ as the time it takes to change an orbit’s angular
momentum by order of its circular angular momentum using a specific torque
over one orbital period of $\tau\sim\frac{GM_{\odot}}{a}$. In practice, we
find that the instability timescale also depends on $N$ due artificially
strong two-body scattering present in our low $N$ simulations (Madigan et al.,
2018). In the limit of infinite $N$, the instability timescale would decrease
by a factor of $\sim 3$. The instability timescales shown in Figure 2 are for
disks with identical $N$ and $M_{\rm disk}$, but different orbital period
distributions, $f(P)\propto P^{-\frac{2\alpha+1}{3}}$. The variation of
instability timescale with $\alpha$ is due to the variation of $t_{\rm sec}$
with $P$.

Figure 3: Box and whisker plot of the instability timescale normalized by the
mean (top) and median (bottom) period within the disk. The median period
properly de-trends the data suggesting the that growth rate at fixed $M_{\rm
disk}$ and $N$ is determined by the median period of an orbit in the disk.

In Figure 3, we normalize the instability timescales shown in Figure 2 with
the mean (top) and median (bottom) orbital periods of the disks. We use box
and whisker plots to summarize the distribution of timescales for each
$\alpha$. The median period successfully de-trends the data to within the
inter-quartile range, suggesting that,

$t_{\rm insta}\propto t_{\rm sec}\propto P_{\rm median},$ (5)

where,

$\displaystyle P_{\rm
median}=2\pi\mu^{-\nicefrac{{1}}{{2}}}\begin{cases}\left(\frac{1}{2}\left(a_{o}^{1-\alpha}+a_{i}^{1-\alpha}\right)\right)^{\frac{3}{2(1-\alpha)}}&\alpha\neq
1\\\ (a_{o}a_{i})^{\nicefrac{{3}}{{4}}}&\alpha=1\\\ \end{cases},$ (6)

for $dN\propto a^{-\alpha}da$ in the range $[a_{i},a_{o}]$.

Previously, we have used $P_{i}$ to define $t_{\rm sec}$ for our disks, but
these results suggest that $P_{\rm median}$ is the correct choice. At
$\alpha=0$, the mean and median periods are equal, and as $\alpha$ increases,
the ratio of the mean to the median period increases. For the power law
distributions explored here, the mean period inherently over-weights the
larger period orbits while the median weights each orbit equally. Because the
orbits in our simulations have equal mass, the median period is the mass-
averaged period.

The variance in $t_{\rm insta}$ at fixed $\alpha$ decreases with increasing
$\alpha$ in Figure 2. This variance is due to sampling noise in the median $P$
from simulation to simulation. The trend in the variance with $\alpha$ is due
to the difference in the slope of the cumulative distribution function (CDF)
of $P$ at the median $P$ for disks with different $\alpha$. We define $F(P)$
as the CDF of the period, $\delta P_{\rm median}=P_{\rm median}-P_{S,{\rm
median}}$ as the difference between the median period of the underlying
distribution and the sample distribution, and $\delta F(P_{\rm
median})=F(P_{\rm median})-F_{S}(P_{\rm median})$ as the difference between
the underlying CDF and sample CDF at the median of the underlying
distribution. We can then approximate,

$\delta P_{\rm median}\approx\frac{dP}{dF(P_{\rm median})}\,\delta F(P_{\rm
median}).$ (7)

$\delta F(P_{\rm median})$ is primarily a function of $N$ which is the same
for all our models. However, $\frac{dP}{dF(P_{\rm median})}$ (the inverse of
the slope of the CDF) decreases with increasing $\alpha$, explaining why the
variance in the median $P$ and therefore $t_{\rm insta}$ decreases with
increasing $\alpha$.

##  4 Discussion

Differential apsidal precession induced by the giant planets in the trans-
Neptunian region suppresses the instability as the coherence time over which
orbits can torque each other is reduced. In Zderic & Madigan (2020), we found
that the inclination instability is suppressed when the timescale for the
instability and the differential precession timescale are approximately equal.
We define the timescale ratio as,

$r=\frac{t_{\rm insta}}{t_{\rm prec}}.$ (8)

We define the differential precession timescale,

$t_{\rm prec}=\frac{1\,{\rm rad}}{\dot{\varpi}_{\rm uqrt}-\dot{\varpi}_{\rm
lqrt}},$ (9)

where $\varpi=\omega+\Omega$ is the longitude of periapsis and
$\dot{\varpi}_{\rm uqrt}$ and $\dot{\varpi}_{\rm lqrt}$ are the upper and
lower quartile precession rates for the disk from our simulations. This
definition emphasizes the importance of differential precession within the
disk in weakening the mutual secular torques between the disk orbits, as
opposed to the absolute value of the precession rate at any point in the disk
(Zderic & Madigan, 2020).

For a scattered disk orbital distribution and the current configuration of
giant planets, we found that the inclination instability was suppressed when
$r\gtrsim 1.2\pm 0.4$. To obtain $t_{\rm prec}$ for a realistic solar system
scattered disk, we scale the measured precession rates in our simulations to
lower disk mass using the secular timescale, $\dot{\varpi}_{\rm disk}\propto
M_{\rm disk}$. As our disk is composed of highly eccentric orbits, we cannot
use a disturbing function approach to get a close-formed expression for the
disk precession rate like those found in Silsbee & Rafikov (2015); Sefilian et
al. (2021), and we must calculate the precession timescale numerically. In
addition, we inject apsidal precession modeling the influences of the giant
planets in the trans-Neptunian region with a quadrupole-level multipole
expansion, i.e. $J_{2}$ precession. For the inclination instability timescale,
we use the $N\rightarrow\infty$ expression for scattered disk found in Zderic
& Madigan (2020) except now including a $P_{\rm median}$ term normalized to
$P_{\rm median}$ when $\alpha=1$,

$t_{\rm insta}=\frac{2.4}{\pi}\frac{M_{\odot}}{M_{\rm disk}}\frac{P_{\rm
median}(\alpha)}{P_{\rm median}(1)}P_{i}.$ (10)

Figure 4: Instability and precession timescale ratio, $r$, as a function of
disk mass, $M_{\rm disk}$, calculated from our simulations. $\alpha$ is
distinguished by color and marker. The critical $r$ value determining the
transition from stability (above) to instability (below) is shown with a grey
band. We find that $\alpha\leq 1$ have similar disk mass requirements needing
more than $20\,M_{\oplus}$ to be unstable, while $\alpha>1$ require less mass
to be unstable.

Using these timescales, we calculate $r$ as a function of $M_{\rm disk}$ for
different $\alpha$ (see Figure 4). We find that the dramatic change in growth
rate for different disk profiles shown in Figure 2 is tempered by the
difference in differential precession in these disks. In fact, we find that
the dependence of $r$ on $M_{\rm disk}$ is nearly identical for $\alpha\leq
1$. At the inner edge of these disks, orbital precession is driven by the
giant planets, $\dot{\varpi}_{J_{2}}\propto a^{-7/2}$. The steeper disk
profiles place more bodies near the inner edge of the disk, and a larger
$t_{\rm prec}$ results. However, we do find that $r$ decreases with $\alpha$
if $\alpha>1$. The differential precession rate in these steeper disks is
similar to the $\alpha=1$ disk while $t_{\rm insta}$ still decreases $\propto
P_{\rm median}$. Therefore, the steepest disk profiles need slightly less mass
to resist the $J_{2}$ precession of the giant planets. The $\alpha=2.5$ disk
profile found by Huang et al. (2022) requires $M_{\rm disk}\gtrsim
10\,M_{\oplus}$ to be unstable, a factor of two lower than the $\alpha=1$
disks used in our prior works. This buckled population will depopulate over
time as the bodies continue to torque each others’ orbits. When torqued to low
orbital angular momenta, disk bodies are vulnerable to scattering and removal
by giant planets. The rate at which the population decreases is left for
future work.

The semi-major axis distribution of the primordial scattered disk effects the
post-instability orbital distributions. In comparison to the $\alpha=1$
distribution explored in our previous work, we find that:

  1. 1.

shallow (small $\alpha$) disks produce large perihelion orbits post-
instability ($q\lesssim 300\,{\rm au}$) in a broad distribution while steeper
disks (large $\alpha$) produce more modestly detached orbits ($q\lesssim
100\,{\rm au}$) in a narrower distribution,

  2. 2.

post-instability, mean perihelion distance increases with increasing semi-
major axis in smaller $\alpha$ disks while the opposite occurs in larger
$\alpha$ disks,

  3. 3.

the post-instability inclinations in the smaller $\alpha$ disks are larger
($\sim 50\pm 40\,{\rm deg}$) compared to larger $\alpha$ disks ($\sim 35\pm
20\,{\rm deg}$),

  4. 4.

all the unstable disks torque orbits to retrograde orientations, but the small
$\alpha$ disks produce them at a higher rate and are capable of producing
retrograde orbits before the disk buckles (preliminary results suggest
$\sim\\!4$ times more retrograde orbits are produced in the $\alpha=0$ disks
compared to the $\alpha=2.5$ disks).

We’ve shown here that the median orbit within the disk determines the
instability’s timescale. The median orbit in a steeper disk also has a lower
initial eccentricity than one in a shallower disk. This means that steeper
disks buckle more weakly and result in lower post-instability inclination and
perihelion distributions. Upcoming observational surveys, such as the Legacy
Survey of Space and Time (LSST) on the Vera Rubin Observatory (Ivezić et al.,
2019), may be able to distinguish between a steep or shallow primordial
scattered disk in the outer solar system using these trends. The discovery of
detached objects with $q\gtrsim 100\,{\rm au}$ at $a\gtrsim 500\,{\rm au}$
would rule out unstable primordial scattered disks with $\alpha\gtrsim 2$, and
the discovery of extremely detached orbits with $q\sim 300\,{\rm au}$ would
require $\alpha\sim 0$. An abundance of highly inclined and retrograde orbits
would indicate a flatter primordial scattered disk semi-major axis
distribution.

## Acknowledgements

We thank the referee for their useful comments. AM gratefully acknowledges
support from the David and Lucile Packard Foundation. This work utilized
resources from the University of Colorado Boulder Research Computing Group,
which is supported by the National Science Foundation (awards ACI-1532235 and
ACI-1532236), the University of Colorado Boulder, and Colorado State
University.

## References

  * Batygin & Brown (2016) Batygin, K., & Brown, M. E. 2016, AJ, 151, 22
  * Becker et al. (2018) Becker, J. C., Khain, T., Hamilton, S. J., et al. 2018, AJ, 156, 81
  * Bernardinelli et al. (2020) Bernardinelli, P. H., Bernstein, G. M., Sako, M., et al. 2020, PSJ, 1, 28
  * Brown & Batygin (2019) Brown, M. E., & Batygin, K. 2019, AJ, 157, 62
  * Duncan & Levison (1997) Duncan, M. J., & Levison, H. F. 1997, Science, 276, 1670
  * Huang et al. (2022) Huang, Y., Gladman, B., Beaudoin, M., & Zhang, K. 2022, ApJ, 938, L23
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Kavelaars et al. (2020) Kavelaars, J. J., Lawler, S. M., Bannister, M. T., & Shankman, C. 2020, Perspectives on the distribution of orbits of distant Trans-Neptunian objects, ed. D. Prialnik, M. A. Barucci, & L. Young (Elsevier), 61–77
  * Lawler et al. (2017) Lawler, S. M., Shankman, C., Kaib, N., et al. 2017, AJ, 153, 33
  * Madigan & McCourt (2016) Madigan, A.-M., & McCourt, M. 2016, MNRAS, 457, L89
  * Madigan et al. (2018) Madigan, A.-M., Zderic, A., McCourt, M., & Fleisig, J. 2018, AJ, 156, 141
  * Napier et al. (2021) Napier, K. J., Gerdes, D. W., Lin, H. W., et al. 2021, The Planetary Science Journal, 2, 59
  * Rein & Liu (2012) Rein, H., & Liu, S.-F. 2012, A&A, 537, A128
  * Sefilian et al. (2021) Sefilian, A. A., Rafikov, R. R., & Wyatt, M. C. 2021, ApJ, 910, 13
  * Shankman et al. (2017) Shankman, C., Kavelaars, J. J., Lawler, S. M., Gladman, B. J., & Bannister, M. T. 2017, AJ, 153, 63
  * Silsbee & Rafikov (2015) Silsbee, K., & Rafikov, R. R. 2015, ApJ, 798, 71
  * Tamayo et al. (2020) Tamayo, D., Rein, H., Shi, P., & Hernand ez, D. M. 2020, MNRAS, 491, 2885
  * Trujillo & Sheppard (2014) Trujillo, C. A., & Sheppard, S. S. 2014, Nature, 507, 471
  * Vokrouhlický & Nesvorný (2019) Vokrouhlický, D., & Nesvorný, D. 2019, Celestial Mechanics and Dynamical Astronomy, 132, 3
  * Zderic et al. (2020) Zderic, A., Collier, A., Tiongco, M., & Madigan, A.-M. 2020, ApJ, 895, L27
  * Zderic & Madigan (2020) Zderic, A., & Madigan, A.-M. 2020, AJ, 160, 50
  * Zderic et al. (2021) Zderic, A., Tiongco, M., Collier, A., et al. 2021, AJ, 162, 278


# Origin of the gamma-ray emission from supernova remnant HB9

Tülün Ergin111Now at Middle East Technical University, Physics Group, Northern
Cyprus Campus, 99738 Kalkanli via Mersin 10, Turkey<EMAIL_ADDRESS>Lab
Saha222Now at Harvard-Smithsonian Center for Astrophysics, Cambridge, MA.
<EMAIL_ADDRESS>Hidetoshi Sano<EMAIL_ADDRESS>Aytap Sezer  Ryo
Yamazaki  Pratik Majumdar  Yasuo Fukui  TUBITAK Space Technologies Research
Institute, ODTU Campus, 06531, Ankara, Turkey  PARCOS and Department of
EMFTEL, Universidad Complutense de Madrid, E-28040 Madrid, Spain  National
Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan  Department
of Electrical-Electronics Engineering, Avrasya University, 61250, Trabzon,
Turkey  Department of Physics and Mathematics, Aoyama Gakuin University,
5-10-1 Fuchinobe, Sagamihara 252-5258, Japan  Saha Institute of Nuclear
Physics, HBNI, 1/AF Bidhannagar, Kolkata 700064, India  Department of Physics,
Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

###### Abstract

HB9 (G160.9+2.6) is a mixed-morphology Galactic supernova remnant (SNR) at a
distance of $\sim$0.6 kpc. Previous analyses revealed recombining plasma
emission in X-rays and an expanding shell structure in H i and CO emission,
which were correlating with the spatial extent of HB9. In GeV energies, HB9
was found to show extended gamma-ray emission with a morphology that is
consistent with the radio continuum emission showing a log-parabola-type
spectrum. The overlap reported between the gas data and the excess gamma-ray
emission at the southern region of the SNR’s shell could indicate a possible
interaction between them. We searched for hadronic gamma-ray emission
signature in the spectrum to uncover possible interaction between the
molecular environment and the SNR. Here we report the results of the gamma-ray
spectral modelling studies of HB9.

###### keywords:

ISM: supernova remnants $-$ HB9 (G160.9+2.6) $-$ gamma-rays: ISM

††journal: Advances in Space Research

##  1 Introduction

Detection of gamma-ray emission from supernova remnants (SNRs) and their
neighbourhood is key for understanding the origin, acceleration, and
propagation of cosmic rays (CRs) with very high energies up to 1015 eV.

CRs (i.e. protons) that escape the shock front of SNRs enter the surrounding
ambient interstellar medium (ISM). Here they interact with protons in the
ambient gas or molecular clouds (MCs) through the hadronic interactions, which
create neutral pions that later decay into gamma rays (Aharonian & Atoyan,
1996; Gabici et al., 2009; Fujita et al., 2009). This hadronic interaction
model is called the illuminated cloud model or runaway CR model (Ptuskin &
Zirakashvili, 2003, 2005; Ohira et al., 2010; Li & Chen, 2010; Ellison &
Bykov, 2011; Malkov et al., 2013; Celli et al., 2019).

Therefore, gamma rays from nearby MCs are crucial to probe the highest energy
protons escaping SNRs. The second well-known hadronic model is the crushed
cloud model or direct interaction model, where the SNR’s shock wave interacts
with MCs (Bykov et al., 2000; Uchiyama et al., 2010; Ohira et al., 2011; Tang
& Chevalier, 2015; Lee et al., 2015; Cardillo et al., 2016). This latter model
is usually observed in middle-aged SNRs.

In the leptonic models (Aharonian & Atoyan, 1999; Porter et al., 2006),
accelerated electrons produce gamma rays through the inverse-Compton (IC) and
Bremsstrahlung processes. This type of gamma-ray emission is usually
accompanied by synchrotron emission observed in radio wavebands. IC emission
is produced by the interaction of the accelerated electrons with the ambient
photons, like cosmic microwave background (CMB) emission. Highly energetic
electrons may be accelerated at the shock fronts of SNRs or they may originate
from pulsars.

HB9 (G160.9+2.6) is a middle-aged mixed-morphology (MM) supernova remnant
(SNR) characterised by centre-filled thermal X-rays and a radio shell (Rho &
Petre, 1998).

HB9 has a large angular size and shell-like morphology with a non-thermal
radio spectral index of $\sim$0.64 (e.g. Leahy & Roger 1991; Leahy &
Aschenbach 1995; Leahy & Tian 2007).

The kinematic distance to HB9 was calculated to be 0.6 $\pm$ 0.3 kpc, which
was derived from the expanding shell structure around the remnant (Sezer et
al., 2019) using H i and 12CO($J$ = 1$-$0) data in the velocity range of
$-$10.5 km s-1 and $+$1.8 km s-1. More recently, Zhao et al. (2020) estimated
the extinction distances to 33 SNRs. Using Gaia stars on the line of sight to
HB9 with infrared/optical photometry, they obtained extinction of each star as
a function of the distance determined by Gaia. Their estimated distance is
0.54 $\pm$ 0.10 kpc, which is fully consistent with results of Sezer et al.
(2019) (0.6 $\pm$ 0.3 kpc) that was calculated using 12CO($J$ = 1$-$0) data.
This result clearly shows the association between HB9 and interstellar gas at
the distance of $\sim$0.6 kpc.

It was also reported that HB9 contained recombining (over-ionised) plasma (RP)
in Suzaku spectra. More recently, Saito et al. (2020) re-analysed Suzaku data
to understand the nature of the hard X-ray component reported by previous
X-ray study (Yamauchi & Koyama, 1993). They found no significant hard X-ray
emission from the central region and no significant feature of an over-ionised
plasma in HB9. Their spectra extracted from the centre of HB9 are well-fitted
with a model consisting of a collisional ionisation equilibrium (CIE) and
ionising plasma (IP).

Fig. 1: Gamma-ray TS maps shown in the energy range of 1$-$300 GeV. Left
Panel: HB9 was not added into the gamma-ray background model while producing
the TS map. Black contours show the TS levels of 30, 59, 89, 119. The white
cross and the dashed circle show the radio centre and extension of HB9; Middle
Panel: The 4850 MHz radio continuum map of HB9 is used as a template in the
gamma-ray background model (shown with green contours); Right Panel: In
addition to the 4850 MHz radio continuum template of HB9, a point-source (PS
J0506.5+4546) was added to the gamma-ray background model. The position of PS
J0506.5+4546 is shown with a yellow cross. For further information see Sezer
et al. (2019). There is an excess of diffuse and extended gamma-ray emission
observable towards the south of the radio-shell of HB9. Fig. 2: Intensity
distributions of H i (left panel) and 12CO($J$ = 1–0) (right panel) in the
gamma-ray analysis region are shown in Equatorial coordinates. The integration
velocity of H i and CO is from $-10.5$ km s-1 to $+1.8$ km s-1. The superposed
white contours correspond to TS values of 25, 49, 64, 81, 100 in the energy
range of 0.2$-$300 GeV. The white cross and the dashed circle show the radio
centre and extension of HB9.

In the GeV gamma-ray energy range, many middle-aged MM SNRs were detected by
the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope
(Fermi) (e.g. Abdo et al. 2009, 2010a, 2010b; Castro & Slane 2010; Ergin et
al. 2017). HB9 has recently been found to be one of these MM SNRs (Araya 2014,
Sezer et al. 2019). Using 5.5 years of Fermi-LAT data, Araya (2014) detected
extended gamma-ray emission centred at the position of HB9 with a significance
of 16$\sigma$ above 0.2 GeV. 10 years of Fermi-LAT data was analysed and
investigated the gamma-ray emission morphology and spectrum within the energy
range of 0.2$-$300 GeV. HB9 was detected as an extended gamma-ray source with
a significance of $\sim$25$\sigma$. The best-fitted spatial model to HB9 in
the 1$-$300 GeV energy range was the 4850 MHz radio continuum map (taken by
the Green Bank Telescope), that has an approximate size of 2∘. In addition, a
new point-source (PS J0506.5+4546) was marginally detected ($\sim$6$\sigma$)
right outside the extended region of HB9.

The spectrum of HB9 was best-fitted to a log-parabola (LP) function with
indices of $\alpha$ = 2.36 $\pm$ 0.05 and $\beta$ = 0.14 $\pm$ 0.05. The total
energy flux was found to be (1.51 $\pm$ 0.08) $\times$ 10-5 MeV cm-2 s-1 in
the energy range of 0.2$-$300 GeV. These results were consistent with previous
studies (Araya, 2014). In the TeV energy range MAGIC observations revealed no
detection with an integral flux upper limit of 1.6 $\times$ 10-11 cm-2 s-1 for
energies above 270 GeV (Carmona, 2009).

This work contains a short summary of our previous paper, Sezer et al. (2019),
and an addition of spectral modelling to it. Here, we aim to investigate
possible gamma-ray emission mechanisms (leptonic/hadronic) by modelling the
gamma-ray spectrum obtained from the whole SNR HB9, which was reported in the
energy range of 0.2$-$300 GeV. The paper is organised as follows: In Section
2, we summarise gamma-ray data analysis results and we estimate the mean ISM
proton density for the remnant (Section 2.1). We then report our preliminary
results about the SED modelling of HB9 in Section 2.2, and we close this work
by conclusions in Section 3.

##  2 Analysis and Results

Figure 1 shows the distribution of the GeV gamma-rays based on the TS maps
produced in the energy range of 1 - 300 GeV, where the left panel shows the TS
map without including the radio continuum template of HB9 and the newly
discovered gamma-ray point source (PS J0506.5+4546) in the gamma-ray
background model. The right panel shows the TS morphology where both PS
J0506.5+4546 and the radio continuum template are included in the background
model. Two sources from the Fourth Fermi-LAT sources (4FGL) catalog (The
Fermi-LAT Collaboration, 2020) are shown by white markers.

PS J0506.5$+$4546 was detected at a significance level of about 6 sigma by
applying the best-fit spatial model (i.e. Radio Template $+$ 1 P. S.), but we
assumed PS J0506.5$+$4546 to be physically not related to HB9, based on the
assumption that it lies outside the shell of the SNR.

However, there are currently no known counterparts corresponding to this
source, apart from the positional overlap with the gas/MCs data, which we have
shown. So, it is not clear whether PS J0506.5+4546 is another source
independent of HB9 or it is part of the SNR’s shell interacting with the
ambient gas/MCs.

There is another region of diffuse gamma-ray emission towards even more south
of the SNR, but neither in our previous paper, nor in this study we analysed
this “rather large/extended region, because we assumed that this region might
not be directly related to HB9, mainly because it lies outside the radio
emission region of the SNR. In addition, this large and diffuse gamma-ray
emission region might be an extension of 4FGL 0503.6+4518, which we also
assumed to be not related to the SNR (but included as a point-source into our
gamma-ray background model).

###  2.1 Estimation of the Mean ISM Proton Density for SNR HB9

In Figure 2, the distribution of the H i (left panel) and 12CO($J$ = 1$-$0)
(right panel) data is shown in the analysis region of HB9, where the 408 MHz
radio continuum intensity contours are overlaid. We calculated a mean ISM
proton density of $\sim$60 cm-3 by assuming following parameters as given by
Sano et al. (2019a):

  * 1.

Distance: 0.6 kpc,

  * 2.

Shell radius: 11 pc ($\sim$1.05 deg),

  * 3.

Shell thickness: 3.5 pc ($\sim$0.33 deg),

  * 4.

CO-to-H2 conversion factor: 2.0 $\times$ 1020 cm-2 (K km s-2)-1,

  * 5.

H i: Optically thin assumption.

The total column density is taken as the summation of twice the molecular
hydrogen column density and the proton column density for atomic hydrogen.

###  2.2 SED Modelling of SNR HB9 and Preliminary Results

We modelled the SED of HB9, which was reported by (Sezer et al., 2019). We
have considered both leptonic, i.e. bremsstrahlung and inverse Compton (IC)
scattering, and hadronic models to explain the observed data (Blumenthal &
Gould, 1970; Kelner et al., 2006). For the relativistic cosmic-ray electrons
and protons we assumed an energy spectrum of the form

$\frac{dN}{dE}=AE^{-\alpha^{\prime}}e^{(-E/E_{\rm cut})},$

where $\alpha^{\prime}$ is the spectral index, A is the amplitude, and $E_{\rm
cut}$ is the cut-off energy of the particle spectrum.

The expression for the total energy (Etotal) of particles is

$E_{total}=\int_{E_{1}}^{E_{2}}E{dN\over dE}dE,$

where $E_{1}$ and $E_{2}$ are limits of the energy of particles.

We used CMB and interstellar dust emission for the calculation of IC fluxes
(Mathis et al., 1983).

Table 1: The estimated spectral parameters for the electron and proton spectra. The total energy (Etotal) calculated for particles for the hadronic ($\pi^{0}$-decay) and leptonic (Bremsstrahlung and IC scattering) spectral models fitted to the gamma-ray SED of HB9 assumes a distance of $\sim$0.6 kpc and nH = 60 cm-3. |  Bremsstrahlung |  $\pi^{0}$-decay |  IC-scattering
---|---|---|---
$\alpha^{\prime}$ |  2.21 |  2.69 |  1.74
Ecut [GeV] |  10 |  1000 |  200
Etotal [erg] |  1.45$\times$1047 |  2.25$\times$1047 |  1.90$\times$1048

For bremsstrahlung and pion-decay processes, we used the mean ISM proton
density (nH) to be 60 cm-3 and the proton-proton cross section from Kelner et
al. (2006) was used for the pion-decay modelling. The preliminary fitted
spectral models are shown in Figure 3. The total energy of particles for each
of the different fitted spectral models is given in Table LABEL:table1.

Fig. 3: Gamma-ray SED and the preliminary leptonic models (IC-scattering on
the upper plot and Bremsstrahlung on the middle plot) and the hadronic (on the
bottom plot) model fitted by assuming nH = 60 cm-3.

##  3 Discussion and Conclusion

We modelled the GeV gamma-ray emission spectrum of HB9, which was obtained
from our previous analysis reported in S19. Using H i and 12CO($J$ = 1$-$0)
data, we estimated the ISM proton density of $\sim$60 cm-3, which is an
average value calculated by assuming the shell radius and thickness to be 11
pc and 3.5 kpc, respectively.

Although the classical signs of shock-cloud interaction (e.g., CO line-
broadening, the high kinetic temperature of shocked clouds) are not present
for HB9, an expanding CO/HI shell related to HB9 was found (shown in Figure 4
Panel (d) of (Sezer et al., 2019)). In general, the expanding shell is thought
to be formed by strong stellar winds from the progenitor of the SNR and the
SNR shock is now likely interacting with the wind-wall (e.g. Sano et al. 2017,
2019b), which are walls of the stellar wind bubble formed in the ISM by swept
up cold and warm gas, in neutral and ionised state (R. Boomsma, 2007).
Therefore, accelerated CRs by the shock fronts of HB9 are mainly interacting
with the inner wall of the stellar wind bubble. The integration velocity
range, therefore, can be determined quite uniquely. We confirmed that the
density contribution is $\sim$70% for atomic hydrogen and $\sim$30% for
molecular hydrogen.

By looking at the gas distribution overlapping HB9, a dominating hadronic
gamma-ray emission was expected as a result of our SED modelling. However, the
observed gamma-ray spectrum of HB9 can be explained by both leptonic and
hadronic models. Within the context of the leptonic model, we found that both
IC and Bremsstrahlung can explain the observed SED individually for two
different sets of model parameters.

The reason for this result is mainly due to the assumption that ISM to be
uniformly distributed within the shell, with a density of $\sim$60 cm-3, which
is not the case. The gaseous interstellar material overlaps only at certain
regions of the SNR’s shell, which has to be taken into account while doing
spectral analysis and modelling the SED. So, if we consider both IC and
Bremsstrahlung processes together to explain the observed SED for a set of
model parameters, we need a robust estimate of the ambient matter density
which is currently lacking.

Although the overall fitting of the models to the gamma-ray SED, does not show
preference of any model over the other (i.e. Leptonic/Hadronic), we note that
the SED of HB9 might give some clues about the dominating gamma-ray emission
model. Namely, when we compare our spectral modelling results for HB9 with
results from Ackermann et al. (2013), we see that interacting SNRs, such as IC
443 and W44, present a spectrum that fits better to a hadronic gamma-ray
emission model. In particular, at gamma-ray energies lower than 200 MeV, there
is a clear separation between Bremsstrahlung models and the Pi-zero Decay
model, where the flux levels of the latter one decrease towards lower
energies.

For the gamma-ray SED of HB9 at energies below 200 MeV, although results of
spectral modelling hint a dominating hadronic emission which is backed up by
CO and HI intensity maps in Figure 2, more data and detailed analyses are
required to conclude whether the hadronic model is really dominating for HB9.

##  4 Acknowledgments

We would like to thank our referees at the Advances in Space Research journal
for their valuable comments.
Facilities: Fermi-LAT, Dominion Radio Astrophysical Observatory, Harvard-
Smithsonian Center for Astrophysics 1.2 m MMW-radio Telescope, Green Bank
Telescope.

## References

  * Abdo et al. (2009) Abdo, A. A., Ackermann, M., Ajello, M. et al., 2009. Fermi LAT Discovery of Extended Gamma-Ray Emission in the Direction of Supernova Remnant W51C. ApJL 706, L1-L6. https://doi.org/10.1088/0004-637X/706/1/L1.
  * Abdo et al. (2010a) Abdo, A. A. et al., 2010a. Observation of Supernova Remnant IC 443 with the Fermi Large Area Telescope. ApJ 712, 459-468. https://doi.org/10.1088/0004-637X/712/1/459.
  * Abdo et al. (2010b) Abdo, A. A. et al., 2010b. Fermi Large Area Telescope Observations of the Supernova Remnant W28 (G6.4-0.1). ApJ 718, 348-356. https://doi.org/10.1088/0004-637X/718/1/348.
  * Ackermann et al. (2013) Ackermann, M. et al. (The Fermi-LAT collaboration), 2013. Detection of the Characteristic Pion-Decay Signature in Supernova Remnants. Science 339, 807-811. https://doi.org/10.1126/science.1231160.
  * Aharonian & Atoyan (1996) Aharonian, F. A., Atoyan, A. M., 1996. On the emissivity of $\pi^{0}$-decay gamma radiation in the vicinity of accelerators of galactic cosmic rays. A&A 309, 917-928.
  * Aharonian & Atoyan (1999) Aharonian, F. A., Atoyan, A. M., 1999. On the origin of TeV radiation of SN 1006. A&A 351, 330-340.
  * Araya (2014) Araya, M., 2014. Fermi LAT observation of supernova remnant HB9. MNRAS 444, 860-865. https://doi.org/10.1093/mnras/stu1484.
  * Blumenthal & Gould (1970) Blumenthal, G.R., Gould, R.J., 1970. Bremsstrahlung, synchrotron radiation, and Compton scattering of high-energy electrons traversing dilute gases. Rev. Mod. Phys. 42, 237–271. http://dx.doi.org/10.1103/RevModPhys.42.237
  * R. Boomsma (2007) Boomsma, R. 2007. The disk-halo connection in NGC 6946 and NGC 253. Ph.D. Thesis, University of Groningen, s.n., ISBNs: 9789036729468.
  * Bykov et al. (2000) Bykov, A. M., Chevalier, R. A., Ellison, D. C., & Uvarov Y. A., 2000. Nonthermal Emission from a Supernova Remnant in a Molecular Cloud. ApJ 538, 203-216. https://doi.org/10.1086/309103.
  * Cardillo et al. (2016) Cardillo, M., Amato, E., & Blasi P., 2016. Supernova remnant W44: a case of cosmic-ray reacceleration. A&A 595, A58. https://doi.org/10.1051/0004-6361/201628669.
  * Carmona (2009) Carmona, E., Costado, M. T., Font, L., Zapatero, J., 2009. Observation of selected SNRs with the MAGIC Cherenkov Telescope. arXiv:0907.1009.
  * Castro & Slane (2010) Castro, D., Slane, P., 2010. Fermi Large Area Telescope Observations of Supernova Remnants Interacting with Molecular Clouds. ApJ 717, 372-378. https://doi.org/10.1088/0004-637X/717/1/372.
  * Celli et al. (2019) Celli, S., Morlino, G., Gabici, S., & Aharonian, F. A., 2019. Exploring particle escape in supernova remnants through gamma rays. MNRAS 490, 4317-4333. https://doi.org/10.1093/mnras/stz2897.
  * Ellison & Bykov (2011) Ellison, D. C. & Bykov, A. M., 2011. Gamma-ray Emission of Accelerated Particles Escaping a Supernova Remnant in a Molecular Cloud. ApJ 731, 87. https://doi.org/10.1088/0004-637X/731/2/87.
  * Ergin et al. (2017) Ergin, T., Sezer, A., Sano, H., Yamazaki, R., Fukui, Y., 2017. Recombining Plasma and Gamma-Ray Emission in the Mixed-morphology Supernova Remnant 3C 400.2. ApJ 842, 22. https://doi.org/10.3847/1538-4357/aa72ee.
  * Fujita et al. (2009) Fujita, Y., Ohira, Y., Tanaka, S. J., Takahara, F., 2009. Molecular Clouds as a Probe of Cosmic-Ray Acceleration in a Supernova Remnant. ApJ 707, L179-L183. https://doi.org/10.1088/0004-637X/707/2/L179.
  * Gabici et al. (2009) Gabici, S., Aharonian, F. A., Casanova, S., 2009. Broad-band non-thermal emission from molecular clouds illuminated by cosmic rays from nearby supernova remnants. MNRAS 396, 1629-1639. https://doi.org/10.1111/j.1365-2966.2009.14832.x.
  * Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., & Bugayov, V. V., 2006. Energy spectra of gamma rays, electrons, and neutrinos produced at proton-proton interactions in the very high energy regime. Phys. Rev. D 74, 034018. https://doi.org/10.1103/PhysRevD.74.034018.
  * Leahy & Roger (1991) Leahy, D. A., Roger, R. S., 1991. Radio Emission from the Supernova Remnant G160.9+2.6 (HB9). AJ 101, 1033. https://doi.org/10.1086/115745.
  * Leahy & Aschenbach (1995) Leahy, D. A., Aschenbach, B., 1995. ROSAT X-ray observations of the supernova remnant HB 9. A&A 293, 853-858.
  * Leahy & Tian (2007) Leahy, D. A., Tian, W. W., 2007. Radio spectrum and distance of the SNR HB9. A&A 461, 1013-1018. https://doi.org/10.1051/0004-6361:20065895.
  * Lee et al. (2015) Lee, S.-H., Patnaude, D. J., Raymond, J. C., Nagataki, S., Slane, P. O., & Ellison, D. C., 2015. Modeling Bright $\gamma$-Ray and Radio Emission at Fast Cloud Shocks. ApJ 806, 71. https://doi.org/10.1088/0004-637X/806/1/71.
  * Li & Chen (2010) Li, H., Chen, Y., 2010. $\gamma$-rays from molecular clouds illuminated by accumulated diffusive protons from supernova remnant W28. MNRAS 409, L35-L38. https://doi.org/10.1111/j.1745-3933.2010.00944.x.
  * Malkov et al. (2013) Malkov, M. A., Diamond, P. H., Sagdeev, R. Z., Aharonian, F. A. & Moskalenko, I. V., 2013. Analytic Solution for Self-regulated Collective Escape of Cosmic Rays from Their Acceleration Sites. ApJ 768, 73. https://doi.org/10.1088/0004-637X/768/1/73.
  * Mathis et al. (1983) Mathis, J. S., Mezger, P. G., & Panagia, N., 1983. Interstellar radiation field and dust temperatures in the diffuse interstellar matter and in giant molecular clouds. A&A 128, 212.
  * Ohira et al. (2010) Ohira, Y., Murase, K., Yamazaki, R., 2010. Escape-limited model of cosmic-ray acceleration revisited. A&A 513, A17. https://doi.org/10.1051/0004-6361/200913495.
  * Ohira et al. (2011) Ohira, Y., Murase, K., Yamazaki, R., 2011. Gamma-rays from molecular clouds illuminated by cosmic rays escaping from interacting supernova remnants. MNRAS 410, 1577-1582. https://doi.org/10.1111/j.1365-2966.2010.17539.x.
  * Porter et al. (2006) Porter, T. A., Mosalenko, I. V., & Strong, A. W., 2006. Inverse Compton Emission from Galactic Supernova Remnants: Effect of the Interstellar Radiation Field. ApJ 648, L29-L32. https://doi.org/10.1086/507770.
  * Ptuskin & Zirakashvili (2003) Ptuskin, V. S.& Zirakashvili, V. N., 2003. Limits on diffusive shock acceleration in supernova remnants in the presence of cosmic-ray streaming instability and wave dissipation. A&A 403, 1-10. https://doi.org/10.1051/0004-6361:20030323.
  * Ptuskin & Zirakashvili (2005) Ptuskin, V. S. & Zirakashvili, V. N., 2005. On the spectrum of high-energy cosmic rays produced by supernova remnants in the presence of strong cosmic-ray streaming instability and wave dissipation. A&A 429, 755-765. https://doi.org/10.1051/0004-6361:20041517.
  * Rho & Petre (1998) Rho, J., Petre, R., 1998. Mixed-Morphology Supernova Remnants. ApJ 503, L167-L170. https://doi.org/10.1086/311538.
  * Saito et al. (2020) Saito, M., Yamauchi, S., Nobukawa, K.K., Bamba, A., Pannuti, T.G., 2020. X-ray emission from the mixed-morphology supernova remnant HB 9. PASJ 72, 65. https://doi.org/10.1093/pasj/psaa042.
  * Sano et al. (2017) Sano, H. et al. 2017. Interstellar gas and X-rays toward the Young supernova remnant RCW 86; pursuit of the origin of the thermal and non-thermal X-ray. JHEAp 15, 1-18. https://doi.org/10.1016/j.jheap.2017.04.002.
  * Sano et al. (2019a) Sano, H. et al. 2019a. Possible Evidence for Cosmic-Ray Acceleration in the Type Ia SNR RCW 86: Spatial Correlation between TeV Gamma-Rays and Interstellar Atomic Protons. ApJ 876, 37. https://doi.org/10.3847/1538-4357/ab108f.
  * Sano et al. (2019b) Sano, H. et al. 2019b. Discovery of Shocked Molecular Clouds Associated with the Shell-type Supernova Remnant RX J0046.5-7308 in the Small Magellanic Cloud. ApJ 881, 85. https://doi.org/10.3847/1538-4357/ab2ade.
  * Sezer et al. (2019) Sezer, A., Ergin, T., Yamazaki, R., Sano, H., Fukui, Y., 2019. Discovery of recombining plasma inside the extended gamma-ray supernova remnant HB9. MNRAS 489, 4300-4310. https://doi.org/10.1093/mnras/stz2461.
  * Tang & Chevalier (2015) Tang, X. & Chevalier, R. A., 2015. Time-dependent Diffusive Shock Acceleration in Slow Supernova Remnant Shocks. ApJ 800, 103. https://doi.org/10.1088/0004-637X/800/2/103.
  * Tashiro et al. (2018) Tashiro, M. et al., 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10699, Proc. SPIE, 1069922.
  * The Fermi-LAT Collaboration (2020) The Fermi-LAT Collaboration, 2020. Fermi Large Area Telescope Fourth Source Catalog. ApJS 247, 33. https://doi.org/10.3847/1538-4365/ab6bcb.
  * Tuohy et al. (1979) Tuohy, I. R., Clark, D. H., Garmire, G. P., 1979. Observation of soft X-ray emission from the SNR HB 9. MNRAS 189, 59-63. https://doi.org/10.1093/mnras/189.1.59P.
  * Yamauchi & Koyama (1993) Yamauchi, S., Koyama, K., 1993. Hard X-Ray Emission from the Region near the Supernova Remnant HB 9 and the Radio Galaxy 4C 46.09. PASJ 45, 545-550.
  * Uchiyama et al. (2010) Uchiyama, Y., Blandford, R. D., Funk, S., Tajima, H., & Tanaka, T., 2010. Gamma-ray Emission from Crushed Clouds in Supernova Remnants. ApJ 723, L122-L126. https://doi.org/10.1088/2041-8205/723/1/L122.
  * Zhao et al. (2020) Zhao, H., Jiang, B., Li J., Chen, B., Yu, B., Wang, Y. 2020. A Systematic Study of the Dust of Galactic Supernova Remnants. I. The Distance and the Extinction. ApJ 891, 137. https://doi.org/10.3847/1538-4357/ab75ef.