# Broad learning system with Takagi-Sugeno fuzzy subsystem for tobacco origin
identification based on near infrared spectroscopy

Di Wang  Simon X. Yang  School of Information Science and Engineering,
Chongqing Jiaotong University, Chongqing 400074, China  School of Engineering,
University of Guelph, Guelph, Ontario, N1G 2W1, Canada
e-mail of corresponding author<EMAIL_ADDRESS>
###### Abstract

Tobacco origin identification is significantly important in tobacco industry.
Modeling analysis for sensor data with near infrared spectroscopy has become a
popular method for rapid detection of internal features. However, for sensor
data analysis using traditional artificial neural network or deep network
models, the training process is extremely time-consuming. In this paper, a
novel broad learning system with Takagi-Sugeno (TS) fuzzy subsystem is
proposed for rapid identification of tobacco origin. Incremental learning is
employed in the proposed method, which obtains the weight matrix of the
network after a very small amount of computation, resulting in much shorter
training time for the model, with only about 3 seconds for the extra step
training. The experimental results show that the TS fuzzy subsystem can
extract features from the near infrared data and effectively improve the
recognition performance. The proposed method can achieve the highest
prediction accuracy (95.59$\%$) in comparison to the traditional
classification algorithms, artificial neural network, and deep convolutional
neural network, and has a great advantage in the training time with only about
128 seconds.

###### keywords:

Sensor data analysis; broad learning system; incremental learning; fuzzy
logic; NIR spectroscopy; tobacco origin identification

##  1 Introduction

Although tobacco is harmful to human health, research on tobacco still has
certain value and significance. Tobacco has been widely known as a commercial
crop, and the quality and flavor of tobacco leaf are affected heavily by the
cultivation geographical region [1, 2]. Therefore, distinguishing the growing
regions of tobacco leaves plays a significant role before putting them into
products [3, 4, 5]. In practical application in tobacco industry, growing
region evaluation is usually conducted by some experts according to the
sensory inspection like sense of smell, taste, and smoke quality, which is
very time-consuming, low efficiency and laborious. The manual evaluation of
growing region depends on the experiences of experts to a great extent, and
the results could be affected by the experts’ emotional state as well as the
environment outside such as the light and temperature of the assessment place.
Such kind of manual evaluation of growing region is not reliable and cannot
meet the requirement of reproducible assessment process for the tobacco
quality control and supervision.

With the development of soft computing technology, many identification
techniques for tobacco origin have been applied. As the distribution of metal
elements in tobacco leaves from different producing areas is different, X-ray
fluorescence spectroscopy was used to rapidly detect the distribution of
various metal elements in tobacco leaves, and discriminant analysis was
combined to identify the tobacco origin [6]. Electronic nose was adopted to
collect data from different sensors and combined principal component analysis,
clustering analysis, and other algorithms to identify tobacco origin [7]. It
has proven that pollen analysis is able to assist with identifying
geographical origin of tobacco [8]. As chemical composition content being
treated as independent variable affecting the producing area of flue-cured
tobacco, Naive Bayes classification algorithm was established for pattern
recognition of tobacco origin [9].

Near infrared (NIR) spectroscopy has been known as a kind of fast,
nondestructive and accurate data analysis technology for origin identification
in recent years [10, 11, 12, 13, 14, 15, 16]. NIR spectroscopy combined with
other soft computing methods has been widely applied for tobacco origin
identification. Zhu et al. used high dimensional feature grouping method in
NIR spectra to identify tobacco growing area [17]. Wang et al. collected
twelve hundred seventy six superior tobacco leaf from four producing areas,
which are Yuxi, Chuxiong, Zhaotong and Dali in Yunnan province, and applied
NIR spectrum projection and correlation methods for tobacco quality analysis
of different producing areas [18].

Although the above soft computing methods obtain better effect than manual
sensory inspection for identifying tobacco growing regions, there is room for
improvement in recognition performance.

Deep learning for deep structural neural network is a popular machine learning
algorithm in recent years [19]. It has achieved breakthrough success in many
fields [19, 20, 21, 22, 23], especially in big data mining. Lu et al.
conducted classification modeling method for NIR spectroscopy of tobacco
leaves based on Convolution Neural Network (CNN) [24]. Although the deep
network shows amazing effects and strong prediction performance, due to the
coupling between layers requires a large number of hyper-parameters, the
complexity not only increases the difficulty of analyzing the deep network
structure, but also consumes a lot of time and resources in the training
process [25]. Wang et al. [26] proposed CNN to discriminate the tobacco
cultivation regions and achieves the prediction precision with 93.16$\%$, but
takes about 6833 seconds for the training process.

In addition, with the arrival of the era of big data, the sensor data
acquisition quantity increases quickly. When the components or property of the
new data are changed, the analysis of new data with the trained model before
will be affected. New data and the previous data need to be retrained for
establishing mathematical model, which is extremely time-consuming and
inconvenient.

Therefore, to address the time-consuming problems happened in deep network
during training for big data and retraining the network due to new sensor
data, the tobacco origin classification model of fuzzy broad learning system
based on Takagi-Sugeno (FBLS-TS) is proposed in this paper. In this soft
computing method, the effective features are extracted adaptively through TS
fuzzy subsystem, and the pseudo-inverse calculation and incremental learning
method are applied for data processing to improve the speed of model training.

##  2 Materials and Methods

In this section, the FBLS-TS model and data source used in this study are
first introduced. After that, the model evaluation criteria are presented in.

###  2.1 Broad learning system

Broad Learning System (BLS) is an efficient incremental learning system
without deep structures, which is proposed by Professor Chen of the University
from Macau university in 2018 [27]. Since there are no multi-layer
connections, BLS does not need to use gradient descent to update the weights,
so the calculation speed is much faster than that of deep learning algorithms
[28, 29]. The principle of the BLS is as follows, it firstly extracts feature
from raw sensor data through linear transformation, and makes the extracted
features to be feature nodes; then the feature nodes are enhanced to be
enhanced nodes; lastly, both of the feature nodes and the enhanced nodes are
taken as the input layer and directly connected to the output layer, and the
connection weight from the input to the output is obtained through the pseudo-
inverse of the network input. Supposing that $A$ is the input matrix which
consists of feature nodes and enhanced nodes, $Y$ is the output, and $W$ is
the weight matrix of the network, the objective of the training is to find the
connection weight, which is given as

$\displaystyle W\text{=}{{A}^{\text{-1}}}Y,$ (1)

where ${{A}^{\text{-1}}}$ is the inverse matrix of the input for the entire
network, and $Y$ is the label of data. The network structure of the BLS is
shown in Figure 1. In the figure, $X$ is the input of the network, $Z_{1}$,
$Z_{2}$, …, $Z_{n}$ are the feature nodes, $H_{1}$, $H_{2}$, …, $H_{m}$ are
the enhanced nodes, $W$ is the weight of network and $Y$ is the output.

Figure 1: Network structure of BLS.

####  2.1.1 Feature nodes

Assuming the NIR data of tobacco is $X_{s\times f}^{{}^{\prime}}$, where $s$
is the number of training sample, and $f$ is the number of features.
Normalization is conducted on $X_{s\times f}^{{}^{\prime}}$ to ensure that the
data was normalized within the same range firstly, and then a column of 1 is
added to the last column of $X_{s\times f}^{{}^{\prime}}$ to make it to be
${{X}_{s\times(f+1)}}$, the purpose of which is to make it convenient to add
the bias in the matrix calculation when generating feature nodes. Supposing
that ${{N}_{1}}$ is the number of feature nodes in per window of the network
structure. Assuming ${{N}_{2}}$ is the number of windows of the network
structure and ${{N}_{3}}$ is the number of enhanced nodes. The process of
generating feature nodes for each window is as follows,

Step 1: Generating a Gaussian random weight matrix of with dimension of
$(f+1)\times{{N}_{1}}$.

Step 2: Putting the value of $w_{e}$ in the variable of
$W_{e}\left\\{i\right\\}$, where $i=1\dots N_{2}$, which represents the number
of iterations and is equal to the number of generated windows.

Step 3: Letting ${{A}_{1}}=Xw_{e}$, the random generated weight matrix is used
to convolve the features of each sample and obtain new features to form
feature nodes of one window.

Step 4: Normalizing ${{A}_{1}}$, then the lasso method is conducted for sparse
autoencoder to get a sparse weight matrix of ${{W}_{z}}$, which is given as

$\displaystyle{{W}_{z}}=\mathrm{argmin}\left(\left\|{{A}_{1}}{{W}_{z}}-X\right\|_{2}^{2}+\lambda{{\left\|{{W}_{z}}\right\|}_{1}}\right),$
(2)

where $\left\|\cdot\right\|_{2}$ is L2 normalization and
$\left\|\cdot\right\|_{1}$ is L1 normalization.

Step 5: Generating feature nodes of ${{Z}_{1}}$ in one window with
normalization method, and the feature nodes of ${{Z}_{1}}$ is given as

$\displaystyle{{Z}_{1}}=\text{norm}(X{{W}_{z}}).$ (3)

The above process is just the calculation for one window, and for ${N_{2}}$
windows, it will generate ${{N}_{1}}$ feature nodes for every window.
Therefore, the matrix of feature node for the whole network is given as

$\displaystyle Z=\left[{{Z}_{1}},{{Z}_{2}},\cdots,{{Z}_{n}}\right],$ (4)

where $Z$ is the matrix with $s\times({{N}_{2}}\times{{N}_{1}})$ dimensions,
and $n={{N}_{2}}\times{{N}_{1}}$.

####  2.1.2 Enhancement nodes

One of characters for BLS is the enhancement node which is introduced to
supplement the feature nodes as the input of the network. Adding enhanced
nodes can increase the nonlinear features of the network and better fit the
nonlinear function. The process of obtaining enhancement nodes is described as
follows,

Step 1: The feature node matrix is standardized and augmented to obtain the
enhancement nodes $H$. The weight matrix of ${{W}_{h}}$ is an orthogonal
normalized random matrix. The purpose of this operation is to map the feature
nodes with linear characteristics to another high-dimensional space, in which
way the network has stronger expression ability.

Step 2: The enhancement nodes are activated to achieve the parameter of
${H}^{\prime}$, which is calculated as

$\displaystyle{H}^{\prime}=\mathrm{tansig}\left(\frac{H{{W}_{h}}s}{\max(H{{W}_{h}})}\right),$
(5)

where $\mathrm{tansig}(\cdot)$ is a hyperbolic tangent activation function,
and $s$ is the zoom factor.

Step 3: The feature nodes and enhancement nodes are combined as input of the
entire network, which is given as

$\displaystyle A=\left[Z,{H}^{\prime}\right],$ (6)

and the dimension of the input is ${{N}_{2}}\times{{N}_{1}}+{{N}_{3}}$.

####  2.1.3 Pseudo-inverse and incremental learning

The purpose of the network is to obtain the connection weight, which is
obtained from

$\displaystyle W={{A}^{-1}}Y.$ (7)

It can be calculated by ridge regression as

$\displaystyle\mathrm{argmin}\left(\left\|AW-Y\right\|_{2}^{2}+\lambda\left\|W\right\|_{2}^{2}\right),$
(8)

then we get

$\displaystyle W={{\left(\lambda I+A{{A}^{T}}\right)}^{-1}}{{A}^{T}}Y,$ (9)

where, the expression of generalized inverse of $A$ is defined as

$\displaystyle{{A}^{+}}=\underset{\lambda\to 0}{\mathop{\lim}}\,{{(\lambda
I+A{{A}^{T}})}^{-1}}{{A}^{T}}.$ (10)

Then, the connection weight is achieved through the above calculations.

In practical applications, more and more sensor data will be collected. When
the sensor data volume reaches the order of billions, it will be extremely
time-consuming to find the pseudo-inverse of such a matrix with hundreds of
millions of rows. At the same time, due to the insufficient fitting ability of
the model, the dimension of sensor data will be increased, that is, the number
of enhanced nodes will be increased. Therefore, incremental learning algorithm
is needed, which is also the essence of BLS.

Incremental learning is a dynamic and progressive updating algorithm. Its
advantage lies in that the updated network connection weight matrix can be
obtained with only a small amount of computation through the previous
calculation results and the obtained new data.

The input of original network is denoted by ${{A}_{n}}$, the new added
enhancement nodes is denoted by $a$. The updated input is given as
${{A}_{n+1}}=[{{A}_{n}}|a]$, and the updated weight is denoted as
${{W}_{new}}={{[{{A}_{n}}|a]}^{-1}}Y$. The generalized inverse of
${{A}_{n+1}}$ is defined as

$\displaystyle{{A}_{n+1}}^{+}={{[{{A}_{n}}|a]}^{+}}=\left[\begin{matrix}A_{n}^{+}-d{{b}^{T}}\\\
{{b}^{T}}\\\ \end{matrix}\right],$ (11)

where $Y$ is the label of training set. From this way, the input and the
weight of the entire network are obtained.

After the training process, the connection weight and normalized parameters
need to be saved. The dimension of weight for the entire network is
${{N}_{2}}\times{{N}_{1}}+{{N}_{3}}\times k$, where $k$ is the number of
classifications.

###  2.2 Fuzzy BLS based on Takagi-Sugeno

The method of FBLS-TS combines the fuzzy logic and neural network, and adjusts
the parameters of the antecedent and subsequent fuzzy system to improve the
reasoning ability of the fuzzy system.

The BLS model generates feature nodes through randomly generating weights and
sparse coefficients, but the FBLS-TS model obtains the feature nodes through
the inference and feature extraction from the original input data by Takagi-
Sugeno (TS) fuzzy subsystem. The network structure is shown in Figure 2.

Figure 2: Network structure of the proposed FBLS-TS.

As seen in the figure, the large dotted box is one adaptive fuzzy subsystem,
and it contains two parts which are antecedent network and subsequent network.
The tobacco NIR data is denoted as
$X={{\left[{{x}_{1}},{{x}_{2}},\cdots,{{x}_{f}}\right]}^{T}}$, and $f$ is the
number of feature variables.

There are four layers in the subsequent network, and the calculation process
is as follows,

1) The first layer is the input layer with $f$ neuronal nodes, and the output
of this layer is the value of every feature variable.

2) The second layer is the membership function layer, which is used to fuzzify
the input data. ${{A}_{1}},{{A}_{2}},\ldots,{{A}_{{{N}_{1}}}}$ are fuzzy sets
obtained by the input data with the Gaussian function, and they denote the
number of the fuzzy rules. The output of this layer is defined as

$\displaystyle\mu_{i}^{j}={{e}^{-\frac{{{\left({{x}_{i}}-{{c}_{ij}}\right)}^{2}}}{{{\sigma}_{ij}}^{2}}}},\qquad{i=1,2,\cdots,f};\quad{j=1,2,\cdots,{k_{i}}},$
(12)

where, the center vector of $c$ is obtained from clustering algorithm of
K-Mean, and denotes the standard deviation.

3) The third layer is the fuzzy rule layer. Every neuron in this layer
represents one fuzzy rule and calculates the fitness. The output of this layer
is given as

$\displaystyle{{w}_{j}}=\mu_{1}^{{{i}_{1}}}\mu_{2}^{{{i}_{2}}}\cdots\mu_{n}^{{{i}_{n}}},$
(13)

where
${i}_{1}\in\\{1,2,\cdots,{{k}_{1}}\\};\,{i}_{2}\in\\{1,2,\cdots,\,{{k}_{2}}\\};\cdots;\,{i}_{f}\in\\{1,2,\cdots,{{k}_{f}}\\};\,j=1,2,\cdots,k;\,k=\underset{i=1}{\overset{f}{\mathop{\prod}}}\,{{k}_{i}}$.

4) The fourth layer is the normalization layer, and its equation is given as

$\displaystyle\overline{{{w}_{j}}}=\frac{{{{w}_{j}}}}{{{\sum\limits_{j=1}^{k}{{{w}_{j}}}}}},\qquad{j=1,2,\cdots,k}.$
(14)

There are three layers in subsequent network, and the calculation process is
as follows,

1) The first layer is the input layer and its function is to convey the value
of feature variable to the next layer.

2) The second layer is the rule layer and every neuronal node represents one
rule. Its output is given as

$\displaystyle{{q}_{j}}={{p}_{j0}}+{{p}_{j1}}{{x}_{1}}+\cdots+{{p}_{jf}}{{x}_{f}}=\sum\limits_{i=0}^{f}{{{p}_{ji}}{{x}_{i}}}.$
(15)

3) The third layer is the output layer of every fuzzy rule, which is used for
generating the feature nodes, and it is calculated as

$\displaystyle{{o}_{j}}={\overline{w}}_{j}{{q}_{j}},\text{ }\quad
j=1,2,\cdots,k.$ (16)

###  2.3 Tobacco database

A total number of 13370 tobacco samples were collected from eight different
areas of Guizhou Province by Guiyang flue-cured tobacco factory in 2018. The
tobacco NIR spectra of the data were recorded by Thermo Antaris 2 (Thermo
Fisher Scientific Inc. Waltham, USA.) with the spectral resolution of 8 cm-1
and 64 scans in the NIR range of 3800 cm-1 to 10000 cm-1, which is shown in
Figure 3.

Figure 3: Raw NIR spectra of the 13370 samples.

NIR spectrum is generated when energy transition occurs due to anharmonic
vibration of molecules and selective absorption of NIR light in specific
frequency bands. It mainly records the absorption of frequency doubling and
frequency combination of vibration of hydrogen groups (C-H, O-H, N-H, etc.).
The absorption intensity and wavelength of NIR light by the same group are
different in different physical and chemical environments, which is the
principle of qualitative analysis.

Data were randomly divided into training set and testing set, and the regional
distribution of training set is shown in Table 1.

Table 1: Regional distribution of training set. Regions |  Sample number
---|---
West |  435
Northwest |  6116
Northeast |  1137
Southeast |  265
North |  2325
South |  114
Middle |  131
Southwest |  173
In total |  10696

###  2.4 Model evaluation

The prediction accuracy is the significant parameter to access the overall
performance in the classification of the tobacco growing area and is defined
as

$\displaystyle{{P}_{\text{a}}}=\frac{{{n}_{r}}}{{{N}_{t}}},$ (17)

where ${{n}_{r}}$ is the number of samples that are predicted correctly and
${{N}_{t}}$ is the number of samples that are used for prediction. In this
study, the number of training set and testing set is 10696 and 2674,
respectively.

##  3 Results and Discussion

The number of feature nodes, windows and enhanced nodes were tested and the
BLS network structure was established according to the experimental results.
The experimental results of incremental learning and the comparative
experimental results before and after model optimization were analyzed. The
results of the model were compared with that of other traditional algorithms.
All simulations were implemented on Windows 10 operating system using Matlab
R2014a, which is running on a laptop with Intel (R) Core (TM) i5-4210U CPU
1.70 GHz and 2.40 GHz RAM.3.1.

###  3.1 Selection of input node number of network

The number of feature nodes, windows and enhancement nodes affects the
prediction performance and computing speed of the model. In order to set
appropriate parameter values, experiments on the training set and testing set
were conducted and results were plotted. The number of windows is denoted by
${{N}_{2}}$, the number of feature nodes in per window is denoted by
${{N}_{1}}$, and the number of enhancement nodes is denoted by ${{N}_{3}}$.
When ${{N}_{3}}$ is fixed as 1000, the model is trained at different values of
${{N}_{1}}$ and ${{N}_{2}}$ for training set. The prediction accuracy and
training time of the model are shown in Figure 4 and Figure 5, respectively.

Figure 4: Prediction accuracy of training set along with ${{N}_{1}}$ and
${{N}_{2}}$ (${{N}_{3}}$=1000).

It can be seen from Figure 4 that, the prediction accuracy of the model
increases along with the increase of the number of windows and feature nodes.
It increases rapidly (from 0.85 to 0.90) when the values of ${{N}_{1}}$ and
${{N}_{2}}$ are relatively small (less than 50). However, it increases slowly
(from 0.90 to 0.91) when the values of ${{N}_{1}}$ and ${{N}_{2}}$ are
relatively large (more than 50). The prediction accuracy achieves as high as
0.9121 when the two parameters are 140 and 150 respectively.

Figure 5: Training time of training set along with ${{N}_{1}}$ and ${{N}_{2}}$
($N_{3}$=1000).

As seen in Figure 5, the training time of the model increases along with the
increase of the parameters. It increases slowly when the values of ${{N}_{1}}$
and ${{N}_{2}}$ are small (less than 110) but increases rapidly when the
values of ${{N}_{1}}$ and ${{N}_{2}}$ are relatively large (more than 110).
The training time gets as high as 131.76 seconds when the two parameters are
140 and 150, respectively.

Set ${{N}_{1}}$=20 and ${{N}_{2}}$=10, the training set was trained with
different number of enhancement nodes, and then the model was tested with the
testing set. The prediction accuracy and training time obtained were shown in
Figure 6. For comparison, the experimental results were plotted while the
parameters were set as ${{N}_{1}}$=140 and ${{N}_{3}}$=1000, which is shown in
Figure 7.

Figure 6: Prediction accuracy and training time of dataset along with
${{N}_{3}}$ (${{N}_{1}}$=20, ${{N}_{2}}$=10). Figure 7: Prediction accuracy
and training time of dataset along with ${{N}_{2}}$ ($N_{1}$=140,
$N_{3}$=1000).

As seen in Figure 6, the curves of prediction accuracy for training set and
testing set increase rapidly (from 0.815 to 0.921 and from 0.820 to 0.895)
when ${{N}_{3}}$ increases from 100 to 2000, then increase slowly along with
the increase of ${{N}_{3}}$, and achieves the highest value with 0.9499 and
0.9088 when ${{N}_{3}}$ is equal to 13000. However, the curve of training time
shows the trend of exponential growth along with the increase of ${{N}_{3}}$,
and the training time is as high as 25.69 seconds when ${{N}_{3}}$ is equal to
13000.

As seen in Figure 7, the curves of prediction accuracy keep the trend of
gradually increase (from 0.924 to 0.941 and from 0.901 to 0.912) when
${{N}_{2}}$ increases from 10 to 150, and the curve of training time for
training set shows the trend of exponential increase (from 5.2 seconds to 131
seconds).

By comparing the results of Figure 6 and Figure 7, it can be seen that the
training time at the highest prediction accuracy obtained in Figure 6 is 25.69
seconds, while it reaches as high as 103.98 seconds at the highest prediction
accuracy obtained in Figure 7, which is 78.29 seconds higher than the former.
The reason lies in that the feature nodes need to be generated by sparse
method and windows need to be iterated at the same time, which are time-
consuming, while the generation of enhanced nodes only needs orthogonal
normalization iteration, which is far less time-consuming than the generation
of feature nodes.

Based on the above analysis, the number of windows in the network is set as
10, the number of feature nodes in each window is set as 20, and the number of
enhanced nodes is set as 13000.

###  3.2 Experimental results and analysis of incremental learning

The incremental learning method can be used to analyze the data when the
sample size of the data or the dimension of the data increases. After
determining the network structure, incremental learning is conducted by
increasing the number of feature nodes and enhancement nodes in the network.
The experimental results are shown in Table 2.

Table 2: Step results of dataset classification with incremental learning method. Number of | Number of |  $P_{a}$ of |  $P_{a}$ of | Extra step | Total training
---|---|---|---|---|---
feature nods | enhancement nods | training set | testing set | time (s) | time (s)
200 | 13000 | 0.94353 | 0.90464 | 11.529 | 11.529
200→210 | 13000→14000 | 0.9469 | 0.90726 | 2.4402 | 13.9692
210→220 | 14000→15000 | 0.94979 | 0.90814 | 2.7496 | 16.7188
220→230 | 15000→16000 | 0.95409 | 0.91090 | 2.9349 | 19.6537
230→240 | 16000→17000 | 0.95559 | 0.91281 | 3.2067 | 22.8604

As seen from Table 2, with the increase of the number of feature nodes and
enhancement nodes each time, the prediction accuracy of the model for training
set and testing set improves (from 0.94353 to 0.95559, and from 0.90464 to
0.91281, respectively). The training time is 11.529 seconds before adding
nodes, and the additional one-step training time after adding nodes is very
small. When 10 feature nodes and 1000 enhancement nodes are added each time,
the additional one-step training time increases from 2.4402 seconds to 3.2067
seconds during four times of incremental learning, and the cumulative training
time is only 22.8604 seconds. The reason lies in that ridge regression method
adopted in the first calculation of node weights, which may take some time to
iterate continuously. However, incremental learning does not need to retrain
the model, and its calculation only involves the calculation of matrix
multiplication of newly added nodes, so the extra one-step training time is
short.

###  3.3 Tobacco comparative experimental results of BLS and FBLS

The BP algorithm is adopted to update the weight of network in the BLS method
in this paper, and two incremental approaches are applied, so three modes are
used for simulation, which are non-incremental learning, increasing the number
of enhancement nodes, and increasing the number of both enhancement nodes and
feature nodes. The experimental results are shown in Table 3. In the table,
BLS_BP represents the BLS method that adopts BP algorithm without incremental
learning, BLS_EN represents the BLS method that increases the number of
enhancement nodes, and BLS_EF represents the BLS method that simultaneously
increases the number of enhancement nodes and feature nodes.

As seen from the Table 3, for all the three approaches, the prediction
accuracy of FBLS method is always higher than that of the BLS method, which
indicates that the feature extraction of TS fuzzy subsystem can effectively
improve the prediction performance of the model. However, the training time of
FBLS method is always larger than that of BLS method. For example, the
training time of FBLS_EN and BLS_EN methods is 128.3030 seconds and 50.8731
seconds respectively. The reason lies in that there are more parameters used
for feature extraction in TS system, and the update of weight process takes
longer time. Therefore, FBLS sacrifices longer time for higher prediction
accuracy.

Table 3: Comparison of experimental results of BLS and FBLS-TS. Method |  $P_{a}$ of |  $P_{a}$ of | Training time | Testing time
---|---|---|---|---
| training set | testing set | (s) | (s)
BLS_BP | 0.8548 | 0.8512 | 57.3322 | 0.2156
FBLS_BP | 0.8840 | 0.8792 | 167.2500 | 0.6922
BLS_EN | 0.9514 | 0.9061 | 50.8731 | 0.6044
FBLS_EN | 0.9872 | 0.9559 | 128.3030 | 0.7028
BLS_EF | 0.9413 | 0.9046 | 91.1280 | 0.9171
FBLS_EF | 0.9854 | 0.9540 | 187\. 8224 | 0.7508

###  3.4 Advantages of the model in training time

BLS and FBLS-TS methods are compared with traditional classification
algorithms (SVM and GA-SVM), artificial neural network (ANN) and convolution
neural network (CNN) in terms of prediction accuracy and training time. The
experimental results are shown in Table 4.

Table 4: Comparison of classification results with different models. Models | $P_{a}$ | Training time (s)
---|---|---
SVM | 0.8637 | 3068.628
GA-SVM | 0.9068 | 4122.037
ANN | 0.8534 | 4091.184
CNN | 0.9316 | 6833.06
BLS | 0.9061 | 50.8731
FBLS-TS | 0.9559 | 128.303

It can be seen from Table 4 that the prediction accuracy of BLS is 0.9061,
which is 0.0424 and 0.0527 higher than SVM and ANN, respectively, and it is
about the same as that of GA-SVM, which proves that the prediction performance
of BLS is superior to the traditional classification algorithms and ANN.
However, the training time of BLS is just 50.8731 seconds, but the training
time of the other three algorithms are all over 3068 seconds, which shows that
BLS has a great advantage in training speed, which is obviously superior to
the existing traditional classification algorithms. Compared with CNN,
although BLS is 0.0255 lower than CNN in terms of prediction accuracy, its
training time is 137 times shorter than CNN. It lies in that CNN is a deep
network structure with a large number of hyper-parameters, which makes the
training process extremely time$-$consuming. However, BLS has no deep network
structure, and expands in the direction of network width with few parameters,
and it updates the model through incremental learning, which saves a lot of
time for training process.

The prediction accuracy of the FBLS-TS model is improved to 0.9559, which is
0.0498 higher than that of BLS, and its training time is 128.303 seconds with
77 seconds longer than that of BLS. This is because FBLS-TS model achieves
higher prediction accuracy by sacrificing training time.

According to the above analysis, the BLS method before and after improvement
has a great advantage over other methods in terms of training time, especially
the improved FBLS-TS method achieves better effect than other algorithms in
terms of both recognition rate and training speed.

##  4 Conclusions

The FBLS-TS model is proposed to discriminate the tobacco growing area based
on the NIR spectra in this study. The experimental results demonstrate that
the FBLS-TS model not only achieves highest prediction accuracy but also has a
great advantage in training time, which can address the time$-$consuming
problems happened in deep network during training process for big data and
retraining the network due to adding more new data. It lies in that the TS
fuzzy subsystem adaptively extracted the effective features from data which
greatly improves the prediction accuracy, and it adopts the pseudo-inverse
calculation and incremental learning method which greatly speed up the
training speed of the model.

The performance of the FBLS-TS model is superior to the other traditional
algorithms on the classification issue of tobacco growing area. This study not
only can be used for classification problems in tobacco industry but also can
be extended for application in agriculture, medicine, food, petrochemical,
environment and other fields.

## Acknowlegments

This research was funded by the Key Program for Science and Technology of
China National Tobacco Corporation (Grant No. 110202102038), National Natural
Science Foundation of China (Grant No. 62103068; Grant No. 51978111), and
Chongqing Municipal Education Commission (Grant No. KJQN202100745).

## References

  * [1] D. Wang, F. Tian, Z. Zhu, W. Pan, Automatic Prediction of Leave Chemical Compositions Based on Nir Spectroscopy with Machine Learning, International Journal of Robotics and Automation 34 (4) (2019) 391–396.
  * [2] D. Wang, F. Tian, S. X. Yang, Z. Zhu, Intelligent Estimate of Chemical Compositions Based on NIR Spectra Analysis, in: 2017 IEEE International Conference on Information and Automation (ICIA), Macau, China, 2017, pp. 472–477.
  * [3] M. S. Rahman, N. F. Ahmed, M. Ali, M. M. Abedin, M. S. Islam, Determinants of Tobacco Cultivation in Bangladesh, Tobacco Control 29 (6) (2020) 692–694.
  * [4] H. Liu, Y. Zhang, X. Zhou, X. You, Y. Shi, J. Xu, Source Identification and Spatial Distribution of Heavy Metals in Tobacco-growing Soils in Shandong Province of China with Multivariate and Geostatistical Analysis, Environmental Science and Pollution Research 24 (6) (2017) 5964–5975.
  * [5] H. Wu, Q. Liu, J. Ma, L. Liu, Y. Qu, Y. Gong, S. Yang, T. Luo, Heavy Metal (loids) in Typical Chinese Tobacco-growing Soils: Concentrations, Influence Factors and Potential Health Risks, Chemosphere 245 (2020) 125591.
  * [6] B. Chen, W. Xing, D. Lu, Application of Tobacco Leaf Origin Identification with X - ray Fluorescent Based on Discriminant Analysis, Journal of Jiangsu University (Natural Science edition) 5 (2015) 545–549.
  * [7] X. Zhu, Y. Zong, Y. Li, J. Xie, Identification of Flue-cured Tobacco from Different Countries with Electronic Noses, Tobacco Science & Technology 3 (2008).
  * [8] S. Williams, S. Hubbard, K. Reinhard, S. Chaves, Establishing Tobacco Origin from Pollen Identification: An Approach to Resolving the Debate, Journal of Forensic Sciences 59 (2014) 1642–1649.
  * [9] S. Wu, T. Liu, J. Ge, Y. Sha, Chemical Composition - Naive Bayesian Classification Algorithm for Tobacco Leaf Origin Pattern Recognition, Journal of Henan Normal University (Natural Science Edition) 01 (2018).
  * [10] L. Wang, J. Zang, Q. Zhang, Z. Niu, G. Hua, N. Zheng, Action Recognition by an Attention-Aware Temporal Weighted Convolutional Neural Network, Sensors 18 (2018) 1979.
  * [11] B. Richter, M. Rurik, S. Gurk, O. Kohlbacher, M. Fischer, Food Monitoring: Screening of the Geographical Origin of White Asparagus Using FT-NIR and Machine Learning, Food Control 104 (2019) 318–325.
  * [12] T. Wu, I. Tung, H. Hsu, C. Kuo, J. Chang, S. Chen, C. Tsai, Y. Chuang, Quantitative Analysis and Discrimination of Partially Fermented Teas from Different Origins Using Visible/Near-Infrared Spectroscopy Coupled with Chemometrics, Sensors 20 (2020) 5451.
  * [13] H. Li, D. Jiang, J.Cao, D. Zhang, Near-Infrared Spectroscopy Coupled Chemometric Algorithms for Rapid Origin Identification and Lipid Content Detection of Pinus Koraiensis Seeds, Sensors 20 (2020) 4905.
  * [14] Q. Xiao, X. Bai, P. Gao, Y. He, Application of Convolutional Neural Network-Based Feature Extraction and Data Fusion for Geographical Origin Identification of Radix Astragali by Visible/Short-Wave Near-Infrared and Near Infrared Hyperspectral Imaging, Sensors 20 (2020) 4940.
  * [15] F. Liu, W. Wang, T. Shen, J. Peng, W. Kong, Rapid Identification of Kudzu Powder of Different Origins Using Laser-induced Breakdown Spectroscopy, Sensors 19 (2019) 1453.
  * [16] D.Wang, L. Xie, S. X. Yang, F. Tian, Support Vector Machine Optimized by Genetic Algorithm for Data Analysis of Near-Infrared Spectroscopy Sensors, Sensors 18 (10) (2018) 3222.
  * [17] C. Zhu, H. Gong, Z. Li, C. Yu, Application of High Dimensional Feature Grouping Method in Near-infrared Spectra of Identification of Tobacco Growing Areas, in: 2016 3rd International Conference on Information Science and Control Engineering (ICISCE), Beijing, China, 2016, pp. 230–234.
  * [18] Y. Wang, X. Ma, Y. Wen, C. Yu, L. Wang, L. Zhao, J. Li, Tobacco Quality Analysis of Producing Areas of Yunnan Tobacco Using Near-infrared (nir) Spectrum, Spectroscopy and Spectral Analysis 33 (2013) 78–80.
  * [19] S. Zhou, B. Tan, Electrocardiogram Soft Computing Using Hybrid Deep Learning CNN-ELM, Applied Soft Computing 86 (2020) 105778.
  * [20] S. J. Fong, G. Li, N. Dey, R. G. Crespo, E. Herrera-Viedma, Composite Monte Carlo Decision Making Under High Uncertainty of Novel Coronavirus Epidemic Using Hybridized Deep Learning and Fuzzy Rule Induction, Applied Soft Computing 93 (2020) 106282.
  * [21] M. Mittal, L. M. Goyal, S. Kaur, I. Kaur, A. Verma, D. J. Hemanth, Deep Learning Based Enhanced Tumor Segmentation Approach for MR Brain Images, Applied Soft Computing 78 (2019) 346–354.
  * [22] Y.-S. Su, C.-F. Ni, W.-C. Li, I.-H. Lee, C.-P. Lin, Applying Deep Learning Algorithms to Enhance Simulations of Large-scale Groundwater Flow in IoTs, Applied Soft Computing 92 (2020) 106298.
  * [23] S. Afrasiabi, M. Afrasiabi, B. Parang, M. Mohammadi, Designing a Composite Deep Learning Based Differential Protection Scheme of Power Transformers, Applied Soft Computing 87 (2020) 105975.
  * [24] M. Lu, K. Yang, P. Song, R. Shu, L. Wang, Y. Yang, H. Liu, J. Li, L. Zhao, Y. Zhang, The Study of Classification Modeling Method for Near Infrared Spectroscopy of Tobacco Leaves Based on Convolution Neural Network, Spectroscopy and Spectral Analysis 38 (2018) 3724–3728.
  * [25] M. Lee, Y. Xiong, G. Yu, G. Li, Deep Neural Networks for Linear Sum Assignment Problems, IEEE Wireless Communications Letters 7 (2018) 962–965.
  * [26] D. Wang, F. Tian, S. Yang, Z. Zhu, D. Jiang, B. Cai, Improved Deep CNN with Parameter Initialization for Data Analysis of Near-Infrared Spectroscopy Sensors, Sensors 20 (2020) 874.
  * [27] C. Chen, Z. Liu, Broad Learning System: An Effective and Efficient Incremental Learning System without the Need for Deep Architecture, IEEE Transactions on Neural Networks and Learning Systems 29 (2017) 10–24.
  * [28] C. Chen, Z. Liu, S. Feng, Universal Approximation Capability of Broad Learning System and Its Structural Variations, IEEE Transactions on Neural Networks and Learning Systems 30 (2018) 1191–1204.
  * [29] S. Luo, C. Zhao, Y. Fu, An Intelligent Human Activity Recognition Method with Incremental Learning Capability for Bedridden Patients, in: 2018 15th International Conference on Control, Automation, Robotics and Vision (ICARCV), IEEE, 2018, pp. 1284–1289.


# A _NICER_ Discovery of a Low-Frequency Quasi-Periodic Oscillation in the
Soft-Intermediate State of MAXI J1535–571

A. L. Stevens Department of Physics & Astronomy, Michigan State University,
567 Wilson Road, East Lansing, MI 48824, USA  Anton Pannekoek Institute for
Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The
Netherlands  Department of Astronomy, University of Michigan, 1085 South
University Avenue, Ann Arbor, MI 48109, USA<EMAIL_ADDRESS>P. Uttley  Anton
Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlands  D. Altamirano Physics and Astronomy,
University of Southampton, Southampton SO17 1BJ, UK  Z. Arzoumanian
Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD
20771, USA  P. Bult  Astrophysics Science Division, NASA Goddard Space Flight
Center, Greenbelt, MD 20771, USA  E. M. Cackett Department of Physics &
Astronomy, Wayne State University, 666 West Hancock Street, Detroit, MI 48201,
USA  A. C. Fabian  Institute of Astronomy, University of Cambridge, Madingley
Road, Cambridge CB3 0HA, UK  K. C. Gendreau  Astrophysics Science Division,
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  K. Q. Ha  Applied
Engineering and Technology Directorate, NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA  J. Homan Eureka Scientific, Inc., 2452 Delmer
Street, Oakland, CA 94602, USA  SRON, Netherlands Institute for Space
Research, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands  A. R. Ingram
Department of Physics, Astrophysics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, UK  E. Kara Department of Astronomy,
University of Maryland, College Park, Maryland 20742, USA  J. Kellogg  Applied
Engineering and Technology Directorate, NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA  R. M. Ludlam Department of Astronomy, University of
Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA  J. M. Miller
Department of Astronomy, University of Michigan, 1085 South University Avenue,
Ann Arbor, MI 48109, USA  J. Neilsen Villanova University, Mendel Hall, Room
263A, 800 East Lancaster Avenue, Villanova, PA 19085, USA  D. R. Pasham MIT
Kavli Institute for Astrophysics and Space Research, MIT, 70 Vassar Street,
Cambridge, MA 02139, USA  R. A. Remillard MIT Kavli Institute for Astrophysics
and Space Research, MIT, 70 Vassar Street, Cambridge, MA 02139, USA  J. F.
Steiner MIT Kavli Institute for Astrophysics and Space Research, MIT, 70
Vassar Street, Cambridge, MA 02139, USA  J. van den Eijnden Anton Pannekoek
Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands

(Received 08/09/2018; Revised 09/13/2018; Accepted 09/14/2018)

###### Abstract

We present the discovery of a low-frequency $\approx$ 5.7 Hz quasi-periodic
oscillation (QPO) feature in observations of the black hole X-ray binary MAXI
J1535–571 in its soft-intermediate state, obtained in 2017 September–October
by the _Neutron Star Interior Composition Explorer_ (_NICER_). The feature is
relatively broad (compared to other low-frequency QPOs; quality factor
$Q\approx 2$) and weak (1.9% rms in 3–10 keV), and is accompanied by a weak
harmonic and low-amplitude broadband noise. These characteristics identify it
as a weak Type A/B QPO, similar to ones previously identified in the soft-
intermediate state of the transient black hole X-ray binary XTE J1550–564. The
lag-energy spectrum of the QPO shows increasing soft lags toward lower
energies, approaching 50 ms at 1 keV (with respect to a 3–10 keV continuum).
This large phase shift has similar amplitude but opposite sign to that seen in
_Rossi X-Ray Timing Explorer_ data for a Type B QPO from the transient black
hole X-ray binary GX 339–4. Previous phase-resolved spectroscopy analysis of
the Type B QPO in GX 339–4 pointed toward a precessing jet-like corona
illuminating the accretion disk as the origin of the QPO signal. We suggest
that this QPO in MAXI J1535–571 may have the same origin, with the different
lag sign depending on the scale height of the emitting region and the observer
inclination angle.

accretion, accretion discs – black hole physics – stars: black holes – X-rays:
binaries – X-rays: individual: MAXI J1535–571

††journal: ApJ Letters††facilities: _NICER_ , ADS, HEASARC††software: Astropy
v2.0.3 (Astropy Collaboration et al., 2013), HEASoft v6.24, Jupyter notebooks
v4.4.0, Matplotlib v2.1.2 (Hunter, 2007), NumPy v1.14.0, Python v3.6.4, SciPy
v1.0.0, XSPEC v12.10.0c (Arnaud, 1996)

##  1 Introduction

MAXI J1535–571 is a newly discovered transient X-ray binary that went into
outburst starting on 2017 September 2. It was independently discovered in the
X-rays by the _Monitor of All-sky X-Ray Image_ (_MAXI_ ; Negoro et al. 2017b)
and the _Neil Gehrels Swift Observatory_ (Markwardt et al., 2017; Kennea et
al., 2017). Its multi-wavelength properties strongly suggested that the source
is an accreting stellar-mass black hole, initially seen in the hard spectral
state (Scaringi & ASTR211 Students, 2017a, b; Negoro et al., 2017a; Russell et
al., 2017b; Dincer, 2017; Britt et al., 2017; Xu et al., 2018).

One week following its discovery, X-ray, radio, and sub-mm detections showed
the source brightening and entering the hard-intermediate spectral state
(Nakahira et al., 2017; Kennea, 2017; Palmer et al., 2017; Tetarenko et al.,
2017; Shidatsu et al., 2017b). In the hard-intermediate state, low-frequency
quasi-periodic oscillations (LF QPOs) were detected in the X-ray light curve
by _Swift_ (Mereminskiy & Grebenev, 2017; Russell et al., 2017a) and the
_Neutron Star Interior Composition Explorer_ (_NICER_ ; Gendreau et al. 2017).
Detailed spectral fits to a _NICER_ observation from 2017 September 13 showed
a narrow Fe K emission line and strong reflection features, pointing toward a
nearly-maximal black hole spin and a possible warp in the accretion disk
(Miller et al., 2018).

The source continued to transition into the soft state in 2017 October and
November (Shidatsu et al., 2017a). Although heavily extinguished (neutral
column density $N_{\rm H}\sim$ few$\times 10^{22}$ cm-2), MAXI J1535–571
reached extremely bright flux levels (up to 5 Crab in 2–20 keV flux; Shidatsu
et al. 2017b). The source faded significantly in the soft state through the
end of 2018 April to a 2–10 keV flux of $<$ 8 mCrab (Negoro et al., 2018),
then showed repeated transitions between soft and hard spectral states at
these low flux levels, beginning in early 2018 May (Russell et al., 2018;
Parikh et al., 2018) and presently ongoing.111For up-to-date information on
MAXI J1535–571 from _MAXI_ , see http://maxi.riken.jp/nakahira/1535monitor/.

LF QPOs are associated with Fourier frequencies of $\sim$ 0.1 to tens of Hz in
the X-ray light curves of accreting black hole X-ray binaries. The three types
of LF QPOs seen in black hole X-ray binaries are classified as Types A, B, and
C (Wijnands et al., 1999; Remillard et al., 2002; Casella et al., 2005; Motta
et al., 2011). Over an entire outburst, Type C QPOs are the most commonly
seen, Type B are sometimes seen, and Type A are rarely seen. Type A LF QPOs
have low-amplitude broadband noise and no harmonic in the average power
spectrum, and they appear in the soft-intermediate state of an outburst when
the light curve has rms $\lesssim$ 3% in the 2–60 keV energy band. Type B LF
QPOs are characterized by weak harmonics and weak broadband noise compared to
the QPO amplitude, and appear in the soft-intermediate state when the rms of
the light curve is $\approx$ 3%–5% in 2–60 keV (see, e.g., Heil et al. 2015b).
Type C LF QPOs are characterized by strong harmonics and strong broadband
noise, and appear in the hard-intermediate state when the rms of the light
curve is $\approx$ 5%–30% in 2–60 keV. Type A QPOs have a weak, broad QPO peak
at $\approx$ 8 Hz, Type B QPOs generally have a strong peak at $\lesssim$ 6 Hz
(Motta et al. 2011 found that the Type B QPO peak frequency was anti-
correlated with the source flux in GX 339–4), and Type C QPOs have peaks that
tend to move or shift in the frequency range $\sim$ 0.1–20 Hz in
(anti-)correlation with evolving spectral parameters (see, e.g., Vignarca et
al. 2003; Stiele et al. 2013). The relative width of the QPO can be quantified
by the coherence, $Q$, computed as a ratio of the QPO’s centroid frequency
($\nu_{\text{centroid}}$) to its FWHM; Type A typically have $Q\lesssim 3$,
Type B typically have $Q\gtrsim 6$, and Type C typically have $7\lesssim
Q\lesssim 12$ (Casella et al., 2005).

Theories of the origin of LF QPO variability encompass several physical
explanations that can broadly be categorized as geometric or intrinsic quasi-
periodic variability. Systematic observational studies suggest that the LF QPO
origin is linked to the geometry of the system due to the correlation of the
QPO amplitude (Schnittman et al., 2006; Motta et al., 2015; Heil et al.,
2015a) and the sign of the energy-dependent lags (van den Eijnden et al.,
2017) with the binary orbit inclination. Type B and C QPOs show different
amplitude-inclination correlations (Motta et al., 2015), which may suggest
distinct origins.

Spectral-timing analysis reveals the nature and scaling of causal
relationships between variations in the distinct emission regions in accreting
black hole systems. These components can be identified from the X-ray spectrum
(e.g., disk blackbody and Comptonized emission). Lags may be computed from the
phase of the Fourier cross spectrum to measure the energy-dependent delays of
variability in a narrow energy band of interest with respect to a (broad)
reference band (Vaughan et al., 1994; Vaughan & Nowak, 1997), and relate them
to the components seen in the X-ray spectrum. For a review and description of
some of the key spectral-timing techniques applied to broadband noise, see
Uttley et al. (2014).

In this Letter, we present the power and lag-energy spectral analysis of a
weak Type A/B QPO from MAXI J1535–571 detected early in its outburst during a
soft-intermediate state. We present spectral-timing data from the _NICER_
payload operating on the International Space Station since 2017 June. _NICER_
’s large collecting area, soft X-ray response, $\approx$ 100 eV energy
resolution, $\approx$ 100 ns time resolution, and high count-rate capability
enable unprecedented spectral-timing studies of QPOs, even of the weak feature
we report here. In Section 2, we describe the data and data reduction
procedure. In Section 3, we present the data selection criteria, average power
spectrum, and lag-energy spectrum. The discussion and conclusions are
presented in Section 4.

##  2 Data

The _NICER_ X-Ray Timing Instrument (XTI) delivers photons onto 56 focal plane
modules (FPMs) that are grouped into seven sets of eight, with each set read
out by one measurement/power unit (MPU; Gendreau et al. 2012; Arzoumanian et
al. 2014). FPMs 11, 20, 22, and 60 are not operational, so they do not
contribute events to the _NICER_ dataset. Furthermore, FPMs 14 and 34 were
removed from our analysis.222FPM 14 was removed because it has been found to
intermittently exhibit faults in its processing that can result in spurious
energy spectral features (K. Hamaguchi 2017, private communication). FPM 34 is
exceptionally sensitive to a solar light-leak that results in excessive noise,
to the extent that noise contribution from this single FPM can exceed the rest
of the XTI ensemble. Our “good” selection consists of 50 operational FPMs,
providing an effective area of approximately 1800 cm2 at 1.5 keV.

In this Letter we present the analysis of _NICER_ XTI observations obtained in
2017 September and October. These data correspond to ObsIDs 1050360101 through
1050361020 and 1130360101 through 1130360114. With a near-daily sampling of
one or multiple few-100 s exposures, these observations provide a total
exposure of $\approx$ 122 ks and $\approx$ $1.6\times 10^{9}$ photon counts.
Based on the source spectral and timing properties (see Section 3.2), these
observations correspond to the hard-intermediate and soft-intermediate states
of accreting black holes.

The data were processed with nicerdas version 2018-03-01$\\_$V003. The data
were cleaned using standard calibration with nicercal and standard screening
with nimaketime. We selected events that were detected outside the South
Atlantic Anomaly, more than 40° away from the bright Earth limb, more that 30°
away from the dark Earth limb, less than 54″ offset in pointing, not flagged
as “overshoot” or “undershoot” resets, and triggered the slow signal chain
(EVENT_FLAGS = x1x000). We also applied a “trumpet filtering” to remove known
background events. These processing steps were automatically carried out in
the standard _NICER_ data reduction pipeline.

For our analysis, we split the _NICER_ event list light curves into
uninterrupted sequential 64 s segments. The XTI has an overall photon time-
stamping resolution of 85 ns (1$\sigma$), but we bin light curves to a time
step of $2^{-8}$ s ($\simeq$ 3.9 ms). For 64 s light curve segments, this
gives 16,384 time bins apiece. This time binning gives a Nyquist frequency of
128 Hz; having a Nyquist frequency much higher than our signal frequency lets
us check for deadtime effects on the Poisson noise level (see Section 4).

Though the background contribution is expected to be negligible ($\sim$ 2
count s-1; Keek et al. 2018; Ludlam et al. 2018), we applied a filter to
remove segments that showed signs of flaring or non-source particle
background. We calculated the average count rate in each FPM for each segment,
and if any FPM count rate was more than four times half the interquartile
range away from the median count rate of all FPMs, that segment was discarded.
Out of 1917 total segment of data, 158 segments were discarded in this way. Of
the 433 segments identified at the end of Section 3.2, 428 segments were kept
and 5 segments were discarded. Here we used the measured interquartile range
instead of the rms deviation to be more robust against bias from outliers.
This filtering method also has the effect of removing data taken during
passage through the polar horns, when there tends to be high particle
radiation (as noted in Ludlam et al. 2018).

Due to the strong extinction seen from MAXI J1535–571, there are few
intrinsically low-energy source counts below $\approx$ 1 keV, and the spectrum
there is dominated by photons redistributed from higher energies. Thus, 1 keV
was chosen as the lower energy bound for the soft band in the analysis. At
high energies, we found that the background was greater than or equal to the
source counts above $\approx$ 10 keV, so we used 10 keV as the higher energy
bound for the hard band in the analysis. The broad energy band used for power
spectral analysis and the lag reference band is 3–10 keV so that the results
can be reasonably compared to previous _RXTE_ results, and because LF QPOs are
typically stronger at higher energies. While the _RXTE_ Proportional Counter
Array (PCA) had a lower limit of 2 keV, its effective area varied dramatically
in the 2–3 keV band: the _NICER_ XTI is $\approx$ 11$\times$ greater in
effective area than the PCA at 2 keV and is about equal in effective area
compared to the PCA at $\approx$ 2.7 keV. Because the effective area curve
weights the energy distribution of the average photon counts, we chose to use
3–10 keV for the broadband analysis.

##  3 Results

###  3.1 Source Evolution

Figure 1:  MAXI J1535–571 spectral and timing properties, for continuous 64 s
segments of the data in the early part of the outburst. Time intervals
containing segments used in our analysis are marked with purple shaded areas
(see Section 3.2 and Figures 2 and 3 for the selection criteria). Top panel:
average 3–10 keV count rate (count s-1). Middle panel: spectral hardness ratio
(ratio of the count rates in 7–10 keV to 1–2 keV bands). Bottom panel:
fractional rms in 3–10 keV power spectra, integrated from 1.5–15 Hz.

In Figure 1 we show the 3–10 keV count rate, spectral hardness ratio, and
fractional rms evolution of MAXI J1535–571 in the early part of its outburst
for sequential 64 s segments of the light curve. Detailed study of the power
spectral evolution throughout the observations showed that in the periods with
high rms ($\gtrsim$ 5%) and a higher spectral hardness ratio, the power
spectrum is dominated by a strong Type C QPO with accompanying broadband noise
and a strong harmonic (P.!Uttley et al. 2018, in preparation). Therefore, we
identify those data with the hard-intermediate state (see, e.g., Nowak 1995;
Homan et al. 2001; Belloni et al. 2005). Here we focus on the observations
showing low rms ($<$ 5%) and relatively soft spectra, which do not show strong
variability, but still show larger fractional variability amplitudes than the
canonical soft states observed in other transient black hole X-ray binaries
(rms $\lesssim$3% in 2–20 keV; Belloni et al. 2005; Heil et al. 2015b). We
therefore identify these observations with the soft-intermediate state. The 5%
rms limit was chosen for our analysis due to the gap in rms between 4% and
5.5% that is visible in the lower panel of Figure 2. The frequency range
1.5–15 Hz was used for the rms computation to encompass some of the broadband
noise and the expected LF QPO frequency range without being potentially
contaminated by very-low-frequency noise or deadtime effects on the high-
frequency Poisson noise. Luckily, as we mention in Section 4, we found that
deadtime effects are very minor thanks to _NICER_ ’s modular design and high
throughput.

###  3.2 Identifying the LF QPO

We looked for Type A and B QPOs in the soft-intermediate state observations of
this source. These QPOs stand out above the broadband noise in the averaged
power spectrum and their frequency tends to stay constant in time for a given
flux, which lets us stack segments easily. However, they are only present in a
small hardness window within the soft-intermediate spectral state (Motta et
al., 2011), and Type B QPOs sometimes appear to sharply turn off and on as a
function of time (Belloni et al., 2005). Because QPOs cannot be readily
identified in 64 s sections of the weakly-varying soft-intermediate state
light curves, we searched for evidence of characteristic timing signatures by
sorting the data according to the spectral hardness ratio (the ratio of the
7–10 keV to 1–2 keV count rate). In data from other X-ray binaries, Type A and
B QPOs have been known to appear only in narrow ranges of hardness as the
system evolves during an outburst (e.g., Belloni 2010). We then selected 64 s
segments in a range of hardness ratios and computed the 3–10 keV average power
spectrum; this step was repeated sliding the hardness ratio window by 30%
increments. This analysis revealed the presence of a fairly broad QPO-like
signal at a frequency of $\approx$ 5.7 Hz, which is present for hardness
ratios $\leq$ 0.031 but not for harder power spectra (plotted with Poisson
noise subtracted in the upper panel of Figure 2).

The presence of an $\approx$ 5.7 Hz QPO-like feature can be excluded at a
$3\sigma$ upper limit of $4.480\times 10^{-5}$ Hz-1 in normalization (0.64% in
rms) in the averaged power spectrum of the harder segments (see Figure 2). The
QPO-like peak seems to shift in frequency and weaken significantly for the
lowest hardness ratios ($<$ 0.02), but appears to be relatively stable in
frequency over the hardness ratio range 0.021–0.031. We also checked the power
spectra in that same hardness range in four different count rate windows that
correspond to days 20$-$23, 23$-$26, 26$-$30, and 35$-$39 (with 2017 September
1 as Day 0, as in Figure 1) to assess if the count rate or elapsed time in
outburst affected the variability. The Poisson-noise-subtracted power spectra
for these four time windows are plotted in Figure 3. The $\approx$ 5.7 Hz QPO-
like feature is evident in the first three windows but not in the last window;
there, its $3\sigma$ upper limit in normalization is $1.080\times 10^{-4}$
Hz-1 (1.0% in rms) and there is a stronger possible feature at $\sim$ 12 Hz.
Therefore, for the rest of our study we combined the data from 64 s sections,
selected in the hardness ratio range 0.021–0.031 with an integrated 1.5–15 Hz
rms $<$ 5% and a 3–10 keV count rate greater than 5000 counts s-1. We note
that the QPO-like feature in Figure 3 appears to be intrinsically broad and
stable in frequency for the first three time windows.

Figure 2:  Top: Poisson-noise-subtracted power spectra in 3–10 keV averaged
from 64 s segments of data with 1.5–15 Hz rms $<$ 5%, for two hardness ratio
ranges: 0.031–0.050 (black) and 0.021–0.031 (purple). A broad QPO-like signal
at $\approx$ 5.7 Hz emerges for the softer hardness ratio range. The error
bars show the 1$\sigma$ level uncertainties. Bottom: the rms plotted with
spectral hardness ratio using the same segment data as in Figure 1. The
segments marked with purple diamonds and black squares represent the segments
used in the purple and black averaged power spectra in the upper plot. The
gray dots represent the rest of the segments from Figure 1.  Figure 3:
Poisson-noise-subtracted power spectra in the hardness ratio range
0.021$-$0.031 in four time windows: days 20$-$23, 23$-$26, 26$-$30, and
35$-$39\. Day 0 is 2017 September 1, as in Figure 1. The $\approx$ 5.7 Hz QPO-
like feature is visible in the power spectra in the first three time windows,
but not in the fourth window. The error bars show the 1$\sigma$ level
uncertainties.

There are 428 good 64 s continuous sections of data that fit our spectral
hardness, rms, and count rate criteria, which gives a total exposure of 27.392
ks. The average count rate in the 3–10 keV energy band is 5921 count s-1. For
the reference energy band used in Section 3.4, which also covers 3–10 keV but
uses MPUs 4–6, the average count rate is 2692 count s-1.

###  3.3 Average Power Spectrum

We computed average power spectra (Figure 4) in the energy band 3–10 keV for
the 428 good segments described above, with all XTI MPUs, using our own
code333The analysis software and processed data products will be available in
the GitHub repository https://github.com/abigailStev/MAXIJ1535_QPO after the
publication of this Letter. (see van der Klis 1989 for an overview of Fourier-
domain methods). The power spectra were geometrically re-binned in increments
of 6% of the frequency. We fitted the power spectra over the range 0.03–128 Hz
with four Lorentzians: two model the broadband noise, one the QPO fundamental,
and one the weak QPO harmonic. From these fits, we found the
$\nu_{\text{centroid}}$, FWHM, and normalization of the Lorentzian components.
The Poisson noise was modeled with a power-law with a slope of 1 for fitting
$\nu P_{\nu}$. The best-fitting parameter values are listed in Table 1. The
errors represent the 90% confidence region computed from the xspec MCMC
routine.

Figure 4:  Power spectrum in the energy band 3–10 keV with the Poisson noise
subtracted, same as the purple power spectrum in Figure 2. The data (purple
points), total model (solid black), and Lorentzian models for the QPO (solid
blue), QPO harmonic (dash-dotted black), and broadband noise (dotted red) are
plotted in the upper panel. The lower panel shows a ratio of the total model
to the data.

The QPO-like feature has a centroid frequency of
$\nu_{\text{centroid}}=5.72^{+0.04}_{-0.06}$ Hz. It is relatively broad, and
the fitted width is influenced by the bump on the left of the peak at $\sim$ 4
Hz and the scatter on the harmonic at $\gtrsim$ 9 Hz. Thus, we heuristically
determined and fixed the quality factor at $Q=2$. The corresponding FWHM of
the feature is 2.86 Hz. 444If the FWHM of the QPO is a free-fitting parameter,
we find $\text{FWHM}=3.0\pm 0.2$ Hz, giving $Q=1.9$. The QPO-like feature
contributes a 3–10 keV rms of 1.9% and a 5–10 keV rms of 2.9%. Although this
feature is quite broad compared to Type C QPOs, it is distinct from the
broadband noise and possesses a clear — if weak — harmonic. Therefore, we will
refer to it as a QPO for the remainder of this Letter. The closest analogy to
this feature appears to be the QPOs originally identified by Wijnands et al.
(1999) and Homan et al. (2001) as Type A-I and Type A-II QPOs in the black
hole transient XTE J1550–564; Type A-I QPOs were later identified by Casella
et al. (2005) as a form of Type B QPO.

The Type A-I (later Type B) QPO in XTE J1550–564 has a frequency of 5.9 Hz, a
$Q$-value of $Q=2.4$, 3% rms in 2–60 keV, a harmonic at 11 Hz, and it was
present in the soft-intermediate state. The Type A-II QPO in XTE J1550–564 has
a frequency of 8.5 Hz, a $Q$-value of $Q=2.2$, 2% rms in 2–60 keV, no
harmonic, and it was also present in the soft-intermediate state. This MAXI
J1535–571 QPO has a similar fundamental frequency and harmonic frequency to
the Type A-I, but Q-value and rms more similar to the Type A-II (though it may
not be inconsistent with the Type A-I when considering differences in the
energy-dependent effective area of the _RXTE_ PCA and _NICER_). Based on these
shared characteristics, the QPO could be a Type A or Type B, so we refer to it
as a Type A/B.

Component | Parameter | Value | Notes
---|---|---|---
Poisson noise | norm. ($\times 10^{-4}$ Hz-1) | $3.357\pm 0.003$ |
BBN1 | FWHM (Hz) | $0.3^{+0.2}_{-0.02}$ |
BBN1 |  $\nu_{\text{centroid}}$ (Hz) | $1\times 10^{-22}$ | [1]
BBN1 | norm. ($\times 10^{-4}$ Hz-1) | $8^{+5}_{-1}$ |
BBN2 | FWHM (Hz) | $2.8^{+0.1}_{-0.2}$ |
BBN2 |  $\nu_{\text{centroid}}$ (Hz) | $0.5^{+0.2}_{-0.1}$ |
BBN2 | norm. ($\times 10^{-3}$ Hz-1) | $1.1\pm 0.1$ |
QPO | FWHM (Hz) | 2.86 | [2]
QPO |  $\nu_{\text{centroid}}$ (Hz) | $5.72^{+0.04}_{-0.06}$ |
QPO | norm. ($\times 10^{-4}$ Hz-1) | $4.1\pm 0.2$ |
Harmonic | FWHM (Hz) | $2^{+2}_{-0.5}$ |
Harmonic |  $\nu_{\text{centroid}}$ (Hz) | $11.1^{+0.3}_{-0.4}$ |
Harmonic | norm. ($\times 10^{-5}$ Hz-1) | $3^{+2}_{-1}$ |
Table 1:  Best-fitting parameter values for the power spectrum with Poisson
noise, with $\chi^{2}=102.25$ for 93 Degrees of Freedom. The errors on the
parameters represent the 90% confidence region. Final column notes: [1]:
Frozen at $1\times 10^{-22}$ to represent the very-low-frequency broadband
noise and eliminate degeneracy between the two broadband noise components.
[2]: Fixed at half the QPO centroid value.

###  3.4 Lag-energy Spectrum

The lag-energy spectrum, computed from the average cross spectrum, measures by
how much the variability in many narrow energy bands leads or lags the
variability in a broad ‘reference’ energy band (see Uttley et al. 2014 for an
overview of lag-energy spectra with examples and recipes). To ensure that the
narrow-band and broad reference band light curves are truly independent, we
use the detector information stored in the event lists as a selection
criterion. We used MPUs 0–3 (inclusive) when extracting narrow bands and MPUs
4–6 (inclusive) for extracting the reference band. This division of MPUs was
chosen to optimize the signal-to-noise in the narrow bands while maintaining
enough signal in the reference band. As noted in Section 2, we removed FPMs
11, 14, 20, 22, 34, and 60, so that there are 50 total FPMs: 27 were used for
the channels of interest, and 23 for the reference band. To further improve
signal-to-noise, we binned coarser than the intrinsic energy resolution of the
data and achieved 50 narrow-band ‘channels of interest’ across the energy
range 1–10 keV. The channels of interest are approximately equally spaced on a
log scale (see Figure 5).

To compare the lag-energy spectra of the broadband noise and the QPO
components, cross-spectra were averaged and resulting phase and time-lags
calculated across three integrated frequency ranges: 0.1–2.0 Hz (intrinsic
broadband noise), 4.29–7.15 Hz (QPO), and 9.6–12.6 Hz (harmonic). The
frequency ranges used for the QPO and harmonic were chosen to be the FWHM
centered on the centroid for the fitted Lorentzian model in the power spectrum
(Table 1; Figure 4). The lag-energy spectra for the broadband noise and QPO
are shown in Figure 5.

Figure 5:  Lag-energy spectra computed over the frequency ranges 0.1–2.0 Hz
(red: broadband noise) and 4.29–7.15 Hz (blue: QPO) plotted in phase
normalized to 1. One “cycle” refers to the characteristic timescale of the
variability, that is, 0.175 s for the QPO quasi-period or 1 s for the middle
of the frequency range used for the broadband noise. We see soft lags for the
broadband noise and QPO.

The broadband noise lags show a “soft lag,” in which the variability at softer
energies lags behind the variability at harder energies. The QPO lags also
show a definitive soft lag shape, with a flattening above 4 keV.

The shape of the QPO lag-energy spectrum is reminiscent of the inflected lag-
energy spectrum shape (like a broken log-linear shape) observed for the Type B
QPO of GX 339–4 (Stevens & Uttley, 2016), but with the opposite sign of lags.
We note for completeness that we also computed the lags over the same
frequency range of the segments that were excluded due to their hardness
values, as a check that this soft lag is due to the presence of the QPO. While
it is not possible to say much about lag features due to signal-to-noise
limitations, the lag-energy spectrum does not appear to follow the same
trends. The harmonic lag-energy spectrum, computed for the frequency range
9.6$-$12.6 Hz, has better signal-to-noise than the power spectrum of the
harder segments (black points in Figure 2), but still has larger error bars
than the QPO lag-energy spectrum. There is significant scatter and it does not
show any trends or features.

##  4 Discussion and Conclusion

We have reported the discovery and initial spectral-timing analysis of a QPO
feature in the soft-intermediate state of the 2017–2018 outburst of the new
candidate black hole transient MAXI J1535–571. This QPO has a centroid
frequency of $5.72^{+0.04}_{-0.06}$ Hz and is relatively broad (Q $\approx$
2). It is also weak, with an integrated fractional rms of 1.9% in 3–10 keV and
2.9% in 5–10 keV. The QPO shows large soft lags with respect to the 3–10 keV
reference band of about 29% of a QPO cycle at 1.3 keV ($\sim$ 50 ms), and the
broadband noise shows 5% phase lags at 1.3 keV (also $\sim$ 50 ms). The
broadband noise shows a smoother lag-energy spectrum compared to the inflected
shape of the QPO lag-energy spectrum. The QPO lags noticeably flatten above
$\sim$ 4 keV.

Soft lags (i.e., positive at softer energies) are also observed for the
similar QPOs in XTE J1550–564, albeit across the harder energy range 2–48 keV
(Wijnands et al., 1999). The lags that we see are approximately of the same
magnitude and the same sign as the lags in Wijnands et al. (1999) when
comparing the overlapping energy range of _NICER_ and _RXTE_.

The QPO analyzed in this Letter has some notable differences with the clear
Type B QPO in GX 339–4 analyzed with lag-energy spectra and phase-resolved
spectroscopy in Stevens & Uttley (2016). There, the QPO had a fractional rms
of $\sim$ 14% in 5–10 keV, and the broadband noise was much weaker relative to
the QPO, but larger in normalization than we see here. The QPO also showed
_hard_ lags of about 10 ms ($\approx$ 30% in normalized phase). It is notable
that the inflected shapes of the lag-energy spectra observed in MAXI J1535–571
and GX 339–4 are similar, but with opposite sign. The inflection may be due to
the phase offset between sinusoidally varying spectral components, as shown in
simulations in Stevens & Uttley (2016) and Stevens (2018). Thus the inflected
shape of the lag-energy spectrum in MAXI J1535–571 would then be linked to a
phase offset between the peaks in power-law emission and the modulation of the
disk blackbody spectrum, where the disk component would _lag_ rather than lead
the power-law emission. In GX 339–4 we inferred that there was a large 30%
phase shift from spectral fitting of the Wien tail, while with _NICER_ data of
MAXI J1535–571 we can probe low-enough energies to measure this lag directly.

The proposed mechanism for the Type B QPO in GX 339–4 was a large-scale-height
precessing power-law emitting region such as the base of the jet (Stevens &
Uttley, 2016). There, the phase lead of the disk relative to power-law
variations in GX 339–4 was attributed to a geometry where the approaching side
of the disk is illuminated by the jet base, leading to enhanced (blue-shifted)
disk emission, before the jet base points toward the observer leading to
enhanced power-law emission (e.g., due to beaming effects, optical depth
effects, and other aspects of the jet emission mechanism). GX 339–4 is
probably a low-inclination system,555The binary inclination angle of GX 339–4
has been constrained to $37^{\circ}<i\lesssim 60^{\circ}$ from X-ray
(Zdziarski et al., 1998) and optical (Heida et al., 2017) observations. Recent
spectral analysis by Wang-Ji et al. (2018) estimates $i\approx 40^{\circ}$
from spectral modeling. while MAXI J1535–571 may be a high-inclination system
(suggested by spectral fits of the iron line in the hard-intermediate state;
Xu et al. 2018; Miller et al. 2018). The opposite lag sign in MAXI J1535–571
could then suggest that the maximum in power-law emission is seen when the jet
is pointing away from the observer (e.g., if solid angle effects dominate over
beaming effects for the expected larger offset of the jet axis to the line of
sight).

We also note that XTE J1550–564, which shows LF QPOs with similar properties
to those seen in MAXI J1535–571, is a high-inclination (more edge-on) source
($74^{\circ}.7\pm 3^{\circ}.8$; Orosz et al. 2011) with the jet and binary
axes in alignment (Steiner & McClintock, 2012). However, XTE J1550–564 also
shows “normal” Type B QPOs (Homan et al., 2001) that show the same lag sign as
in GX 339–4 (i.e., opposite to what the same source shows for the different
Type A QPOs; Remillard et al. 2002). This suggests that the sign of the lag
may instead (or also) be due to the scale height of the emitting region, which
then affects whether optical depth effects are important when considering
whether the power-law flux variation lags or leads the soft flux variation.
The changing sign of the lag in the same source is then linked to a change in
the geometry (scale height) of the emitting region (as it is not possible for
the inclination of a source to rapidly change). The sign of the lag in MAXI
J1535–571, and the _time_ lag being $\sim$ 50 ms for both the QPO and the
broadband noise, could therefore mean that there is a compact, self-
interacting emitting region that strongly couples the broadband noise lag to
the QPO lag. Our hypothesis on the emission region and mechanism for this QPO
will be rigorously tested in a follow-up paper using phase-resolved
spectroscopy analysis as in Stevens & Uttley (2016) and Stevens (2018). It is
possible that MAXI J1535–571 also exhibited a “normal” Type B QPO with a
higher $Q$-value and higher rms, and it was missed by _NICER_ due to a lengthy
Sun-avoidance data gap in winter 2017–2018.

It is interesting to note that the fitted Poisson noise level in Section 3.3
is only 1% lower than that expected from the mean count rate (noise level of
$2/\text{countrate}=2/5921=3.378\times 10^{-4}$ in units of fractional rms2
Hz-1). Poisson noise levels are reduced by instrumental deadtime effects,
which suppress the observed noise variance due to the resulting
(anti-)correlations between successive photon counts. The data show that the
fraction of photons lost to deadtime in the _NICER_ detectors is indeed
remarkably small for such a bright source ($\sim$ 16,000 counts s-1 for the
whole bandpass).

Finally, we note that these observations highlight the enormous potential of
_NICER_ for transforming our understanding of accreting compact objects via
spectral-timing methods. At the peak flux levels of the soft-intermediate
state of MAXI J1535–571, the total _NICER_ count rate exceeded 16,000 count
s-1, with very little deadtime and hence minimal spectral distortion.
Furthermore, these large phase lags would not have been observable by the
_RXTE_ PCA due to the cut off in response below $\approx$ 3 keV. The
combination of such large count rates with good energy resolution and a soft
X-ray response is a revolutionary capability. _NICER_ points the way to a
bright future for our understanding of the innermost regions of accreting
compact objects.

The authors appreciate the _NICER_ data reduction pipeline development led by
C.B. Markwardt. A.L.S. is supported by an NSF Astronomy and Astrophysics
Postdoctoral Fellowship under award AST-1801792. D.A. and A.I. acknowledge
support from the Royal Society. E.M.C. gratefully acknowledges support from
NSF CAREER award AST-1351222. A.C.F. acknowledges support from ERC Advanced
Grant 340442. R.M.L. acknowledges funding through a NASA Earth and Space
Science Fellowship. D.R.P. is supported by NASA through an Einstein fellowship
(PF6-170156). J.F.S. has been supported by NASA Einstein Fellowship grant No.
PF5-160144. J.v.d.E. acknowledges support from NWO. This work was supported by
NASA through the _NICER_ mission and the Astrophysics Explorers Program. This
research made use of data and software provided by the High Energy
Astrophysics Science Archive Research Center (HEASARC) and NASA’s Astrophysics
Data System Bibliographic Services.

## References

  * 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
  * Arzoumanian et al. (2014) Arzoumanian, Z., Gendreau, K. C., Baker, C. L., et al. 2014, in Proc. SPIE, Vol. 9144, Space Telescopes and Instrumentation 2014: Ultraviolet to Gamma Ray, 914420
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,  A&A, 558, A33
  * Belloni et al. (2005) Belloni, T., Homan, J., Casella, P., et al. 2005,  A&A, 440, 207
  * Belloni (2010) Belloni, T. M. 2010, States and Transitions in Black Hole Binaries, 53
  * Britt et al. (2017) Britt, C. T., Bahramian, A., & Strader, J. 2017, ATel , 10816, ads
  * Casella et al. (2005) Casella, P., Belloni, T., & Stella, L. 2005,  ApJ, 629, 403
  * Dincer (2017) Dincer, T. 2017, ATel , 10716, ads
  * Gendreau et al. (2017) Gendreau, K., Arzoumanian, Z., Markwardt, C., et al. 2017, ATel , 10768, ads
  * Gendreau et al. (2012) Gendreau, K. C., Arzoumanian, Z., & Okajima, T. 2012, in Proc. SPIE, Vol. 8443, Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, 844313
  * Heida et al. (2017) Heida, M., Jonker, P. G., Torres, M. A. P., & Chiavassa, A. 2017,  ApJ, 846, 132
  * Heil et al. (2015a) Heil, L. M., Uttley, P., & Klein-Wolt, M. 2015a,  MNRAS, 448, 3348
  * Heil et al. (2015b) Heil, L. M., Uttley, P., & Klein-Wolt, M. 2015b,  MNRAS, 448, 3339
  * Homan et al. (2001) Homan, J., Wijnands, R., van der Klis, M., et al. 2001,  ApJS, 132, 377
  * Hunter (2007) Hunter, J. D. 2007,  CSE , 9, 90
  * Keek et al. (2018) Keek, L., Arzoumanian, Z., Bult, P., et al. 2018,  ApJ, 855, L4
  * Kennea (2017) Kennea, J. A. 2017, ATel , 10731, ads
  * Kennea et al. (2017) Kennea, J. A., Evans, P. A., Beardmore, A. P., et al. 2017, ATel , 10700, ads
  * Ludlam et al. (2018) Ludlam, R. M., Miller, J. M., Arzoumanian, Z., et al. 2018,  ApJ, 858, L5
  * Markwardt et al. (2017) Markwardt, C. B., Burrows, D. N., Cummings, J. R., et al. 2017, GRB Coordinates Network , 21788, ads
  * Mereminskiy & Grebenev (2017) Mereminskiy, I. A., & Grebenev, S. A. 2017, ATel , 10734, ads
  * Miller et al. (2018) Miller, J. M., Gendreau, K., Ludlam, R. M., et al. 2018,  ApJ, 860, L28
  * Motta et al. (2011) Motta, S., Muñoz-Darias, T., Casella, P., Belloni, T., & Homan, J. 2011,  MNRAS, 418, 2292
  * Motta et al. (2015) Motta, S. E., Casella, P., Henze, M., et al. 2015,  MNRAS, 447, 2059
  * Nakahira et al. (2017) Nakahira, S., Negoro, H., Mihara, T., et al. 2017, ATel , 10729, ads
  * Negoro et al. (2017a) Negoro, H., Kawase, T., Sugizaki, M., et al. 2017a, ATel , 10708, ads
  * Negoro et al. (2017b) Negoro, H., Ishikawa, M., Ueno, S., et al. 2017b, ATel , 10699, ads
  * Negoro et al. (2018) Negoro, H., Mihara, T., Nakahira, S., et al. 2018, ATel, 11568, ads
  * Nowak (1995) Nowak, M. A. 1995,  PASP, 107, 1207
  * Orosz et al. (2011) Orosz, J. A., Steiner, J. F., McClintock, J. E., et al. 2011,  ApJ, 730, 75
  * Palmer et al. (2017) Palmer, D. M., Krimm, H. A., & Swift/BAT Team. 2017, ATel , 10733, ads
  * Parikh et al. (2018) Parikh, A. S., Russell, T. D., Wijnands, R., et al. 2018, ATek s , 11652, ads
  * Remillard et al. (2002) Remillard, R. A., Sobczak, G. J., Muno, M. P., & McClintock, J. E. 2002,  ApJ, 564, 962
  * Russell et al. (2017a) Russell, T. D., Altamirano, D., Tetarenko, A. J., et al. 2017a, ATel , 10899, ads
  * Russell et al. (2018) Russell, T. D., Altamirano, S. R. D., Miller-Jones, J. C. A., et al. 2018, The Astronomer’s Telegram , 11611, ads
  * Russell et al. (2017b) Russell, T. D., Miller-Jones, J. C. A., Sivakoff, G. R., Tetarenko, A. J., & Jacpot Xrb Collaboration. 2017b, ATel , 10711, ads
  * Scaringi & ASTR211 Students (2017a) Scaringi, S., & ASTR211 Students. 2017a, ATel , 10702, ads
  * Scaringi & ASTR211 Students (2017b) Scaringi, S., & ASTR211 Students. 2017b, ATel , 10704, ads
  * Schnittman et al. (2006) Schnittman, J. D., Homan, J., & Miller, J. M. 2006,  ApJ, 642, 420
  * Shidatsu et al. (2017a) Shidatsu, M., Nakahira, S., Negoro, H., et al. 2017a, ATel , 11020, ads
  * Shidatsu et al. (2017b) Shidatsu, M., Nakahira, S., Negoro, H., et al. 2017b, ATel , 10761, ads
  * Steiner & McClintock (2012) Steiner, J. F., & McClintock, J. E. 2012,  ApJ, 745, 136
  * Stevens (2018) Stevens, A. L. 2018, PhD thesis, Universiteit van Amsterdam (Netherlands)
  * Stevens & Uttley (2016) Stevens, A. L., & Uttley, P. 2016,  MNRAS, 460, 2796
  * Stiele et al. (2013) Stiele, H., Belloni, T. M., Kalemci, E., & Motta, S. 2013,  MNRAS, 429, 2655
  * Tetarenko et al. (2017) Tetarenko, A. J., Russell, T. D., Miller-Jones, J. C. A., Sivakoff, G. R., & Jacpot Xrb Collaboration. 2017, ATel , 10745, ads
  * Uttley et al. (2014) Uttley, P., Cackett, E. M., Fabian, A. C., Kara, E., & Wilkins, D. R. 2014,  A&A Rev., 22, 72
  * van den Eijnden et al. (2017) van den Eijnden, J., Ingram, A., Uttley, P., et al. 2017,  MNRAS, 464, 2643
  * van der Klis (1989) van der Klis, M. 1989, in NATO Advanced Science Institutes (ASI) Series C, Vol. 262, NATO Advanced Science Institutes (ASI) Series C, ed. H. Ögelman & E. P. J. van den Heuvel, 27
  * Vaughan et al. (1994) Vaughan, B., van der Klis, M., Lewin, W. H. G., et al. 1994,  ApJ, 421, 738
  * Vaughan & Nowak (1997) Vaughan, B. A., & Nowak, M. A. 1997,  ApJ, 474, L43
  * Vignarca et al. (2003) Vignarca, F., Migliari, S., Belloni, T., Psaltis, D., & van der Klis, M. 2003,  A&A, 397, 729
  * Wang-Ji et al. (2018) Wang-Ji, J., García, J. A., Steiner, J. F., et al. 2018,  ApJ, 855, 61
  * Wijnands et al. (1999) Wijnands, R., Homan, J., & van der Klis, M. 1999,  ApJ, 526, L33
  * Xu et al. (2018) Xu, Y., Harrison, F. A., García, J. A., et al. 2018,  ApJ, 852, L34
  * Zdziarski et al. (1998) Zdziarski, A. A., Poutanen, J., Mikolajewska, J., et al. 1998,  MNRAS, 301, 435