# Continuous Prediction of Lower-Limb Kinematics From Multi-Modal Biomedical
Signals

Chunzhi Yi, Feng Jiang∗, Shengping Zhang, Hao Guo, Chifu Yang, Zhen Ding,
Baichun Wei, Xiangyuan Lan∗, and Huiyu Zhou Chunzhi Yi, Hao Guo, Zhen Ding and
Chifu Yang are with the School of Mechatronics Engineering, Harbin Institute
of Technology, Harbin, Heilongjiang, 150001 China e-mail:
chunzhiyi@hit.edu.cn(C.Yi).Chunzhi Yi, Hao Guo, Baichun Wei and Feng Jiang are
with the School of Computer, Harbin Institute of Technology, Harbin,
Heilongjiang, 150001 China and Pengcheng Laboratory, Shenzhen, Guangdong,
China e-mail: fjiang@hit.edu.cn(F.Jiang).Shengping Zhang is with the School of
Computer Science and Technology, Harbin Institute of Technology, Weihai
264209, China.Xiangyuan Lan is with the Department of Computer Science, Hong
Kong Baptist University, Hong Kong email<EMAIL_ADDRESS>(Corresponding author: Feng Jiang, Xiangyuan Lan). Huiyu Zhou is with the
School of Informatics, University of Leicester, Leicester LE1 7RH, U.K.

###### Abstract

The fast-growing techniques of measuring and fusing multi-modal biomedical
signals enable advanced motor intent decoding schemes of lower-limb
exoskeletons, meeting the increasing demand for rehabilitative or assistive
applications of take-home healthcare. Challenges of exoskeletons’ motor intent
decoding schemes remain in making a continuous prediction to compensate for
the hysteretic response caused by mechanical transmission. In this paper, we
solve this problem by proposing an ahead-of-time continuous prediction of
lower-limb kinematics, with the prediction of knee angles during level walking
as a case study. Firstly, an end-to-end kinematics prediction
network(KinPreNet)111The demo is avalable on https://youtu.be/CMwtrVaeLZ4, and
the source code is attached to the IEEE Code Ocean., consisting of a feature
extractor and an angle predictor, is proposed and experimentally compared with
features and methods traditionally used in ahead-of-time prediction of gait
phases. Secondly, inspired by the electromechanical delay(EMD), we further
explore our algorithm’s capability of compensating response delay of
mechanical transmission by validating the performance of the different
sections of prediction time. And we experimentally reveal the time boundary of
compensating the hysteretic response. Thirdly, a comparison of employing EMG
signals or not is performed to reveal the EMG and kinematic signals’
collaborated contributions to the continuous prediction. During the
experiments, EMG signals of nine muscles and knee angles calculated from
inertial measurement unit (IMU) signals are recorded from ten healthy
subjects. Our algorithm can predict knee angles with the averaged RMSE of 3.98
deg which is better than the 15.95-deg averaged RMSE of utilizing the
traditional methods of ahead-of-time prediction. The best prediction time is
in the interval of 27ms and 108ms. To the best of our knowledge, this is the
first study of continuously predicting lower-limb kinematics in an ahead-of-
time manner based on the electromechanical delay (EMD).

###### Index Terms:

Electromechanical Delay, Reponse Delay Compensation, Kinematics Prediction
Electromyography, Long Short-Term Memory.

##  I Introduction

THE increasing number of gait-impaired patients, led by recent advancement in
medical treatment and expanded life expectancy [1, 2], casts a rising demand
on rehabilitative and assistive devices for take-home healthcare [3, 4, 5].
Lower-limb exoskeletons, functioned as a home-healthcare device for both
rehabilitation and assistive applications, have been enabled by the fast-
growing techniques of measuring and fusing multi-modal biomedical signals [4,
6]. Fusing multi-modal biomedical signals sought to decode human motor intent,
i.e. to perceive subject-specific gait characteristics, thus can result in
improved assistive and/or rehabilitative performance of assistance. Decoding
human motor intent, which highly relies on kinematics-related information, is
still a challenging topic of exoskeletons especially for assisting subjects
still with mobility. Traditionally, pressure insoles allow the robot to
automatically detect users’ gait phases according to the different pressure
distribution under feet during one gait cycle. However, such movement intent
decoding might result in a severe intent misjudgment and thus a potential
injury if stepping on a protuberance, while cumbersome customized shoes have
to be embedded into the robot [7]. Alternatively, some movement intent
decoding schemes utilized kinematics-related characteristics to detect key
timings of gait [8, 9, 10], demonstrating the superiority of employing
kinematics information. However, on one hand, the key timings provide discrete
information of gait and movement intent, thus might limit the further
improvement of assistive performances. On the other hand, as revealed by the
simulation study [11], the response delay caused by the transmission of the
mechanical structure of exoskeletons will impede the close coordination
between human and exoskeleton and thus greatly degrade the assistive
performance. How to compensate the transmission delay is ignored by currently
existing control and movement intent decoding schemes of lower-limb
exoskeletons.

The acquisition of continuous kinematics could be enabled by techniques like
adaptive oscillator (AO) or Bayesian filter. AO [12, 6] learns periodic
characteristics of locomotion by extracting frequency components from
kinematics and then using the components to reconstruct the reference curves.
According to such reference curves, AO estimates the gait percentage of
stride, the information of which is limited to reflect kinematics in detail.
In [13], iterative extended Kalman filter was employed to make a one-step-
ahead prediction of richer gait-related information, including gait events and
trajectories of lower-limb joints. However, the methods in [12, 6, 13], based
on learning from the previous reference curves, rendered strong dependence on
the periodic characteristics of gait. Additionally, the one-step-ahead
prediction might be too short to fully compensate for the response delay of
the mechanical transmission.

Neural signals, which are generated prior to the corresponding movements,
provide a promising solution to the compensation of response delay.
Myoelectric signals, one of the main sources for obtaining neural activation
have been widely used in the control of exoskeletons [14], powered prostheses
[15, 16, 17] and rehabilitation robots [18, 19, 20, 21, 22, 23, 24]. In
decoding EMG signals, a pioneer work of H. Huang [25], which investigated the
EMG-based pattern recognition (PR) methods to identify subjects’ locomotion
modes, demonstrated the feasibility of using EMG signals to make an ahead-of-
time prediction. Following works like [26] and [27] applied such EMG-based
pattern recognition methods on different locomotion modes, which provided
additional demonstrations. However, given that recent assistive strategies of
lower-limb exoskeletons are developed to require continuous information of
lower-limb kinematics, such PR-based methods just focus on the qualitative and
discrete presentation of lower limb motions, making it difficult to provide
further information.

One solution to continuous EMG-to-kinematics mapping is the musculoskeletal
model combined with dynamic models. A study on upper limbs [28], which
included musculoskeletal and dynamic models into Kalman filter, demonstrated
the feasibility of calculating kinematics based on EMG signals. However, if we
apply this method on lower limbs, the ground-feet wrench, which is intractable
to estimate using wearable sensors, is unavoidably employed in dynamic models.
This issue impedes musculoskeletal and dynamic models’ combined usage on lower
limbs, thus results in either time-costing offline parameter identification of
musculoskeletal models [29] or a real-time framework enabled by force plates
[30]. Alternatively, Brantley et al. [31] directly mapped EMG into knee angles
using an unscented Kalman filter, which was limited by the drawbacks of the
Kalman filter itself. All the above-mentioned EMG-to-movement mapping methods,
albeit capable of continuously acquiring kinematics, did not demonstrate the
capability of making ahead-of-time predictions.

The goal of this study is to make a continuous prediction of kinematics ahead
of time. Particularly, a case study is performed on knee angle prediction. We
propose to explore and exploit the EMG’s characteristics of its onset before
the onset of its corresponding movements (denoted as electromechanical delay
(EMD)) by the means of an end-to-end LSTM network, in order to make the
prediction. The architecture of the algorithm is made comparisons with
traditionally used EMG features and prediction algorithms. Furthermore, based
on experimental validations, we explore the time boundary of how much
transmission delay can be compensated by the continuous kinematics prediction
and reveal how EMG and kinematics signals can contribute to the final
prediction performance. The main contributions of this paper are summarized as
follow:

  * •

To the best of our knowledge, this is the first study of continuously
predicting lower-limb kinematics in an ahead-of-time manner based on the
electromechanical delay (EMD).

  * •

We experimentally revealed EMG and kinematics signals’ collaborated
contribution to the continuous kinematics prediction.

  * •

We further explore the EMD-inspired prediction time through different trials
of prediction time and experimentally revealed our algorithm’s time boundary
of compensating mechanical response delay.

This paper has been organized as follows. Related works are presented in
Section II. Section III details the methodology and experiments. Experimental
results have been explained in Section IV. Discussion and conclusion of the
entire research work have been given in Sections V and VI, respectively.

##  II Related Work

Challenges of constructing the algorithm include solving the continuous
prediction, determining prediction time and extracting effective features of
EMG signals.

###  II-A Continuous and Ahead-Of-Time Prediction

The continuous and ahead-of-time prediction can be inspired by previous works
from two aspects. Firstly, the ahead-of-time prediction of discrete gait
phases, was demonstrated with feasibility by the pioneer works of Huang et.
al. [25]. In Huang’s works, support vector machine (SVM) classifier was
leveraged to continuously classify locomotion modes and predict the
transmission of locomotion modes ahead of time, which presented better
performances than linear discriminant analysis (LDA). Secondly, methods of
upper-limb movement regression demonstrated the robustness and accuracy of
Recurrent Neural Network (RNN) on such kinematics regression tasks [32, 33].
Due to the long-term discrepancy problem of RNN, Long Short Term Memory model
(LSTM), which was developed for processing sequence information with multi-
layer neuron cells [34], could be proper to decode the continuity correlation
between non-static EMG signals and kinematics under noisy measurements. Thus,
in this study, a comparison was made between LSTM and SVM.

Figure 1: Our teacher-student paradigm used in our method provided a cross-
modal supervision for the . The upper pipeline provided a training
supervision, whereas the bottom learned to predict joint angles using EMG and
IMU signals.

###  II-B Prediction Time

The prediction time, which is denoted by the interval between the timings of
the prediction and the actual movement, is related to the data labelling of
the algorithm, as will be shown in Section III-D. Determining the prediction
time can give a precise reference on how much the mechanical delay can be
compensated. Huang et. al. [26] proposed to make a prediction based on
identifying the critical timing of two stable patterns’ transition. That is,
the proposed prediction time was related to the fixed consequence of gait
phases. However, in this study, it might be difficult to leverage the discrete
information of gait phases to perform the task of continuous kinematics
prediction. Inspired by the fact that EMG signal is generated 25-125 ms before
its corresponding actual motion [35, 36, 37], denoted as the electromechanical
delay (EMD), our algorithm utilizes this characteristic to perform continuous
kinematics prediction and we further explore the maximum of transmission delay
that can be compensated by the prediction time of the algorithm.

###  II-C Feature extraction of EMG signals

Feature extraction of EMG signals is also a challenge for the continuous
prediction due to the lack of features with demonstrated effectiveness and
robustness on continuously predicting lower-limb kinematics. Huang et. al.
made ahead-of-time prediction of discrete information using four time-domain
features [26]. Ngeo et. al. [38] investigated simultaneous and multiple finger
kinematics estimation based on comparing four time-domain features and
envelope of multi-channel surface EMG signals. The result showed that
traditional EMG time-domain features outperformed filtered EMG envelop.
Traditional time-domain features allow us to effectively recognize patterns
and joint continuous parameters. However, they ignored the influence of EMG
fluctuation and the antagonism of muscles, therefore cannot reflect the
magnitude matchup between EMG signals and relative movements. Thus, in this
paper, a feature extractor that could extract in-depth information of multi-
channel EMG signals is enabled by the end-to-end training of the LSTM network
and compared with traditionally used time-domain features.

##  III Methodology

###  III-A Overview of the Proposed Model

Our method, as shown in Fig. 1, follows a teacher-student paradigm. The top
pipeline calculated joint angles using data from the stereo-photogrammetric
framework, which provided a cross-modal supervision; the bottom pipeline
predicted angles using IMU and EMG signals, which performed an ahead-of-time
prediction. To be specific, we developed an end-to-end LSTM network to predict
kinematics with data labelled according to different prediction time. The end-
to-end network was designed that it can be divided into feature extractor and
motion predictor. Other than the LSTM feature extractor and LSTM motion
predictor, we also employed traditionally-used time-domain features and SVM,
in order to contribute to a comparison study.

###  III-B Cross-Modal Supervision

The reason of employing the stereo-photogrammetric framework was to provide a
”gold standard” of joint angles so as to have the angle prediction framework
supervised with accurate labels. As shown in the top pipeline of Fig. 1, the
measured angles were obtained from videos. The measured 3-D locations of
markers were firstly filtered by a 4th order Butterworth low-pass filter
(cutoff frequency 6 Hz). The coordinate system of each segment was set
according to the definition of anatomical orientations and the placement of
markers [39, 40]. Particularly, the markers attached to analytical landmarks
were recognized and tracked by the stereo-photogrammetric system. And the
bone-embedded frames, i.e., the coordinate system of each segment, were
determined according to its definition [41] and at least three markers of each
segment. In this way, a multi-link kinematics model can be built. Then, pose
estimation was applied to extract positions and orientations of segments from
markers by the means of a global optimization process, which was reported with
the advantage of minimizing the soft tissue and measurement error [42, 43].
The pose estimation can be formulated as

$\displaystyle\min\limits_{q}(\sum\limits_{i}^{N}(\omega_{i}p_{i}^{measured}-T(q)p_{i}))$
(1)

where $q$ denoted the generalized coordinates of the multi-link model,
$\omega_{i}$ denoted the weight of the ith landmark, $T(q)$ denoted the
transmission from the local segment-fixed coordinate frame to the laboratory
frame, $p_{i}^{measured}$ denoted the measured position vector of the $i$th
landmark with respect to the laboratory frame and $p_{i}$ denoted the position
vector of the $i$th landmark with respect to the segment-fixed coordinate
frame. Finally, a 7-link inverse kinematics modeling of subjects (the
conventional gait model) was performed after pose estimation, while
anthropometric measurements of subjects were used to scale the model. Joint
angles were calculated using the joint axes of the inverse kinematics model
and the generalized coordinates of the adjacent segments [41, 44].

The stereo-photogrammetric framework took images as input and calculated joint
angles as measured angles $\hat{\theta}$. The measured angles $\hat{\theta}$
provided cross-modal supervision for the angle prediction framework, which
learned to make an ahead-of-time prediction of angles as stated in the
following sections.

###  III-C Signal Preprocessing and Data Windowing

Figure 2: Data windowing scheme. The feature vector was formed by combining
EMG features of a sliding window and the current IMU-calculated knee angle,
where $\theta$ denoted the IMU-calculated knee angle.

As presented in Fig.1, knee angles were firstly calculated by the method
proposed in [45] using IMU measurements. The calculation accuracy of the
method was presented to be 1.67 deg during dynamic motions. Then, the
calculated angles combined with EMG signals segmented by sliding windows were
used to form feature vectors, as shown in Fig.2. To be specific, following the
experience of literatures [38, 26], the length of the sliding window was
148.5ms, and the increment was 13.5ms. The feature vector was constructed by
combining EMG features extracted from nine channels with the IMU-calculated
knee angle at the end of the sliding window, which is given by

$\displaystyle\bm{x}=[\bm{f_{1}},\bm{f_{2}},\ldots,\bm{f_{9}},\theta]$ (2)

where $\bm{f_{i}},i=1,\cdots,9$ denoted the extracted features from each EMG
channel, $\theta$ denoted the calculated knee angles from IMU measurements.

###  III-D Feature Extraction

Features in the time domain (FT): Past studies have demonstrated the
effectiveness of time-domain features traditionally used in PR-based
prediction for EMG-controlled lower-limb prostheses [26, 46]. In order to test
their performance in continuous motion prediction, we first extracted such
traditional features from EMG signals. Four traditionally-used EMG time-domain
features [46] (mean absolute value, number of slope sign changes, number of
zero crossings and waveform length) were simultaneously extracted in this
experiment. As presented in Eq. (2), EMG features, $\bm{f_{n}}$, each
containing four time-domain features were extracted from nine channels to
construct a 37-element feature vector with a knee angle $\theta$.

Features from LSTM (FL): Given that lower-limb movements are driven by the
coordinative and antagonistic activations of muscles, FT, which was extracted
from each single channel, cannot fully reflect the coordinative and
antagonistic efforts of muscle activations, i.e. inhibition and excitation
mechanisms [47], during level walking. Thus, information that represented such
mechanism and/or muscle activation should be included in the extracted
features in order to potentially improve the performance of the continuous
prediction. However, currently used methods cannot perform this task well. The
envelope of EMG signals, although was mostly regarded as muscle activations,
cannot give a full insight of the coordinative and antagonistic efforts of
muscles. To this end, we employed an artificial neural network with the aim of
automatically extracting EMG features, expecting to break the potential
bottleneck. Herein, LSTM, with the capability of forgetting useless
information, jointly processed the EMG signals of all the channels, thus could
effectively extract the inhibition and excitation mechanisms among EMG signals
and highlight the correlation between joint efforts of muscle activations and
movements. Particularly, a four-layer LSTM was used to extract features from
the nine-channel EMG signals. The four-layer extractor consisted of three LSTM
layers (40 as hidden size and 60 time steps) and a 40-by-9 fully connected
layer. The topology was determined by multiple trials with the aim of
realizing the simplest architecture while preserving enough accuracy. In order
to maximally remain the basic information into EMG signals, the input and
output size of the extractor was set to be equal. Then, the nine EMG features
extracted by the feature extractor were combined with the calculated knee
angles and concatenated into the 10-element feature vectors through Eq. (2),
following the data windowing scheme presented in Fig. 2. During the training
session, both the LSTM predictor and the LSTM extractor were trained together
as an end-to-end network. During the working session, we employed the well-
trained feature extractor to work with other predictors.

Features from LSTM and Time domain (FLT): When a comparison was made between
FT and FL, one question of interest was whether they incorporate complementary
information to each other. If so, a combination of FT and FL could provide
more global insight of EMG signals. To answer this question, FT and FL were
combined together, then used to construct a 46-element feature vector for
testing their performance on each predictor.

###  III-E Angle Predictor

Due to its successful usage on processing sequence signal in translation [46],
LSTM was employed to perform the many-to-one mapping between EMG and its
consequent movement. A 5-layer LSTM(40-40-40-80-1, the last two layers of
which were fully connected layers), which incorporated a forget gate and
remember gate, was utilized as a predictor. Particularly, the 5-layer
predictor was consisted of three LSTM layers (40 as hidden size and 60 as time
step) and two fully connected layers with the size of 40*80 and 80 *1.

Following SVM’s demonstrated performance in predicting locomotion modes and
real-time application of EMG-controlled hand prostheses, SVM was selected to
contribute to a comparison study. The regression version of SVM, support
vector regression (SVR) was employed as the other angle predictor. The applied
kernel function was the radial basis kernel function(RBF).

###  III-F Labeling and Prediction Time

Figure 3: Data labeling. Feature vectors were labelled by angles measured by
the stereo-photogrammetric framework. The prediction time was set according to
EMD. $\hat{\theta}$ denoted the angle measured by the framework

As shown in Fig.3, each feature vector was labeled with the knee angle
$\hat{\theta}$ measured by the stereo-photogrammetric framework after
prediction time. The length of the prediction time $P_{1}$ was set inspired by
the physiological time interval of EMD.

To evaluate the influence of prediction time and to explore the time boundary,
six time sections (27 ms, 54 ms, 81 ms, 108 ms, 135ms and 162ms) were selected
considering the requirement of data synchronization and the inspiration of the
normal physiological time interval of EMD. All data were collected and labeled
before the training session. Cure parameters of each well-trained model,
including the LSTM extractor, the SVM predictor and the LSTM predictor, varied
with different sets of prediction time.

##  IV Experiments

###  IV-A Experimental Setup and Design

Ten healthy subjects (eight men and two women, age = 25$\pm$5 years, height =
1.75$\pm$0.05 m, weight = 67.6$\pm$12 kg) are asked to walk with self-selected
speeds. As shown in Fig. 4, nine muscles from one leg were selected to
incorporate most functional muscles relative to normal walking, including:
rectus femoris (RF), vastus lateralis muscle (VL), vastus medialis muscle
(VM), tibialis anterior muscle (TA), soleus (SL), biceps femoris muscle (BF),
semitendinosus (ST), gastrocnemius muscle medial head (GM) and gastrocnemius
muscle lateral head (GL).

Figure 4: Selected muscles and sensor attachment on subjects.

Surface electrodes (Delsys Trigno, IM type $\&$ Avanti type) are placed on
target muscles through palpation and skin preparation. In addition, data from
nine-axis inertial measurement units (IMUs) is provided by the inertial
sensors embedded in electrodes. Sixteen retro-reflective markers are attached
to subjects’ pelvis and lower limbs. The markers are attached to the
analytical landmarks of body [48] according to the instructions and principles
of [48, 40], and guaranteed by palpation [25]. The 3-D locations of the
markers are recorded (100 Hz) using a 8-camera vedio system (Vicon, Oxford,
UK). The signals from EMG sensors (1111.11 Hz), IMUs (74Hz) and the vedio
system (100 Hz) are recorded and synchronized by time stamps.

In the experiment, three minutes of standing still are provided for
initializing the joint angle calculation method. Five trials are performed on
each subject, each trial lasts at least three minutes. Rest periods are
allowed between trials to avoid fatigue. Before starting the experiment,
anthropometric measurements like the height, weight and lengths of each lower-
limb segments are measured in order to scale the model of the stereo-
photogrammetric framework.

The experiment protocol is approved by Chinese Ethics Committee of Registering
Clinical Trials and all participants have been informed of the content and
their right to withdraw from the study at any time, without giving any
explanation.

###  IV-B Implementation Details

The training parameters, i.e. batch size, optimizer, Epoch, and the learning
rate of extractor and predictor is set as batch size= 1, optimizer=Adam,
Epoch=30, the initial learning rate of extractor $\&$ predcitor=0.001 and
0.0001, respectively. Every 20 times of training, the learning rates reduce
20$\%$.

Experiments are constructed to test the performance of different motion
predictors (P), extracted features sets (F) and prediction time (T), during
which each set of the three indicators are evaluated. During such experiments,
10-fold cross-validation is applied to evaluate our algorithms without the
loss of generality. In the cross-validation procedure, data of a subject are
used as the testing database, while data of the remaining subjects are used as
the training database. This procedure is repeated for each subject so that
data from each subject could be used as a testing database once. All the data
collected during experiments form the dataset, which consists of over
90,000,000 sample points of EMG signals, over 670,000 sample points of IMU-
calculated knee angles and over 900,000 sample points of measured knee angles.
A comparison is also made to evaluate the prediction’s reliance on periodic
gait characteristics versus EMG, in order to reveal the collaborated
contribution of EMG and kinematic signals.

With the well-trained LSTM extractor, SVR predictor and LSTM predictor, the
predicted knee angles are analyzed to evaluate the performance of different
sets of (P, F, T). In order to distinguish the separate influence of the
extracted features (F) and prediction time (T), we present our results by
averaging each of them, which is

$\displaystyle
I_{F}=\frac{\sum_{T}I_{F,T}}{n_{T}},I_{T}=\frac{\sum_{F}I_{F,T}}{n_{F}}$ (3)

where $I$ denotes any evaluation index (e.g. RMSE, SNR, R-value or adjusted
R), F denotes a feature extractor, T denotes a time section and $n_{T}$
,$n_{F}$ denote the amount of time sections and the number of feature
extractors, which is equal to 4 and 3 respectively. The one-way ANOVA is
performed on the results to depict the repeated measures analysis of variance.
A significant level is set to $\alpha$ = 0.05 for all testing.

Figure 5: The performance plots of different prediction time.

(a) RMSE (b) SNR (c) R-value

Figure 6: (a) RMSE, (b) SNR and (c) R-value of the performance of different
prediction time over ten able-bodied subjects. $*$ indicates a statistically
significant difference (one-way ANOVA, P$<$0.05).

###  IV-C Evaluation of The Prediction Time

As discussed above, six time sections are labeled as prediction time to access
the interval that could be used to compensate for mechanical delay with the
best performance. To do this, the Pearsons correlation coefficient (R-value),
the root mean square error (RMSE), the signal-to-noise ratio (SNR) and the
adjusted R square are calculated with different prediction time.

Fig. 6 presents the average performance of the six sections of prediction
time. It can be seen from the figure that the accuracy of the last two time
sections (135ms and 162ms) is obviously worse than the rest four sections of
prediction time, regardless of the metrics. The relatively inferior
performance of the last two time sections indicates the boundary of the
prediction time.

Fig. 6 presents the comparison of different prediction time sections depicted
by RMSE, SNR and R-value. Note that in this figure, only four time sections
are presented. Different grey levels denote different prediction time. In Fig.
6(a), there is no significant difference of SVR predictor among all the
prediction time (P$>$0.05). The best results are 12.99 deg for 27ms, 14.31 deg
for 54ms, 14.34 deg for 81ms and 14.36 deg for 108ms, respectively. Also, the
statistical analysis shows there is no significant difference (P$>$0.05) among
the RMSE of the LSTM predictor. The best results of 27ms, 54ms, 81ms and 108ms
are 0.66 deg, 0.69 deg, 1.48 deg, 1.36 deg, respectively.

As shown Fig. 6(b), the SNRs of the SVM predictor among different prediction
time are statistically equal (P$>$0.05). The best results are 2.49dB for 27ms,
2.52dB for 54ms, 2.52dB for 81ms and 2.43dB for 108ms, respectively. Averaged
SNR of the LSTM predictor generally decreases with the increment of prediction
time. The statistical analysis shows there is no significant difference
between 54ms and 108ms (P$>$0.05). The best results are 24.47dB for 27ms,
25.3dB for 54ms, 24.87dB for 81ms and 23.98dB for 108ms, respectively.

In Fig. 6(c), the averaged R-value of both predictors presents a generally
decreasing trend as the prediction time increases. No significant difference
is found in the results of the SVM predictor between 27ms and 108ms
(P$>$0.05). And the best results are 0.89 for 27ms, 0.84 for 54ms, 0.76 for
81ms and 0.67 for 108ms, respectively. For the LSTM predictor, a significant
difference exists between 27ms and 108ms (P=0.009) and between 54ms and 108ms
(P=0.0248). The best results for the LSTM predictor are 0.91 for 27ms, 0.91
for 54ms, 0.908 for 81ms and 0.908 for 108ms, respectively.

It is in a statistic paradigm that the prediction time sections are set and
the performance of different prediction time sections is evaluated. Accuracy
improvement can be expected if subject-specific EMD-inspired prediction time
is employed for predicting joint angles. However, EMD varies with muscle
fatigue and across the continuous repetitive motions. The intractable online
detection of EMD time would make such a prediction time setting paradigm
impossibly achievable for exoskeletons. Thus, we employ time sections, rather
than a specific value, in order to study the effects of EMD-inspired
prediction time on angle prediction and to provide a quantified reference for
the compensation of exoskeletons’ transmission delay. According to the
performance of all the six time sections, the boundary of prediction time is
experimentally explored. In addition, the prediction time influences the
prediction accuracy. Compared with [13] that predicted one time-step ahead,
our RMSEs are slightly larger, which might be attributed to our significantly
larger prediction time. And the performance of prediction time between 27ms
and 108ms does not show a significant difference and the performance of
prediction time out of this interval presents obvious deterioration. This
result indicates that any value in the time interval of 27ms and 108ms can be
determined as the prediction time. The prediction time section for controlling
exoskeletons is recommended to be initially set between 54ms and 81ms and
tuned according to the predicting performance. The inference time of executing
our algorithm is 5-7ms, achieved by Delsys Trigno SDK module and NVIDIA TX2
computing module. In this way, our prediction time is large enough to cover
the execution of our algorithm and still leaves enough room for compensating
the mechanical transmission delay.

###  IV-D Effectiveness of The Predictors

For the purpose of evaluating the performance of SVR and LSTM on the
continuous kinematics prediction, an accuracy comparison needs to be made
between the two predictors. We use R-value, SNR and RMSE to quantify the
performance.

As presented in Fig. 6, the overall performance of the predictors depicted by
all the three indexes present that the LSTM predictor significantly
outperforms the SVR predictor. (P$<$0.001). It is also shown in Fig. 9 and
Fig.9 that all the four indexes of the LSTM predictor are obviously better
than those of the SVM predictor. It should be noted that the relatively larger
RMSE of LSTM shown in Fig. 6 is due to the averaging over all the feature
sets. It can be seen in Fig. 9 (a) that the LSTM predictor with the FL feature
set is of good accuracy, with the RMSE of 3.98 deg.

The devised prediction method, which conducts an end-to-end kinematics
prediction network (KinPreNet), achieves the best accuracy of predicting knee
angles. Moreover, the accuracy of the KinPreNet is better than that of the
IMU-based angle calculation method we use [45], which indicates the
predictor’s ability of improving the errors of the input. And our better
results, compared with the results of RCNN for predicting upper-limb motions
[32], suggest the benefit of our architecture. Based on the results, it can be
concluded that LSTM is more suitable for making a continuous and ahead-of-time
prediction of knee angles. LSTM model benefits from its capability of learning
to extract information deeply hidden in the input features through the
regression process. Deep hierarchical representations of input features can be
explored by the devised LSTM predictor. The SVR predictor, with the lesser
capability of exploring the hidden information of features, could not explore
the unobvious characters of the biomedical signals for this study.

Figure 7: Statistic adjusted R-square results averaged over ten able-bodied
subjects. The adjusted R-square results refer to all designed EMG features
using the SVM and LSTM predictors were shown for individual prediction times:
(a) 27ms, (b) 54ms, (c) 81ms, (d) 108ms. $*$ indicates a statistically
significant difference (one-way ANOVA, P$<$0.05). FT, FL and FTL denote EMG
features from time domain, EMG features from LSTM and EMG features from time
domain and LSTM, respectively.

Figure 8: Statistic results of prediction comparison based on different
features over ten able-bodied subjects.$*$ indicates a statistically
significant difference (one-way ANOVA, P$<$0.05). (a) FT, (b) FL and (c) FTL
denote EMG features from time domain, EMG features from LSTM and EMG features
from time domain and LSTM, respectively.

Figure 9: Statistic results of prediction comparison based on soly kinematics
signals and fusion of kinematics and EMG signals over ten able-bodied
subjects. $*$ indicates a statistically significant difference (one-way ANOVA,
P$<$0.05). (a) RMSE, (b) R-value and (c) SNR denote different indexes for
evaluation.

###  IV-E Validity of The Extracted Feature Sets

To validate the completeness of information incorporated in each feature set
and their consequent performance on angle prediction, R value, RMSE and SNR
are estimated to present prediction accuracy. In addition, the adjusted R
square is employed to evaluate to what degree the accuracy variation of
predicted angles is correlated with different feature sets.

The average adjusted R-square for each prediction time of ten subjects is
shown in Fig. 9. The performance of FL outperforms the performance of FT and
FTL. FL significantly outperforms FT with all prediction time in the results
of the SVR predictor (P$<$0.05). Although the adjusted R square of FL is shown
to be better than FTL in both predictors, the difference of FL and FTL are not
statistically significant (P$>$0.05). The FL-based prediction achieves better
performance excepted for the prediction from prediction time of 54ms and 108ms
with the LSTM predictor. In addition, with FL as the extracted feature, the
average adjusted R-square of the LSTM predictor is 27$\%$ \- 54$\%$ higher
than that of the SVR predictor. Statistical comparison of adjusted R-square
among prediction time shows that significant difference is not found when
prediction time changes for any EMG feature set (P$>$0.05).

Fig.9 depicts the performance with different EMG feature sets. The results of
different EMG feature sets are represented by different grey levels. In Fig.
9(a), averaged RMSE of SVR predictor are performed with three different EMG
feature sets (FL $>$ FTL $>$ FT). The best results of FT, FL and FTL are 15.95
deg, 15.07 deg and 12.99 deg, respectively. No significant difference is found
among features ( P(FT, FL)=0.077 ). Averaged results with the LSTM predictor
with FT, FL and FTL are 10.71$\pm$4.30 deg, 3.98$\pm$3.26 deg and
7.13$\pm$4.60 deg, respectively. Significant differences are found in the
results of the LSTM predictor among all the EMG feature sets (P$<$0.01).

In Fig. 9(b), SNR of the LSTM predictor significantly outperforms that of the
SVR predictor. There is no significant difference in the results among
different EMG feature sets of the SVR predictor (P$>$0.05). The best results
of FT, FL and FTL are 2.09dB, 2.52dB and 1.95dB, respectively. Average SNR of
the LSTM predictor are depicted with different EMG feature sets (FL $>$ FTL
$>$ FT). Significant differences are found in the results of the LSTM
predictor among all the features (P$<$0.01).

In Fig. 9(c), there are obviously significant differences existed in R-value
of the SVM predictor between FT $\&$ FL and FL $\&$ FTL (P$<$0.001).
Significant differences also exist in pairs of FT and FTL (P=0.0277). The
average R-value of FL can up to 0.7$\pm$0.13 which is about 34$\%$ \- 40$\%$
higher than that of other EMG feature sets. Similar to SVR predictor,
significant differences are found in R-value results of the LSTM predictor
between FT $\&$ FL and FL $\&$ FTL (P$<$0.005). The average R-value of FL with
LSTM predictor is 0.88$\pm$0.04.

The LSTM feature extractor is established with the help of the LSTM predictor.
Obtaining model parameters of the LSTM prediction system in the training
session, the former part of this system, defined as the LSTM extractor, can be
used to extract features from EMG signals. As shown in Fig. 9, the R-value,
RMSE and SNR of FL outperform those of FT and FLT with significant
differences, which validates the effectiveness of the LSTM extractor on
prediction accuracy. The adjusted R square value of FL demonstrates the LSTM
extractor’s high correlation with the performance, which represents its
capability of explaining knee angle variations. In addition, no significant
difference exists among different sets of predictors and prediction time,
indicating the stability of the LSTM extractor.

The reason of such outstanding performance is twofold. Firstly, the joint
training of the LSTM extractor and predictor, which propagates the error of
angle prediction back to the extractor during the training session, reinforces
the correlation between the extracted EMG features and knee angles. The
reinforced correlation contributes to well-decoded information from multiple
EMG channels. The results of [38], without presenting the capability of the
ahead-of-time prediction, presented relatively lower accuracy, which also
gives a side proof of the benefit of joint training. Secondly, the LSTM
extractor explores deep information inside the multi-channel EMG signals.
LSTM, developed for processing sequence information with multi-layer neuron
cells, is proper to decode the continuity correlation between non-static EMG
signals and kinematics under noisy measurements. Meanwhile, rather than
extracting features from every single channel, the LSTM extractor extracts
features from multiple channels simultaneously contributing to a comprehensive
metric among muscles, which is related to muscles antagonism during level
walking.

Surprisingly, the performance of FTL is presented to be just slightly better
than that of FT, while no improvement is observed comparing with the
performance of FL. FT, representing the overall information of the signal in a
sliding analysis window, is just general time-domain features of EMG signals.
In contrast, features extracted from the LSTM extractor, due to the function
of the remember gate and forget gate in the LSTM model, magnify the effect of
dynamic components and minify the effect of common resting components of EMG.
Thus, FTL, regarded as FL corrupted by FT, still contributes to a better
prediction than FT. Each FT is extracted from a single channel of EMG
electrodes, the effectiveness of which does not depend on FT from other
channels. On the contrary, FL that is extracted with joint information of
multi-channel EMG signals works as a whole in predicting knee angles. Hence,
if compared with FL, mixing FT with FL brings some redundancy information into
the extracted features, which contributes to a worsened performance. But
compared with FT, such a mixture adds some deep characters into features,
which improves the accuracy of the prediction.

###  IV-F Effectiveness of Employing EMG Signals

The ahead-of-motion prediction might result from two factors: the pseudo-
periodic characteristics of gait and the electromechanical delay. In order to
distinguish the influence of employing EMG signals, a comparison between
predictions solely from kinematic signals and the fusion of kinematic and EMG
signals is performed. Due to the outstanding performance of FL combining with
the LSTM predictor, the comparison under other sets of the feature extractors
and predictors is with little value. Thus, the comparison is just performed
using FL and the LSTM predictor.

Fig. 9 presents the prediction performance with different prediction time
using solely kinematics signals and the fusion of EMG and kinematics signals.
There are significant differences existing in RMSE and SNR between the
prediction from different signals, regardless of prediction time (P $<$ 0.05).
And it can be seen from Fig.9 that the performance of the prediction from EMG
and kinematics signals generally outperforms that from solely kinematics
signals.

The performance comparison between prediction from kinematics signal and the
fusion of EMG and kinematics signals demonstrates the benefits of employing
EMG signals in the regression-based motion predicting. It can be concluded
that although the pseudo-periodic characteristics of gait could make a
contribution to the ahead-of-motion prediction, employing EMG signals plays a
necessary role given the significantly improved accuracy. This phenomenon
meets the significantly smaller prediction time reported in [13] that solely
exploited the pseudo-periodic characteristics of gait. Two factors might
result in performance improvement. Firstly, EMG signals from antagonistic
muscle pairs around the knee relate to the moment and angular acceleration of
the knee, which incorporate vital information of knee motion changes. Thus,
the intra-gait knee angle changes, which can hardly be predicted by the
periodic gait characteristics, could be covered by features from EMG signals.
Secondly, EMD provides ahead-of-motion information, which could consequently
improve the prediction accuracy.

## Acknowledgment

This work is partly funded by National Natural Science Foundation of China
(No. 61872112) and National Key Research and Development Program of China
(Nos. 2018YFC0806802 and 2018YFC0832105).

## References

  * [1] A. E. Emery, “Population frequencies of inherited neuromuscular diseases—a world survey,” _Neuromuscular Disorders_ , vol. 1, no. 1, pp. 19 – 29, 1991\.
  * [2] K. M. Steele, A. Rozumalski, and M. H. Schwartz, “Muscle synergies and complexity of neuromuscular control during gait in cerebral palsy,” _Developmental Medicine & Child Neurology_, vol. 57, no. 12, pp. 1176–1182, 2015.
  * [3] E. Martini, S. Crea, A. Parri, L. Bastiani, U. Faraguna, Z. McKinney, R. Molino-Lova, L. Pratali, and N. Vitiello, “Gait training using a robotic hip exoskeleton improves metabolic gait efficiency in the elderly,” _Scientific reports_ , vol. 9, no. 1, pp. 1–12, 2019.
  * [4] G. Durandau, D. Farina, G. Prieto, I. Dimbwadyo-Terrer, S. Lerma Lara, J. Pons, J. Moreno, and M. Sartori, “Voluntary control of wearable robotic exoskeletons by patients with paresis via neuromechanical modeling,” _Journal of NeuroEngineering and Rehabilitation_ , vol. 16, 07 2019.
  * [5] L. Awad, P. Kudzia, D. Arumukhom Revi, T. Ellis, and C. Walsh, “Walking faster and farther with a soft robotic exosuit: Implications for post-stroke gait assistance and rehabilitation,” _IEEE Open Journal of Engineering in Medicine and Biology_ , vol. PP, pp. 1–1, 04 2020.
  * [6] V. R. Garate, A. Parri, T. Yan, M. Munih, R. M. Lova, N. Vitiello, and R. Ronsse, “Walking assistance using artificial primitives: a novel bioinspired framework using motor primitives for locomotion assistance through a wearable cooperative exoskeleton,” _IEEE Robotics & Automation Magazine_, vol. 23, no. 1, pp. 83–95, 2016.
  * [7] A. Razak, A. Hadi, A. Zayegh, R. K. Begg, and Y. Wahab, “Foot plantar pressure measurement system: A review,” _Sensors_ , vol. 12, no. 7, pp. 9884–9912, 2012.
  * [8] L. M. Mooney, E. J. Rouse, and H. M. Herr, “Autonomous exoskeleton reduces metabolic cost of human walking during load carriage,” _Journal of neuroengineering and rehabilitation_ , vol. 11, no. 1, p. 80, 2014.
  * [9] Y. Ding, I. Galiana, A. T. Asbeck, S. M. M. De Rossi, J. Bae, T. R. T. Santos, V. L. de Araujo, S. Lee, K. G. Holt, and C. Walsh, “Biomechanical and physiological evaluation of multi-joint assistance with soft exosuits,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 25, no. 2, pp. 119–130, Feb 2017.
  * [10] Y. Ding, F. A. Panizzolo, C. Siviy, P. Malcolm, I. Galiana, K. G. Holt, and C. J. Walsh, “Effect of timing of hip extension assistance during loaded walking with a soft exosuit,” _Journal of neuroengineering and rehabilitation_ , vol. 13, no. 1, p. 87, 2016.
  * [11] S. Fang, A. L. Kinney, M. E. Reissman, and T. Reissman, “Simulation on the effect of gait variability, delays, and inertia with respect to wearer energy savings with exoskeleton assistance,” in _2019 IEEE 16th International Conference on Rehabilitation Robotics (ICORR)_. IEEE, 2019, pp. 506–511.
  * [12] G. Chen, P. Qi, Z. Guo, and H. Yu, “Gait-event-based synchronization method for gait rehabilitation robots via a bioinspired adaptive oscillator,” _IEEE Transactions on Biomedical Engineering_ , vol. 64, no. 6, pp. 1345–1356, 2016.
  * [13] K. Tanghe, F. De Groote, D. Lefeber, J. De Schutter, and E. Aertbeliën, “Gait trajectory and event prediction from state estimation for exoskeletons during gait,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , pp. 1–1, 2019.
  * [14] Y. Sankai, “HAL: Hybrid assistive limb based on cybernics,” in _Robotics research_. Springer, 2010, pp. 25–34.
  * [15] Y. Huang, K. B. Englehart, B. Hudgins, and A. D. Chan, “A Gaussian mixture model based classification scheme for myoelectric control of powered upper limb prostheses,” _IEEE Transactions on Biomedical Engineering_ , vol. 52, no. 11, pp. 1801–1811, 2005.
  * [16] S. Amsüss, P. M. Goebel, N. Jiang, B. Graimann, L. Paredes, and D. Farina, “Self-correcting pattern recognition system of surface emg signals for upper limb prosthesis control,” _IEEE Transactions on Biomedical Engineering_ , vol. 61, no. 4, pp. 1167–1176, 2013.
  * [17] A. D. Chan and K. B. Englehart, “Continuous myoelectric control for powered prostheses using hidden Markov models,” _IEEE Transactions on Biomedical Engineering_ , vol. 52, no. 1, pp. 121–124, 2004.
  * [18] O. Van Der Niet Otr, H. A. Reinders-Messelink, R. M. Bongers, H. Bouwsema, and C. K. Van Der Sluis, “The i-LIMB hand and the DMC plus hand compared: a case report,” _Prosthetics and orthotics international_ , vol. 34, no. 2, pp. 216–220, 2010.
  * [19] A. M. Grabowski and S. D’Andrea, “Effects of a powered ankle-foot prosthesis on kinetic loading of the unaffected leg during level-ground walking,” _Journal of neuroengineering and rehabilitation_ , vol. 10, no. 1, p. 49, 2013\.
  * [20] J. C. Perry and J. Rosen, “Design of a 7 Degree-of-Freedom upper-limb powered exoskeleton,” in _The First IEEE/RAS-EMBS International Conference on Biomedical Robotics and Biomechatronics, 2006. BioRob 2006._ IEEE, 2006, pp. 805–810.
  * [21] C. Connolly, “Prosthetic hands from touch bionics,” _Industrial Robot: An International Journal_ , vol. 35, no. 4, pp. 290–293, 2008.
  * [22] J. Liu and P. Zhou, “A novel myoelectric pattern recognition strategy for hand function restoration after incomplete cervical spinal cord injury,” _IEEE transactions on neural systems and rehabilitation engineering_ , vol. 21, no. 1, pp. 96–103, 2012.
  * [23] T. Matsubara and J. Morimoto, “Bilinear modeling of EMG signals to extract user-independent features for multiuser myoelectric interface,” _IEEE Transactions on Biomedical Engineering_ , vol. 60, no. 8, pp. 2205–2213, 2013\.
  * [24] Z. Li, B. Wang, C. Yang, Q. Xie, and C.-Y. Su, “Boosting-based EMG patterns classification scheme for robustness enhancement,” _IEEE Journal of Biomedical and Health Informatics_ , vol. 17, no. 3, pp. 545–552, 2013.
  * [25] H. Huang, T. A. Kuiken, R. D. Lipschutz _et al._ , “A strategy for identifying locomotion modes using surface electromyography,” _IEEE Transactions on Biomedical Engineering_ , vol. 56, no. 1, pp. 65–73, 2008.
  * [26] H. Huang, F. Zhang, L. J. Hargrove, Z. Dou, D. R. Rogers, and K. B. Englehart, “Continuous locomotion-mode identification for prosthetic legs based on neuromuscular–mechanical fusion,” _IEEE Transactions on Biomedical Engineering_ , vol. 58, no. 10, pp. 2867–2875, 2011.
  * [27] G. R. Naik, S. E. Selvan, S. P. Arjunan, A. Acharyya, D. K. Kumar, A. Ramanujam, and H. T. Nguyen, “An ICA-EBM-based sEMG classifier for recognizing lower limb movements in individuals with and without knee pathology,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 26, no. 3, pp. 675–686, 2018.
  * [28] J. Han, Q. Ding, A. Xiong, and X. Zhao, “A state-space EMG model for the estimation of continuous joint movements,” _IEEE Transactions on Industrial Electronics_ , vol. 62, no. 7, pp. 4267–4275, 2015.
  * [29] D. Ao, R. Song, and J. Gao, “Movement performance of human–robot cooperation control based on emg-driven hill-type and proportional models for an ankle power-assist exoskeleton robot,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 25, no. 8, pp. 1125–1134, 2017\.
  * [30] G. Durandau, D. Farina, and M. Sartori, “Robust real-time musculoskeletal modeling driven by electromyograms,” _IEEE Transactions on Biomedical Engineering_ , vol. 65, no. 3, pp. 556–564, March 2018.
  * [31] J. A. Brantley, T. P. Luu, S. Nakagome, and J. L. Contreras-Vidal, “Prediction of lower-limb joint kinematics from surface EMG during overground locomotion,” in _2017 IEEE International Conference on Systems, Man, and Cybernetics (SMC)_. IEEE, 2017, pp. 1705–1709.
  * [32] P. Xia, J. Hu, and Y. Peng, “EMG-based estimation of limb movement using deep learning with recurrent convolutional neural networks,” _Artificial organs_ , vol. 42, no. 5.
  * [33] L. Wang and T. S. Buchanan, “Prediction of joint moments using a neural network model of muscle activations from EMG signals,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 10, no. 1, pp. 30–37, 2002.
  * [34] M. Karimi, D. Wu, Z. Wang, and Y. Shen, “DeepAffinity: interpretable deep learning of compound–protein affinity through unified recurrent and convolutional neural networks,” _Bioinformatics_ , vol. 35, no. 18, pp. 3329–3338, 02 2019.
  * [35] J. T. Blackburn, D. R. Bell, M. F. Norcross, J. D. Hudson, and L. A. Engstrom, “Comparison of hamstring neuromechanical properties between healthy males and females and the influence of musculotendinous stiffness,” _Journal of Electromyography and Kinesiology_ , vol. 19, no. 5, pp. e362–e369, 2009.
  * [36] L. Li and B. S. Baum, “Electromechanical delay estimated by using electromyography during cycling at different pedaling frequencies,” _Journal of Electromyography and Kinesiology_ , vol. 14, no. 6, pp. 647 – 652, 2004.
  * [37] H. Begovic, G. Zhou, T. Li, Y. Wang, and Y.-P. Zheng, “Detection of the electromechanical delay and its components during voluntary isometric contraction of the quadriceps femoris muscle,” _Frontiers in physiology_ , vol. 5, p. 494, 12 2014.
  * [38] J. G. Ngeo, T. Tamei, and T. Shibata, “Continuous and simultaneous estimation of finger kinematics using inputs from an EMG-to-muscle activation model,” _Journal of neuroengineering and rehabilitation_ , vol. 11, no. 1, p. 122, 2014.
  * [39] R. B. Davis, S. Õunpuu, D. Tyburski, and J. R. Gage, “A gait analysis data collection and reduction technique,” _Human Movement Science_ , vol. 10, no. 5, pp. 575 – 587, 1991.
  * [40] D. Sainsbury, “Serge van sint jan, color atlas of skeletal landmark definitions , churchill livingstone/elsevier (2007) isbn 9-78-0443-10315-5 181,” _Physiotherapy_ , vol. 95, pp. 234–235, 09 2009.
  * [41] A. Cappozzo, F. Catani, U. Della Croce, and A. Leardini, “Position and orientation in space of bones during movement: anatomical frame definition and determination,” _Clinical Biomechanics_ , vol. 10, no. 4, pp. 171 – 178, 1995.
  * [42] T.-W. Lu and J. O’connor, “Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints,” _Journal of biomechanics_ , vol. 32, no. 2, pp. 129–134, 1999.
  * [43] A. Leardini, C. Belvedere, F. Nardini, N. Sancisi, M. Conconi, and V. Parenti-Castelli, “Kinematic models of lower limb joints for musculo-skeletal modelling and optimization in gait analysis,” _Journal of biomechanics_ , vol. 62, pp. 77–86, 2017.
  * [44] C. Charalambous, “Measurement of lower extremity kinematics during level walking,” _Classic Papers in Orthopaedics_ , pp. 397–398, 01 2014.
  * [45] C. Yi, J. Ma, H. Guo, J. Han, H. Gao, F. Jiang, and C. Yang, “Estimating three-dimensional body orientation based on an improved complementary filter for human motion tracking,” _Sensors_ , vol. 18, no. 11, p. 3765, 2018.
  * [46] I. Sutskever, O. Vinyals, and Q. V. Le, “Sequence to sequence learning with neural networks,” in _Advances in neural information processing systems_ , 2014, pp. 3104–3112.
  * [47] H. Geyer and H. Herr, “A muscle-reflex model that encodes principles of legged mechanics produces human walking dynamics and muscle activities,” _IEEE transactions on neural systems and rehabilitation engineering : a publication of the IEEE Engineering in Medicine and Biology Society_ , vol. 18, pp. 263–73, 04 2010.
  * [48] A. Leardini, F. Biagi, A. Merlo, C. Belvedere, and M. G. Benedetti, “Multi-segment trunk kinematics during locomotion and elementary exercises,” _Clinical Biomechanics_ , vol. 26, no. 6, pp. 562 – 571, 2011\.


# Images of magnetospheric reconnection-powered radiation around supermassive
black holes

Benjamin Crinquand<EMAIL_ADDRESS>Univ. Grenoble Alpes, CNRS, IPAG,
38000 Grenoble, France
Department of Astrophysical Sciences, Peyton Hall, Princeton University,
Princeton, NJ 08544, USA  Benoît Cerutti  Univ. Grenoble Alpes, CNRS, IPAG,
38000 Grenoble, France
Guillaume Dubus  Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
Kyle Parfrey  School of Mathematics, Trinity College Dublin, Dublin 2, Ireland
Alexander Philippov  Center for Computational Astrophysics, Flatiron
Institute, 162 Fifth Avenue, New York, NY 10010, USA  Department of Physics,
University of Maryland, College Park, MD 20742, USA

(September 3, 2024)

###### Abstract

Accreting supermassive black holes can now be observed at the event-horizon
scale at millimeter wavelengths. Current predictions for the image rely on
hypotheses (fluid modeling, thermal electrons) which might not always hold in
the vicinity of the black hole, so that a full kinetic treatment is in order.
In this letter, we describe the first 3D global general-relativistic particle-
in-cell simulation of a black-hole magnetosphere. The system displays a
persistent equatorial current sheet. Synthetic radio images are computed by
ray-tracing synchrotron emission from nonthermal particles accelerated in this
current sheet by magnetic reconnection. We identify several time-dependent
features of the image at moderate viewing angles: a variable radius of the
ring, and hot spots moving along it. In this regime, our model predicts that
most of the flux of the image lies inside the critical curve. These results
could help understand future observations of black-hole magnetospheres at
improved temporal and spatial resolution.

#### Introduction

At least part of the formidable energy released by active galactic nuclei
(AGN) may originate from the rotational energy of the central black hole,
extracted by the means of strong magnetic fields threading its event horizon
[1]. These magnetic fields are sustained by currents flowing in an accretion
flow, which in the case of a low-luminosity AGN is hot, geometrically thick,
and collisionless. Recently, for the first time, the Event Horizon Telescope
(EHT) collaboration was able to spatially resolve the immediate vicinity of
the low-luminosity AGN M87*. The image shows a circular asymmetric ring
encompassing a depression in brightness [2]. From the high polarization
fraction and brightness temperature of the radiation, we know that the light
captured by the EHT is optically thin synchrotron emission emitted by mildly
relativistic leptons. These measurements are consistent with the scenario of
dynamically important poloidal magnetic fields at event horizon scales, thus
favoring the “magnetically arrested disk” (MAD) scenario for the accretion
flow [3, 4].

However, the location of the main emission sites, as well as the mechanisms
causing plasma heating and particle acceleration, are still poorly
constrained. A widely adopted approach has been to interpret the emission as
synchrotron radiation by thermal [e.g. 5, 6] or non-thermal [7, 8, 9]
electrons, energized for example by subgrid kinetic turbulent dissipation or
magnetic reconnection. General-relativistic magnetohydrodynamic (GRMHD)
simulations can be educated with prescriptions for electron heating used to
model the local synchrotron emissivity, which allows one to synthesize images
when coupled with GR ray-tracing algorithms. In this framework, the image
should prominently display a black-hole shadow determined by the properties of
the spacetime [10, 11, 12]. It is still unclear, though, to what extent the
features of the image are related to the astrophysical details of the emission
[13, 14, 15, 16].

Recently, several GRMHD studies of MAD disks (Ripperda  _et al._ [17],
Chashkina  _et al._ [18] in 2D and Ripperda  _et al._ [19] in 3D) have
reported dynamical transitions from an accreting state, during which magnetic
flux is accumulated onto the horizon, to a highly magnetized quasi-force-free
state, which is characterized by a stark drop in the accretion rate and the
existence of a thin equatorial current sheet. This configuration favors fast
dissipation of the magnetic energy accumulated in the magnetosphere through
magnetic reconnection in the current sheet, triggering efficient particle
acceleration and high-energy nonthermal radiation [20, 21]. Thus, it has been
suggested that gamma-ray flares from radio galaxies (such as M87*) could
naturally occur when the central supermassive black hole is in this quasi-
force-free state, hereafter labeled “high-energy flaring state”.

In this regime, with the innermost zone almost depleted of electron-ion
plasma, predictions for the image based on GRMHD do not apply, whereas kinetic
simulations can accurately capture the physics and dynamics of the highly
magnetized and collisionless plasma. Global kinetic particle-in-cell (PIC)
simulations of black-hole magnetospheres have only recently begun to be
carried out [22, 23, 24, 25, 26]. Crinquand  _et al._ [24] (hereafter, C21)
previously studied the dynamics and high-energy gamma emission of the
equatorial current sheet that develops generically within the ergosphere from
a large class of ordered magnetic configurations [27]. Conversely, in this
letter, we characterize the low-energy emission at radio wavelengths from a
black hole in a high-energy flaring state, embedded in a quasi-force-free
magnetosphere, and discard any contribution from the accretion flow.

#### Simulations

We use the GRPIC code GRZeltron to simulate magnetospheric pair plasma. The
background spacetime is described by the Kerr metric with a dimensionless spin
parameter $a=0.99$. We use spherical horizon-penetrating Kerr-Schild
coordinates $\left(r,\theta,\varphi\right)$, and model self-consistent plasma
injection by including inverse Compton scattering and photon-photon
annihilation [23]. We have run both a 2D and a 3D simulation, using the setup
described in C21 with initially paraboloidal magnetic field lines switching
polarity at the equator. The 3D simulation has a resolution of
$1024\left(r\right)\times 256\left(\theta\right)\times
512\left(\varphi\right)$, whereas the 2D one has $1536\left(r\right)\times
1024\left(\theta\right)$. The simulation domain is
$r\in[0.9\,r_{\mathrm{h}},10\,r_{\mathrm{g}}],\theta\in[0.1,\pi-0.1],\varphi\in[0,2\pi]$,
with $r_{\mathrm{g}}$ the gravitational radius of the black hole and
$r_{\mathrm{h}}$ the radius of the event horizon. Our fiducial 2D and 3D runs
do not have synchrotron cooling. Because of a higher numerical cost, the
values of the dimensionless input parameters are slightly less realistic in 3D
than in 2D: the fiducial Larmor radius $r_{\mathrm{L}}$ of particles is
$r_{\mathrm{L}}=10^{-4}\,r_{\mathrm{g}}$ in 3D instead of $2\times
10^{-6}\,r_{\mathrm{g}}$ in 2D, whereas the background soft photon energy is
$\varepsilon_{0}=10^{-2}\,m_{\mathrm{e}}c^{2}$ instead of $10^{-3}$. In both
cases, the scale separation is necessarily smaller than that of a realistic
AGN. M87*, for example, is characterized by a ratio
$r_{\mathrm{L}}/r_{\mathrm{g}}\approx 10^{-14}$, which is currently unfeasible
numerically. We have also confirmed that 2D simulations are not sensitive to
the precise initial distribution of plasma and photons.

The general behavior of the 3D simulation is very similar to that of the 2D
one. The magnetic field develops a strong toroidal component, with opposite
signs above and below the equator. Thus, a current sheet is needed to support
the associated discontinuity. Pairs are accelerated by nonideal electric
fields in the equatorial plane, triggering high-energy photon upscattering and
even stronger pair creation. The densities are significantly larger in the
current sheet than in the polar regions. The initial current sheet quickly
fragments into flux ropes (see Fig. 1). The 3D evolution of the reconnecting
current sheet is governed by the tearing instability, which produces flux
ropes elongated in the direction orthogonal to the reconnecting field, and the
drift-kink instability, which can corrugate the sheet along the reconnecting
field [20, 28]. Besides, the flux ropes are orbiting around the black hole and
are constantly being sheared by differential rotation, and plasma can be
pushed away from the reconnecting sheet and ejected along the separatrix (the
last open field line connected to the black hole). These competing effects
drive the current sheet toward a state of twisted and disordered flux ropes.
It should be pointed out that our simulation ran for around
$40\,r_{\mathrm{g}}/c$, which is too short to observe the full cyclic behavior
observed in C21, during which a giant magnetic island, formed between the
current sheet and the separatrix, replenishes the black hole with magnetic
flux. Although these flux ropes could lose coherence for longer integration
times, we note that the simulation has run in steady state for about
$20\,r_{\mathrm{g}}/c$, which is longer than the characteristic growth time of
the tearing and drift-kink instabilities which could disrupt it [29] and
longer than the orbital period at the event horizon
($4\pi/\omega_{\mathrm{h}}\sim 14\,r_{\mathrm{g}}/c$, with
$\omega_{\mathrm{h}}=ac/2r_{\mathrm{h}}$).

Figure 1: Snapshot of the total plasma density, normalized by the fiducial
Goldreich-Julian density, at $t=40\,r_{\mathrm{g}}/c$. Magnetic field lines
are represented as blue solid lines. The black surface marks the event
horizon.

#### Ray-tracing

Figure 2: Time-averaged synchrotron images from the 2D simulation. The flux is
in arbitrary units. The three columns show three different viewing angles
$\alpha_{\mathrm{obs}}$, measured with respect to the spin axis of the black
hole. The white cross marks the position of the black hole. The blue dashed
line denotes the critical curve for $a=0.99$ at each viewing angle.

We use the public ray-tracing code geokerr [30] to perform geodesic
integration of synchrotron photons from the emission point to a distant
observer’s screen, in post-processing (see the Supplemental Material, which
includes Ref. [31]). When analyzing the 2D simulation, we focus on the time
dependence on the image. By virtue of axisymmetry, deposition of photons on
the observer’s screen is independent of the observer’s azimuthal viewing
angle. By contrast, when processing the 3D simulation, we effectively take a
snapshot of the simulation but allow for different azimuthal viewing angles
$\varphi_{\mathrm{obs}}$. This way, we can separately assess the effects of
time variability (inaccessible to the short 3D simulation) and nonaxisymmetry
of the current sheet (inaccessible to the 2D axisymmetric simulation) on the
image.

We only show the contribution of the current sheet. Indeed, synchrotron
emission by particles accelerated in the polar gaps is unlikely to be
observable. Only a small fraction of the Blandford-Znajek electromagnetic
power is dissipated in these gaps, which probably achieve plasma densities too
low to yield significant emission. Besides, unlike particles accelerated in
the current sheet, particles emitted in the polar regions have low pitch
angles, implying a comparatively low synchrotron emission. We have also run a
2D simulation with synchrotron cooling turned on, where the synchrotron
radiation-reaction force is rescaled, so that the characteristic cooling time
of an electron should be a few simulation time steps [32]. In that case, the
morphology of the image is practically unchanged whereas the polar cap
emission is quenched. This is consistent with previous PIC studies of magnetic
reconnection with synchrotron cooling [33, 34, 35]: particles are accelerated
near X points where the magnetic field vanishes, allowing particles to reach
high energies before escaping the reconnecting layer and cooling.

#### Averaged images

Fig. 2 shows images from our 2D simulation averaged over
$100\,r_{\mathrm{g}}/c$, which are very similar to the images from our 3D
simulation averaged over the azimuthal viewing angle. The leftmost panel
represents a view close to face-on, which is the most relevant to the M87*
system, assuming that the spin of the black hole is aligned with the large-
scale jet [36]. In general, the image is made up of two distinct components:
an inner and an outer ring. The rings’ brightness is higher on the left side
of the image, as a result of relativistic beaming: this is mainly a
geometrical effect due to the Kerr metric, rather than Doppler boosting. The
inner ring is the direct image of the current sheet, whereas the outer ring is
the lensed image. On each image, we have plotted the critical curve
corresponding to each inclination, i.e. the set of points on the screen hit by
photons which have orbited an arbitrarily large number of times around the
black hole close to bound, marginally stable spherical orbits [37, 38]. This
curve matches well the outer ring emission, confirming its interpretation as
the lensed image of the current sheet. However, the lensed image accounts for
only a small portion (around $10\%$) of the flux of the image. Our model hence
predicts that in a high-energy flaring state, most of the flux of the image
lies inside the critical curve. In this configuration, the black-hole shadow
prediction breaks down because most of the emission comes from within the
ergosphere and is not spherically distributed. We find that a major fraction
of the current sheet radiation originates from a very compact zone, within
$2-3\,r_{\mathrm{g}}$. Both rings become thinner as the observing frequency
increases, since the emission originates from the most energetic particles
which have been accelerated deeper within the current sheet.

Figure 3: Top row: Snapshot of the image from the 2D simulation (leftmost
panel) and differential maps illustrating the time evolution. The three maps
show the difference between successive snapshots at later times and the first
one. Bottom row: Image from the 3D simulation at a given azimuthal viewing
angle $\varphi_{\mathrm{obs}}$ (leftmost panel) and differential maps
illustrating the dependence on $\varphi_{\mathrm{obs}}$. The three maps show
the difference between images at higher $\varphi_{\mathrm{obs}}$ and the first
one.

#### Variability

The first row of Fig. 3 illustrates the time variability in the 2D
axisymmetric simulation, for a viewing angle
$\alpha_{\mathrm{obs}}=22.5^{\circ}$. Images are labeled using the observer’s
Boyer-Lindquist coordinate time. At a certain time during the simulation, a
bright circular feature appears between the inner and outer rings (top left
panel). As time elapses, this feature shrinks in size until it overlaps the
inner ring. Its typical lifetime is about $10\,r_{\mathrm{g}}/c$. By computing
the synchrotron emissivity of the plasma, we can relate this emission to a
giant magnetic island accretion event (see C21). The second row of Fig. 3
illustrates the dependence on $\varphi_{\mathrm{obs}}$ of images from the 3D
simulation, for the same viewing angle. Both the inner and outer rings show
nonaxisymmetric features in the form of brighter spots. They are more
pronounced on the inner ring, which is less impacted by gravitational lensing.
The presence of hot spots remains visible even after blurring the image at a
resolution $\simeq 5\,r_{\mathrm{g}}$. These hot spots can be traced back to
the different equatorial orbiting flux ropes. They radiate relatively
isotropically, so that every flux rope is visible regardless of
$\varphi_{\mathrm{obs}}$. This is why the positions of the hot spots along the
ring shift continuously clockwise with increasing $\varphi_{\mathrm{obs}}$ in
Fig. 3. The hot spots appear to move counterclockwise along the ring, with a
period $\sim 4\pi/\omega_{\mathrm{h}}\sim 14\,r_{\mathrm{g}}/c$, to a fixed
observer at a given $\varphi_{\mathrm{obs}}$.

Although the ratio between $r_{\mathrm{g}}$ and the electron plasma skin depth
is unrealistically small in these simulations, the hierarchical merger of
magnetic islands should proceed until plasmoids with macroscopic sizes are
formed. The final number of synchrotron-emitting flux ropes, and hence of hot
spots on the image, is determined by their characteristic escape time from the
merging region [39]. The flux ropes are constantly reformed within the current
sheet to compensate for their escape, along the separatrix or into the black
hole. Cerutti and Giacinti [39] predict the final number of macroscopic
plasmoids to be $\sim\pi/\beta_{\mathrm{rec}}$, independently of the initial
number of plasmoids in the current sheet and the separation of scales. This is
consistent with the number of hot spots we observe in the inner ring of the
image.

#### Observability in the radio band

The typical flux from the inner $10\,r_{\mathrm{g}}$ in M87* at the frequency
$\nu_{0}=$230\text{\,}\mathrm{GHz}$$ is $1\text{\,}\mathrm{Jy}$ [2],
translating into a radiated power $\nu F_{\nu}\approx
10^{40}\,$\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{s}}^{-1}$$.
The total jet power $L_{\mathrm{jet}}$ is likely to lie within the range
$10^{42}-10^{45}\,$\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{s}}^{-1}$$
[40, 41]. Assuming reconnection in the equatorial plane powers the emission,
an upper bound on the total dissipated energy is
$L_{\mathrm{dis}}\sim\beta_{\mathrm{rec}}L_{\mathrm{jet}}$, with
$\beta_{\mathrm{rec}}\approx 0.1$ the collisionless reconnection rate [42].
The synchrotron cooling time scale $\tau_{\mathrm{s}}$ is related to the
inverse Compton cooling time $\tau_{\mathrm{IC}}$ in the Thomson regime by
$\tau_{\mathrm{s}}/\tau_{\mathrm{IC}}=U_{\mathrm{IC}}/\left(B_{0}^{2}/8\pi\right)\approx
10^{-2}$, with $B_{0}\approx$10\text{\,}$-$100\text{\,}\mathrm{G}$$ the
typical magnetic field strength at horizon scale [4] and
$U_{\mathrm{IC}}\approx
10^{-2}\,$\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{cm}}^{-3}$$
the energy density of low-energy photons in the inner zone ($\sim
10\,r_{\mathrm{g}}$) [43, 44]. Since
$\tau_{\mathrm{s}}\ll\tau_{\mathrm{IC}}\ll r_{\mathrm{g}}/c$, most of the
magnetic energy dissipated by magnetic reconnection, initially carried by
nonthermal particles, is converted into synchrotron emission.

The spectrum of reconnection-accelerated particles in pair plasma with zero
guide field can largely be described as a power-law
$\mathop{}\mathopen{}\mathrm{d}{N}/\mathop{}\mathopen{}\mathrm{d}{\gamma}\propto\gamma^{-p}$
[20, 21, 45]. The index $p$ decreases with increasing upstream magnetization
$\sigma$, going from $p\approx 2$ for $\sigma\gtrsim 10$ to $1$ at $\sigma\gg
10$, in the absence of synchrotron cooling [45]. As $p$ gets closer to $1$,
more and more energy is carried by high-energy particles. The synchrotron
spectrum emitted by such a population of leptons also follows a power-law,
with a spectral energy density depending on frequency as $\nu
F_{\nu}\propto\nu^{(3-p)/2}$. The normalization is such that the total
radiated power matches $L_{\mathrm{dis}}$. The peak of the total synchrotron
spectrum cannot extend significantly beyond the burn-off limit at the photon
energy
$h\nu_{\mathrm{rad}}\approx$20\text{\,}$\left(\beta_{\mathrm{rec}}/0.1\right)$\text{\,}\mathrm{MeV}$$,
with $h$ the Planck constant ($\nu_{\mathrm{rad}}\sim 10^{10}\nu_{0}$). The
characteristic synchrotron frequency of the photons emitted by the highest-
energy particles $\nu_{\mathrm{s}}\sim\sigma^{2}eB_{0}/m_{\mathrm{e}}c$ lies
above $\nu_{\mathrm{rad}}$ for $\sigma\gtrsim 10^{6}$. The value of $\sigma$
in M87* is likely to be much larger than $1$, but its exact value is
uncertain.

If $\sigma\gtrsim\approx 10^{6}$, reconnection in the bare current sheet
occurs in the strong cooling regime. Assuming $p\lesssim 2$, the radiated
power $\nu F_{\nu}$ at the radio frequency
$\nu_{0}=$230\text{\,}\mathrm{GHz}$$ approximately reads
$L_{\mathrm{dis}}(\nu_{0}/\nu_{\mathrm{rad}})^{(3-p)/2}$, yielding

$\nu F_{\nu}\approx
10^{33}\left(\dfrac{L_{\mathrm{jet}}}{10^{44}\,$\text{\,}\mathrm{erg}\text{{\hbox{{}\\!\cdot\\!{}}}}{\mathrm{s}}^{-1}$}\right)\left(\dfrac{h\nu_{\mathrm{rad}}}{$20\text{\,}\mathrm{MeV}$}\dfrac{$230\text{\,}\mathrm{GHz}$}{\nu_{0}}\right)^{\textstyle\frac{p-1}{2}}\,$\text{\,}\mathrm{erg}\text{{\hbox{{}\\!\cdot\\!{}}}}{\mathrm{s}}^{-1}$,$
(1)

where we have neglected gravitational redshift. Here, the most energetic
photons have energies in the $\text{\,}\mathrm{MeV}$ range, so that they can
produce pairs in the upstream [35, 46]. We show in the Supplemental Material
that these secondary pairs also cool efficiently through synchrotron
radiation, and emit in the radio band. Their total luminosity is
${10}^{39}\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{s}}^{-1}$,
which is much larger than the estimate from Eq. (1), but is still outshone by
the accretion disk, though marginally. Hence, the radio image would be dimmer
in the flaring state.

On the other hand, $\sigma$ could be smaller than $10^{6}$ but still large
enough for reconnection to occur in the relativistic regime. Although the
multiplicities resulting from electromagnetic cascades in the polar spark gaps
are not very high [47, 23], additional mass loading from the accretion flow
could reduce $\sigma$ with respect to the values defined by the plasma
densities produced in these gaps. In that case, the particle spectrum retains
its non-radiative shape approximately up to
$\gamma\sim\mathrm{few}\times\sigma$. Beyond, the spectrum steepens, the exact
shape depending on the relative strengths of cooling and secondary
acceleration mechanisms. To give a simple estimate of the power at radio
frequencies, we assume that most of the power is carried by particles below
$\sigma$. The energetic constraint on the synchrotron radio flux is then
relaxed with respect to the $\sigma>10^{6}$ case: effectively,
$\nu_{\mathrm{rad}}$ must be replaced by
$\nu_{\mathrm{s}}\sim\sigma^{2}eB_{0}/m_{\mathrm{e}}c$ in Eq. (1). Fixing the
power-law index at $p_{0}=1$, a typical value in the case of $e^{\pm}$ non-
radiative reconnection up to $\sigma$, we obtain $\nu F_{\nu}\approx
10^{40}\left(L_{\mathrm{jet}}/10^{44}\,$\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{s}}^{-1}$\right)\left(\sigma/10^{3}\right)^{p_{0}-3}\,$\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{s}}^{-1}$.$
The radiation emitted by the current sheet in the high-energy flaring state
could then reach levels similar to the quiescent state and be observable by
the EHT. To obtain this estimate, we have neglected the contribution of
additional acceleration mechanisms, beyond the impulsive acceleration of
particles by the reconnected electric field, to the energy budget [48, 49,
50].

#### Conclusion

In this letter, we have provided a physically motivated, first-principles
model for the image of a low-luminosity AGN. The validity of our results is
limited to describing a high-energy flaring state, in which the accretion flow
has receded drastically, rather than the quiescent state observed by the EHT.
However, even if the current sheet that we obtain from our GRPIC simulations
results from idealized initial conditions, very similar configurations are
reached by GRMHD simulations, which also include accretion physics. As long as
the upstream magnetization is much larger than $1$, the morphology of the
image should only depend on the black-hole spin.

Since Sgr A* is expected to vary on time scales of several minutes, future
observations will be able to observe changes in the morphology of the image
within the EHT observation window. From our simulations, we expect the radius
of the bright ring to change with time, and possible hotspots to move along
the ring, thereby providing a test of our model. In M87*, these hot spots
should fully rotate along the ring in $5$ days, whereas in Sgr A* it would
take about $5$ minutes.

###### Acknowledgements.

The authors would like to thank the anonymous referees for insightful
comments, as well as Bart Ripperda and Fabio Bacchini for helpful discussions.
This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 863412). Computing resources were provided by TGCC and
CINES under the allocation A0090407669 made by GENCI, and by Princeton
Research Computing, a consortium of groups including the Princeton Institute
for Computational Science and Engineering (PICSciE) and the Office of
Information Technology’s High Performance Computing Center and Visualization
Laboratory at Princeton University.

Multi-dimensional simulations of ergospheric pair discharges around black
holes:
Supplemental Material

##  I Numerical techniques

###  I.1 Field equations and charge deposition

In GRZeltron, we solve the electromagnetic field equations derived by
Komissarov [27] in the $3+1$ formalism:

$\displaystyle\partial_{t}\boldsymbol{B}$ $\displaystyle=$
$\displaystyle-\boldsymbol{\nabla}\times\boldsymbol{E},$ (2)
$\displaystyle\partial_{t}\boldsymbol{D}$ $\displaystyle=$
$\displaystyle\boldsymbol{\nabla}\times\boldsymbol{H}-4\pi\boldsymbol{J},$ (3)

where $\boldsymbol{H}$ and $\boldsymbol{E}$ are auxiliary fields defined by

$\displaystyle\boldsymbol{H}$ $\displaystyle=$
$\displaystyle\alpha\boldsymbol{B}-\boldsymbol{\beta}\times\boldsymbol{D},$
(4) $\displaystyle\boldsymbol{E}$ $\displaystyle=$
$\displaystyle\alpha\boldsymbol{D}+\boldsymbol{\beta}\times\boldsymbol{B}.$
(5)

We have taken $c=1$ by convention. In these equations, $\alpha$ is the lapse
function and $\boldsymbol{\beta}$ is the shift vector associated with a
particular foliation of spacetime. We apply axial symmetry at the $\theta=0$
and $\theta=\pi$ boundaries: $D^{\varphi}=0$, $B^{\theta}=0$ and $\partial
D^{r}/\partial\theta=0$. We use outflowing boundary conditions at the outer
radial boundary, and zero-gradient conditions for $D^{\theta}$ and $B^{r}$ at
the inner radial boundary, within the event horizon.

The current density $\boldsymbol{J}$ acts as a source term for the
electromagnetic fields. It is obtained by depositing the current densities of
macro-particles on the grid. We use a first-order interpolation of the
particle coordinates. If a macroparticle with charge $q$ and weight $w$ is
located between the nodes $r\in[r_{i},r_{i+1}]$ and
$\theta\in[\theta_{j},\theta_{j+1}]$, then the charge deposited on each
neighboring node is

$\displaystyle q_{i,j}$ $\displaystyle=\dfrac{V_{i+1,j+1}}{V}\,q\,w,$
$\displaystyle q_{i,j+1}$ $\displaystyle=\dfrac{V_{i+1,j}}{V}\,q\,w,$ (6)
$\displaystyle q_{i+1,j}$ $\displaystyle=\dfrac{V_{i,j+1}}{V}\,q\,w,$
$\displaystyle q_{i+1,j+1}$ $\displaystyle=\dfrac{V_{i,j}}{V}\,q\,w,$ (7)

where we have defined $V=(r_{i+1}-r_{i})(\theta_{j+1}-\theta_{j})$, and
$V_{i,j}=\absolutevalue{(r-r_{i})(\theta-\theta_{j})}$.

###  I.2 Divergence cleaning

This choice of deposition is conceptually simple and robust, but does not
maintain the Maxwell-Gauss equation
$\boldsymbol{\nabla}\cdot\boldsymbol{D}=4\pi\rho$ to machine precision (with
$\rho$ the charge density deposited by the particle distribution). This lack
of charge conservation could disrupt the simulation on long time scales. To
alleviate this problem, we perform elliptic divergence cleaning every $N_{c}$
time steps, with $N_{c}=25$ in our case. We solve Poisson’s equation
$\Delta\Phi=-(4\pi\rho-\boldsymbol{\nabla}\cdot\boldsymbol{D})$ for the
correction $\Phi$ to the total electric potential (with $\boldsymbol{D}$ the
electric field before the cleaning), and then correct the electric field by
the amount $\delta\boldsymbol{D}=-\boldsymbol{\nabla}\Phi$. For example, a 2D
discretization of the general-relativistic Maxwell-Gauss equation reads

$\displaystyle\left(\Delta\phi\right)_{i,j}$
$\displaystyle=\dfrac{1}{\sqrt{h}_{i,j}}\left(C_{1}\left(\phi_{i+1,j}-\phi_{i,j}\right)+C_{2}\left(\phi_{i-1,j}-\phi_{i,j}\right)+C_{3}\left(\phi_{i,j+1}-\phi_{i,j}\right)+C_{4}\left(\phi_{i,j-1}-\phi_{i,j}\right)\right)$
(8)
$\displaystyle=\dfrac{\left(D^{r}\sqrt{h}\right)_{i+1/2,j}-\left(D^{r}\sqrt{h}\right)_{i-1/2,j}}{{\Delta
r}\sqrt{h}_{i,j}}+\dfrac{\left(D^{\theta}\sqrt{h}\right)_{i,j+1/2}-\left(D^{\theta}\sqrt{h}\right)_{i,j-1/2}}{{\Delta\theta}\sqrt{h}_{i,j}}-4\pi\rho_{i,j},$
(9)

where the weight coefficients area are

$\displaystyle C_{1}$
$\displaystyle=\dfrac{\sqrt{h}_{i+1/2,j}h^{rr}_{i+1/2,j}}{{\Delta
r}^{2}},\quad C_{2}=\dfrac{\sqrt{h}_{i-1/2,j}h^{rr}_{i-1/2,j}}{{\Delta
r}^{2}},$ (10) $\displaystyle C_{3}$
$\displaystyle=\dfrac{\sqrt{h}_{i,j+1/2}h^{\theta\theta}_{i,j+1/2}}{{\Delta\theta}^{2}},\quad
C_{4}=\dfrac{\sqrt{h}_{i,j-1/2}h^{\theta\theta}_{i,j-1/2}}{{\Delta\theta}^{2}}.$
(11)

##  II Ray-tracing

Let us consider a photon with $4$-momentum $p_{\mu}$. Geodesic motion in Kerr
spacetime is completely integrable, since three conserved quantities can be
attributed to each geodesic: the energy $E=-p_{t}$, the angular momentum
$L=p_{\varphi}$, and the Carter constant $Q$, defined as [31]

$Q={p_{\theta}}^{2}-\cos^{2}{\theta}\left(a^{2}\dfrac{E^{2}r_{\mathrm{g}}^{2}}{c^{2}}-\dfrac{L^{2}}{\sin^{2}{\theta}}\right).$
(12)

Null geodesics are independent of the photon energy $E$, so we define the
dimensionless parameters $\ell=cL/Er_{\mathrm{g}}$ and
$q^{2}=Qc^{2}/E^{2}r_{\mathrm{g}}^{2}$. If a photon reaches spatial infinity,
$E$ is its redshifted energy as measured by an observer at infinity. A photon
trajectory is fully determined by a given $\left(\ell,q^{2}\right)$,
supplemented by the initial signs of $\dot{r}$ and $\dot{\theta}$. Two
geodesics are associated with these parameters (determined by the choice of
one of these two signs), whereas the other choice of sign determines the
direction of travel along that geodesic.

Let us consider a distant observer equipped with a screen, the center of which
lies on the black hole. The screen is perpendicular to the line of sight of
the distant observer to the black hole, which has an inclination
$\alpha_{\mathrm{obs}}$ with respect to the spin axis. Let us assume that a
geodesic with parameters $(\ell,q^{2})$ reaches this screen at the position
$(X,Y)$. The impact parameter $Y$ is the apparent displacement of the photon
in the direction parallel to the spin axis of the black hole, whereas $X$ is
the apparent displacement in the direction perpendicular to the projected
axis. They are given by [37]

$\displaystyle\dfrac{X}{r_{\mathrm{g}}}$
$\displaystyle=-\dfrac{\ell}{\sin{\alpha_{\mathrm{obs}}}},$ (13)
$\displaystyle\dfrac{Y}{r_{\mathrm{g}}}$
$\displaystyle=\pm\sqrt{q^{2}+a^{2}\cos^{2}{\alpha_{\mathrm{obs}}}-\ell^{2}\,\mathrm{cotan}^{2}{\alpha_{\mathrm{obs}}}}.$
(14)

The sign of $Y$ depends on the initial sign of $\dot{\theta}$.

At every time step, each macro-particle can emit a macro-photon along its
direction of motion, by virtue of relativistic beaming. To simplify matters,
and because we do not claim to model realistic spectra given our limited scale
separation, we assume that the power spectrum of any macro-photon is
monoenergetic at the characteristic local synchrotron frequency
$\nu_{\mathrm{s}}$. The total distribution of synchrotron photons is
downsampled to maintain a reasonable computational cost. We assume that the
plasma is optically thin everywhere and neglect annihilation between two
synchrotron photons, so that once emitted they no longer interact with the
plasma. As a result, their information can be saved and post-processed. Rather
than the locally emitted frequency $\nu_{\mathrm{s}}$, we actually store the
redshifted frequency $\nu_{\infty}$ as measured by an observer at infinity. It
is computed as $h\nu_{\infty}=\alpha h\nu_{\mathrm{s}}-\beta^{i}p_{i}$, where
$h$ is the Planck constant, $p_{i}$ the photon $3$-momentum, $\alpha$ the
lapse function and $\beta^{i}$ the shift vector of the Kerr spacetime.

The screen has a resolution $200\times 200$ and a field of view
$8\,r_{\mathrm{g}}\times 8\,r_{\mathrm{g}}$. The angular resolution is
$9^{\circ}$ in the polar direction and $15^{\circ}$ in the azimuthal
direction. To reduce shot noise, we apply a median filter on all images with a
window of $3$ pixels. In all displayed images, we only select photons with
frequencies between $eB_{0}/m_{\mathrm{e}}c$ and $10\,eB_{0}/m_{\mathrm{e}}c$,
a band that lies roughly in the middle of the simulated synchrotron spectrum.
This allows us to focus on photons emitted by accelerated particles and to get
rid of the diffuse background emission, while retaining a pass-band wide
enough to have good statistics.

##  III Luminosity of secondary pairs if $\sigma>10^{6}$

We define $\gamma_{\mathrm{rad}}$ as the Lorentz factor of particles at which
the synchrotron drag matches the accelerating electric force.
$\gamma_{\mathrm{rad}}$ can therefore be expressed as
$\beta_{\mathrm{rec}}B_{0}=\sigma_{\mathrm{T}}U_{\mathrm{B}}\gamma_{\mathrm{rad}}^{2}$
(with $\sigma_{\mathrm{T}}$ the Thomson cross section and $U_{\mathrm{B}}$ the
magnetic energy density), yielding $\gamma_{\mathrm{rad}}\approx 10^{6}$ in
the case of M87*. In the regime where $\sigma>\gamma_{\mathrm{rad}}$, most of
the dissipated energy is radiated close to the broad peak at the synchrotron
burnoff limit at $\epsilon_{\gamma}\approx$20\text{\,}\mathrm{MeV}$$. Photons
with such energy can annihilate with each other and produce secondary pairs
with characteristic Lorentz factors
$\gamma_{\mathrm{s}}\sim\epsilon_{\gamma}/m_{\mathrm{e}}c^{2}\approx 100$. The
synchrotron cooling time $t_{\mathrm{syn}}$ of these secondary pairs is given
by

$\dfrac{ct_{\mathrm{syn}}}{r_{\mathrm{g}}}\sim\dfrac{m_{\mathrm{e}}c^{2}}{\gamma_{\mathrm{s}}\sigma_{\mathrm{T}}r_{\mathrm{g}}U_{B}},$
(15)

with $U_{B}$ the magnetic energy density. Even for secondary pairs, we have
$ct_{\mathrm{syn}}/r_{\mathrm{g}}\sim 10^{-3}\ll 1$, for a magnetic strength
of $100\text{\,}\mathrm{G}$. Consequently, secondary pairs radiate almost
instantly their energy before they have time to escape. Most of this energy is
emitted around
$h\gamma_{\mathrm{s}}^{2}eB_{0}/m_{\mathrm{e}}c\approx$10\text{\,}\mathrm{meV}$$,
so that these pairs contribute to the $\text{\,}\mathrm{mm}$ wavelengths
considered in this paper.

The optical depth to pair production for photons at energy $\epsilon_{\gamma}$
can be estimated as

$\tau_{\gamma\gamma}\sim\sigma_{\gamma\gamma}l\dfrac{U_{\mathrm{s}}\left(\epsilon_{\mathrm{s}}\right)}{\epsilon_{\mathrm{s}}},$
(16)

where $l\sim 10\,r_{\mathrm{g}}$ is the typical size where high-energy photons
are produced, $\sigma_{\gamma\gamma}\sim 0.2\,\sigma_{\mathrm{T}}$ is the peak
cross section for photon annihilation, and
$U_{\mathrm{IC}}\left(\epsilon_{\mathrm{s}}\right)$ is the energy density of
photons at the target energy
$\epsilon_{\mathrm{s}}=\left(m_{\mathrm{e}}c^{2}\right)/\epsilon_{\gamma}$
(with $\epsilon_{\mathrm{s}}\approx\epsilon_{\gamma}$ in the case of photon-
photon annihilation in the $\text{\,}\mathrm{MeV}$ range). Estimating
$U_{\mathrm{s}}\left(\epsilon_{\mathrm{s}}\right)=L_{\mathrm{dis}}/4\pi
l^{2}c$, we obtain $\tau_{\gamma\gamma}\sim
10^{-4}-10^{-5}\left(L_{\mathrm{jet}}/10^{44}\,$\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{s}}^{-1}$\right)\left(10\,r_{\mathrm{g}}/l\right)$
[19]. In this optically thin regime, the pair creation rate is then given by
$\dot{N}\simeq\tau_{\gamma\gamma}L_{\mathrm{dis}}/\epsilon_{\gamma}$. Because
the energy in the secondary pairs is quickly radiated away, the secondary
pairs synchrotron luminosity is given by
$L_{\mathrm{sec}}\sim\gamma_{\mathrm{s}}m_{\mathrm{e}}c^{2}\dot{N}\sim\tau_{\gamma\gamma}L_{\mathrm{dis}}$.
For
$L_{\mathrm{jet}}=10^{44}\,$\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{s}}^{-1}$$,
we obtain $L_{\mathrm{sec}}\approx
10^{39}-10^{40}\,$\text{\,}\mathrm{erg}\text{${}\\!\cdot\\!{}$}{\mathrm{s}}^{-1}$$.

## References

  * Blandford and Znajek [1977] R. D. Blandford and R. L. Znajek, MNRAS 179, 433 (1977).
  * Event Horizon Telescope Collaboration  _et al._ [2019a] Event Horizon Telescope Collaboration _et al._ , ApJ 875, L1 (2019a), arXiv:1906.11238 [astro-ph.GA] .
  * Event Horizon Telescope Collaboration  _et al._ [2021a] Event Horizon Telescope Collaboration _et al._ , ApJ 910, L12 (2021a).
  * Event Horizon Telescope Collaboration  _et al._ [2021b] Event Horizon Telescope Collaboration _et al._ , ApJ 910, L13 (2021b).
  * Mościbrodzka  _et al._ [2016] M. Mościbrodzka, H. Falcke, and H. Shiokawa, A&A 586, A38 (2016), arXiv:1510.07243 [astro-ph.HE] .
  * Event Horizon Telescope Collaboration  _et al._ [2019b] Event Horizon Telescope Collaboration _et al._ , ApJ 875, L5 (2019b), arXiv:1906.11242 [astro-ph.GA] .
  * Davelaar  _et al._ [2019] J. Davelaar, H. Olivares, O. Porth, T. Bronzwaer, M. Janssen, F. Roelofs, Y. Mizuno, C. M. Fromm, H. Falcke, and L. Rezzolla, A&A 632, A2 (2019), arXiv:1906.10065 [astro-ph.HE] .
  * Scepi  _et al._ [2021] N. Scepi, J. Dexter, and M. C. Begelman, arXiv e-prints , arXiv:2107.08056 (2021), arXiv:2107.08056 [astro-ph.HE] .
  * Fromm  _et al._ [2022] C. M. Fromm, A. Cruz-Osorio, Y. Mizuno, A. Nathanail, Z. Younsi, O. Porth, H. Olivares, J. Davelaar, H. Falcke, M. Kramer, and L. Rezzolla, A&A 660, A107 (2022), arXiv:2111.02518 [astro-ph.HE] .
  * Falcke  _et al._ [2000] H. Falcke, F. Melia, and E. Agol, ApJ 528, L13 (2000), arXiv:astro-ph/9912263 [astro-ph] .
  * Narayan  _et al._ [2019] R. Narayan, M. D. Johnson, and C. F. Gammie, ApJ 885, L33 (2019), arXiv:1910.02957 [astro-ph.HE] .
  * Bronzwaer  _et al._ [2021] T. Bronzwaer, J. Davelaar, Z. Younsi, M. Mościbrodzka, H. Olivares, Y. Mizuno, J. Vos, and H. Falcke, MNRAS 501, 4722 (2021), arXiv:2011.00069 [astro-ph.HE] .
  * Gralla  _et al._ [2019] S. E. Gralla, D. E. Holz, and R. M. Wald, Phys. Rev. D 100, 024018 (2019), arXiv:1906.00873 [astro-ph.HE] .
  * Nalewajko  _et al._ [2020] K. Nalewajko, M. Sikora, and A. Różańska, A&A 634, A38 (2020).
  * Vincent  _et al._ [2021] F. H. Vincent, M. Wielgus, M. A. Abramowicz, E. Gourgoulhon, J. P. Lasota, T. Paumard, and G. Perrin, A&A 646, A37 (2021), arXiv:2002.09226 [gr-qc] .
  * Lockhart and Gralla [2022] W. Lockhart and S. E. Gralla, MNRAS 509, 3643 (2022), arXiv:2107.06948 [astro-ph.HE] .
  * Ripperda  _et al._ [2020] B. Ripperda, F. Bacchini, and A. A. Philippov, ApJ 900, 100 (2020), arXiv:2003.04330 [astro-ph.HE] .
  * Chashkina  _et al._ [2021] A. Chashkina, O. Bromberg, and A. Levinson, MNRAS 508, 1241 (2021), arXiv:2106.15738 [astro-ph.HE] .
  * Ripperda  _et al._ [2022] B. Ripperda, M. Liska, K. Chatterjee, G. Musoke, A. A. Philippov, S. B. Markoff, A. Tchekhovskoy, and Z. Younsi, ApJ 924, L32 (2022), arXiv:2109.15115 [astro-ph.HE] .
  * Sironi and Spitkovsky [2014] L. Sironi and A. Spitkovsky, ApJ 783, L21 (2014), arXiv:1401.5471 [astro-ph.HE] .
  * Guo  _et al._ [2014] F. Guo, H. Li, W. Daughton, and Y.-H. Liu, Phys. Rev. Lett. 113, 155005 (2014), arXiv:1405.4040 [astro-ph.HE] .
  * Parfrey  _et al._ [2019] K. Parfrey, A. Philippov, and B. Cerutti, Phys. Rev. Lett. 122, 035101 (2019), arXiv:1810.03613 [astro-ph.HE] .
  * Crinquand  _et al._ [2020] B. Crinquand, B. Cerutti, A. Philippov, K. Parfrey, and G. Dubus, Phys. Rev. Lett. 124, 145101 (2020), arXiv:2003.03548 [astro-ph.HE] .
  * Crinquand  _et al._ [2021] B. Crinquand, B. Cerutti, G. Dubus, K. Parfrey, and A. Philippov, A&A 650, A163 (2021), arXiv:2012.09733 [astro-ph.HE] .
  * Bransgrove  _et al._ [2021] A. Bransgrove, B. Ripperda, and A. Philippov, Phys. Rev. Lett. 127, 055101 (2021), arXiv:2109.14620 [astro-ph.HE] .
  * El Mellah  _et al._ [2021] I. El Mellah, B. Cerutti, B. Crinquand, and K. Parfrey, arXiv e-prints , arXiv:2112.03933 (2021), arXiv:2112.03933 [astro-ph.HE] .
  * Komissarov [2004] S. S. Komissarov, MNRAS 350, 427 (2004).
  * Cerutti  _et al._ [2014] B. Cerutti, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, ApJ 782, 104 (2014), arXiv:1311.2605 [astro-ph.HE] .
  * Zenitani and Hoshino [2007] S. Zenitani and M. Hoshino, ApJ 670, 702 (2007), arXiv:0708.1000 [astro-ph] .
  * Dexter and Agol [2009] J. Dexter and E. Agol, ApJ 696, 1616 (2009), arXiv:0903.0620 [astro-ph.HE] .
  * Carter [1968] B. Carter, Physical Review 174, 1559 (1968).
  * Cerutti  _et al._ [2016] B. Cerutti, A. A. Philippov, and A. Spitkovsky, MNRAS 457, 2401 (2016), arXiv:1511.01785 [astro-ph.HE] .
  * Cerutti  _et al._ [2013] B. Cerutti, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, ApJ 770, 147 (2013), arXiv:1302.6247 [astro-ph.HE] .
  * Kagan  _et al._ [2016] D. Kagan, E. Nakar, and T. Piran, ApJ 833, 155 (2016), arXiv:1610.05219 [astro-ph.HE] .
  * Hakobyan  _et al._ [2019] H. Hakobyan, A. Philippov, and A. Spitkovsky, ApJ 877, 53 (2019), arXiv:1809.10772 [astro-ph.HE] .
  * Walker  _et al._ [2018] R. C. Walker, P. E. Hardee, F. B. Davies, C. Ly, and W. Junor, ApJ 855, 128 (2018), arXiv:1802.06166 [astro-ph.HE] .
  * Cunningham and Bardeen [1973] C. T. Cunningham and J. M. Bardeen, ApJ 183, 237 (1973).
  * Teo [2003] E. Teo, General Relativity and Gravitation 35, 1909 (2003).
  * Cerutti and Giacinti [2021] B. Cerutti and G. Giacinti, arXiv e-prints , arXiv:2111.04337 (2021), arXiv:2111.04337 [astro-ph.HE] .
  * Broderick  _et al._ [2015] A. E. Broderick, R. Narayan, J. Kormendy, E. S. Perlman, M. J. Rieke, and S. S. Doeleman, ApJ 805, 179 (2015), arXiv:1503.03873 [astro-ph.HE] .
  * Prieto  _et al._ [2016] M. A. Prieto, J. A. Fernández-Ontiveros, S. Markoff, D. Espada, and O. González-Martín, MNRAS 457, 3801 (2016), arXiv:1508.02302 [astro-ph.GA] .
  * Sironi  _et al._ [2016] L. Sironi, D. Giannios, and M. Petropoulou, MNRAS 462, 48 (2016), arXiv:1605.02071 [astro-ph.HE] .
  * Broderick and Tchekhovskoy [2015] A. E. Broderick and A. Tchekhovskoy, ApJ 809, 97 (2015), arXiv:1506.04754 [astro-ph.HE] .
  * EHT MWL Science Working Group  _et al._ [2021] EHT MWL Science Working Group _et al._ , ApJ 911, L11 (2021), arXiv:2104.06855 [astro-ph.HE] .
  * Werner  _et al._ [2016] G. R. Werner, D. A. Uzdensky, B. Cerutti, K. Nalewajko, and M. C. Begelman, ApJ 816, L8 (2016), arXiv:1409.8262 [astro-ph.HE] .
  * Mehlhaff  _et al._ [2021] J. M. Mehlhaff, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, MNRAS 508, 4532 (2021), arXiv:2107.01297 [astro-ph.HE] .
  * Chen and Yuan [2020] A. Y. Chen and Y. Yuan, ApJ 895, 121 (2020), arXiv:1908.06919 [astro-ph.HE] .
  * Petropoulou and Sironi [2018] M. Petropoulou and L. Sironi, MNRAS 481, 5687 (2018), arXiv:1808.00966 [astro-ph.HE] .
  * Hakobyan  _et al._ [2021] H. Hakobyan, M. Petropoulou, A. Spitkovsky, and L. Sironi, ApJ 912, 48 (2021), arXiv:2006.12530 [astro-ph.HE] .
  * Zhang  _et al._ [2021] H. Zhang, L. Sironi, and D. Giannios, ApJ 922, 261 (2021), arXiv:2105.00009 [astro-ph.HE] .