# The evolution of the corona in MAXI J1535$-$571 through type-C quasi-
periodic oscillations with Insight-HXMT

Yuexin Zhang,1,2 Mariano Méndez,1 Federico García,1,8 Shuang-Nan Zhang,2,5
Konstantinos Karpouzas,1,3 Diego Altamirano,3 Tomaso M. Belloni,4 Jinlu Qu,2,5
Shu Zhang,2 Lian Tao,2 Liang Zhang,3,2 Yue Huang,2 Lingda Kong,2,5 Ruican
Ma,2,5 Wei Yu,2,5 Divya Rawat,6,7 and Candela Bellavita9
1Kapteyn Astronomical Institute, University of Groningen, P.O. BOX 800, 9700
AV Groningen, The Netherlands
2Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, People’s Republic of China
3School of Physics and Astronomy, University of Southampton, Southampton, SO17
1BJ, UK
4INAF-Osservatorio Astronomico di Brera, via E. Bianchi 46, I-23807 Merate,
Italy
5University of Chinese Academy of Sciences, Chinese Academy of Sciences,
Beijing 100049, People’s Republic of China
6Department of physics, IIT Kanpur, Kanpur, Uttar Pradesh 208016, India
7Inter-University Center for Astronomy and Astrophysics, Ganeshkhind, Pune
411007, India
8Instituto Argentino de Radioastronomía (CCT La Plata, CONICET; CICPBA; UNLP),
C.C.5, (1894) Villa Elisa, Buenos Aires, Argentina
9Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La
Plata, Paseo del Bosque, B1900FWA La Plata, Argentina
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Type-C quasi-periodic oscillations (QPOs) in black hole X-ray transients can
appear when the source is in the low-hard and hard-intermediate states. The
spectral-timing evolution of the type-C QPO in MAXI J1535$-$571 has been
recently studied with Insight-HXMT. Here we fit simultaneously the time-
averaged energy spectrum, using a relativistic reflection model, and the
fractional rms and phase-lag spectra of the type-C QPOs, using a recently
developed time-dependent Comptonization model when the source was in the
intermediate state. We show, for the first time, that the time-dependent
Comptonization model can successfully explain the X-ray data up to 100 keV. We
find that in the hard-intermediate state the frequency of the type-C QPO
decreases from 2.6 Hz to 2.1 Hz, then increases to 3.3 Hz, and finally
increases to $\sim$ 9 Hz. Simultaneously with this, the evolution of corona
size and the feedback fraction (the fraction of photons up-scattered in the
corona that return to the disc) indicates the change of the morphology of the
corona. Comparing with contemporaneous radio observations, this evolution
suggests a possible connection between the corona and the jet when the system
is in the hard-intermediate state and about to transit into the soft-
intermediate state.

###### keywords:

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

††pagerange: The evolution of the corona in MAXI J1535$-$571 through type-C
quasi-periodic oscillations with Insight-HXMT–The evolution of the corona in
MAXI J1535$-$571 through type-C quasi-periodic oscillations with Insight-HXMT

##  1 Introduction

Black hole low-mass X-ray binaries (LMXBs) consist of a low-mass star and a
central black hole which accretes materials from the star via Roche-lobe
overflow. The inner part of an accreting black hole system includes an
accretion disc and a hot Comptonization region, which is called the corona
(for a review, see Gilfanov, 2010). As the accretion disc is close to the
black hole, the infalling gas releases strong gravitational energy, producing
a multi-temperature blackbody emission that can be detected with a temperature
around 0.3–2.0 keV in the soft X-ray band. Some of the soft photons are
inverse Comptonized into hard X-ray photons in the corona, producing roughly a
power-law continuum in the energy spectrum (see Done et al., 2007, for a
review). A fraction of the Comptonized photons illuminate the inner accretion
disc and can be Compton back-scattered in the disc, producing a relativistic
reflection spectrum with characteristic emission lines, among which the most
prominent feature is a broad iron $K_{\alpha}$ line around 6.4 keV, and a
Compton hump at around 20 keV (Fabian et al., 1989; García et al., 2014). Due
to the relativistic property of the spacetime, the X-ray reflection spectrum
has recently been developed as a powerful tool to test General Relativity (for
a review, see Bambi et al., 2021).

A typical outburst of a black hole X-ray transient follows a ‘q’ path in the
hardness intensity diagram (HID; Homan et al., 2001; Fender et al., 2004),
showing a hysteresis phenomenon between the soft and hard states (Maccarone &
Coppi, 2003). Before the outburst, the source spends time in quiescence with
the X-ray luminosity being several orders of magnitude lower than during
outburst. As the outburst starts, the source enters the low-hard state with a
dominant hard-photon emission; as mass accretion rate from the disc onto the
black hole increases the source at some point quickly (on the order of days)
transits to the high-soft state, during which the X-ray spectrum is disc
dominated and the power-law emission becomes steeper; finally the source
returns back to the low-hard state and eventually to quiescence as mass
accretion rate decreases; a relativistic jet can appear in the transitions
between the hard and the soft states (Fender et al., 2004; Remillard &
McClintock, 2006). Usually two types of relativistic jets are observed in
these sources; the jets can be classified by the radio spectral index and the
jet morphology observed: a small-scale, optically thick, steady jet and an
extended, optically thin, transient jet (for a review, see Fender, 2006). In
the hard state, or even the hard-intermediate state, a steady jet is observed
(e.g., Fender, 2001; Stirling et al., 2001; Russell et al., 2019), while
during the transition from the the hard-intermediate to the soft-intermediate
state, the emission of the steady jet is quenched (Fender et al., 2004;
Russell et al., 2011). Around the transition to the soft state, no steady jet
is observed but a transient jet is launched, which is bright and consists of
discrete relativistic ejecta from the black hole (Mirabel & Rodríguez, 1994;
Tingay et al., 1995; Corbel et al., 2004; Miller-Jones et al., 2012; Russell
et al., 2019).

X-ray emission from black hole binaries (BHBs) shows variability over a broad
range of frequencies, and at well defined frequencies (for a recent review,
see Ingram & Motta, 2019). Quasi-periodic oscillations (QPOs; Van der Klis,
1989; Nowak, 2000; Belloni et al., 2002) are the narrow peaks in the power
density spectrum (PDS). The low-frequency QPOs (LFQPOs) in BHBs can be
classified into three classes, namely type A, B, and C, based on the shape of
the broadband noise in the PDS, the root mean square (rms) amplitude and phase
lags of the QPOs, and the spectral state of the source (Casella et al., 2005;
Huang et al., 2018). High-frequency QPOs, with frequency up to $\sim$ 350 Hz,
in black holes X-ray binaries are less common (e.g., Morgan et al., 1997;
Remillard et al., 1999; Méndez et al., 2013). The dynamic origin of LFQPOs is
explained mainly by geometric effects or instabilities in the accretion flow
(Stella & Vietri, 1998; Tagger & Pellat, 1999; Ingram et al., 2009; You et
al., 2018; Ma et al., 2021). The radiative properties of QPOs, rms amplitude
and phase lags, provide extra information. The 0.1–10 Hz integrated rms
amplitude, for instance, is an indicator of accretion regimes in black hole
transients (Muñoz-Darias et al., 2011). Type-C QPOs often appear when the
source is in the low-hard state and the hard-intermediate state, with central
frequency in the range $\sim 0.1$–$15$ Hz and associated strong broadband
noise in the PDS (Casella et al., 2005). The rms amplitude of the type-C QPOs
can increase with energy, and reach $\sim$ 15% at 100 keV (e.g., Huang et al.,
2018; Bu et al., 2021), indicating that the radiative mechanism of the type-C
QPO has to be related to the corona.

The phase lags are the phase difference of two light curves at different
energies in the same frequency range (Nowak et al., 1999). Hard (positive)
lags are believed to be due to disc fluctuations that propagate through the
corona (Miyamoto et al., 1988), while soft (negative) lags may come from the
reprocessing of the hard photons in the disc when the corona illuminates the
disc (Lee et al., 2001). The rms spectrum and the phase-lag spectrum of the
LFQPOs can indicate how the disc and the Comptonization region influence the
evolution of the radiative properties of the LFQPOs (e.g., Reig et al., 2000;
Belloni et al., 2005; Zhang et al., 2020a; Kong et al., 2020).

Although the corona around a black hole is widely believed to consist of hot
electrons with temperatures up to $\sim$ 100 keV that can give rise to the
Comptonized spectrum (Zdziarski et al., 1996; Życki et al., 1999), the
geometry of the corona is still under debate (Galeev et al., 1979; Haardt &
Maraschi, 1991; Poutanen et al., 2018). It is therefore of great significance
to determine the geometry of the Comptonization region. Models proposed to
constrain the properties of the corona describe either time-averaged energy
spectrum or the X-ray variability. For instance, García et al. (2014) proposed
a time-averaged reflection model called relxill that can measure the height of
the corona assuming a lamppost geometry. Using this model You et al. (2021)
provided information about how the corona moves. Ingram et al. (2019)
presented a relativistic transfer function model, reltrans, that calculates
the X-ray time delays and energy shifts based on reverberation mapping. Both
these models have been applied to observations (e.g., Xu et al., 2018; Wang et
al., 2021), however, neither the time-averaged reflection model nor the
reverberation mapping can explain the radiative properties of the QPOs.

Recently, Karpouzas et al. (2020) developed a time-dependent Comptonization
model based on the idea proposed by Lee & Miller (1998) and Kumar & Misra
(2014). The model of Karpouzas et al. (2020) assumes the QPO to be a
sinusoidal coherent oscillation of the Comptonized X-ray spectrum. The
oscillation comes from the temperature of the corona, the inner disc, and the
external heating source which is necessary to keep the thermal equilibrium of
the system. In this model, the corona is assumed to be spherically symmetric
and partially covering the accretion disc, while blackbody seed photons from
the accretion disc are Compton up-scattered in the corona, leading to hard
lags. Since the corona covers the disc to some extent, a fraction of the hard
photons in the corona can impinge back onto the accretion disc, resulting in
reprocessing of the hard photons and soft lags. The fraction of the feedback
photons to the Comptonized photons, ranging from 0 to 1, is called the
feedback fraction. The steady-state version of this model describes the
Comptonized continuum in the time-averaged energy spectrum, similar, for
instance, to nthcomp (Zdziarski et al., 1996; Życki et al., 1999), while the
time-dependent model can fit simultaneously the fractional rms and phase-lag
spectra. The latest version (Bellavita et al., 2022) of the model (Karpouzas
et al., 2020) incorporates a disc blackbody as the soft-photon source so the
steady-state version of this model is the same as nthcomp. Karpouzas et al.
(2020) first applied the model to the kilohertz QPOs in the neutron star 4U
1636$-$53\. Later on this model was further applied to the LFQPOs of BHBs:
Karpouzas et al. (2021) measured a variable corona in GRS 1915+105 as the
frequency of the type-C QPO changes using RXTE data; García et al. (2021)
applied a dual Comptonization model to the type-B QPO of MAXI J1348+603 with
NICER, showing that there may be two Comptonization regions near the black
hole where the outer region can be the jet. Here we use this time-dependent
Comptonization model to study the evolution of the Comptonization region of
MAXI J1535$-$571 through type-C QPOs with Insight-HXMT.

MAXI J1535$-$571 is a new X-ray transient discovered independently by MAXI/GSC
(Negoro et al., 2017) and Swift/BAT (Kennea et al., 2017) on September 2 2017.
Follow-up X-ray and radio observations suggested that this new source is a
black hole X-ray binary system (Negoro et al., 2017; Russell et al., 2017).
Miller et al. (2018) analyzed the NICER data and proposed that the black hole
spin is $0.994\pm 0.002$ and inclination angle is $67.4^{\circ}\pm
0.8^{\circ}$. A broadband X-ray spectral study also indicated that MAXI
J1535$-$571 has a high spin ($>0.84$) and high inclination angle ($\sim
60^{\circ}$) (Xu et al., 2018). Type-C QPOs were discovered in the low-hard
state of MAXI J1535$-$571, based on Swift, XMM-Newton, and NICER observations
(Stiele & Kong, 2018). Bhargava et al. (2019) showed a tight correlation
between the type-C QPOs and the power-law spectral index, $\Gamma$, using
AstroSat data. Huang et al. (2018) presented a systematic timing study of MAXI
J1535$-$571 with Insight-HXMT, classifying different types of LFQPOs and
confirming the high inclination angle nature of this system. Using part of the
same Insight-HXMT dataset where type-C QPOs were detected, Kong et al. (2020)
analyzed the energy and fractional rms spectra, and gave a picture of a
shrinking and hardening corona. The dataset of MAXI J1535$-$571 with Insight-
HXMT provides broadband energy range as well as strong QPO signals that we can
utilize to fit with the time-dependent Comptonization model of Karpouzas et
al. (2020).

In this paper, we further explore the Insight-HXMT observations of MAXI
J1535$-$571 during the outburst in 2017/2018. We fit simultaneously the
fractional rms, phase-lag, and time-averaged energy spectra of each
observation using the time-dependent Comptonization model of Karpouzas et al.
(2020). This paper is organized as follows: In Section 2 we describe the
reduction of the Insight-HXMT data of MAXI J1535$-$571, and explain how we
generate the power density spectra and the cross spectra for different energy
bands. We also explain how we simultaneously fit the fractional rms, lags, and
time-averaged energy spectra. In Section 3 we show the temporal evolution of
MAXI J1535$-$571, the dependence of the spectral parameters with QPO frequency
and, most importantly, the evolution of the corona. In Section 4 we discuss
our results and propose a physical picture to explain the evolution of corona
and the possible connection between the corona and the jet.

##  2 Observations and Data Analysis

All observations presented here were carried out with Insight-HXMT. The
payload of Insight-HXMT (Zhang et al., 2014) consists of three instruments:
the Low Energy X-ray Telescope (LE), the Medium Energy X-ray Telescope (ME),
and the High Energy X-ray Telescope (HE). The LE, ME, and HE cover,
respectively, the 1–15 keV energy range with 1 ms time resolution, the 5–30
keV energy range with 240 $\mu$s time resolution, and the 20–250 keV energy
range with 4 $\mu$s time resolution.

Insight-HXMT triggered Target of Opportunity (ToO) observations of MAXI
J1535$-$571 in the period 2017 September 6–23. We use the Insight-HXMT Data
Analysis Software (HXMTDAS) v2.04 to extract and reduce the data. We use the
following criteria to establish the good time intervals (GTIs) to filter the
data: we require that the offset of the pointing to the source is less than
0.04∘ and the source is at least 10∘ above Earth horizon. To filter out high-
background intervals, we make a cut on the geometric cutoff rigidity with COR
$>8$ GeV, and we select events that are detected at least 300 s before and
after the passage of the satellite through the South Atlantic Anomaly (SAA).
We finally select all non-blind detectors to generate the products for LE, ME,
and HE, respectively.

We use GHATS
111http://www.brera.inaf.it/utenti/belloni/GHATS_Package/Home.html to compute
the Fourier power density spectrum (PDS) of photons in the LE in 1–10 keV
energy band, ME in 10–35 keV energy band, and HE in 26–100 keV energy band. We
set the time length of each Fast Fourier Transformation (FFT) segment to 25 s
and the time resolution to 3 ms so that the lowest frequency of the PDS is
0.04 Hz and the Nyquist frequency is 167 Hz. We average the PDS within one
observation, subtract the Poisson level by calculating the average power in
the frequency range of $>$ 100 Hz, and normalize the PDS to units of rms2 per
Hz (Belloni & Hasinger, 1990). The background count rate used to convert to
rms units is estimated with the tools LE(ME/HE)BKGMAP in HXMTDAS. Finally, we
apply a logarithmic rebin in frequency to the PDS data such that the size of a
bin increases by $\exp(1/100)$ compared to that of the previous one.

We use XSPEC v12.11.1 (Arnaud, 1996) to fit the PDS with a model consisting of
several Lorentzian functions in the frequency range of 0.5–30.0 Hz. Three
Lorentzians represent the low-frequency broadband noise, high-frequency
broadband noise, and the QPO. Two extra Lorentzians may be needed to represent
the harmonic and the subharmonic of the QPO. Following Huang et al. (2018) and
Kong et al. (2020), the eight observations with type-C QPOs in which all the
three instruments are active, are selected for further study. Our observations
1–8 correspond to observations in Table 2 of Huang et al. (2018) and Table 3
of Kong et al. (2020); all these observations are in the HIMS. Dynamical power
spectra of these observations show that the QPO frequency does not change
significantly during an observation. We regard the best-fitting model for each
observation as a baseline to perform further fitting on the separate energy
bands.

To study the rms spectrum of the QPO, we select the following energy bands,
LE: 1.0–2.5 keV, 2.5–3.5 keV, 3.5–5.0 keV, and 5.0–10.0 keV, ME: 10.0–12.0
keV, 12.0–14.0 keV, 14.0–17.0 keV, and 17.0–35.0 keV, HE: 26.0–30.3 keV,
30.3–35.5 keV, 35.5–55.0 keV, and 55.0–100.0 keV. We use the background in the
corresponding energy band and use the square root of the normalization of the
Lorentzian representing the type-C QPOs to calculate the fractional rms
amplitude, hereafter rms (Belloni & Hasinger, 1990). To calculate the error of
the rms we include the uncertainties of the background count rates.

For each observation, we compute Fourier cross spectra for different energy
bands to calculate the phase-lag spectrum of the QPOs (Vaughan & Nowak, 1997;
Nowak et al., 1999) using the overlapping Good Time Intervals (GTIs) for all
the three instruments. The energy bands are the same as for the rms spectra.
As we do for the power density spectrum, we use the same time resolution of 3
ms and same time segment of 25 s. We regard the lowest energy band (1.0–2.5
keV) as the reference band. The phase lags of the QPOs are calculated from the
average of the real and imaginary parts of the cross spectrum in the frequency
range of $\nu_{0}\pm\text{FWHM}/2$, where $\nu_{0}$ is the centroid frequency
of the QPOs. A positive lag means that the hard photons lag the soft ones.

###  2.1 Fitting the time-averaged, rms, and phase-lag spectra separately

We use XSPEC v12.11.1 (Arnaud, 1996) to fit the time-averaged, rms, and phase-
lag spectra. We adopt the energy bands 2–10 keV (LE), 10–35 keV (ME), and
26-100 keV (HE) for the fitting of the time-averaged spectrum, adding a
systematic error of 1%. Since we assume that the variability at the QPO
frequency comes from the corona, we fit the time-averaged and the time-
dependent spectra using the same model in which we consider the Comptonization
of the seed photons from the accretion disc as the component that drives the
variability. We start by fitting the time-averaged spectrum in each
observation with the model const*tbabs*(diskbb+nthcomp), where the
multiplicative constant is used to account for cross-calibration uncertainties
between three instruments. In the absorption model tbabs, the parameter
$N_{\text{H}}$ represents the hydrogen column density, with the cross-section
and solar-abundance tables of Verner et al. (1996) and Wilms et al. (2000),
respectively. To fit the emission of the accretion disc, we use diskbb
(Mitsuda et al., 1984), which has two parameters: the inner disc temperature,
$kT_{\text{in}}$, and a normalization. To fit the hard Comptonized continuum,
we use nthcomp (Zdziarski et al., 1996; Życki et al., 1999), which has the
following parameters: the corona temperature, $kT_{\text{e}}$, the X-ray
photon index, $\Gamma$, the seed photon temperature, $kT_{\text{s}}$, and a
normalization. We assume that the seed photon source is the disc so when
performing the fitting the $kT_{\text{in}}$ in diskbb and $kT_{\text{s}}$ in
nthcomp are linked.

Fits with the model const*tbabs*(diskbb+nthcomp) show residuals at $\sim$ 6–7
keV and $\sim$ 20–30 keV, suggesting that there may exist a reflection
component (e.g., Xu et al., 2018; Jiang et al., 2020). We therefore add a
relativistic reflection component, relxillCp (García et al., 2014), to the
model that now is const*tbabs*(diskbb+nthcomp+relxillCp). We link
$kT_{\text{e}}$ and $\Gamma$ in relxillCp to $kT_{\text{e}}$ and $\Gamma$ in
nthcomp and, as before, $kT_{\text{s}}$ in nthcomp to $kT_{\text{in}}$ in
diskbb. In relxillCp, the black hole spin, $a_{*}$, inclination angle, $i$,
disc inner radius, $R_{\text{in}}$, logarithm of the ionization parameter,
$\log\xi$, iron abundance, $A_{\text{Fe}}$, and normalization of the model are
free during the fits, whereas the reflection fraction is fixed at $-1$ since
we only consider the reflected emission (the direct emission is given by
nthcomp); we also link the emissivity indices $q_{1}$ and $q_{2}$ during the
fits.

To fit the rms and phase-lag spectra, we utilize the model of Karpouzas et al.
(2020) with a disc blackbody spectrum for the seed-photon source (Bellavita et
al. 2022), which we call vkompthdk in this paper. Since both vkompthdk and
nthcomp come from the same calculation of the Comptonization spectrum, the
time-averaged version of vkompthdk is the same as nthcomp, and both models
share the same parameters. However, the time-dependent version of vkompthdk
has extra parameters, namely the size of the corona, $L$, feedback fraction,
$\eta$, amplitude of the variability of the external heating rate,
$\delta\dot{H}_{\text{ext}}$, QPO frequency and a reference lag. The reference
lag is an additive parameter that accounts for the fact that the reference
band of the lags is arbitrary, and hence the lag spectrum can be shifted up
and down arbitrarily (e.g., García et al., 2021).

###  2.2 A joint fitting

Since our aim is to fit the rms, phase-lag, and time-averaged spectra
simultaneously, we use vkompthdk instead of nthcomp to fit the three spectra
given that the time-averaged version of vkompthdk is the same as nthcomp. We
therefore fit jointly the three spectra in each observation where the energy
ranges we use per instrument are 1–10 keV, 10–35 keV, and 26–100 keV for the
rms and phase-lag spectra, and 2–10 keV, 10–35 keV, and 26–100 keV for the
time-averaged spectra. The model that we use is written as
const*tbabs*(diskbb+vkompthdk+relxillCp). We let $N_{\text{H}}$ and the
normalizations of diskbb and relxillCp free to fit the time-averaged spectrum,
while we fix all of them at 0 to fit the rms and phase-lag spectra. For
vkompthdk the normalization is free when fitting the time-averaged spectrum
and fixed at 1 for the rms and phase-lag spectra (see García et al., 2021).
Besides these parameters, the other parameters that are always free are:
$kT_{\text{in}}$ in diskbb, $kT_{\text{e}}$ and $\Gamma$ in the time-averaged
version of vkompthdk, where $kT_{\text{s}}$ in vkompthdk is linked to
$kT_{\text{in}}$ in diskbb, $q_{1}$, $a_{*}$, $i$, $R_{\text{in}}$, $\log\xi$,
and $A_{\text{Fe}}$ in relxillCp, where $kT_{\text{e}}$ and $\Gamma$ in
relxillCp are linked to the same parameters in vkompthdk. The reflection
fraction is fixed to $-1$ so that relxillCp only provides the reflected
component, while vkompthdk provides the direct emission from the hard
Comptonization component. When we use the time-dependent version of vkompthdk
to fit the rms and phase-lag spectra, the QPO frequency is fixed at the
frequency of the QPO in the PDS in each observation, while $L$, $\eta$,
$\delta\dot{H}_{\text{ext}}$ and the reference lag, which act only on the
phase-lag spectrum as mentioned in 2.1, are free.

We fit the data from the three instruments, i.e. LE, ME, and HE,
simultaneously, so in each observation we have three time-averaged, three rms
and three phase-lag spectra, respectively, and the parameters in the model of
each component are linked across different instruments. Initially the spin,
$a_{*}$, and the inclination angle, $i$, in the reflection model are free to
change from one observation to the other. However, we notice that these
parameters are high and more or less constant when the QPO frequencies are
low, but drop when the QPO frequencies are high. Miller et al. (2018) and Xu
et al. (2018) suggested a high spin ($>0.84$) and a high inclination angle
($>60^{\circ}$) for MAXI J1535$-$571\. Since in principle these parameters
should not change from one observation to the other, we fix the spin and
inclination angle at the average values of all our fits, i.e. 0.998 and 60.2∘,
respectively. In the fits the hydrogen column density $N_{\text{H}}$ is also
fixed at the average value obtained from the fits of the eight observations.

##  3 Results

###  3.1 The light curve and HID of MAXI J1535$-$571

Figure 1: Left panel: MAXI light curve of MAXI J1535$-$571 in black with the
time of all the Insight-HXMT observations indicated by the vertical dashed
lines. The selected MJD range of MAXI observations is from 57918 to 58312. The
red dashed lines represent the observations that we analyze with type-C QPOs
in the intermediate state. Right panel: Hardness intensity diagram (HID) of
MAXI J1535$-$571 using MAXI data in black. The selected MJD range of MAXI
observations is from 57997 to 58174. The blue and red points indicate the
simultaneous Insight-HXMT observations. The red points represent the
observations that we analyze with type-C QPOs in the intermediate state. Note
that some of the red points overlap each other. There are 3 red points in the
left group and 5 red points in the right group.

In the left panel of Figure 1 we show the MAXI light curve of MAXI J1535$-$571
in the 2–20 keV band during the 2017/2018 outburst, from MJD 57918 to MJD
58312. The X-ray emission of MAXI J1535$-$571 is undetectable until $\sim$ MJD
58000. From that date onwards the source displays a rise up to a maximum
intensity of $\sim$ 7 counts cm-2 s-1 on MJD 58016, and after that it decays
gradually. The vertical dashed lines in Figure 1 show that the Insight-HXMT
observations of MAXI J1535$-$571 only cover part of the full outburst,
especially the rising phase and the end of the outburst. Huang et al. (2018)
reported the type-C QPOs detected by Insight-HXMT in the intermediate state of
the source. Those observations are marked with red dashed lines in Figure 1.
Radio observations of MAXI J1535$-$571 start from MJD 58000 until just after
MJD 58050, covering the whole intermediate states of the source (Russell et
al., 2019; Chauhan et al., 2021).

The MAXI hardness intensity diagram (HID) of MAXI J1535$-$571 in the right
panel of Figure 1 shows that the source traces the characteristic ‘q’ path,
which is typical for other black hole transients (Fender et al., 2004). During
the transition to the high-soft state, the hardness ratio drops as the count
rate increases. Insight-HXMT observations after the intermediate state show
that the source goes into a low-luminosity state. We also plot in red the
eight observations with type-C QPOs that we analyze in detail in this paper.
As shown in Table 1, according to the classification by Huang et al. (2018),
observations 1–5 were in the hard-intermediate state (HIMS), while
observations 6–8 were in the soft-intermediate state (SIMS). We note, however,
that the presence of a type-C QPO in those observations, also classified as
such by Huang et al. (2018), indicates that the source was instead in the
HIMS. Huang et al. (2018) found a type-B QPO, which is part of the definition
of the SIMS, in observation 701, between the first five and the last three
observations in our sample. This indicates that during the first five
observations in our sample the source was in the HIMS, it went briefly into
the SIMS in observation 701 of Huang et al. (2018, this observation is not
included in our sample), and then went back to the HIMS in observations
901–903 in our sample. Note that in the right panel of Figure 1 a blue point
(obs ID P011453500601) between the red points corresponds also to an
observation in the intermediate state with type-C QPOs. However we exclude
this observation since the GTIs of LE, ME, and HE do not overlap and thus we
cannot compute the phase-lag spectrum of this observation.

###  3.2 Joint fits of the rms, phase-lag, and time-averaged spectra

Figure 2: Two representative PDS of MAXI J1535$-$571 using the 1–10 keV LE
data of Insight-HXMT . Left panel: The PDS (top) and the best-fit residuals
(bottom) of Obs 1. Right panel: The PDS (top) and the best-fit residuals
(bottom) of Obs 8.

Figure 3: Two representative rms (top panel), phase-lag (middle panel), and time-averaged (bottom panel) spectra with joint fitting of MAXI J1535$-$571 using Insight-HXMT. The first panel shows the spectra for Obs 1 (QPO frequency 2.6 Hz), while the second panel shows the residuals with respect to the best-fitting model. The third and fourth panels show the same for Obs 8 (QPO frequency 8.1 Hz). The black, red, and green colors indicate the LE, ME, and HE data of Insight-HXMT, respectively. Table 1: Observations and QPOs of MAXI J1535$-$571\. The error bar indicates the 90% confidence level. Each observation ID has a prefix of P011453500. | Obs IDs | MJD | QPO Frequency (Hz) | FWHM
---|---|---|---|---
Obs 1 | 144 | 58008.44 | $2.61\pm 0.02$ | $0.32\pm 0.06$
Obs 2 | 145 | 58008.58 | $2.58\pm 0.02$ | $0.63\pm 0.05$
Obs 3 | 301 | 58011.20 | $2.12\pm 0.02$ | $0.35\pm 0.10$
Obs 4 | 401 | 58012.26 | $2.73\pm 0.02$ | $0.29\pm 0.06$
Obs 5 | 501 | 58013.26 | $3.34\pm 0.02$ | $0.39\pm 0.05$
Obs 6 | 901 | 58017.10 | $8.88\pm 0.08$ | $0.93\pm 0.19$
Obs 7 | 902 | 58017.25 | $9.02\pm 0.06$ | $1.04\pm 0.19$
Obs 8 | 903 | 58017.39 | $8.08\pm 0.11$ | $1.04\pm 0.34$
Table 2: Best-fitting spectral parameters for the model
tbabs*(diskbb+vkompthdk+relxillCp) used to fit the data of MAXI J1535$-$571\.
The error bars indicate the 90% confidence level.

| Obs 1 | Obs 2 | Obs 3 | Obs 4 | Obs 5 | Obs 6 | Obs 7 | Obs 8
---|---|---|---|---|---|---|---|---
$N_{\text{H}}$ (1022) | $5.553^{*}$ | $5.553^{*}$ | $5.553^{*}$ | $5.553^{*}$ | $5.553^{*}$ | $5.553^{*}$ | $5.553^{*}$ | $5.553^{*}$
$kT_{\text{s}}$ (keV) | $0.27\pm 0.03$ | $0.23\pm 0.02$ | $0.24\pm 0.01$ | $0.21\pm 0.05$ | $0.23\pm 0.05$ | $0.78\pm 0.06$ | $0.70\pm 0.05$ | $0.71\pm 0.04$
Norm${}_{\text{disc}}$ ($10^{3}$) | $1700^{+2500}_{-1000}$ | $5900^{+9800}_{-3100}$ | $5500^{+1700}_{-1300}$ | $16000^{+192000}_{-10000}$ | $3600^{+33100}_{-3100}$ | $7.4^{+2.1}_{-1.6}$ | $13\pm 3$ | $8.3^{+1.5}_{-2.6}$
$\Gamma$ | $2.49\pm 0.03$ | $2.49\pm 0.01$ | $2.43\pm 0.02$ | $2.59\pm 0.02$ | $2.60\pm 0.02$ | $2.83^{+0.02}_{-0.05}$ | $2.88\pm 0.02$ | $2.80\pm 0.02$
$kT_{\text{e}}$ (keV) | $32\pm 3$ | $34\pm 2$ | $34\pm 2$ | $44^{+8}_{-4}$ | $39\pm 3$ | $202^{+44}_{-53}$ | $250_{-40}$ | $250_{-48}$
$q_{1}$ | $8.1^{+1.8}_{-1.3}$ | $9.7^{+(P)}_{-1.1}$ | $6.0\pm 0.3$ | $6.1\pm 0.2$ | $5.9\pm 0.2$ | $8.4\pm 0.4$ | $8.9\pm 0.2$ | $8.7\pm 0.3$
$a_{*}$ | $0.998^{*}$ | $0.998^{*}$ | $0.998^{*}$ | $0.998^{*}$ | $0.998^{*}$ | $0.998^{*}$ | $0.998^{*}$ | $0.998^{*}$
$i$ (∘) | $60.2^{*}$ | $60.2^{*}$ | $60.2^{*}$ | $60.2^{*}$ | $60.2^{*}$ | $60.2^{*}$ | $60.2^{*}$ | $60.2^{*}$
$R_{\text{in}}$ (ISCO) | $1.64\pm 0.06$ | $1.66\pm 0.03$ | $1.00^{+0.06}$ | $1.00^{+0.28}$ | $1.00^{+0.06}$ | $1.08\pm 0.02$ | $1.08\pm 0.01$ | $1.07\pm 0.02$
$\log\xi$ | $2.83^{+0.26}_{-0.07}$ | $2.89\pm 0.07$ | $2.86^{+0.07}_{-0.11}$ | $2.7\pm 0.1$ | $3.09^{+0.08}_{-0.05}$ | $4.70_{-0.04}$ | $4.70_{-0.03}$ | $4.6^{+(P)}_{-0.1}$
$A_{\text{Fe}}$ | $0.95^{+0.75}_{-0.12}$ | $0.58\pm 0.06$ | $0.51^{+0.11}_{-0.01}$ | $0.72\pm 0.09$ | $0.88^{+0.08}_{-0.11}$ | $2.1\pm 0.3$ | $3.2^{+1.0}_{-0.3}$ | $2.0\pm 0.2$
Norm${}_{\text{refl}}$ | $0.53\pm 0.09$ | $0.60\pm 0.05$ | $1.1\pm 0.1$ | $1.8\pm 0.3$ | $1.7\pm 0.1$ | $3.3\pm 0.2$ | $3.0\pm 0.2$ | $3.2\pm 0.2$
$L$ ($10^{3}$ km) | $5.1\pm 2.0$ | $4.4\pm 1.2$ | $3.4\pm 0.4$ | $3.0\pm 0.9$ | $3.2^{+0.3}_{-0.6}$ | $5.5^{+1.2}_{-0.9}$ | $4.6\pm 0.7$ | $6.9^{+1.5}_{-1.2}$
$\eta$ | $0.85^{+0.13}_{-0.10}$ | $0.74\pm 0.06$ | $0.77\pm 0.05$ | $0.90\pm 0.05$ | $0.96^{+(P)}_{-0.04}$ | $0.54^{+0.18}_{-0.15}$ | $0.48\pm 0.16$ | $0.78^{+(P)}_{-0.26}$
$\delta\dot{H}_{\text{ext}}$ | $0.085\pm 0.007$ | $0.097\pm 0.004$ | $0.0817\pm 0.003$ | $0.072\pm 0.004$ | $0.077\pm 0.007$ | $0.089\pm 0.015$ | $0.095\pm 0.006$ | $0.13\pm 0.02$
Norm${}_{\text{Comp}}$ | $0.23^{+0.20}_{-0.12}$ | $0.53^{+0.35}_{-0.19}$ | $0.51\pm 0.06$ | $2.5^{+1.9}_{-1.3}$ | $1.5\pm 1.0$ | $0.014^{+0.013}_{-0.006}$ | $0.045\pm 0.021$ | $0.025^{+0.015}_{-0.009}$
$\chi^{2}/\nu$ | $1194.8/1426$ | $1284.8/1426$ | $1415.4/1426$ | $1358.0/1426$ | $1484.7/1426$ | $1404.8/1426$ | $1381.6/1426$ | $1330.6/1426$

  * 1

The letter $(P)$ indicates that the positive 90% error of the parameter pegged
at the upper boundary.

In Table 1 we show the Insight-HXMT observation IDs of MAXI J1535$-$571 with
type-C QPOs; the QPO frequency is in the range of 2.12–9.02 Hz, and the source
is in the intermediate state. Our classification of the QPOs as type-C QPOs is
consistent with the classification in these observations by Huang et al.
(2018) although, as we mentioned above, Huang et al. (2018) misclassified the
state of the source in the last three observations presented here. From MJD
58008 to MJD 58013, when the type-C QPO frequency decreases from 2.61 Hz to
2.12 Hz and then increases to 3.34 Hz, the source is in the HIMS. According to
Huang et al. (2018), on MJD 58015 a type-B QPO is detected, which indicates
that on this date the source is in the SIMS. (Notice that we do not include
this observation and this QPO in our analysis.) Finally, the three
observations on MJD 58017 contain again a type-C QPO with a frequency that
shows a fast increase from 3.34 Hz to around 9 Hz. The presence of the type-C
QPO indicates that, after a short excursion, the source is again in the HIMS.
The variable frequency of the type-C QPOs in the HIMS is consistent with the
measurements of Bhargava et al. (2019) who had a much better time sampling.
The left panel of Figure 2 shows the PDS of observation 1, while the right
panel shows the PDS of observation 8. Both panels display a strong signal of
the type-C QPOs and the harmonics. In the first and third panels of Figure 3
we plot the rms, phase-lag, and time-averaged spectra when the QPO frequencies
are 2.61 Hz and 8.08 Hz, while the second and fourth panels show the residuals
to the best-fitting model. When the QPO frequency is 2.61 Hz, the rms first
increases with energy and gradually flattens above 10 keV, while the phase
lags are always negative (soft) and their magnitude increases with energy.
When the QPO frequency is 8.08 Hz, the rms increases with energy, gradually
flattens above 10 keV, and drops slightly above 30 keV, while the phase lags
are negative and more or less constant with energy, but with larger error bars
than in the case when the QPO frequency is at 2.61 Hz. The fact that the lags
are soft indicates that soft photons emitted by the disc might be largely due
to hard photons from the illuminating corona after being reprocessed (Lee et
al., 2001; Karpouzas et al., 2020).

We fit the rms, phase-lag, time-averaged spectra jointly using the method
described in Section 2.2. As an example, Figure 3 shows the best-fitting model
to the rms, phase-lag, and time-averaged spectra when the QPO frequency is
2.61 Hz (first panel, with residuals in the second panel) and 8.08 Hz (third
panel, with residuals in the fourth panel).

In observations 1–5, when the QPO frequency is in the range of 2.12–3.34 Hz,
the power-law index of the Comptonized component in the time-averaged spectrum
is in the range of 2.4–2.6 and the disc temperature is around 0.3 keV, while
in the observations 6–8, when the QPO frequency is about 9 Hz the power law is
steeper, with $\Gamma\sim$ 2.8, and the temperature of the disc increases to
$\sim$ 1.0 keV (see Figure 4 and Table 2). In observations 1–5 the temperature
of the corona is $kT_{\text{e}}\sim$ 40 keV, whereas in observations 6–8
$kT_{\text{e}}$ changes quickly and reaches the upper limit of 250 keV that we
set in the fits. The change of the values of spectral parameters in the time-
averaged spectrum indicates that MAXI J1535$-$571 may be about to undergo a
state transition from the HIMS to the SIMS. The normalization of the disc
component, Norm${}_{\text{disc}}$, decreases as the temperature of the disc,
$kT_{\text{s}}$, increases; this anti-correlation is consistent with the
findings of Tao et al. (2018). The emissivity index, $q_{1}$, in the relxillCp
component is usually larger than 6, indicating that the corona is close to the
black hole (García et al., 2014, 2018) or the radial dependence of the disc
ionization state (Svoboda et al., 2012). The best-fitting values of the inner
radius of the accretion disc, $R_{\text{in}}$, are consistent with the
innermost stable circular orbit (ISCO), which suggests that the inner disc is
stable from MJD 58008 to MJD 58017. In the HIMS, the logarithm of the disc
ionization parameter is around 3 and the iron abundance is less than that of
the Sun; in the late stage of the HIMS, the disc is more ionized with the
ionization parameter increasing by an order of the magnitude and the best-
fitting iron abundance is around 2 times the solar abundance.

###  3.3 The evolution of the geometry of the corona

Figure 4: The evolution of the parameters of MAXI J1535$-$571 vs. QPO
frequency. The parameters are, respectively, (a) disc temperature,
$kT_{\text{s}}$, (b) photon index, $\Gamma$, (c) corona temperature,
$kT_{\text{e}}$, (d) corona size, $L$, (e) feedback fraction, $\eta$, and (f)
variability of the external heating rate, $\delta\dot{H}_{\text{ext}}$. The
error bars represent the 90% confidence level, while red points show 95% lower
limits. Black dashed broken line is the best fit to the data by a broken line
with a break frequency at 3.25 Hz. The numbers on each panel indicate the
observations in Table 1. The horizontal black dashed lines in panels (c) and
(e) are the upper boundary of 250 keV and 1, respectively. Figure 5: Temporal
evolution of the corona size $L$ (black) and the feedback fraction $\eta$
(red) of MAXI J1535$-$571\. The arrow in the point means the lower limit of
the measurement. The red horizontal dashed line is the upper limit of 1,
according to the definition of the time-dependent Comptonization model.

In Figure 4 we show the evolution of the size of the corona, $L$, the feedback
fraction, $\eta$, and the amplitude of the variability of the external heating
rate, $\delta\dot{H}_{\text{ext}}$, obtained from the time-dependent version
of the model vkompthdk that fits the rms and phase-lag spectra. As the
frequency of the type-C QPO increases from 2.12 Hz to 3.34 Hz, the size of the
corona decreases from 5100 km to 3200 km. At the same time the feedback
fraction remains broadly constant in the range of 0.7–1. When the QPO
frequency increases further, from 3.34 Hz to 9.02 Hz, the feedback fraction
decreases to around 0.5 as the size of the corona increases from 3200 km to
around 5500 km. The amplitude of the variability of the external heating rate
shows the same evolution trend as the size of the corona. The Lense-Thirring
radius (Ingram et al., 2009), calculated from the QPO frequency assuming a
black hole mass of 10 $M_{\odot}$ (Sridhar et al., 2019) and a spin of 0.998,
goes from 200 km to 100 km as the QPO frequency evolves from 2 Hz to 9 Hz.
These radii are well below the size of the corona, which suggests that the
corona always covers a large fraction of the inner disc, consistent with the
relatively large values of $\eta$.

If we order the size of the corona and feedback fraction according to time
(see Table 2 and Figure 5), in observations 1–5 the size of the corona starts
at 5100 km, then decreases to 3200 km while, at the same time, the feedback
fraction first decreases slightly from 0.85 to 0.74, and then increases
marginally and stays at values larger than 0.9. In observations 5–8 the corona
size increases from 3200 km to around 5500 km, while the feedback fraction
decreases from 0.96 to around 0.5.

Using NICER data, Rawat et al. (2022) found that the relation of the time lags
of the QPO between two broad energy bands, the spectral parameters of the
source, $kT_{\text{s}}$, and $\Gamma$, the corona size $L$ and the feedback
fraction $\eta$ vs. QPO frequency show a statistically significant break at a
QPO frequency $\nu_{\text{c}}=3.5$ Hz. (Because of NICER’s lack of high-energy
coverage, different from us they could not measure $kT_{\text{e}}$ in their
data.) To explore this in our observations, we first fit each relation in
Figure 4 independently with either a straight or a broken line. Since in the
latter fits the values of the break frequency, $\nu_{\text{c}}$, of the
individual fits are consistent with each other within the 3-$\sigma$ errors,
we fit all the relation again simultaneously, this time linking the break
frequency of all fits. Differently from Rawat et al. (2022), we find that
since the F-test probability is 0.03, the fit with a broken line is only
marginally better than that with a straight line, possibly due to the fact
that we have no observation with QPO frequency in the range of 4–8 Hz. Our
best-fitting broken line gives $\nu_{\text{c}}=3.25^{+0.60}_{-0.21}$ Hz,
consistent with the value found by Rawat el al. (2022). We also fit a broken
power-law instead of broken line to the size of the corona and the external
heating rate vs. QPO frequency, however, the reduced $\chi^{2}$ and the break
frequency do not change significantly compared to the fits with a broken line.

##  4 Discussion

Figure 6: Schematic figure illustrating the time evolution of the corona and
the jet of MAXI J1535$-$571 (see text for a more detailed explanation). From
MJD 58008 to MJD 58011 the size of the corona contracts; from MJD 58011 to MJD
58013 the corona contracts vertically but expands horizontally. In the period
before MJD 58013 the radio spectrum of the source is consistent with optically
thick synchrotron emission from a compact steady jet and on MJD 58013 the jet
is quenched. From MJD 58013 to MJD 58017 the corona expands vertically but
contracts horizontally and after MJD 58017 the radio data are consistent with
optically thin synchrotron emission from a transient jet.

We have analyzed all the observations of the 2017/2018 outburst of MAXI
J1535$-$571 with Insight-HXMT up to 100 keV. The light curve and ‘q’ path in
the HID (Figure 1) show that MAXI J1535$-$571 behaves like a typical black
hole X-ray transient. We find that the temperature of the disc and the photon
index of the hard Comptonized component increase with time, indicating a
softening of the source spectrum as the source moves along the HIMS. At the
same time the inner radius of the accretion disc remains rather stable. We
jointly fit the time-averaged, rms and phase-lag spectra of the QPO of eight
Insight-HXMT observations of MAXI J1535$-$571 with type-C QPOs; the frequency
of the QPO ranges from 2.12 Hz to 9.02 Hz during the HIMS. As shown in Figure
4, as the frequency of the type-C QPO increases, the size of the corona, $L$,
and the amplitude of the variability of the external heating rate,
$\delta\dot{H}_{\text{ext}}$, decrease and then increase slightly. At the same
time the feedback fraction, $\eta$, does the opposite.

Previous studies have reported the evolution of the physical parameters of
MAXI J1535$-$571, e.g., spin of the black hole, inclination angle of the
accretion disc, and height of the corona. For instance, Huang et al. (2018)
analyzed the Insight-HXMT observations, some of which we present in this
paper, and concluded that the system inclination angle must be high, judging
from the soft lags of the type-C QPOs. From fits to the time-averaged spectra
with the reflection model rexilllpCp, Xu et al. (2018) reported a high spin
and a corona that is very close to the black hole in the low hard state.
Miller et al. (2018) used relxill to deduce a high spin and a high emissivity
index, which also indicates that the corona is close to the black hole. All of
these authors concluded that the inclination angle of MAXI J1535$-$571 is
high. However, the spectral analysis mentioned before just focus on one
observation of NuSTAR or NICER. In the hard/soft intermediate state, we use
the relxillCp model to fit eight observations, showing an average spin of
0.998 and inclination angle of $60.2^{\circ}$. Our fits show that the
emissivity index is large ($>6$), indicating a corona close to black hole
(Wilkins & Fabian, 2012) or a radial dependence of the disc ionization during
these observations (Svoboda et al., 2012). The time-dependent Comptonization
model vkompthdk requires a large corona to explain the lags of the photons. If
the corona is large, the radial change of the disc ionization profile could be
responsible for the high value of the emissivity index in MAXI J1535$-$571.

In an outburst of a BHB, as the disc emission increases during the softening
of the spectrum in the LHS and HIMS (e.g., Méndez & van der Klis, 1997; Huang
et al., 2018; Zhang et al., 2020b), the rate at which the corona cools down
due to inverse-Compton scattering would increase. The temperature of the
corona, and hence the energy at which the spectrum cuts off at high energies,
may either decrease or increase, depending on the balance between the inverse-
Compton cooling and the mechanism that heats up the corona (Merloni & Fabian,
2001; Karpouzas et al., 2020). Our results show that during the HIMS, the
high-energy cutoff increases from 32 keV to at least 250 keV. This trend is
the same as the ones shown in Joinet et al. (2008), Del Santo et al. (2008)
and Motta et al. (2009). Del Santo et al. (2008) showed this phenomenology in
the HIMS data of GRO J1655$-$40 with INTEGRAL, which displays a non-thermal
hard tail at high energy. This tail becomes more dominant and the high-energy
cutoff increases when the source is about to transit into the soft state.
Similarly, Motta et al. (2009) showed that in the case of GX 339$-$4 the
power-law cutoff energy first decreases as the source moves in the LHS and
HIMS, and increases again just before the transition to the SIMS. The
evolution of the power-law cutoff energy in these two sources is consistent
with the behavior of the corona temperature that we report here for MAXI
J1535$-$571 in the HIMS.

In the framework of Lense-Thirring (LT) precession as the explanation of the
dynamical origin of the type-C QPOs, Ingram et al. (2009) assumed that a
torus-like corona precesses within the inner radius of a truncated disc (see
also, You et al., 2018). Assuming that the mass of the black hole in MAXI
J1535$-$571 is $10M_{\odot}$ (Sridhar et al., 2019) and the spin is 0.998 from
our fits, the LT radius decreases from 200 km to 100 km as the QPO frequency
increases from 2 Hz to 9 Hz. The LT radius is much smaller than the corona
size we obtain from our fits, the minimum of which is around 3000 km. The
time-dependent Comptonization model (Karpouzas et al., 2020) that we use does
not explain the dynamical origin of the QPO, but only requires that there is a
coupled oscillation of the physical quantities: corona temperature,
$kT_{\text{e}}$, temperature of the inner accretion disc, $kT_{\text{s}}$,
external heating rate, $\delta\dot{H}_{\text{ext}}$, and the time-dependent
photon number density that is dynamically scattering in the plasma. We favor
the scenario in which the inner disc oscillates since, if the source of soft
photons precesses, the line-of-sight soft photons oscillate and would
potentially cause a perturbation in the corona temperature.

Assuming that the time lags of the broadband noise component are due to
reverberation, Kara et al. (2019) concluded that the corona of MAXI J1820+070
contracts as the source is in the low-hard state. Wang et al. (2021) proposed
that the source expands again when evolving further in the intermediate state.
Using the time-dependent Comptonization model, Karpouzas et al. (2021) showed
that in GRS 1915+105 the size of the corona first decreases from 10000 km to
100 km and then increases to 1000 km as the QPO frequency increases from 0.5
Hz to 6 Hz, while the feedback fraction increases monotonically (see also
Méndez et al., 2022). Our result of the size of the corona shows the same
trend in MAXI J1535$-$571 as that in GRS 1915+105, but the size of the corona
in the case of MAXI J1535$-$571 is about one order of magnitude larger than in
GRS 1915+105. In MAXI J1535$-$571, different from GRS 1915+105 the feedback
fraction first remains broadly constant in the range of 0.7–1 and then
decreases as the QPO frequency increases, providing a different picture of the
evolution of the corona.

Similar to the case of GRS 1915+105 (Méndez et al., 2022), the large size of
the corona in MAXI J1535$-$571 indicates that the corona would cover the disc
beyond the inner radius. However, at the same time that the corona is large,
the feedback fraction is small. This suggests that the corona is not
spherical, but extends perpendicularly to the accretion disc. In these
conditions the assumption of a spherical corona in the model of Karpouzas et
al. (2020) would not be satisfied. In this case $L$ must therefore represent a
characteristic size of the corona such that the magnitude of the predicted
lags can match the observed ones (see Méndez et al., 2022).

Our results of the temporal evolution of the corona may be compared to the
evolution of the radio jet detected in this source. Radio measurements
(Russell et al., 2019, 2020; Chauhan et al., 2021) show that during the
intermediate state of MAXI J1535$-$571 a steady, optically thick jet gradually
quenches and then changes into a transient jet. More specifically, on MJD
58008, the radio observations detect a steady, optically thick jet. On MJD
58013 the radio emission weakens and the jet is quenched. On MJD 58017, the
radio emission reappears, but this time the spectrum is optically thin,
indicating a transient jet with discrete ejecta.

Based on our results and the radio measurements of Russell et al. (2019) and
Chauhan et al. (2021), in Figure 6, we present schematically a possible
scenario of the evolution of the X-ray corona and the simultaneous radio jet.
The parameters from the reflection suggest a rather stable inner radius of the
disc that is very close to the ISCO, so the disc inner radius does not change
significantly in Figure 6. Initially, the QPO frequency is 2.61 Hz, and the
size of the corona is 5100 km. When the QPO frequency first decreases from
2.61 Hz to 2.12 Hz, the corona size decreases from 5100 km to 3400 km and the
feedback fraction decreases slightly from 0.85 to 0.77 in the time period from
MJD 58008 to MJD 58011 (Figure 5). As the QPO frequency increases from 2.12 Hz
to 3.34 Hz, the corona size remains constant at around 3000–3400 km. Since the
feedback fraction increases from 0.77 to 0.96, suggesting that the corona
covers a larger part of the disc, the corona should expand along the
horizontal direction but contract along the vertical direction. This phase
corresponds to the time period from MJD 58011 to MJD 58013. In the period MJD
58013–58017, the QPO changes from 3.34 Hz to 8.08 Hz and the corona expands
again. Since in this phase the feedback fraction decreases from 0.96 to $\sim$
0.5, indicating that a smaller part of the disc is covered by the corona, the
expansion of the corona must be along the vertical direction.

Using different approaches, a similar evolution of the corona has been
reported. Using the reflection model relxillCp (García et al., 2014) to fit
the Insight-HXMT data, You et al. (2021) found that the corona in MAXI
J1820+070 outflows faster when it moves closer to the black hole, which
suggests a jet-like corona and that the jet gains energy as the corona
outflows. Utilizing the reltrans model (Ingram et al., 2019) to fit the
broadband time-lag spectra, Wang et al. (2021) found a connection between the
corona and the jet in the intermediate state of MAXI J1820+070. The picture
proposed by Wang et al. (2021) is in agreement with ours: The quenching of the
steady jet follows the contraction of the corona, while the transient jet
appears after the expansion of the corona. More observations and spectral
modelling using the time-dependent Comptonization model in the future would
provide a complete picture of the corona-jet evolution, and the role of the
disc during this evolution.

## Acknowledgements

We would like to thank the referee for constructive comments that helped
improve this paper. This work is supported by the National Key R&D Program of
China (2021YFA0718500). We thank Professor Bei You for helpful discussions.
Y.Z. acknowledges support from China Scholarship Council (CSC 201906100030).
M.M. and F.G. acknowledge the research programme Athena with project number
184.034.002, which is (partly) financed by the Dutch Research Council (NWO).
F.G. is a CONICET researcher. F.G. acknowledges support by PIP 0113 (CONICET).
This work received financial support from PICT-2017-2865 (ANPCyT). T.M.B.
acknowledges financial contribution from the agreement ASI-INAF n. 2017-14-H.0
and from PRIN INAF 2019 n.15. S.Z. acknowledges support of the National
Natural Science Foundation of China under grant U1838201. L.T. acknowledges
funding support from the National Natural Science Foundation of China (NSFC)
under grant No. 12122306. D.A. acknowledges support from the Royal Society.
L.Z. acknowledges support from the Royal Society Newton Funds. C.B. is a CIN
fellow.

## Data Availability

The data used in this article are accessible at the Insight-HXMT website
http://hxmtweb.ihep.ac.cn.

## References

  * Arnaud (1996) Arnaud K. A., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 17
  * Bambi et al. (2021) Bambi C., et al., 2021, Space Sci. Rev., 217, 65
  * Bellavita & et al., (2022) Bellavita C., et al., 2022, in prep
  * Belloni & Hasinger (1990) Belloni T., Hasinger G., 1990, A&A, 230, 103
  * Belloni et al. (2002) Belloni T., Psaltis D., van der Klis M., 2002, ApJ, 572, 392
  * Belloni et al. (2005) Belloni T., Homan J., Casella P., van der Klis M., Nespoli E., Lewin W. H. G., Miller J. M., Méndez M., 2005, A&A, 440, 207
  * Bhargava et al. (2019) Bhargava Y., Belloni T., Bhattacharya D., Misra R., 2019, MNRAS, 488, 720
  * Bu et al. (2021) Bu Q. C., et al., 2021, ApJ, 919, 92
  * Casella et al. (2005) Casella P., Belloni T., Stella L., 2005, ApJ, 629, 403
  * Chauhan et al. (2021) Chauhan J., et al., 2021, Publ. Astron. Soc. Australia, 38, e045
  * Corbel et al. (2004) Corbel S., Fender R. P., Tomsick J. A., Tzioumis A. K., Tingay S., 2004, ApJ, 617, 1272
  * Del Santo et al. (2008) Del Santo M., Malzac J., Jourdain E., Belloni T., Ubertini P., 2008, MNRAS, 390, 227
  * Done et al. (2007) Done C., Gierliński M., Kubota A., 2007, A&ARv, 15, 1
  * Fabian et al. (1989) Fabian A. C., Rees M. J., Stella L., White N. E., 1989, MNRAS, 238, 729
  * Fender (2001) Fender R. P., 2001, MNRAS, 322, 31
  * Fender (2006) Fender R., 2006, Jets from X-ray binaries. pp 381–419
  * Fender et al. (2004) Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105
  * Galeev et al. (1979) Galeev A. A., Rosner R., Vaiana G. S., 1979, ApJ, 229, 318
  * García et al. (2014) García J., et al., 2014, ApJ, 782, 76
  * García et al. (2018) García J. A., et al., 2018, ApJ, 864, 25
  * García et al. (2021) García F., Méndez M., Karpouzas K., Belloni T., Zhang L., Altamirano D., 2021, MNRAS, 501, 3173
  * Gilfanov (2010) Gilfanov M., 2010, X-Ray Emission from Black-Hole Binaries. p. 17, doi:10.1007/978-3-540-76937-8_2
  * Haardt & Maraschi (1991) Haardt F., Maraschi L., 1991, ApJ, 380, L51
  * Homan et al. (2001) Homan J., Wijnands R., van der Klis M., Belloni T., van Paradijs J., Klein-Wolt M., Fender R., Méndez M., 2001, ApJS, 132, 377
  * Huang et al. (2018) Huang Y., et al., 2018, ApJ, 866, 122
  * Ingram & Motta (2019) Ingram A. R., Motta S. E., 2019, New Astron. Rev., 85, 101524
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Ingram et al. (2019) Ingram A., Mastroserio G., Dauser T., Hovenkamp P., van der Klis M., García J. A., 2019, MNRAS, 488, 324
  * Jiang et al. (2020) Jiang J., Fürst F., Walton D. J., Parker M. L., Fabian A. C., 2020, MNRAS, 492, 1947
  * Joinet et al. (2008) Joinet A., Kalemci E., Senziani F., 2008, ApJ, 679, 655
  * Kara et al. (2019) Kara E., et al., 2019, Nature, 565, 198
  * Karpouzas et al. (2020) Karpouzas K., Méndez M., Ribeiro E. M., Altamirano D., Blaes O., García F., 2020, MNRAS, 492, 1399
  * Karpouzas et al. (2021) Karpouzas K., Méndez M., García F., Zhang L., Altamirano D., Belloni T., Zhang Y., 2021, MNRAS, 503, 5522
  * Kennea et al. (2017) Kennea J. A., Evans P. A., Beardmore A. P., Krimm H. A., Romano P., Yamaoka K., Serino M., Negoro H., 2017, The Astronomer’s Telegram, 10700, 1
  * Kong et al. (2020) Kong L. D., et al., 2020, Journal of High Energy Astrophysics, 25, 29
  * Kumar & Misra (2014) Kumar N., Misra R., 2014, MNRAS, 445, 2818
  * Lee & Miller (1998) Lee H. C., Miller G. S., 1998, MNRAS, 299, 479
  * Lee et al. (2001) Lee H. C., Misra R., Taam R. E., 2001, ApJ, 549, L229
  * Ma et al. (2021) Ma X., et al., 2021, Nature Astronomy, 5, 94
  * Maccarone & Coppi (2003) Maccarone T. J., Coppi P. S., 2003, MNRAS, 338, 189
  * Méndez & van der Klis (1997) Méndez M., van der Klis M., 1997, ApJ, 479, 926
  * Méndez et al. (2013) Méndez M., Altamirano D., Belloni T., Sanna A., 2013, MNRAS, 435, 2132
  * Méndez et al. (2022) Méndez M., Karpouzas K., García F., Zhang L., Zhang Y., Belloni T., Altamirano D., 2022, Nature Astronomy, in press
  * Merloni & Fabian (2001) Merloni A., Fabian A. C., 2001, MNRAS, 321, 549
  * Miller-Jones et al. (2012) Miller-Jones J. C. A., et al., 2012, MNRAS, 421, 468
  * Miller et al. (2018) Miller J. M., et al., 2018, ApJ, 860, L28
  * Mirabel & Rodríguez (1994) Mirabel I. F., Rodríguez L. F., 1994, Nature, 371, 46
  * Mitsuda et al. (1984) Mitsuda K., et al., 1984, PASJ, 36, 741
  * Miyamoto et al. (1988) Miyamoto S., Kitamoto S., Mitsuda K., Dotani T., 1988, Nature, 336, 450
  * Morgan et al. (1997) Morgan E. H., Remillard R. A., Greiner J., 1997, ApJ, 482, 993
  * Motta et al. (2009) Motta S., Belloni T., Homan J., 2009, MNRAS, 400, 1603
  * Muñoz-Darias et al. (2011) Muñoz-Darias T., Motta S., Belloni T. M., 2011, MNRAS, 410, 679
  * Negoro et al. (2017) Negoro H., et al., 2017, The Astronomer’s Telegram, 10708, 1
  * Nowak (2000) Nowak M. A., 2000, MNRAS, 318, 361
  * Nowak et al. (1999) Nowak M. A., Vaughan B. A., Wilms J., Dove J. B., Begelman M. C., 1999, ApJ, 510, 874
  * Poutanen et al. (2018) Poutanen J., Veledina A., Zdziarski A. A., 2018, A&A, 614, A79
  * Rawat & et al., (2022) Rawat D., et al., 2022, submitted to MNRAS
  * Reig et al. (2000) Reig P., Belloni T., van der Klis M., Méndez M., Kylafis N. D., Ford E. C., 2000, ApJ, 541, 883
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Remillard et al. (1999) Remillard R. A., Morgan E. H., McClintock J. E., Bailyn C. D., Orosz J. A., 1999, ApJ, 522, 397
  * Russell et al. (2011) Russell D. M., Miller-Jones J. C. A., Maccarone T. J., Yang Y. J., Fender R. P., Lewis F., 2011, ApJ, 739, L19
  * Russell et al. (2017) Russell T. D., Miller-Jones J. C. A., Sivakoff G. R., Tetarenko A. J., Jacpot Xrb Collaboration 2017, The Astronomer’s Telegram, 10711, 1
  * Russell et al. (2019) Russell T. D., et al., 2019, ApJ, 883, 198
  * Russell et al. (2020) Russell T. D., et al., 2020, MNRAS, 498, 5772
  * Sridhar et al. (2019) Sridhar N., Bhattacharyya S., Chandra S., Antia H. M., 2019, MNRAS, 487, 4221
  * Stella & Vietri (1998) Stella L., Vietri M., 1998, ApJ, 492, L59
  * Stiele & Kong (2018) Stiele H., Kong A. K. H., 2018, ApJ, 868, 71
  * Stirling et al. (2001) Stirling A. M., Spencer R. E., de la Force C. J., Garrett M. A., Fender R. P., Ogley R. N., 2001, MNRAS, 327, 1273
  * Svoboda et al. (2012) Svoboda J., Dovčiak M., Goosmann R. W., Jethwa P., Karas V., Miniutti G., Guainazzi M., 2012, A&A, 545, A106
  * Tagger & Pellat (1999) Tagger M., Pellat R., 1999, A&A, 349, 1003
  * Tao et al. (2018) Tao L., et al., 2018, MNRAS, 480, 4443
  * Tingay et al. (1995) Tingay S. J., et al., 1995, Nature, 374, 141
  * Van der Klis (1989) Van der Klis M., 1989, in Ögelman H., van den Heuvel E. P. J., eds, NATO Advanced Study Institute (ASI) Series C Vol. 262, Timing Neutron Stars. p. 27
  * Vaughan & Nowak (1997) Vaughan B. A., Nowak M. A., 1997, ApJ, 474, L43
  * Verner et al. (1996) Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 1996, ApJ, 465, 487
  * Wang et al. (2021) Wang J., et al., 2021, ApJ, 910, L3
  * Wilkins & Fabian (2012) Wilkins D. R., Fabian A. C., 2012, MNRAS, 424, 1284
  * Wilms et al. (2000) Wilms J., Allen A., McCray R., 2000, ApJ, 542, 914
  * Xu et al. (2018) Xu Y., et al., 2018, ApJ, 852, L34
  * You et al. (2018) You B., Bursa M., Życki P. T., 2018, ApJ, 858, 82
  * You et al. (2021) You B., et al., 2021, Nature Communications, 12, 1025
  * Zdziarski et al. (1996) Zdziarski A. A., Johnson W. N., Magdziarz P., 1996, MNRAS, 283, 193
  * Zhang et al. (2014) Zhang S., Lu F. J., Zhang S. N., Li T. P., 2014, in Takahashi T., den Herder J.-W. A., Bautz M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9144, Space Telescopes and Instrumentation 2014: Ultraviolet to Gamma Ray. p. 914421, doi:10.1117/12.2054144
  * Zhang et al. (2020a) Zhang L., et al., 2020a, MNRAS, 494, 1375
  * Zhang et al. (2020b) Zhang L., et al., 2020b, MNRAS, 499, 851
  * Życki et al. (1999) Życki P. T., Done C., Smith D. A., 1999, MNRAS, 309, 561


This work has been submitted to the IEEE for possible publication. Copyright
may be transferred without notice, after which this version may no longer be
accessible.

# Detecting Drill Failure in the Small Short-sound Drill Dataset

Thanh Tran  Nhat Truong Pham  and Jan Lundgren This research was supported by
the EU Regional Fund, the MiLo Project (No. 20201888), and the Acoustic sensor
set for the AI monitoring systems (AISound ) project. (Thanh Tran and Nhat
Truong Pham are co-first authors.)Thanh Tran and Jan Lundgren are with the
Department of Electronics Design, Mid Sweden University, Sundsvall, Sweden
(e-mails<EMAIL_ADDRESS>jan.lundgren@miun.se).Nhat Truong Pham is with
the Institute for Computational Science, Ton Duc Thang University, Ho Chi Minh
City, 72915, Vietnam (e-mail: phamnhattruong.st@tdtu.edu.vn).

###### Abstract

Monitoring the conditions of machines is vital in the manufacturing industry.
Early detection of faulty components in machines for stopping and repairing
the failed components can minimize the downtime of the machine. This article
presents an approach to detect the failure occurring in drill machines based
on drill sounds from Valmet AB. The drill dataset includes three classes:
anomalous sounds, normal sounds, and irrelevant sounds, which are also labeled
as “Broken”, “Normal”, and “Other”, respectively. Detecting drill failure
effectively remains a challenge due to the following reasons. The waveform of
drill sound is complex and short for detection. Additionally, in realistic
soundscapes, there are sounds and noise in the context at the same time.
Moreover, the balanced dataset is small to apply state-of-the-art deep
learning techniques. To overcome these aforementioned difficulties, we
augmented sounds to increase the number of sounds in the dataset. We then
proposed a convolutional neural network (CNN) combined with a long short-term
memory (LSTM) to extract features from log-Mel spectrograms and learn global
high-level feature representation for the classification of three classes. A
leaky rectified linear unit (Leaky ReLU) was utilized as the activation
function for our proposed CNN instead of the rectified linear unit (ReLU).
Moreover, we deployed an attention mechanism at the frame level after the LSTM
layer to learn long-term global feature representations. As a result, the
proposed method reached an overall accuracy of 92.35% for the drill failure
detection system.

Attention mechanism, Convolutional neural network, Drill fault detection,
Leaky ReLU, Log-Mel spectrogram, Long short-term memory

##  1 Introduction

Drill fault detection systems are widely used in factories to prevent machine
failure. The drilling machine is included 90 or 120 drill bits to drill
thousand of small holes on the surface of the metal [1]. Broken drill bits can
seriously damage a product because they lack holes on the metal surface. The
maintenance technician stops the machine every 10 min to identify any broken
drill bits and change them before re-initiating the drilling machine.
Therefore, a fault detection system for the drilling machine is very crucial
to minimize the downtime of the machine as well as the maintenance cost.

Many studies have been conducted on detecting and diagnosing drill failure in
the past decade. Choi  _et al._ [2] proposed a time domain and frequency
domain feature extraction method named Characteristic Parameters of the Drill
Failure (CPDF). In the second step, a multilayer perceptron (MLP) was used to
predict drill failure based on the drill state index threshold. This would be
lower the error rate. In order to solve this issue as well as improve the
accuracy of diagnosis for drill failure, Skalle  _et al._ [3] proposed a
method based on symptom detection (e.g., Soft Fm, Cutting Accumulation, Local
Dog Leg). In [4], Kumar  _et al._ used vibration signals to detect and
classify drill failures using three different classifier schemes: Artificial
Neural Network (ANN), support vector machines (SVM), and Bayesian classifiers.
Because vibration-based signals often contain noise, several techniques were
required to remove noise and separate the sources to improve fault detection
accuracy.

Researchers have used sound and vibration analysis in recent years to detect
and classify faults [5]. As a result of the development of deep learning, CNNs
(Convolutional Neural Networks) were used to automatically extract features to
diagnose and classify faults on machines, especially drill machines. Due to
the advantage of acoustic analysis over vibrations, Glowacz [6] proposed an
acoustic-based fault detection method for electric impact drills and coffee
grinders. These acoustic features, including the root mean square (RMS) and a
method selection of amplitude using a multi-expanded filter
(MSAF-17-MULTIEXPANDED-FILTER-14), were used to classify fault status by the
nearest neighbor classifier. Additionally, to detect the fault in electric
impact drills, it is necessary to determine the fault of the gearbox device of
the drill, since the gears are the main component of the power transmission.
Jing  _et al._ [7] proposed a method for detecting electric impact drill
failure by using logistic regression from time-varying loudness and acoustic
signals.

Recently, a number of techniques have been investigated in the field of fault
detection and machine condition monitoring. Hou  _et al._ [8] used wavelet
packet energy to extract features from acoustic signals, then applied a
feature selection method based on the Pearson correlation coefficient to
select features. The selected features were used to classify the fault status
with a neural network classifier. In addition to synchronous hydraulic motors,
this approach can be applied to other rotating machines as well. In another
approach, Wang  _et al._ [9] proposed a multimodal method to detect bearing
faults by fusing acoustic and vibration signals collected from the
accelerometer and the microphone using the 1 dimensional CNN.

In recent years, deep learning has had a great deal of success in the
detection and diagnosis of mechanical faults by using vibration and acoustics
signals [10], [11], [12], [13], [14], [15], [16]. Besides, recent studies have
demonstrated that image representations of sound signals can be used to train
the deep learning architecture for sound classification tasks. Researchers
have proposed a lot of image representations for sounds, such as MFCC (Mel
Frequency Cepstral Coefficients) [17], [18], spectrogram [19], Mel spectrogram
[20]. Additionally, many state-of-the-art deep learning models have been used
for sound classification. Boddapati  _et al._ [19] compared the classification
accuracy of AlexNet and GoogleLNet on three different feature representations
of sound (spectrogram, MFCC, and cross recurrence plot). A variant of
conditional neural networks, called masked conditionaL neural network (MCLNN)
has been proposed by Medhat  _et al._ [21] for classifying sounds. Researchers
have used dilated CNNs with dilated filters and leaky ReLU activation
functions [17], [22]. The effect of modulating the dilation rate in dilated
CNN on sound classification was compared in Chen  _et al._ [22]. Recent
studies have shown that recurrent neural networks (RNN) produce excellent
results for variable-length sound sequences. Zhang  _et al._ [23] proposed a
CNN architecture to learn spectro-temporal features and a bidirectional gated
recurrent unit (Bi-GRU) with a frame-level attention mechanism for sound
classification. Wang  _et al._ [24] proposed a CNN architecture with a
parallel temporal-spectral attention mechanism to capture certain frames where
sound events occur and pay attention to varying frequency bands.

Our article proposed an approach to detect drill machine failures based on
drill sounds from Valmet AB. Most drilling fault detection studies, however,
used a large, balanced dataset. In fact, the anomalous sounds that were
captured while drill machines were damaged are only a small fraction of the
overall dataset. Deep learning models might not do well in the prediction of
minority classes due to imbalanced datasets in real-world applications. Also,
the extracted features from raw sound signals are insufficient for
classification since the sample durations for sounds in our dataset are around
20.83 ms and 41.67 ms. This makes it more challenging to compare our results
to those of previous research in the field of drilling sounds classification.
As a result, an end-to-end deep learning system faces many challenges when it
comes to detecting drill faults. To overcome these difficulties, we re-sampled
the dataset (undersampling) as described in [1]. Then we applied data
augmentation to generate more samples of the dataset. The augmentation methods
were shifting the sound by 5 ms and increasing the volume by 2. These sounds
in the augmented dataset were converted into log-Mel spectrograms. In
addition, we proposed a CNN combined with an attention-based LSTM for
classifying drill sounds. Feature maps were extracted from the log-Mel
spectrograms using CNN, and then an LSTM layer was used to learn high-level
global feature representation from extracted features. We used Leaky ReLU in
CNN instead of using ReLU to alleviate the potential problem that CNN stops
learning when the ReLU has a value of less than zero. Leaky ReLU helps CNN
continue learning when input values are negative. To focus on the important
parts of drill sounds and discard the unnecessary parts, an attention layer
was added after the LSTM.

The rest of the article is organized as follows. Our proposed method is
described in Section II. Section III presents experiment results and
discussion. Finally, section IV is the conclusion.

##  2 Proposed Methodology

Our proposed architecture is described as shown in Fig. 1. Firstly, original
sounds were applied audio augmentation methods to increase the number of
samples in the dataset. In the next step, we proposed a small CNN architecture
that includes five layers to generate features from the Mel spectrograms of
sounds. Finally, these features were used as the input of the LSTM with the
attention mechanism to learn high-level feature representation. The details of
the layers in our proposed model are described in Table 1, where $nC$ is the
number of classes.

Figure 1: Our proposed methodology. Table 1: The layers of the proposed model. Layer | Kernel/Size | Output shape
---|---|---
Input | — | 100$\times$96$\times$1
Conv 1 | 3$\times$3 | 100$\times$96$\times$128
Max_pooling 1 | 2$\times$4 | 50$\times$24$\times$128
Conv 2 | 3$\times$3 | 50$\times$24$\times$128
Max_pooling 2 | 2$\times$4 | 25$\times$6$\times$128
Conv 3 | 3$\times$3 | 25$\times$6$\times$256
LSTM | 256 | —
Attention | — | —
FC 1 | — | 64
FC 2 | — | nC

###  2.1 Data Augmentation

This study is a continuation of the study presented in [1]. Therefore, we used
the same dataset. Valmet AB drills small holes in metal plates with multiple
machines. There are two type of drilling machines in a factory that are 90 and
120 bits. Fig. 2 shows a healthy drill bit and a broken drill bit [1].

Figure 2: A healthy drill bit (on the left side) and a broken drill bit (on
the right side) [1].

In this dataset, sound from a drill machine in Sundsvall, Sweden was recorded
with four AudioBox iTwo Studio microphones. For capturing drill sounds, 96 kHz
was used as the sampling rate. The dataset has two parts. There are 833 files
in one part. In this part, the length of each sound is 20.83 ms. The other
part contains 40,417 files. In this part, the length of each sound is 41.67
ms. On average, broken drill bits account for just 0.16% of the total number
of drill sounds recorded (67 anomalous sounds out of 41,250 total drill
sounds). Therefore, broken drills are hard to detect. To proceed further, it
is therefore important to create a balanced dataset. In the majority classes
(normal sounds and unrelated sound classes), we selected randomly the same
amount of sounds as in the minority classes (anomalous sounds). We have 201
sounds in the balanced dataset after re-sampling. There is a limited amount of
data to work within the balanced dataset.

Although we tuned hyper-parameters of our model to adapt to the drill fault
detection task, however, the lack of sound when the drills were broken is
still a big challenge. To tackle this challenge, we applied data augmentation
to the original sounds. Thus, the number of sounds in the dataset increased.
Besides, data augmentation helps improve the generalization capability of our
proposed model. There are many methods of audio augmentation such as time-
stretching, pitch-shifting, volume control, noise addition, etc. It is not
appropriate to apply some augmentation methods to the sounds in our dataset
since they are very short at only 20.83 ms or 41.67 ms. Experimentation
revealed that only time-shifting and volume control data augmentation methods
are effective for our data set.

In this article, we applied time-shifting and volume control to generate
syntactical sounds. We did not add noise to the sound as an augmentation
method because the sound in our dataset is very short. Noise makes it
difficult to classify sounds. Matlab provides a simple function,
audioDataAugmenter to augment the sound.

####  2.1.1 Time-shifting

A time-shifting is the process of shifting the sound backward or forward at
random. We shifted the starting point of the sound by 5 ms to the right, then
padded it back to its original length. Fig. 3a shows the time representation
of the original fault sound and augmented sound using time-shifting.

(a) Time representation of the original fault sound and augmented sound using
time-shifting.

(b) Time representation of the original fault sound and the augmented sound
using volume control.

Figure 3: The time representation of the original fault sound and the
augmented sound using time-shifting and volume control.

(a) An anomalous sound (a broken drill).

(b) A time-shifted sound.

(c) A volume gained sound.

Figure 4: Log-Mel spectrograms of an original anomalous sound, the augmented
sound using time-shifting, the augmented sound using volume control.

####  2.1.2 Volume control

We increased the volume by multiplying the audio by a random amplitude factor.
The volume gain was set as 2 dB. Using this technique, we can gain some in-
variance concerning the audio input gain. In Fig. 3b, we show the time
representation of the original fault sound and augmented sound using the
volume control.

###  2.2 Convert Sounds into Log-Mel Spectrograms

Recent advances in the field of image classification using CNN for multiple
classes with high accuracy motivated us to investigate the ability of image
representation of sounds to detect drill failures. In this paper, drill sounds
were converted into log-Mel spectrograms to feed into our proposed CNN. The
log-Mel spectrogram was generated as follows. Given a raw drill sound, we
computed the Mel spectrogram using short-time Fourier transform (STFT) with
Hamming windows of 100 ms and the hop length of 50 ms, the length of FFT was
2,048, the sampling rate was 96 kHz, and the number of Mel-filter bank was 96.
Since the authors in [25] found that the log-scaled Mel spectrogram improves
the classification accuracy compared to the Mel spectrogram. Therefore, we
took the logarithm of the Mel spectrogram as the input of our proposed CNN
architecture. Fig. 4 shows log-Mel spectrograms of an original anomalous sound
and its augmented sound using volume control and time-shifting.

###  2.3 Extract Features using CNN with Leaky ReLU

We proposed a CNN architecture only used for extracting features from log-Mel
spectrograms. As a result, we used the third convolutional layer for
extracting features instead of adding a dense layer at the end. Additionally,
instead of using ReLU as the activation function, we used the Leaky ReLU. The
experiment results show that using Leaky ReLU can improve the classification
accuracy of our dataset. To learn global high-level feature representation,
the extracted features were fed into LSTM with an attention mechanism.

Our proposed CNN architecture consisted of three convolutional layers and two
max-pooling layers, and six batch normalization layers with the Leaky ReLU
activation functions, as shown in Fig. 5. Log-Mel spectrograms were fed into
our proposed CNN to extract high-level features for the classification task.
Firstly, we utilized three convolutional layers with 3$\times$3 filter kernel
sizes. Three convolutional layers have 128, 128, and 256 feature maps,
respectively. Secondly, we added a max-pooling layer with 2$\times$4 filter
kernel sizes after the first two convolutional layers. We also added a pair of
batch normalization (BN) layers with Leaky ReLU before and after the
convolutional layers to normalize the features and reduce over-fitting.

Figure 5: Our proposed CNN with Leaky ReLU.

The equation for ReLU is $f(x)=max(0,x)$. When the input of the layer is
negative, the ReLU is equal to zero. Consequently, gradient descents reach the
value of zero and cannot converge to the local minimum. For Leaky ReLU, there
is always a small slope to allow the weight update of the accumulated
gradient. Therefore, although ReLU can compute faster, we used Leaky ReLU
instead of the ReLU so that the layers did not stop learning when the slope of
the ReLU is zero. The Leaky ReLU activation function [26] is described by the
equation 1:

$f(x)=\begin{cases}x,&\text{if}\ x>0\\\ \alpha
x,&\text{otherwise},\end{cases}$ (1)

where $\alpha$ was set to 0.3 in this research.

###  2.4 Global Feature Learning using LSTM and Attention Mechanism

In this article, we utilized LSTM [27] to learn sequential feature maps that
are extracted from our proposed CNN. Each LSTM unit can be updated as in
equations 2 – 7:

$f_{t}=\sigma\big{(}W_{f}X_{t}+U_{f}h_{t-1}+b_{f}\big{)},$ (2)
$i_{t}=\sigma\big{(}W_{i}X_{t}+U_{i}h_{t-1}+b_{i}\big{)},$ (3)
$o_{t}=\sigma\big{(}W_{o}X_{t}+U_{o}h_{t-1}+b_{o}\big{)},$ (4)
$\tilde{c}_{t}=\tau\big{(}W_{c}X_{t}+U_{c}h_{t-1}+b_{c}\big{)},$ (5)
$c_{t}=f_{t}\odot\tilde{c}_{t-1}+i_{t}\odot\tilde{c}_{t},$ (6)
$h_{t}=o_{t}\odot\tau(c_{t}),$ (7)

where all variables are defined as follows:

  * •

$X_{t}$ is the mini-batch input;

  * •

$i_{t}$ is the input gate;

  * •

$f_{t}$ is the forget gate;

  * •

$o_{t}$ is the output gate;

  * •

$\tilde{c}_{t}$ is the input cell;

  * •

$c_{t}$ is the cell state;

  * •

$h_{t}$ is the hidden state;

  * •

$\sigma$ is the _sigmoid_ function;

  * •

$\tau$ is the _tanh_ function;

  * •

$W$, $U$ are the weight matrices;

  * •

$b$ is the bias parameter;

  * •

$t$ is the time step.

Due to different frame-level features contribute unequally to classifying
event sound classes. Hence, we added an attention mechanism layer [28] after
LSTM to specific points in a sequence when computing its output. Additionally,
during the transition from the normal state of the drill bit to the broken
state, the pitch of the audio changes. Therefore, the features extracted from
the log-Mel spectrogram right at the moment the drill bit cracks will have an
abnormality. The purpose of the attention layer is to focus on that anomaly.
The attention mechanism has been widely used in the sequence-to-sequence
model. For the LSTM, the output of attention $att$ can be defined as below:

$att=\sum_{t=1}^{T}\alpha_{t}h_{t},$ (8)

where $h_{t}$ denotes the $t_{th}$ hidden output from the LSTM at time step
$t$, $T$ represents the total number of time steps in the input sequence, and
the $\alpha_{t}$ is the attention weight can be computed as follows:

$\alpha_{t}=\frac{\exp{\big{(}W\cdot
h_{t}\big{)}}}{\sum_{k=1}^{T}\exp{\big{(}W\cdot h_{t}\big{)}}}.$ (9)

##  3 Experimental Results

###  3.1 Experimental Setup

####  3.1.1 Datasets

After applying time-shifting and volume control augmentation methods to 201
original sounds from three categories (anomaly, normal, and irrelevant
sounds), we obtained 603 sounds. To train our end-to-end model, we converted
these sounds in the augmented dataset into log-Mel spectrograms. We used
around 70% of the dataset (420 log-Mel spectrograms) for training and 30% (183
log-Mel spectrograms) for testing, in which when training the model, 420 was
divided by the ratio of 70/30 again for training and validation sets.

####  3.1.2 Hyper-parameters and training setup

The proposed deep learning model was trained on Intel CORE i5 8th Gen with
NVIDIA Graphics Card 1050Ti. We used Keras library [29] with TensorFlow
toolkit [30] that are popular deep learning frameworks to implement and deploy
our proposed deep learning model. Additionally, the Librosa library [31] was
used to generate log-Mel spectrograms from original drill sounds.

For hyper-parameters optimization, we used the Adam optimizer [32] with a
learning rate of 0.001, a batch size of 4, a momentum of 0.9, and 100 epochs.
During training, categorical cross-entropy was used as the loss function
$L_{f}$ to update the network weights. It is defined as follows:

$L_{f}=-\sum_{n=1}^{nC}y_{n}\log\big{(}\hat{y}_{n}\big{)},$ (10)

where $nC$ is the number of classes, $y_{n}$ is the ground truth, and
$\hat{y}_{n}$ is the predicted class probabilities for the $n_{th}$ element of
model predictions. Furthermore, to avoid over-fitting and to improve the
generalized model, we applied early stopping to train the network with the
patience of 5.

###  3.2 Results

As shown in Table III, the overall accuracy of our proposed method, CNN using
the Leaky ReLU activation function combined with attention-based LSTM (CNN -
LSTM - Attention - Leaky ReLU), was 92.35%. The confusion matrix for our
proposed method is in Fig. 6. Table II shows the F1-score, precision, and
recall for each class in the augmented dataset.

Our previous work in [1] used the Continuous wavelet transform (CWT) and the
Mel spectrogram as the input of machine learning classifiers (linear
discriminant, SVM, and KNN). In contrast, in this paper, we proposed an end-
to-end model that utilized log-Mel spectrograms as inputs and produced an
accuracy of 92.35%, whereas the previous work produced an accuracy of 80.25%.

Figure 6: The confusion matrix for the proposed model (CNN - LSTM - Attention
- Leaky ReLU) on the augmented dataset. Table 2: The classification results of
CNN with Leaky ReLU in conjunction with attention-based LSTM

.

| Precision | Recall | F1-score | No.
---|---|---|---|---
Anomaly sound (“Broken”) | 0.96 | 0.87 | 0.91 | 61
Normal sound (“Normal”) | 1.00 | 0.93 | 0.97 | 61
Irrelevant sound (“Other”) | 0.83 | 0.97 | 0.89 | 61
Accuracy |  |  | 0.92 | 183
Macro avg | 0.92 | 0.93 | 0.92 | 183
Weighted avg | 0.92 | 0.93 | 0.92 | 183

###  3.3 Discussion

This part aims to demonstrate that data augmentation, the Leaky ReLU and the
attention mechanism can affect the overall accuracy of the proposed method
when combined with CNN and LSTM. The mean accuracy of all the experiments is
shown in Table III for comparison. Using CNN with Leaky ReLU activation
function in conjunction with attention-based LSTM achieves the highest
accuracy of 92.35%. The following experiments were conducted to validate the
effectiveness of our proposed method:

Table 3: The comparison of different models. Model | Mean Accuracy (%)
---|---
CNN - Leaky ReLU | 69.95
CNN - LSTM - Leaky ReLU | 86.89
CNN - LSTM - Attention - Leaky ReLU | 92.35
CNN - LSTM - Attention - ReLU | 90.16

####  3.3.1 CNN - Leaky ReLU

In CNN architecture, we run experiments with the Leaky ReLU activation
function, which has layers similar to Fig. 5. The experiment parameters were
identical to the CNN architecture in our proposed method. However, we used two
fully connected layers at the end of CNN for classification. According to
Table V, the overall accuracy for this method was only 69.95%, which is lower
than the accuracy of our proposed method (92.35%). Fig. 7 shows the confusion
matrix for this method.

Figure 7: CNN - Leaky ReLU.

####  3.3.2 CNN - LSTM - Leaky ReLU

This part experimented with CNN uses the Leaky ReLU activation function in
conjunction with LSTM. This experiment tests whether incorporating an
attention layer into the model is effective. In this method, accuracy achieved
86.89%, which is less than what we propose in section IIIB (accuracy of
92.35%). It is clear that the accuracy of the model was improved by including
the attention layer. In theory, with the attention layer, the LSTM is supposed
to invest more computing power of that small but important part of the input,
so the network enhances these parts and fades out the rest. The confusion
matrix for this method is shown in Fig. 8.

Figure 8: CNN - LSTM - Leaky ReLU.

####  3.3.3 CNN - LSTM - Attention - ReLU

This part experimented with CNN architecture uses the ReLU activation function
in conjunction with attention-based LSTM. The confusion matrix for this method
is shown in Fig. 9. In this experiment, we run the model with ReLU activation
to prove it is less effective than Leaky ReLU activation on our dataset. When
using the ReLU activation function, the accuracy was 90.16%, while using Leaky
ReLU, the accuracy was higher (92.35%). As Leaky ReLU has a slope of 0.3
instead of 0, CNN can train faster and avoid the “dying ReLU” problem on our
dataset.

Figure 9: CNN - LSTM - Attention - ReLU.

####  3.3.4 The comparison between using the proposed method on the original
dataset and the augmented dataset

To test the efficiency of the data augmentation process, the proposed model in
section II was run on both the original and augmented datasets. The results
are shown in Table IV, the accuracy on the augmented dataset reached 92.35%
whereas the accuracy on the original dataset only reached 84.13%. The accuracy
of our proposed method on the augmented dataset (603 sounds) is clearly higher
than on the original dataset (201 sounds). The confusion matrix for our
proposed method on the original dataset is shown in Fig. 10.

Figure 10: The proposed method (CNN - LSTM - Attention - Leaky ReLU) on the original dataset. Table 4: The comparison of the original and augmented datasets using the same proposed method (CNN - LSTM - Attention - Leaky ReLU). Dataset | Mean Accuracy (%)
---|---
Augmented dataset (603 sounds) | 92.35
Original dataset (201 sounds) | 84.13

##  4 Conclusion

In this article, we proposed a deep learning model for drill fault detection.
We applied time-shifting and volume control augmentation methods to increase
the number of sounds in the small dataset. The sounds in the augmented dataset
were converted into log-Mel spectrograms. These log-Mel spectrograms were used
to train our proposed CNN architecture with the Leaky ReLU activation function
in conjunction with attention-based BiLSTM for detecting drill failure. It was
found that the overall accuracy of our proposed system reached 92.35%. This
method has a huge potential to be used to diagnose faults in industrial
machines. It is a non-invasive method of diagnosing machine failure using
short sounds or small datasets. For future work, we consider using both sounds
and images to detect the drill failure to obtain better drill fault detection
results.

## References

  * [1] T. Tran and J. Lundgren, “Drill fault diagnosis based on the scalogram and mel spectrogram of sound signals using artificial intelligence,” _IEEE Access_ , vol. 8, pp. 203 655–203 666, 2020.
  * [2] Y. J. Choi, M. S. Park, and C. N. Chu, “Prediction of drill failure using features extraction in time and frequency domains of feed motor current,” _International Journal of Machine Tools and Manufacture_ , vol. 48, no. 1, pp. 29–39, 2008.
  * [3] P. Skalle, A. Aamodt, and O. E. Gundersen, “Detection of symptoms for revealing causes leading to drilling failures,” _SPE Drilling & Completion_, vol. 28, no. 02, pp. 182–193, 2013.
  * [4] A. Kumar, J. Ramkumar, N. K. Verma, and S. Dixit, “Detection and classification for faults in drilling process using vibration analysis,” in _2014 International Conference on Prognostics and Health Management_. IEEE, 2014, pp. 1–6.
  * [5] P. Henriquez, J. B. Alonso, M. A. Ferrer, and C. M. Travieso, “Review of automatic fault diagnosis systems using audio and vibration signals,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 44, no. 5, pp. 642–652, may 2014. [Online]. Available: http://ieeexplore.ieee.org/document/6670088/
  * [6] A. Glowacz, “Fault detection of electric impact drills and coffee grinders using acoustic signals,” _Sensors_ , vol. 19, no. 2, p. 269, 2019.
  * [7] Y. Jing, H. Su, S. Wang, W. Gui, and Q. Guo, “Fault diagnosis of electric impact drills based on time-varying loudness and logistic regression,” _Shock and Vibration_ , vol. 2021, 2021.
  * [8] J. Hou, H. Sun, A. Xu, Y. Gong, and D. Ning, “Fault diagnosis of synchronous hydraulic motor based on acoustic signals,” _Advances in Mechanical Engineering_ , vol. 12, no. 4, p. 1687814020916107, 2020.
  * [9] X. Wang, D. Mao, and X. Li, “Bearing fault diagnosis based on vibro-acoustic data fusion and 1d-cnn network,” _Measurement_ , vol. 173, p. 108518, 2021\.
  * [10] K. Polat, “The fault diagnosis based on deep long short-term memory model from the vibration signals in the computer numerical control machines,” _Journal of the Institute of Electronics and Computer_ , vol. 2, no. 1, pp. 72–92, 2020.
  * [11] A. González-Muñiz, I. Díaz, and A. A. Cuadrado, “Dcnn for condition monitoring and fault detection in rotating machines and its contribution to the understanding of machine nature,” _Heliyon_ , vol. 6, no. 2, p. e03395, 2020.
  * [12] D. Verstraete, A. Ferrada, E. L. Droguett, V. Meruane, and M. Modarres, “Deep learning enabled fault diagnosis using time-frequency image analysis of rolling element bearings,” _Shock and Vibration_ , vol. 2017, 2017.
  * [13] S. Zhang, S. Zhang, B. Wang, and T. G. Habetler, “Deep learning algorithms for bearing fault diagnostics—a comprehensive review,” _IEEE Access_ , vol. 8, pp. 29 857–29 881, 2020.
  * [14] Z. Chen, X. Chen, C. Li, R.-V. Sanchez, and H. Qin, “Vibration-based gearbox fault diagnosis using deep neural networks,” _Journal of Vibroengineering_ , vol. 19, no. 4, pp. 2475–2496, 2017.
  * [15] M. M. Islam and J.-M. Kim, “Motor bearing fault diagnosis using deep convolutional neural networks with 2d analysis of vibration signal,” in _Canadian Conference on Artificial Intelligence_. Springer, 2018, pp. 144–155.
  * [16] L. Xueyi, L. Jialin, Q. Yongzhi, and H. David, “Semi-supervised gear fault diagnosis using raw vibration signal based on deep learning,” _Chinese Journal of Aeronautics_ , vol. 33, no. 2, pp. 418–426, 2020.
  * [17] X. Zhang, Y. Zou, and W. Shi, “Dilated convolution neural network with LeakyReLU for environmental sound classification,” _International Conference on Digital Signal Processing, DSP_ , vol. 2017-Augus, 2017.
  * [18] N. Davis and K. Suresh, “Environmental sound classification using deep convolutional neural networks and data augmentation,” in _2018 IEEE Recent Advances in Intelligent Computational Systems (RAICS)_ , 2018, pp. 41–45.
  * [19] V. Boddapati, A. Petef, J. Rasmusson, and L. Lundberg, “Classifying environmental sounds using image recognition networks,” _Procedia Computer Science_ , vol. 112, pp. 2048–2056, 2017. [Online]. Available: http://dx.doi.org/10.1016/j.procs.2017.08.250
  * [20] Z. Mushtaq, S. F. Su, and Q. V. Tran, “Spectral images based environmental sound classification using CNN with meaningful data augmentation,” _Applied Acoustics_ , vol. 172, p. 107581, 2021. [Online]. Available: https://doi.org/10.1016/j.apacoust.2020.107581
  * [21] F. Medhat, D. Chesmore, and J. Robinson, “Masked conditional neural networks for environmental sound classification,” _arXiv_ , vol. 608014, no. 608014, pp. 21–33, 2018.
  * [22] Y. Chen, Q. Guo, X. Liang, J. Wang, and Y. Qian, “Environmental sound classification with dilated convolutions,” _Applied Acoustics_ , vol. 148, pp. 123–132, 2019. [Online]. Available: https://doi.org/10.1016/j.apacoust.2018.12.019
  * [23] Z. Zhang, S. Xu, S. Zhang, T. Qiao, and S. Cao, “Attention based convolutional recurrent neural network for environmental sound classification,” _Neurocomputing_ , no. xxxx, 2020. [Online]. Available: https://doi.org/10.1016/j.neucom.2020.08.069
  * [24] H. Wang, Y. Zou, D. Chong, and W. Wang, “Environmental Sound Classification with Parallel Temporal-spectral Attention,” 2019. [Online]. Available: http://arxiv.org/abs/1912.06808
  * [25] K. Choi, G. Fazekas, M. Sandler, and K. Cho, “A comparison of audio signal preprocessing methods for deep neural networks on music tagging,” in _2018 26th European Signal Processing Conference (EUSIPCO)_. IEEE, 2018, pp. 1870–1874.
  * [26] A. L. Maas, A. Y. Hannun, and A. Y. Ng, “Rectifier nonlinearities improve neural network acoustic models,” _in ICML Workshop on Deep Learning for Audio, Speech and Language Processing_ , vol. 28, 2013.
  * [27] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [28] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Advances in neural information processing systems_ , 2017, pp. 5998–6008.
  * [29] F. Chollet _et al._ , “Keras,” https://github.com/fchollet/keras, 2015\.
  * [30] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015, software available from tensorflow.org. [Online]. Available: https://www.tensorflow.org/
  * [31] B. McFee, A. Metsai, M. McVicar, S. Balke, C. Thomé, C. Raffel, F. Zalkow, A. Malek, Dana, K. Lee, O. Nieto, D. Ellis, J. Mason, E. Battenberg, S. Seyfarth, R. Yamamoto, viktorandreevichmorozov, K. Choi, J. Moore, R. Bittner, S. Hidaka, Z. Wei, nullmightybofo, D. Hereñú, F.-R. Stöter, P. Friesch, A. Weiss, M. Vollrath, T. Kim, and Thassilo, “librosa/librosa: 0.8.1rc2,” May 2021. [Online]. Available: https://doi.org/10.5281/zenodo.4792298
  * [32] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2015. [Online]. Available: http://arxiv.org/abs/1412.6980