# Continual Learning-Based MIMO Channel Estimation: A Benchmarking Study
††thanks: This work was partially supported by the Discovery Grants Program of
the Natural Sciences and Engineering Research Council of Canada (NSERC).

Mohamed Akrout∗, Amal Feriani∗, Faouzi Bellili, Amine Mezghani, and Ekram
Hossain  Department of Electrical and Computer Engineering, University of
Manitoba, Winnipeg, Canada
{akroutm<EMAIL_ADDRESS>{Faouzi.Bellili, Amine.Mezghani,
<EMAIL_ADDRESS>
###### Abstract

With the proliferation of deep learning techniques for wireless communication,
several works have adopted learning-based approaches to solve the channel
estimation problem. While these methods are usually promoted for their
computational efficiency at inference time, their use is restricted to
specific stationary training settings in terms of communication system
parameters, e.g., signal-to-noise ratio (SNR) and coherence time. Therefore,
the performance of these learning-based solutions will degrade when the models
are tested on different settings than the ones used for training. This
motivates our work in which we investigate continual supervised learning (CL)
to mitigate the shortcomings of the current approaches. In particular, we
design a set of channel estimation tasks wherein we vary different parameters
of the channel model. We focus on Gauss-Markov Rayleigh fading channel
estimation to assess the impact of non-stationarity on performance in terms of
the mean square error (MSE) criterion. We study a selection of state-of-the-
art CL methods and we showcase empirically the importance of catastrophic
forgetting in continuously evolving channel settings. Our results demonstrate
that the CL algorithms can improve the interference performance in two channel
estimation tasks governed by changes in the SNR level and coherence time.

###### Index Terms:

Continual learning, task transferability, channel estimation, MIMO, Gauss-
Markov channel model.

**footnotetext: equal contribution, order determined randomly.

##  I Introduction

###  I-A Background and motivation

Deep neural networks (DNN) are widely applied to solve wireless communication
problems, such as the channel estimation problem, where they played important
roles in improving the system performance and reducing computational
complexity [1, 2]. However, these learning-based methods are usually trained
and evaluated on _fixed datasets_ and/or _stationary propagation environments_
[3]. Indeed, critical parameters of communication system models such as the
signal-to-noise ratio (SNR) and coherence time ($T_{c}$) are often assumed to
be fixed for these learning-based methods. Since wireless communications
systems are _inherently non-stationary_ , trained models need to be tested on
continuously and dynamically changing datasets. This shift in data
distribution results in the degradation or even the failure of the trained
DNNs, and this problem is known in the literature as _catastrophic forgetting_
[4]. One solution to avoid catastrophic forgetting is to train an ensemble or
a pool of models for different training scenarios (a.k.a., tasks) and use the
appropriate model at inference time. Such approaches are inefficient, because
$i)$ they do not scale well since the number of models will increase as a
function of the number of tasks, and $ii)$ the appropriate model selection at
inference time is challenging. Thus, there is an urge to design solutions that
can learn incrementally and adapt to changes in the environment.

Continual learning (CL) methods [5, 6, 7] were introduced to train DNNs on
sequentially changing tasks without forgetting the previously learned ones. CL
training strategies sidestep the impact of averaging conflicting gradients on
the learning dynamics of DNNs [6]. As one example, consider training a DNN
sequentially on two different tasks $\mathcal{T}_{\textrm{SNR,high}}$ then
$\mathcal{T}_{\textrm{SNR,low}}$ at high and low SNR levels, respectively.
Once the training on the first task $\mathcal{T}_{\textrm{SNR,high}}$ is done,
an equilibrium in the space of DNN parameters is reached ensuring good
performance on the task $\mathcal{T}_{\textrm{SNR,high}}$. The subsequent
training on the second task $\mathcal{T}_{\textrm{SNR,low}}$, however, will
pull all the DNN parameters toward the direction of its gradients, thereby
leading to forgetting the first task $\mathcal{T}_{\textrm{SNR,high}}$. For
this reason, gradient-based training algorithms require that, in expectation,
all tasks are present to keep the conflicting tension between gradients
stemming from multiple tasks [6]. With CL techniques, the DNN is expected to
adapt to task changes by not forgetting the learned skill used in the previous
task.

There are different CL scenarios depending on the level of difficulty [7]. The
easiest one, called task-incremental learning (Task-IL), assumes the identity
of the task to be always provided which makes task-specific training possible.
The typical approach for this scenario is to train a “multi-headed” DNN where
each task has its output units and the rest of the network is shared between
tasks. The second scenario is domain-incremental learning (Domain-IL) where
the task identity is not known at inference time. A common example in this
scenario is where the input distribution is changing but the nature of the
task remains the same.

In this paper, we formulate the wireless channel estimation problem in a
multiple-input multiple-output (MIMO) wireless system as a domain-IL problem
for two main reasons: $i)$ the distribution shift in the DNN inputs due to
changes in the propagation environment properties, and $ii)$ the DNN is
expected to seamlessly estimate the channel, i.e., without knowing when the
environment has changed.

###  I-B Prior Work

####  I-B1 Overview of continual learning methods

Continual learning is an active research field within the machine learning
community [8]. To tackle catastrophic forgetting, CL-based methods can be
divided into three categories based on how task information is stored and used
throughout the sequential learning process (see [5] for a comprehensive
review):

Replay methods. These methods store a subset of the training task samples in a
memory buffer and revisit them while learning a new task. The task samples are
either reused as DNN inputs for rehearsal [9, 10, 11], or to constrain the
optimization of the new task loss to prevent previous task interference [12,
13].

Regularization-based methods. Unlike replay methods, regularization-based
methods do not store task samples to respect the data privacy, which also
eases the memory requirements [14, 15, 16]. Alternatively, they introduce an
extra regularization term in the loss function of the DNN, thereby
consolidating previous knowledge when learning on new tasks.

Parameter isolation methods. This family of methods append new model
parameters to each task. Therefore, DNNs can grow new layer branches dedicated
for the training for new tasks, while freezing [17] or masking out [18, 19]
branches trained on previous tasks.

For real-world problems, the aforementioned CL categories should satisfy the
following three constraints [5]:

  * •

No task boundaries: tasks do not arrive sequentially and are not delineated by
clear boundaries during training,

  * •

No task identifiers: the inference step does not identify the task difficulty
to the DNN using task identifiers,

  * •

Constant memory: the task dataset collected during the training should fit a
buffer with a bounded memory footprint.

####  I-B2 Continual learning for wireless communications

Little effort has been devoted to investigate CL methods in wireless
communications problems. A CL formulation for wireless system design was
adopted in [20] to optimize wireless resource allocation in a dynamic
environment. In this work, the DNN makes use of a continuously updated data
buffer containing representative historical data from different dynamic
settings of the environment. Another line of work proposed a CL-based minimum
mean-square error (CL-MMSE) technique for mmWave channel estimation to address
high frequency signal characteristics [21]. However, CL-MMSE ignores pilot
transmission overhead and only varies the number of receive antennas between 8
and 128 to generate the CL tasks.

###  I-C Contributions

Different from the aforementioned work, we present an exhaustive benchmark of
CL approaches for MIMO channel estimation. The contributions of this paper are
as follows:

  * •

We formulate the channel estimation problem as a domain-IL problem and propose
two different benchmarks based on the coherence time and SNR considerations;

  * •

We conduct an exhaustive empirical comparison of seven CL methods with a
detailed analysis of their performance;

  * •

For the first time, we provide a benchmark for CL for wireless communications.

##  II System Model and Assumptions

We consider a MIMO communication system under Gauss-Markov Rayleigh fading
channels between $N$ transmit antennas and $M$ receive antennas. The
transmission lasts for $B$ time blocks each of which is used to transmit $K$
symbols. The channel estimation method is illustrated in Fig. 1 where the
block-fading channel is assumed to remain constant over the entire $b$-th time
block ($b\in\\{1,\dots,B\\}$), while the total duration of $KB$ symbol periods
is used for channel estimation within the overestimated channel coherence time
$T_{c}$.

Figure 1:  An illustration of the considered channel estimation method. The
training time is divided into $B$ time-blocks of duration $K$ symbol periods
each. The channel is constant during the $b$-th time-block and varies between
different time-blocks while the same pilot symbols $\mathbf{x}[1]$, . . . ,
$\mathbf{x}[K]$ are reused from one block to another.

For every time-block $b$, the received signal is given by:

$\bm{y}_{b}=\bm{H}_{b}\,\bm{x}_{b}+\bm{n}_{b},\,\quad\quad b=1,\dots,B.$ (1)

In (1), $\bm{x}_{b}\in\mathbb{C}^{N}$ represents the transmitted symbol vector
and $\bm{n}_{b}\in\mathbb{C}^{M}$ is an additive complex white Gaussian noise
vector, i.e.,
$\bm{n}_{b}\sim\mathcal{CN}(\mathbf{0}_{M},\sigma^{2}_{\bm{n}_{b}}\,\mathbf{I}_{M})$.
Here, at a given SNR level $\rho$ [dB], the noise variance
$\sigma^{2}_{\bm{n}_{b}}$ is defined as:

$\sigma^{2}_{\bm{n}_{b}}=\mathbb{E}\Big{[}\left\lVert\bm{H}_{b}\,\bm{x}_{b}\right\rVert_{2}^{2}\Big{]}\,10^{-\rho/10}.$
(2)

Moreover, $\bm{H}_{b}$ is the time-varying Rayleigh fading channel evolving
according to the Gauss-Markov model [22]:

$\bm{H}_{b}=\sqrt{\alpha}\,\bm{H}_{b-1}+\sqrt{1-\alpha}\,\bm{W}_{b},\,\quad\quad
b=2,\dots,B.$ (3)

where $\alpha\in[0,1[$ is the channel memory factor. Moreover, $\bm{H}_{1}$
and $\bm{W}_{b}$ are independent and identically distributed realizations
drawn from the complex Gaussian distribution
$\mathcal{CN}(\mathbf{0}_{NM},\mathbf{I}_{NM})$. By recursively injecting the
expression of $\bm{H}_{b-1}$ in (3) for $b=1,\dots,B$, we obtain:

$\bm{H}_{b}=\sqrt{\alpha}^{b-1}\,\bm{H}_{1}+\sqrt{1-\alpha}\,\sum\limits_{p=2}^{b}\,\sqrt{\alpha}^{b-p}\,\bm{W}_{p}.$
(4)

Therefore, the sequence $\\{\bm{H}_{b}\\}_{b}$ in (4) forms a discrete-time
Markov chain that evolves as a random walk process starting from the initial
channel $\bm{H}_{1}$ scaled by the $b$-th block-dependent factor
$\sqrt{\alpha}^{b-1}$. This Markov chain models the stochastic channel aging
which induces significant shifts in channel distribution. The channel aging
phenomenon results from the variation of the channel between the instant when
it is learned using pilots and the instant when it is used for signal
processing. This time variation is due to user mobility and processing delay.

(a) varying $T_{c}$ at $\rho=10\,\textrm{dB}$

(b) varying $\rho$ at $T_{c}=20~{}\textrm{symbol periods}$

Figure 2: Illustration of the distribution shift in the received signal due to
changes in (a) the time coherence $T_{c}$ and (b) the SNR level $\rho$.

##  III Continual Learning-Based Channel Estimation

###  III-A Motivating Example

As mentioned before, the wireless channel is non-stationary, thus in this
example, we will study the impact of the non-stationarity of the channel
$\bm{H}_{b}$ on the distribution of the received signal $\bm{y}_{b}$ given in
(1). This is because $\bm{y}_{b}$ is the input vector to the neural network
for channel estimation. In this context, we consider two main changing factors
impacting the received signal $\bm{y}_{b}$: $i)$ the SNR level $\rho$ in (2)
and the coherence time $T_{c}=K$ as depicted in Fig. 1.

For better illustration, we consider the MIMO communication system described
in Section II with $N=8$ and $M=5$. We then examine the following two
scenarios:

  * •

In Fig. 2(a), we fix the SNR level to $\rho=10$ dB and vary the coherence time
$T_{c}$ between $10$ and $250$ transmit symbol periods. We observe a positive
shift (i.e., translation) in the empirical density of the received signal as
the coherence time increases***Here the density has higher values than 1 which
does not contradict basic probability theory because the density corresponds
to the probability per unit value (here $\sim 10^{-3}$) of the random variable
representing the received power..

  * •

In Fig. 2(b), we fix the coherence time $T_{c}$ to $20$ symbol periods and
vary the SNR between $0$ and $25$ dB. It is seen that the empirical density of
the received power* ‣ • ‣ III-A is unimodal at lower SNR levels (i.e., 0, 5
and 10 dB) and becomes multimodal at higher SNR levels (i.e., 15 and 20 dB).

This example showcases that the change in the characteristics of the wireless
channel will lead to a shift in the distribution of the received signal. This
motivates the design of channel estimation algorithms based on CL to avoid
catastrophic forgetting resulting from the non-stationarity of wireless
channels. In the next section, we explain the design of these CL tasks for
channel estimation.

###  III-B Task protocols

To comprehensively compare the performance of various CL approaches for
wireless channel estimation, we consider two sets of domain-IL tasks,
$\mathcal{T}_{\textrm{SNR}}$ and $\mathcal{T}_{T_{c}}$, defined as follows:

  * •

$\mathcal{T}_{\textrm{SNR}}\triangleq\bigcup\limits_{i}\big{\\{}\mathcal{T}_{\textrm{SNR}=i[\textrm{dB}]}\big{\\}}$
where different channel estimation tasks are generated by varying the SNR
level $\rho$. In this work, we assume the SNR ranges from $0\,\textrm{dB}$ to
$20\,\textrm{dB}$ with a step of $2$ dB, resulting in $11$ tasks in total;

  * •

$\mathcal{T}_{T_{c}}\triangleq\bigcup\limits_{j}\big{\\{}\mathcal{T}_{T_{c}=j}\big{\\}}$
where a sequence of $20$ tasks are generated by varying the index $j$ between
$5$ and $100$ symbol periods with a step size of $5$ symbol periods.

Note that we consider these tasks as domain-IL since an _implicit_ task-
dependent transformation will be applied on the inputs of the channel
estimation model (i.e., received signal). However, it is important to
highlight that this setting is slightly different from the common assumption
in the CL literature where domain-IL tasks are defined using _direct_
transformation on the inputs (e.g., rotation, permutation).

Furthermore, we scrutinize how changing the task order affects the channel
estimation. We do so by considering two different task orderings of the sets
$\mathcal{T}_{\textrm{SNR}}$ and $\mathcal{T}_{T_{c}}$:

  * •

Curriculum: the tasks are ordered from easy to hard which imitates the
meaningful learning order in human curricula. For instance, tasks are
presented to the DNN from high to low SNR and $T_{c}$ values;

  * •

Random: the tasks are presented to the DNN in a random order.

The task ordering hypothesis suggests that starting from the curriculum order
of tasks would lead to a better result as compared to the random order of
tasks, and some CL methods have shown little sensitivity to the task order in
the context of image classification [5]. Whether the same behavior is observed
for channel estimation problems will be discussed in Section V.

##  IV Experimental Setup

###  IV-A Benchmarking Methods

We consider the following CL methods:

  * •

Synaptic Intelligence (SI) [23]: this is a regularization-based method in
which a regularization term is added to penalize changes in parameters when
they move away from their values after training on previous tasks;

  * •

Experience Replay (ER) [9, 24, 25]: this replay-based method concatenates the
training samples from the current task with old samples from previous tasks to
construct training batches;

  * •

Averaged-GEM (A-GEM) [26]: this is an lightweight version of the Gradient
Episodic Memory (GEM) [27] that uses an episodic memory to leverage old
training data and introduce constraints to maintain similar predictions on
previous tasks while learning the new one. This method uses task boundaries to
store the task samples in the buffer;

  * •

Averaged-GEM with a reservoir buffer (A-GEM-R) [28]: this method is similar to
A-GEM. The main difference is in how the replay buffer is populated. This
method uses reservoir memory buffer to avoid using task boundaries;

  * •

Function Distance Regularization (FDR) [29]: this is another replay-based
method in which the network outputs and training samples are stored to align
previous and current outputs using an $L_{2}$-regularization loss. FDR relies
on task boundaries to populate the replay buffer;

  * •

Dark Experience Replay (DER) [28]: this is also a replay-based method in which
the outputs of the model’s final layer are stored instead of the ground truth
for the previously learned tasks. This method adopts reservoir sampling to
avoid using task boundaries. A penalty term is added to the loss to minimize
the difference between the current and the previous model outputs;

  * •

Dark Experience Replay++ (DER++) [28]: this method is an extension of the DER
algorithm where, in addition to the model outputs, the ground truth of
previous tasks is stored and another penalty term is added to reduce the
difference between the output of the model on the buffer datapoints and their
ground truth labels.

We further include the following two baselines:

  * •

No-CL: the DNN is trained sequentially on all the tasks without taking into
consideration the catastrophic forgetting. This can be seen as an upper bound
on the MSE;

  * •

Joint training: the DNN is trained on the joint data from all the tasks, which
leads to a lower bound on the MSE.

The implementation of the different CL methods is based on the one published
in [30]. For a broader reproducible benchmark of CL methods in wireless
communication, we will publish on Github the SNR and $T_{c}$ tasks as well as
instructions needed to reproduce the main experimental results presented in
Section V.

###  IV-B Training Details

####  IV-B1 Hyperparameters

We simulate a MIMO system as described in Section II with $N=8$ and $M=5$. We
set the carrier frequency $f_{c}$ to $3$ GHz, the path-loss exponent to $2.5$,
the channel memory factor $\alpha$ to $0.01$. and the distance between the
transmitter and the receiver to $20$ meters. We also consider a 2-PAM pilot
sequence of $10$ symbols.

####  IV-B2 Dataset generation

A dataset of $20000$ samples is generated following the system dynamics
described in Section II. The generated data is shuffled and divided such that
$90\%$ of the generated data is used for training and the remaining is
reserved for testing.

####  IV-B3 Architecture and training

For all the tasks, we use a fully-connected neural network with five hidden
layers of $3\times M$ ReLU units followed by an output layer with $N$ units.
All models are trained using stochastic gradient descent for a fair comparison
between all CL methods. For each task, the model is trained for $100$ epochs
to reach a satisfactory performance. Due to computational complexity, we avoid
the online scenario (i.e., train the model for one epoch) since it will not be
possible to distinguish between the effects of forgetting and under-fitting
[30].

###  IV-C Evaluation Metrics

While monitoring the MSE performance of the model across tasks, we assess its
ability to transfer knowledge by measuring:

  * •

The backward transfer which quantifies the influence of learning a task
$t_{i}$ on the performance of a previous task $t_{k}$ for $k<i$. A positive
(resp. negative) backward transfer is observed when learning a particular task
$t_{i}$ increases (resp. decreases) the performance on some preceding task
$t_{k}$. Large negative backward transfer is also known as catastrophic
forgetting;

  * •

The forward transfer which describes the impact of learning a task $t_{i}$ on
the performance of future tasks $t_{k}$ for $k>i$. Positive forward transfer
is observed when the model is able to perform “zero-shot” learning by
exploiting the feature correlation across tasks.

These two transfer measures are the common evaluation metrics adopted in the
CL literature [31, 30, 5]. More specifically, the model is evaluated on all
previously learned tasks after it has finished learning a task $t_{i}$. Hence,
after learning all the $T$ tasks, we get an error matrix
$E\in\mathbb{R}^{T\times T}$, where each element $e_{i,j}$ is the MSE error of
the task $t_{j}$ after learning the task $t_{i}$ (e.g., $e_{T,i}$ stands for
the error on the task $i$ after learning all the $T$ tasks). Further, we
denote by $\bar{b}$ the error vector of a randomly-initialized model. The
considered evaluation metrics are defined as follows [30]:

$\displaystyle\textbf{Average MSE}~{}(\texttt{AVG\\_MSE})$
$\displaystyle=\frac{1}{T}\sum_{i=1}^{T}e_{T,i},$ (5a)
$\displaystyle\textbf{Forward transfer}~{}(\texttt{FWT})$
$\displaystyle=\frac{1}{T-1}\sum_{i=2}^{T}\left(\bar{b}_{i}-e_{i-1,i}\right),$
(5b) $\displaystyle\textbf{Backward transfer}~{}(\texttt{BWT})$
$\displaystyle=\frac{1}{T-1}\sum_{i=1}^{T-1}\left(e_{T,i}-e_{i,i}\right).$
(5c)

##  V Numerical Results and Analysis

(a) $\mathcal{T}_{\textrm{SNR}}$ – Curriculum

(b) $\mathcal{T}_{\textrm{SNR}}$ – Random

(c) $\mathcal{T}_{T_{c}}$ – Curriculum

(d) $\mathcal{T}_{T_{c}}$ – Random

(e) $\mathcal{T}_{\textrm{SNR}}$ – Curriculum

(f) $\mathcal{T}_{\textrm{SNR}}$ – Random

(g) $\mathcal{T}_{T_{c}}$ – Curriculum

(h) $\mathcal{T}_{T_{c}}$ – Random

Figure 3: Top: The final MSE after learning all the tasks; Bottom: the
evolution of the MSE as more tasks are learned

TABLE I: Forward and backward transfer for all tasks Task | Metric | CL Methods
---|---|---
NO-CL | SI | A-GEM | A-GEM-R | ER | FDR | DER | DER++
$\mathcal{T}_{\text{SNR}}$ – Curriculum |  FWT $\Uparrow$ | 0.14 | 0.16 | 0.15 | 0.14 | 0.16 | 0.18 | 0.2 | 0.2
BWT $\Downarrow$ | 0.2 | 0.15 | 0.15 | 0.16 | 0.1 | 0.12 | 0.09 | 0.07
$\mathcal{T}_{\text{SNR}}$ – Random |  FWT $\Uparrow$ | 0.12 | 0.15 | 0.14 | 0.14 | 0.16 | 0.18 | 0.2 | 0.19
BWT $\Downarrow$ | 0.25 | 0.18 | 0.18 | 0.23 | 0.13 | 0.14 | 0.11 | 0.08
$\mathcal{T}_{T_{C}}$ – Curriculum |  FWT $\Uparrow$ | 0.12 | 0.14 | 0.13 | 0.14 | 0.16 | 0.17 | 0.2 | 0.19
BWT $\Downarrow$ | 0.15 | 0.12 | 0.14 | 0.14 | 0.12 | 0.11 | 0.09 | 0.07
$\mathcal{T}_{T_{C}}$ – Random |  FWT $\Uparrow$ | 0.13 | 0.17 | 0.14 | 0.16 | 0.19 | 0.19 | 0.22 | 0.2
BWT $\Downarrow$ | 0.23 | 0.16 | 0.21 | 0.21 | 0.13 | 0.13 | 0.1 | 0.08

All the results reported in this section are averaged across 5 runs, each one
involving a different random seed. We benchmark the performance of all the
algorithms described in Section IV-A with respect to the metrics detailed in
(5). We set the buffer size for all the replay-based methods to $200$.

###  V-A Estimation Accuracy

Figs. 3(a)-(d) depict the performance in terms of average MSE at the end of
all tasks while Figures 3(e)-(h) depict the evolution of the MSE. First, we
observe that the baseline NO-CL achieves the lowest MSE as expected which
emphasizes the inevitable impact of catastrophic forgetting when continual
learning is not taken into account. Besides, training on all tasks
simultaneously (i.e., joint training method) leads to a lower bound on the MSE
performance. Overall, the experience replay methods (i.e., ER, DER, DER++)
outperform the other CL algorithms. It is also seen that they achieve the
lowest MSE performance on both SNR and time coherence tasks. Moreover, we
observe that the MSE values are lower in the curriculum order compared to the
random one for both SNR and time coherence tasks. These results confirm the
task ordering hypothesis as defined in Section III-B which assumes that
forgetting becomes less severe when the curriculum order is adopted.

###  V-B Task Transferability

We present the results of the forward and backward transfer values in Table I.
It is observed that DER and DER++ achieve the best (i.e., highest) forward and
(i.e., lowest) backward performance on all SNR and time coherence tasks with
both curriculum and random task orders. Again, the random task ordering harms
both forward and backward transfers as compared to the curriculum ordering.
Furthermore, it is interesting to note that the gap between random and
curriculum orders for SNR tasks is smaller than the one for time coherence
tasks. This is justified by the following two arguments. First, the ability to
continuously learn is highly dependent on the number of total tasks. Hence,
the fact that the number of time coherence tasks is higher than the number of
SNR tasks makes learning with the former more challenging. Second, the SNR
step size (2 dB) is more granular than the time coherence step size (5 symbol
periods). This means that the obtained SNR tasks have a smoother increase in
task difficulty (or equivalently a high correlation across tasks) as compared
to the time coherence tasks, thereby improving the overall task
transferability.

##  VI Conclusion

We have benchmarked the state-of-the-art continual learning methods on the
Gauss-Markov Rayleigh fading channel estimation problem. After motivating how
the SNR level and the time coherence lead to a significant shift in the
distribution of the received power, we have varied these two parameters to
define SNR and time coherence tasks. We have also examined the learning of
these two tasks in both curriculum and random orders. The simulation results
suggest that experience replay methods outperform the other continual learning
approaches. This benchmark can be further extended to more diversified
wireless communications tasks by varying the communication distance and/or the
path-loss exponent, to determine whether or not the same experimental
conclusions hold for them.

## References

  * [1] A. Le Ha, T. Van Chien, T. H. Nguyen, W. Choi _et al._ , “Deep learning-aided 5g channel estimation,” in _2021 15th international conference on ubiquitous information management and communication (IMCOM)_. IEEE, 2021, pp. 1–7.
  * [2] Q. Hu, F. Gao, H. Zhang, S. Jin, and G. Y. Li, “Deep learning for channel estimation: Interpretation, performance, and comparison,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 4, pp. 2398–2412, 2020\.
  * [3] Y. Shen, Y. Shi, J. Zhang, and K. B. Letaief, “Lorm: Learning to optimize for resource management in wireless networks with few training samples,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 1, pp. 665–679, 2019.
  * [4] R. M. French, “Catastrophic forgetting in connectionist networks,” _Trends in cognitive sciences_ , vol. 3, no. 4, pp. 128–135, 1999.
  * [5] M. D. Lange, R. Aljundi, M. Masana, S. Parisot, X. Jia, A. Leonardis, G. G. Slabaugh, and T. Tuytelaars, “A continual learning survey: Defying forgetting in classification tasks,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 44, no. 7, pp. 3366–3385, 2022.
  * [6] R. Hadsell, D. Rao, A. A. Rusu, and R. Pascanu, “Embracing change: Continual learning in deep neural networks,” _Trends in cognitive sciences_ , vol. 24, no. 12, pp. 1028–1040, 2020.
  * [7] G. M. van de Ven and A. S. Tolias, “Three scenarios for continual learning,” _CoRR_ , vol. abs/1904.07734, 2019.
  * [8] S. Dohare, A. R. Mahmood, and R. S. Sutton, “Continual backprop: Stochastic gradient descent with persistent randomness,” _arXiv preprint arXiv:2108.06325_ , 2021.
  * [9] S.-A. Rebuffi, A. Kolesnikov, G. Sperl, and C. H. Lampert, “icarl: Incremental classifier and representation learning,” in _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , 2017, pp. 2001–2010.
  * [10] D. Rolnick, A. Ahuja, J. Schwarz, T. Lillicrap, and G. Wayne, “Experience replay for continual learning,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [11] D. Isele and A. Cosgun, “Selective experience replay for lifelong learning,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 32, no. 1, 2018.
  * [12] D. Lopez-Paz and M. Ranzato, “Gradient episodic memory for continual learning,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [13] A. Chaudhry, M. Ranzato, M. Rohrbach, and M. Elhoseiny, “Efficient lifelong learning with a-gem,” _arXiv preprint arXiv:1812.00420_ , 2018.
  * [14] Z. Li and D. Hoiem, “Learning without forgetting,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 40, no. 12, pp. 2935–2947, 2017\.
  * [15] J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A. A. Rusu, K. Milan, J. Quan, T. Ramalho, A. Grabska-Barwinska _et al._ , “Overcoming catastrophic forgetting in neural networks,” _Proceedings of the national academy of sciences_ , vol. 114, no. 13, pp. 3521–3526, 2017.
  * [16] C. V. Nguyen, Y. Li, T. D. Bui, and R. E. Turner, “Variational continual learning,” _arXiv preprint arXiv:1710.10628_ , 2017.
  * [17] J. Xu and Z. Zhu, “Reinforced continual learning,” _Advances in Neural Information Processing Systems_ , vol. 31, 2018.
  * [18] A. Mallya and S. Lazebnik, “Packnet: Adding multiple tasks to a single network by iterative pruning,” in _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , 2018, pp. 7765–7773.
  * [19] J. Serra, D. Suris, M. Miron, and A. Karatzoglou, “Overcoming catastrophic forgetting with hard attention to the task,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 4548–4557.
  * [20] H. Sun, W. Pu, X. Fu, T.-H. Chang, and M. Hong, “Learning to continuously optimize wireless resource in a dynamic environment: A bilevel optimization perspective,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 1900–1917, 2022.
  * [21] S. Kumar, S. K. Vankayala, B. S. Sahoo, and S. Yoon, “Continual learning-based channel estimation for 5g millimeter-wave systems,” in _2021 IEEE 18th Annual Consumer Communications & Networking Conference (CCNC)_. IEEE, 2021, pp. 1–6.
  * [22] W. C. Jakes and D. C. Cox, _Microwave mobile communications_. Wiley-IEEE press, 1994.
  * [23] F. Zenke, B. Poole, and S. Ganguli, “Continual learning through synaptic intelligence,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 3987–3995.
  * [24] M. Riemer, I. Cases, R. Ajemian, M. Liu, I. Rish, Y. Tu, and G. Tesauro, “Learning to learn without forgetting by maximizing transfer and minimizing interference,” _arXiv preprint arXiv:1810.11910_ , 2018.
  * [25] A. Robins, “Catastrophic forgetting, rehearsal and pseudorehearsal,” _Connection Science_ , vol. 7, no. 2, pp. 123–146, 1995.
  * [26] A. Chaudhry, M. Ranzato, M. Rohrbach, and M. Elhoseiny, “Efficient lifelong learning with a-gem,” _arXiv preprint arXiv:1812.00420_ , 2018.
  * [27] D. Lopez-Paz and M. Ranzato, “Gradient episodic memory for continual learning,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [28] P. Buzzega, M. Boschini, A. Porrello, D. Abati, and S. Calderara, “Dark experience for general continual learning: a strong, simple baseline,” _Advances in neural information processing systems_ , vol. 33, pp. 15 920–15 930, 2020.
  * [29] A. S. Benjamin, D. Rolnick, and K. Kording, “Measuring and regularizing networks in function space,” _arXiv preprint arXiv:1805.08289_ , 2018.
  * [30] P. Buzzega, M. Boschini, A. Porrello, D. Abati, and S. Calderara, “Dark experience for general continual learning: a strong, simple baseline,” in _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, Eds., 2020.
  * [31] D. Lopez-Paz and M. Ranzato, “Gradient episodic memory for continual learning,” in _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , I. Guyon, U. von Luxburg, S. Bengio, H. M. Wallach, R. Fergus, S. V. N. Vishwanathan, and R. Garnett, Eds., 2017, pp. 6467–6476.


Synchrotron cutoff in Ultraluminous X-ray sources
Ghosh et al.

Synchrotron cutoff in Ultraluminous X-ray sources

Tanuman Ghosh

Tanuman Ghosh
Astronomy and Astrophysics, Raman Research Institute, C. V. Raman Avenue, Sadashivanagar, Bangalore 560080, India

Shiv Sethi
Astronomy and Astrophysics, Raman Research Institute, C. V. Raman Avenue, Sadashivanagar, Bangalore 560080, India

Vikram Rana
Astronomy and Astrophysics, Raman Research Institute, C. V. Raman Avenue, Sadashivanagar, Bangalore 560080, India

The origin of spectral curvature at energies $E\simeq 10$ keV in ultraluminous X-ray sources is not well understood.
In this paper, we propose a novel mechanism based on synchrotron radiation to explain this cutoff. We show that relativistic plasma can give rise to observed spectral curvature for neutron star magnetic fields due to the variation in the latitude of synchrotron radiation. We analyze the  data of two bright pulsar ULXs, NGC 5907 ULX1 and NGC 7793 P13, and provide estimates of the physical parameters of these sources. We fit the data for synchrotron emission at various latitudes and show that the spectral cutoff in these cases can be explained for a large range of acceptable physical parameters, e.g., a semi-relativistic plasma with $\gamma \simeq 20$ for high latitudes or a highly relativistic plasma ($\gamma \simeq 10^5$) for emission close to the electron's orbital plane in a typical magnetic field of $B\simeq 10^{12} \, \rm G$. We also discuss how such an emission mechanism can be distinguished from other proposed models. A corollary to our study is that most ULXs might be neutron stars as they display such a spectral cutoff.

§ INTRODUCTION

Ultraluminous X-ray sources (ULXs) are some of the brightest known  X-ray sources ($\rm L_x > 10^{39} \lumcgs$). Their luminosities exceed the classical Eddington limit of a $10$ black hole (see Kaaret et al., 2017 for a recent review). In addition, many  ULXs display a unique spectral curvature at energies $\simeq 10$ keV as shown  by broadband X-ray data (e.g., Bachetti et al., 2013, Walton et al., 2013). This spectral feature is one of the distinctive characteristics of ULXs compared to the hard state of Galactic X-ray binaries (XRBs) and active galactic nuclei (AGNs). The discovery of a neutron star ULX [Bachetti et al., 2014] changed the perception of these sources, and a foremost conjecture is that a large fraction of the ULX population is neutron stars (e.g., King & Lasota, 2016, King et al., 2017, King & Lasota, 2020). Many theoretical models have been studied to explain emission mechanisms that generate such high luminosity from neutron stars (e.g., Mushtukov et al., 2015, Mushtukov et al., 2017, Mushtukov et al., 2018, Mushtukov et al., 2019). The origin of spectral cutoff, however, has no compelling theoretical model. Recent observational studies provide phenomenological models which invoke physical scenarios like compton scattering in the coronal region in low magnetic sources like black holes or comptonization in the accretion column in highly magnetized neutron stars (see e.g., West et al., 2018, Walton et al., 2018, Walton et al., 2020). In this paper, we propose an alternative model based on synchrotron radiation from different latitudes to explain the observed spectral cutoff. We explore the possible physical scenarios of this phenomenon in the context of ULXs and estimate physical parameters related to both the luminosity and the spectral cutoff in ULXs.

Synchrotron radiation is one of the most prevalent radiative processes in astrophysics [Rybicki & Lightman, 1986]. While the non-relativistic synchrotron radiation, the cyclotron radiation, provides a discrete spectrum, the emission by relativistic particles yields near-continuum spectrum owing to the higher harmonics contributing more predominantly to the observed spectrum (see Landau & Lifshitz, 1975, Rybicki & Lightman, 1986 for a review). The astrophysical implications of synchrotron radiation are well studied in multiple wavelengths, including soft to hard X-rays (e.g., Longair, 2011, Heinz, 2004, Maccarone, 2005, Markoff et al., 2005, Kisaka & Tanaka, 2017, Kisaka & Tanaka, 2017, Riegler, 1970). In this paper, we explore the impact of high-latitude, optically-thin, classical  synchrotron emission on the  spectrum of the radiation for a range of speeds  encompassing a broad range from semi-relativistic  to ultra-relativistic electrons.

In the next section, we briefly review the physics of synchrotron radiation relevant to our work. We also provide
approximate analytical expressions that allow one to study the emission from semi-relativistic to highly relativistic electrons for a range of latitudes. In  3, we provide details of the data we use and its pre-processing. The main results are presented in  4. In  5, we summarize our findings and discuss how our proposed method can be distinguished from other models.

§ SYNCHROTRON RADIATION: SEMI-RELATIVISTIC TO ULTRA RELATIVISTIC TRANSITION

We assume a geometric construct in which the incoherent synchrotron radiation originates close to the surface of a neutron star. As the length scale of magnetic fields is much larger than the curvature of the gravitating body, we can assume that the magnetic field lines are straight on scales from which the observed synchrotron emission occurs. Without loss of generality, we assume that the magnetic field is in the $z$-direction of cartesian geometry and the charged particles move in a circular motion around the uniform magnetic field lines in the xy plane.

The angular distribution of the radiated power in $n$th harmonic (or an angular frequency of observation, $\omega$) for a single electron   ($\rm erg \, sec^{-1}$)  per unit solid angle ($d\Omega$) can be expressed as [Landau & Lifshitz, 1975]:

\begin{equation}
dI_n = \frac{e^2\omega^2}{2\pi c} \left[\tan^2\theta J^2_n(n\beta\cos\theta)  +\beta^2 J'^2_n(n\beta\cos\theta) \right] {d\Omega}
\label{eq:synrad_ang}
\end{equation}

Here $\beta = v/c$,   $B$ is the magnetic field strength, and $\theta$ is the angle between radiated emission and the particle's orbital plane. $J_n(x)$ is the Bessel function  and $J_n'(x)$ is  its derivative. The integer  $n$ denotes the discrete energy levels of electron's energy with  $\omega = n \omega_B$. $\omega_B = eB/\gamma m_e c$,    $\gamma = 1/\sqrt{1-\beta^2}$ is the relativistic boost.

The aim of our study is to analyze   the  emission from both semi-relativistic and  ultra-relativistic plasmas. Eq. (<ref>) allows for the transition from cyclotron to synchrotron radiation. If the argument of the Bessel functions is small, $\beta \ll 1$, the emission is dominated by low multipoles, $n \simeq 1$ (cyclotron radiation with most of the radiation occurring at $\omega \simeq \omega_B$). As the argument of Bessel functions approaches unity, the
contribution of higher multipoles  increases. In the ultra-relativistic case with $\theta \simeq 0$ (emission close to the plane of rotation), the emission is dominated by multipoles $n \lesssim \gamma^3$, with an exponential cut-off at large frequencies. For $\gamma \gg 1$, the  spectral gap between successive multipoles $\Delta\omega =\omega_B \ll \omega$, and the emission spectrum
is near-continuum (synchrotron radiation). We discuss the case of non-zero  $\theta$ below.

Eq. (<ref>) gives the synchrotron spectrum for a single electron of energy $E = m_e c^2 \gamma$. We consider a range of electron energies and model the electron energy distribution using an exponential cutoff power law $f(\gamma) = N\gamma^{-p} \exp(-\gamma/\gamma_{\rm max})$ in the range $\gamma_{\rm min}$ and $\gamma_{\rm max}$ (e.g., Reynolds & Keohane, 1999). $N$ gives the overall normalization. For our work, we treat $\gamma_{\rm min}$ as a free parameter and $\gamma_{\rm max} = 1000 \gamma_{\rm min}$. We use the energy spectral index $p =2.2$, which is consistent with the shock acceleration mechanism (e.g., Allen et al., 2001). For this case, if $\gamma_{\rm max}$ is larger than $\gamma_{\rm min}$ by more than a few factors of 10, its impact on our results is found to be negligible.   The factors needed for conversion to flux units for comparison with the data are absorbed in the definition of $N$: $N = \rho_N V /D^2$, where $\rho_N$ is the number density
of relativistic electrons, $V$ is the volume of the emitting region, and $D$ is the luminosity distance to the source.

As noted above, Eq. (<ref>) allows one to analyze the transition from cyclotron to synchrotron radiation. For large $\gamma$, the emission is dominated by large $n$ and is restricted to an angle $\theta \simeq 1/\gamma$ centered on the plane of the orbit. For intermediate
$\gamma$ or semi-relativistic electrons ($\gamma \lesssim 10$), it is possible to have substantial emissions from higher
In this paper, we explore the possibility that the observed radiation could emanate from high latitudes with respect to the plane of the orbit.  In this case,  $\beta' = \beta \cos\theta$ acts as the effective velocity parameter in  Eq. (<ref>) and determines the frequency at which the synchrotron spectrum begins to fall exponentially.  In Figure <ref>, we
show synchrotron spectra for different values of $\theta$. As expected, for a fixed $\beta$ and $B$,  the spectral cut-off shifts to smaller harmonics $n$ for
larger $\theta$.

Figure <ref> is based on the numerical evaluation of Bessel functions in
Eq. (<ref>). One can gain more direct insight into the
relevant physics with  analytic approximations.  In the literature,  such analytic expressions have been computed for  angle-averaged emission for $\beta \simeq 1$ (e.g., Schwinger et al., 1998). However, such approximations are not valid here as the relevant parameter for us is $\beta \cos\theta$, which can deviate significantly from unity for large angles even for $\beta \simeq 1$. We find that   it is possible to approximate the Bessel function and its derivative in Eq. (<ref>) using the stationary phase approximation even when $\beta'$ deviates significantly from unity  (Schwinger et al., 1998 employs this method in  the angle-averaged case for $\beta \simeq 1$). This allows us to obtain the following approximate expressions for the Bessel function and its derivative (see Appendix <ref> for details):

\begin{align}
J_n(n\beta') & \simeq  0.447 n^{-1/3} \beta'^{-1/3} \> \> {\rm for}\> \> n \le n_c \nonumber \\
J_n(n\beta') & \simeq 0.335 n^{-1/2} (1-\beta')^{-1/4}\beta'^{-1/4} \nonumber \\
&\times \exp\left(-2\frac{\sqrt{2}}{3} n \beta'^{-1/2}(1-\beta')^{3/2} \right ) \> \> {\rm for}\> \> n \ge n_c \nonumber \\
J'_n(n\beta') & \simeq  0.411 n^{-2/3} \beta'^{-2/3} \> \> {\rm for}\> \> n \le n_c \nonumber \\
J'_n(n\beta') & \simeq 0.474 n^{-1/2} (1-\beta')^{1/4}\beta'^{-3/4}  \nonumber \\
&\times \exp\left (-2 \frac{\sqrt{2}}{3} n \beta'^{-1/2}(1-\beta')^{3/2} \right ) \> \> {\rm for}\> \> n \ge n_c
\label{eq:bess_appro}
\end{align}

\begin{equation}
n_c \simeq \frac{\beta'^{1/2}}{(1-\beta')^{3/2}}.
\label{eq:cutoffnu}
\end{equation}

$n_c$ denotes the harmonic at which spectral cutoff occurs. The analytic expressions
given in Eq. (<ref>) agree with numerical results to better
than 10% in the acceptable range of $\beta' \gtrsim 0.3$.  Also, for
$\beta' \simeq 1$ ($\beta \simeq 1$ and $\theta \simeq 0$), $1/(1-\beta') \simeq 2\gamma^2$, and  $n_c \simeq \gamma^3$ which agrees with the angle-averaged case  (Schwinger et al., 1998).

Our aim in this paper is to explain  spectral cutoff as observed in ULXs. For
fitting the  X-ray continuum spectral data, we require $n_c \gg 1$, which constrains the latitude $\theta \lesssim 70^{\circ}$. One can consider the intriguing
possibility that the observed spectrum could arise from a set of discrete lines
(though it is unlikely, as we argue below). This requires:
$n_c \ge 1$, which gives $0.3 \lesssim \beta \lesssim 1$ and $0.3 \lesssim \cos\theta \lesssim 1$ [<https://www.wolfram.com/mathematica/>]. We note that the exponential terms in Eq. (<ref>) adequately capture the cut-off frequency in the entire parameter range of interest, which is key to modeling the ULX cut-off frequency.   While we compute Bessel functions numerically for data analysis, these analytic expressions help us interpret our results.

In Figure <ref>, we display synchrotron spectra for emission from different latitudes. Eq. (<ref>) allows us to understand the spectral shapes seen in the figure. The spectral cut-off occurs
at an angular frequency $\omega \simeq n_c \omega_B$. For emission close to the plane of the rotation ($\theta \simeq 0$), $n_c \simeq \gamma^3$. However, for larger angles $n_c < \gamma^3$, as Eq. (<ref>) shows, and the spectral cut-off shifts to smaller frequencies. As we discuss later, the spectral cut-off in the data we analyze occurs at $E \simeq 10~\, \rm keV$, which is possible for a range of $\gamma$, $B$, and $\theta$ as will be discussed below in more detail. Even though we only assume electron motion
in the plane perpendicular to the magnetic field, our results do not qualitatively change if the electron has a $z$-component of velocity. This case can be incorporated into our analysis by altering $B$ to $B_\perp \equiv B \cos\chi$, where $\chi$ is the angle between the velocity vector and the magnetic field (e.g., Landau & Lifshitz, 1975).

For $B=5 \times 10^{11}$ Gauss and $\gamma=10$, the synchrotron spectra are shown as a function of multipoles  (Eq. (<ref>))  for different latitudes. The figure shows the role of high-latitude emission  in  introducing spectral curvature.

§ DATA

We utilize the  observations of two bright pulsar ULXs, NGC 5907 ULX1 (RA:15 15 58.62, DEC: +56 18 10.3; Israel et al., 2017, Walton et al., 2015, Fürst et al., 2017) and NGC 7793 P13 (RA: 23 57 50.9, DEC: -32 37 26.6; Fürst et al., 2016, Israel et al., 2017, Walton et al., 2018) to compare  our  theoretical model against  data. The distances to the host galaxies are $\simeq 17.1$ Mpc (e.g.,  Fürst et al., 2017) and $\simeq 3.5$ Mpc (e.g., Walton et al., 2018),  respectively. These two sources were observed by  several times in the past decade, which provide us an opportunity to verify the consistency of the theory over long-term spectral evolution of the sources.  In particular, for NGC 7793 P13, we detect two distinct flux states. The choice of the instrument is motivated by its energy coverage, which allows us to model the spectral curvature of the source. The broadband spectra of ULXs are generally fitted with multiple components: neutral absorption, accretion disk (geometrically thin or slim), and a phenomenological model of either magnetic or non-magnetic Comptonization processes (see e.g., Kaaret et al., 2017). The thermal disk component and neutral absorption mostly play a role in the soft energy regime ($E \lesssim 5 \, \rm keV$). Our aim in this paper is to explain the spectral cutoff in ULXs, which occurs in a higher energy range ($E \simeq \rm 10 keV$). To minimize contamination from soft components and to adequately model the spectral break, we study  the energy  range $\simeq 5\hbox{--}25$ keV in this paper.  In our study, we consider all the available  data sets for both the sources. However, for NGC 5907 ULX1, there are a few observations for which  the  signal-to-noise ratio (S/N) is poor owing to the faintness of the source. We do not utilize these data for our analysis.

§.§ Data reduction process

The  data are extracted using the  version 6.29 [<https://heasarc.gsfc.nasa.gov/docs/software/heasoft/>]. We use  tool to extract cleaned products and  tool to extract the source and background spectra and the response files from both FPMA and FPMB modules. In general we follow the method outlined in previous works (e.g., Israel et al., 2017, Walton et al., 2015, Fürst et al., 2017, Fürst et al., 2016, Israel et al., 2017, Walton et al., 2018, Lin et al., 2022) for data reduction of these two sources. We choose the source photon extraction region as $50$ arcsec radius circle for both sources. The background regions are selected as $100$ arcsec radius circle in all the cases. The number of counts per energy bin for grouping the spectra are a minimum of $30$ counts per energy bin for all NGC 5907 ULX1 spectra and $50$ counts per bin for NGC 7793 P13 spectra where the source is in high flux state, and $20$ counts per energy bin for low flux state (observation IDs - 30502019002, 30502019004, 50401003002, 90601327002) of the source.

After we obtain the spectra, we use  [Arnaud, 1996] spectral analysis package to convert the spectra into flux units for further analysis.  spectra beyond $\simeq 25$ keV are dominated by the background for both the sources, and therefore spectral data above this energy are not utilized in our analysis. In some low flux state observations, the background starts dominating well below $\simeq 20$ keV; however, to provide similar treatment to all observations, we take spectra up to $\simeq 25$ keV for all cases. We fit the  spectra with a cutoff power-law model (in  the syntax is ). The  model represents the instrumental cross-calibration differences, and  is the continuum representing an exponentially cutoff power-law spectrum. This model for $5.0\hbox{--}25.0$ keV spectra give statistically good fit for both the sources. For the cutoffpl model, we fix the index to $0.59$, a typical value for ULX pulsars (see e.g., Walton et al., 2020). We do not consider neutral absorption, since it plays a role  only in softer regime of the spectra. We then convert the spectral counts into flux $\rm \nu F_\nu$ (ergs/cm$^2$/sec) by the  command and take the data points in the energy range $5.0\hbox{--}25.0$ keV ($\simeq 1.2\hbox{--}6.0 \times 10^{18}$ Hz) to perform further analysis described in the  <ref>. We have further verified the robustness of this data extraction procedure with another model such as a simple  of photon index $0$, instead of the  model, and extract the spectra in flux unit using . We find that our results are insensitive to the  choice of  models used to generate the flux data points.

§ ANALYSIS AND RESULTS

From Eq. <ref>, one can verify that $\gamma$ and $B$ are degenerate with each other if $\beta \simeq 1$. As we wish to explore a range of electron speeds from semi-relativistic to ultra-relativistic, this degeneracy cannot be removed. [$\beta \simeq 1$ approximation is appropriate for ultra-relativistic case or semi-relativistic case in higher latitudes.] Thus we choose $\gamma_{\rm min}/B$ as one of the parameters in the analysis. We use three parameters—$\gamma_{\rm min}/B$, $\theta$ and $N$—in our analysis and later fix one more parameter to deal with residual degeneracies. We also convert the model to $\rm \nu F_\nu$ unit by appropriately scaling Eq. <ref> by a multiplicative factor $n = \omega/\omega_B$, where $\omega = 2\pi\nu$.

For our model, the observed flux can be written  as:

\begin{align}
\nu F_\nu = N'
\int_{x_{\rm min}}^{x_{\rm max}} S x^{-p} \exp\left(-x/x_{\rm max}\right) dx
\label{eq:total_spec}
\end{align}

where $x =\gamma/B$,

\begin{equation}
N' = \frac{N}{\int_{x_{\rm min}}^{x_{\rm max}} x^{-p} \exp\left(-x/x_{\rm max}\right) dx}
\end{equation}

\begin{equation}
S = \frac{2\pi e^2\nu^3}{c\nu_B} \left[\tan^2\theta J^2_\frac{\nu}{\nu_B}\left(\frac{\nu}{\nu_B}\cos\theta\right)  + J'^2_\frac{\nu}{\nu_B}\left(\frac{\nu}{\nu_B}\cos\theta\right) \right]
\end{equation}

To explore the congruence of the data and the model, we choose two different statistical methods– frequentist approach and Bayesian analysis. Given that it is hard to determine the best-fit and the errors on the three parameters simultaneously, we fix the angle $\theta$ and keep the other two parameters $\gamma_{\rm min}/B=x_{\rm min}$ and $N$ free to vary. For the frequentist approach, we first carry out a minimum $\chi^2$ analysis using the  [Virtanen et al., 2020] “curve_fit" tool [<https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html>]. The best-fit parameters, $\chi^2$ values, and 1$\sigma$ errors (computed using covariance matrix) are given in Table <ref>. For each data set, we consider three values of $\theta = 1^{\circ}, 15^{\circ}, 30^{\circ}$ for our analysis.  As discussed in the foregoing, this choice is based on the acceptable range of latitudes to ensure $n_c \gg 1$. In Figure <ref> we plot the data, the bestfit curve and   spectral residuals for one observation for each source.

In the Bayesian analysis, we utilize the Markov chain Monte Carlo (MCMC) method using    package [<https://emcee.readthedocs.io/en/stable/>] [Foreman-Mackey et al., 2013]. We find convergence in each case and the computed posterior probabilities agree with the results obtained using the frequentist method.

We next discuss the physical implications of the parameter range suggested by statistical analyses. We find that for fixed emission angles, the spectrum is scaled by the parameter $x_{\rm min}$ and the overall amplitude is scaled by $N$. For both the sources, we do not find significant long-term spectral variability in the high-energy band as the parameter $x_{\rm min}$ remains nearly the same for all observations (Table <ref>). However, due to flux changes in different epochs of observation, the parameter $N$ varies significantly.

For the range of $\theta$ shown in Table <ref>, the estimated range of
$x_{\rm min} = \gamma_{\rm min}/B$
varies from $10^{-6} \, \rm G^{-1}$ to $10^{-11} \, \rm G^{-1}$. It can be  verified from the expression for $\nu_B = \omega_B/(2\pi)$ and Eq. (<ref>) that the results given in Table <ref> correspond
to a cutoff frequency of around 10 keV. The allowed range of $x_{\rm  min}$ encompasses a large range of
particle speeds and magnetic field strengths.  At higher latitudes, our results are consistent with  semi-relativistic electrons, $\gamma_{\rm min} \simeq 20$ and $B \simeq 10^{12} \, \rm G$. This magnetic field strength is expected on the surface of neutron stars (e.g., Caballero & Wilms, 2012, Pétri, 2016). The overall normalization $N = \rho_N V/D^2$ is highly uncertain as both the relativistic electron density $\rho_N$ and the volume of emission region $V$ are very poorly determined even theoretically. Typically, in neutron star magnetosphere, the lower limit of the plasma density is given by the Goldreich-Julian limit [Goldreich & Julian, 1969] which depends on pulsar spin period, magnetic field strength, and alignment of spinning axis with magnetic field line. Depending on the volume of emission region, we find that the estimated number density could be comparable to or higher than the Goldreich-Julian limit [Goldreich & Julian, 1969] for a $1$ sec spinning pulsar with $B \simeq 10^{12}$ G , i.e., $7 \times 10^{10}$ particles $\rm cm^{-3}$. In reality, the plasma density can be significantly higher than the Goldreich-Julian limit (see Lyutikov & Gavriil, 2006 and references therein). This means our results are consistent with this theoretical expectation. We also determine  that, for a range of acceptable parameters, the emitting region is optically thin to synchrotron self-absorption and compton scattering.

If we restrict the maximum limit of magnetic field on the NS surface to the Schwinger limit of $4.4 \times 10^{13}$ G, the magnetic field strength at which the quantum effects become important,  then the maximum value of $\gamma_{\rm min}$ can be estimated. For lower $x_{\rm min}$ (i.e., lower $\gamma_{\rm min}$ or higher $B$), we get emission at higher latitudes which requires lower value of $N$ to explain the observed ULX flux. On the other hand, when $x_{\rm min}$ is higher, we get emission closer to the plane of orbit and higher $N$ is required to generate such high flux in these sources. It would be possible in the future to constrain all the parameters adequately if we can have at least one parameter determined  from other data. Our results point to the possibility that the spectral curvature in ULXs might have a common origin and all the ULXs are possibly highly magnetized neutron stars. This theoretical model can also be employed to explain  the high energy cutoff in  X-ray binary pulsar sources. Essentially, this model suggest that the curvature in the spectrum is governed by plasma velocity, magnetic field strength, electron number density, and emission angle.

The best fit parameters and $\chi^2$ from python curve fit  for 7  observations of NGC 5907 ULX1 and 10  observations of NGC 7793 P13 are listed. The errors are calculated for  1-$\sigma$ confidence from the covariance matrix using the parameter .

3c|$\theta = 1^{\circ}$   3c|$\theta = 15^{\circ}$   3c|$\theta = 30^{\circ}$

$\rm Observations$     x_min       N     χ^2 /d.o.f     x_min       N      χ^2 /d.o.f      x_min        N      χ^2 /d.o.f

(10^-7  G^-1)      (10^-17  cm^-2)          (10^-10  G^-1)      (10^-18  cm^-2)            (10^-11  G^-1)      (10^-19  cm^-2)

NGC 5907 ULX1

30002039005     (9.54 ±0.56)     (1.10 ±0.18)     57/49      (2.81 ±0.16)     (0.70 ±0.12)     57/49      (3.45 ±0.20)     (2.97 ±0.49)     57/49

30302004006     (6.38 ±0.36)     (1.00 ±0.12)     66/52    (1.88 ±0.11)      (0.63 ±0.08)     66/52     (2.31 ±0.13)     (2.72 ±0.34)     66/52

30302004008     (5.94 ±0.35)     (0.81 ±0.10)     62/49    (1.75 ±0.10)     (0.51 ±0.06)     62/49     (2.15 ±0.12)     (2.21 ±0.27)     62/49

80001042002     (7.69 ±0.27)      (2.52 ±0.23)     80/78     (2.27 ±0.08)      (1.60 ±0.14)      80/78     (2.78 ±0.10)      (6.84 ±0.61)      80/78

80001042004     (6.57 ±0.26)      (1.69 ±0.15)     82/77     (1.94 ±0.08)     (1.07 ±0.10)      82/77     (2.38 ±0.09)     (4.58 ±0.41)      83/77

90501331002     (5.97 ±0.46)     (0.47 ±0.08)     34/39     (1.76 ±0.14)      (0.30 ±0.05)      34/39     (2.16 ±0.16)      (1.28 ±0.21)     34/39

90601323002     (7.23 ±0.29)     (1.66 ±0.16)      89/88     (2.13 ±0.09)      (1.05 ±0.10)      90/88     (2.62 ±0.10)      (4.50 ±0.43)     90/88

NGC 7793 P13

30302005002     (7.66 ±0.20)       (3.07 ±0.20)      77/82     (2.26 ±0.06)     (1.95 ±0.13)      77/82     (2.77 ±0.07)     (8.30 ±0.55)     77/82

30302005004     (7.46 ±0.15)     (4.58 ±0.23)      111/117     (2.20 ±0.04)      (2.90 ±0.15)      111/117     (2.70 ±0.05)      (12.40 ±0.62)     111/117

30302015002      (7.53 ±0.14)      (5.95 ±0.27)      147/133     (2.22 ±0.04)      (3.77 ±0.17)      147/133     (2.73 ±0.05)      (16.10 ±0.73)      148/133

30302015004     (7.53 ±0.17)      (4.75 ±0.26)      130/123     (2.22 ±0.05)      (3.02 ±0.17)      130/123     (2.72 ±0.06)      (12.88 ±0.70)      130/123

30502019002     (7.41 ±0.53)      (0.65 ±0.11)     42/61     (2.18 ±0.16)      (0.41 ±0.07)     42/61     (2.68 ±0.19)      (1.77 ±0.30)      42/61

30502019004     (6.93 ±0.62)      (0.57 ±0.12)      32/42     (2.04 ±0.18)     (0.36 ±0.08)      32/42     (2.51 ±0.22)      (1.55 ±0.32)      32/42

50401003002     (7.37 ±0.79)     (0.59 ±0.15)      25/31     (2.17 ±0.23)      (0.38 ±0.10)      25/31     (2.67 ±0.28)      (1.60 ±0.41)      25/31

80201010002     (7.26 ±0.10)      (6.20 ±0.22)      156/204     (2.14 ±0.03)      (3.94 ±0.14)      156/204     (2.63 ±0.04)      (16.81 ±0.58)     156/204

90301326002     (7.20 ±0.14)      (7.11 ±0.33)      104/123     (2.12 ±0.04)     (4.51 ±0.21)     104/123     (2.61 ±0.05)     (19.26 ±0.89)      105/123

90601327002     (9.65 ±1.01)     (1.00 ±0.30)      35/34     (2.84 ±0.30)      (0.63 ±0.19)      35/34     (3.49 ±0.36)      (2.70 ±0.80)     35/34

Example plot showing spectral fitting of the theoretical model on one observation for each source. The top panels show the data with model over plotted and corresponding residuals are shown in bottom panels. The model shown here are for $1^{\circ}$ angle in case of NGC 5907 ULX1 (left; Observation - 80001042002) and $30^{\circ}$ angle for NGC 7793 P13 (right; Observation - 90301326002).

§ DISCUSSION

In this paper, we seek to explain the X-ray spectra of two known pulsar ULXs. In particular, we focus on the  cut-off in such spectra at $E \simeq 10$ keV, which is a generic feature of  many pulsar ULXs.
We propose synchrotron radiation at a range of latitudes  as the possible physical process to explain this observed spectral shape. Our main results, based on the analysis  of 17  spectra of the two sources, are summarized in Table <ref>. Figure <ref> shows the fit and its residual for one spectrum for each source.

Other  models that have been explored  to explain the dominant emission in hard  X-ray range in ULXs invoke the comptonization from coronal region of a non-magnetic sources and comptonization from magnetized column in neutron stars (see e.g., West et al., 2018, Walton et al., 2020). Compton scattering in presence of high magnetic field in neutron star is a possible candidate to explain high luminosity in these sources [Mushtukov et al., 2015]. In principle, there could be two possible ways to distinguish our proposed scenario from these models.

Table <ref> shows that the best-fit value $x_{\rm min} = \gamma_{\rm min}/B \simeq 2 \times 10^{-11}$ for $\theta \simeq 30^\circ$. This could correspond, for instance, to a semi-relativistic electron ($\gamma_{\rm min} \simeq 20$) along with a magnetic field $B\simeq 10^{12} \, \rm G$. In such cases, the fundamental mode of emission, $\nu_B = eB/(2\pi m_ec\gamma_{\rm min}) \simeq  600\, \rm eV$. As this is larger than the spectral resolution of  in energy range of interest, the observed spectrum could be a set of discrete cyclotron lines. In practice, such an interpretation could be difficult owing to mixing with larger $\gamma$ values and the width of spectral lines, which are  difficult to ascertain. This would also require either re-analysis of the data or new data which is beyond the scope of the paper. Our analysis raises the intriguing possibility that the discreteness of the spectrum could probe the latitude of the emission. As we have already discussed in  <ref>, an  upper limit on the  latitude  of emission can be obtained by  requiring $n_c \gg 1$; this yields a stringent upper bound $\theta \simeq 70{^\circ}$. For fitting  continuum X-ray spectral data  in the energy range of interest,  this requirement    motivates the upper limit of   $\theta \simeq 30^{\circ}$ we use in this paper.

Another possible probe of our model could be the polarization of received photons. The photons emerging from higher latitudes would be elliptically polarized while those from closer to the plane of the orbit would be linearly polarized. While the non-magnetic comptonization will not show polarized emission, the magnetic comptonization and the synchrotron radiation could display different degrees of polarization. Modern X-ray polarimeters such as IXPE [Weisskopf et al., 2016] and upcoming mission POLIX [Paul et al., 2010] might be able to address these questions.

We provide a brief summary of our main results and perspectives below:

1. We propose that spectral cutoff in  ULXs arise from  classical, high-latitude, and  optically-thin synchrotron radiation. For classical radiation, the cutoff occurs at energies $\gamma^3 \nu_B$ for radiation close to the plane of the orbit, but the cutoff frequency shifts to much smaller frequencies for high-latitude emission. Quantum effects only dominate the cutoff for energies close to the electron rest mass  and hence cannot be responsible for the observed  cutoff at $E \simeq 10~\, \rm keV$.

2. This model is compared with the 17  spectra, corresponding to different flux states of two ULXs. The observed fluxes are  modeled using four theoretical parameters. Given the degeneracy between these parameters, only two parameters can be  estimated from the data.  To test  the robustness of our statistical analysis, we  carry out both frequentist and bayesian analysis (using MCMC). While our analysis yields a large range of possible theoretical models, the most interesting case
corresponds to high-latitude  emission  ($\theta \simeq 30^\circ$) from a  semi-relativistic plasma from the surface of  the neutron star ( $B\simeq 10^{12} \, \rm G$ and $\gamma \simeq 20$). The statistical analysis also allows us
to establish that the plasma is optically thin for a plausible range of parameters.

3. It is possible to verify the model using X-ray polarization data which might be available in the near future. Another possible probe of the semi-relativistic plasma could be the discreteness of the observed spectrum, which we have briefly discussed.

We would like to thank the referee for the valuable suggestions that helped further improve the manuscript. We would like to thank Keith Arnaud from HEASARC helpdesk for his valuable suggestions regarding eeufspec tool in XSPEC. This research has utilized archival data (available at the High Energy Astrophysics Science Archive Research Center (HEASARC)) obtained with , a project led by Caltech, funded by NASA, and managed by the NASA Jet Propulsion Laboratory (JPL), and has made use of the NuSTAR Data Analysis Software () jointly developed by the ASI Space Science Data Centre (SSDC, Italy) and the California Institute of Technology (Caltech, USA).

; Harrison et al., 2013

HEASOFT (<https://heasarc.gsfc.nasa.gov/docs/software/heasoft/>; Nasa High Energy Astrophysics Science Archive Research Center (Heasarc), 2014), Mathematica (<https://www.wolfram.com/mathematica/>; Wolfram Research, Inc., 2022)

§ ANALYTIC APPROXIMATION OF BESSEL FUNCTION

The starting point of approximating the Bessel function and its derivatives is
the integral representation of these functions (e.g. [Schwinger et al., 1998, Landau & Lifshitz, 1975]),

\begin{align}
J_n(z) = \int_0^\pi d\phi \frac{1}{\pi} \cos(z\sin\phi-n\phi)\\
J'_n(z) = -\int_0^\pi d\phi \frac{1}{\pi} \sin\phi \sin(z\sin\phi-n\phi)
\end{align}

Here $z = n\beta\cos\theta = n \beta'$. As the integrads are highly oscillatory, the main contribution to the integrals arise from  regions near $\phi=0$ when the phase is large (stationary phase approximation). This is ensured by the condition $n\beta' \gg 1$. Expanding the phase  of $J_n(z)$ around $\phi = 0$, we get:

\begin{equation}
z\sin\phi-n\phi = n \beta'\sin\phi-n\phi = - n \left[\phi(1-\beta')+\frac{\beta'\phi^3}{3!}\right ]
\end{equation}

Making the substitution,  $\phi = (1-\beta')^{1/2}x/\beta'^{1/2} $ yields:

\begin{equation}
\phi(1-\beta')+\frac{\beta'\phi^3}{3!} = \frac{(1-\beta')^{3/2}}{\beta'^{1/2}} \left (x+\frac{x^3}{6}\right )
\end{equation}

In this case, the stationary phase points are located at:

\begin{equation}
x = \pm \sqrt{2}i
\end{equation}

Following the procedure outlined in Schwinger et al., 1998,  in the neighborhood of the stationary phase point, we can write:

\begin{equation}
x = \sqrt{2}i+\xi,
\end{equation}

where $\xi$ is real and  small, which gives:

\begin{equation}
x+\frac{x^3}{6} = \sqrt{2} i \left(\frac{2}{3}+\frac{\xi^2}{2}\right)
\end{equation}

This allows us to write:

\begin{equation}
J_n(n\beta') = \int_0^\infty dx \frac{1}{\pi} \frac{(1-\beta')^{1/2}}{\beta'^{1/2}} \cos\left(n\left[\frac{(1-\beta')^{3/2}}{\beta'^{1/2}} \left(x+\frac{x^3}{6}\right)\right]\right)
\label{eq:bessaprro1}
\end{equation}

First, we deal with the case
when $n (1-\beta')^{3/2}/\beta'^{1/2} \ll  1$. In this case,  the main contribution to the integral  come from the region where $x$ is large.
Given that most of the contribution to the integral comes from regions,
where the phase is close to unity, the integration limit can be extendend to  infinity (e.g. [Schwinger et al., 1998]). Solving the resultant integral,  we get:

\begin{equation}
J_n(n\beta')  \simeq  0.447 n^{-1/3} \beta'^{-1/3} \> \> {\rm for}\> \> n \frac{(1-\beta')^{3/2}}{\beta'^{1/2}} \ll 1
\end{equation}

For $n (1-\beta')^{3/2}/\beta'^{1/2} \gg  1$,  the integral can be written as:

\begin{equation}
J_n(n\beta') = Re \int_0^\infty dx \frac{1}{\pi} \frac{(1-\beta')^{1/2}}{\beta'^{1/2}} \exp\left(i n(1-\beta')^{3/2} (x+x^3/6)/\beta'^{1/2}\right)
\end{equation}

This can readily be integrated as most of contribution arises from regions
close to $x\simeq 0$:

\begin{equation}
J_n(n\beta') \simeq 0.335 n^{-1/2} (1-\beta')^{-1/4}\beta'^{-1/4} \times \exp\left(-2\frac{\sqrt{2}}{3} n \beta'^{-1/2}(1-\beta')^{3/2} \right )
\end{equation}

Following a similar calculation procedure, we obtain the approximate forms of $J'_n(n\beta')$ in Eq. <ref>.

Our analysis extends the procedure outlined by Schwinger et al., 1998 for the extreme relativistic case $\beta' \simeq 1$ to  arbitrary $\beta'$.  We note that the analytic expressions we derive yield excellent fits to  numerical results for $\beta' \gtrsim 0.3$.

[Allen et al., 2001]
Allen, G. E., Petre, R., & Gotthelf, E. V. 2001, , 558, 739,

[Arnaud, 1996]
Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series,
Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H.
Jacoby & J. Barnes, 17

[Bachetti et al., 2013]
Bachetti, M., Rana, V., Walton, D. J., et al. 2013, , 778, 163,

[Bachetti et al., 2014]
Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, , 514,
202, 10.1038/nature13791

[Caballero & Wilms, 2012]
Caballero, I., & Wilms, J. 2012, , 83, 230.

[Foreman-Mackey et al., 2013]
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, ,
125, 306, 10.1086/670067

[Fürst et al., 2017]
Fürst, F., Walton, D. J., Stern, D., et al. 2017, , 834, 77,

[Fürst et al., 2016]
Fürst, F., Walton, D. J., Harrison, F. A., et al. 2016, , 831,
L14, 10.3847/2041-8205/831/2/L14

[Goldreich & Julian, 1969]
Goldreich, P., & Julian, W. H. 1969, , 157, 869,

[Harrison et al., 2013]
Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ,
770, 103, 10.1088/0004-637X/770/2/103

[Heinz, 2004]
Heinz, S. 2004, , 355, 835, 10.1111/j.1365-2966.2004.08361.x

[Israel et al., 2017]
Israel, G. L., Belfiore, A., Stella, L., et al. 2017a,
Science, 355, 817, 10.1126/science.aai8635

[Israel et al., 2017]
Israel, G. L., Papitto, A., Esposito, P., et al. 2017b,
, 466, L48, 10.1093/mnrasl/slw218

[Kaaret et al., 2017]
Kaaret, P., Feng, H., & Roberts, T. P. 2017, , 55, 303,

[King & Lasota, 2016]
King, A., & Lasota, J.-P. 2016, , 458, L10,

[King & Lasota, 2020]
—. 2020, , 494, 3611, 10.1093/mnras/staa930

[King et al., 2017]
King, A., Lasota, J.-P., & Kluźniak, W. 2017, , 468, L59,

[Kisaka & Tanaka, 2017]
Kisaka, S., & Tanaka, S. J. 2017a, , 837, 76,

[Kisaka & Tanaka, 2017]
Kisaka, S., & Tanaka, S. J. 2017b, in Journal of Physics
Conference Series, Vol. 932, Journal of Physics Conference Series, 012015,

[Landau & Lifshitz, 1975]
Landau, L. D., & Lifshitz, E. M. 1975, The classical theory of fields

[Lin et al., 2022]
Lin, L. C.-C., Hu, C.-P., Takata, J., et al. 2022, , 924, 65,

[Longair, 2011]
Longair, M. S. 2011, High Energy Astrophysics

[Lyutikov & Gavriil, 2006]
Lyutikov, M., & Gavriil, F. P. 2006, , 368, 690,

[Maccarone, 2005]
Maccarone, T. J. 2005, , 360, L68,

[Markoff et al., 2005]
Markoff, S., Nowak, M. A., & Wilms, J. 2005, , 635, 1203,

[Mushtukov et al., 2019]
Mushtukov, A. A., Ingram, A., Middleton, M., Nagirner, D. I., & van
der Klis, M. 2019, , 484, 687, 10.1093/mnras/sty3525

[Mushtukov et al., 2017]
Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Ingram, A.
2017, , 467, 1202, 10.1093/mnras/stx141

[Mushtukov et al., 2015]
Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J.
2015, , 454, 2539, 10.1093/mnras/stv2087

[Mushtukov et al., 2018]
Mushtukov, A. A., Tsygankov, S. S., Suleimanov, V. F., & Poutanen, J.
2018, , 476, 2867, 10.1093/mnras/sty379

[Nasa High Energy Astrophysics Science Archive Research Center (Heasarc), 2014]
Nasa High Energy Astrophysics Science Archive Research Center (Heasarc).
2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source
Code Library, record ascl:1408.004.

[Paul et al., 2010]
Paul, B., Rishin, P. V., Maitra, C., et al. 2010, in The First Year of
MAXI: Monitoring Variable X-ray Sources, 68

[Pétri, 2016]
Pétri, J. 2016, Journal of Plasma Physics, 82, 635820502,

[Reynolds & Keohane, 1999]
Reynolds, S. P., & Keohane, J. W. 1999, , 525, 368,

[Riegler, 1970]
Riegler, G. R. 1970, , 226, 1041, 10.1038/2261041a0

[Rybicki & Lightman, 1986]
Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in

[Schwinger et al., 1998]
Schwinger, J., DeRaad, L., Milton, K., & Tsai, W.-Y. 1998, Classical
Electrodynamics (1st ed.) (CRC Press)

[Virtanen et al., 2020]
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17,
261, 10.1038/s41592-019-0686-2

[Walton et al., 2013]
Walton, D. J., Fuerst, F., Harrison, F., et al. 2013, , 779, 148,

[Walton et al., 2015]
Walton, D. J., Harrison, F. A., Bachetti, M., et al. 2015, , 799,
122, 10.1088/0004-637X/799/2/122

[Walton et al., 2018]
Walton, D. J., Fürst, F., Heida, M., et al. 2018a,
, 856, 128, 10.3847/1538-4357/aab610

[Walton et al., 2018]
Walton, D. J., Fürst, F., Harrison, F. A., et al.
2018b, , 473, 4360, 10.1093/mnras/stx2650

[Walton et al., 2020]
Walton, D. J., Pinto, C., Nowak, M., et al. 2020, , 494, 6012,

[Weisskopf et al., 2016]
Weisskopf, M. C., Ramsey, B., O'Dell, S., et al. 2016, in Society of
Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905,
Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed.
J.-W. A. den Herder, T. Takahashi, & M. Bautz, 990517,

[West et al., 2018]
West, L. A., Lehmer, B. D., Wik, D., et al. 2018, , 869, 111,

[Wolfram Research, Inc., 2022]
Wolfram Research, Inc. 2022, Mathematica, Version 13.1.