# An X-ray study for white dwarf binary AR Scorpii

J. Takata11affiliation: Department of Astronomy, School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China , X.F.
Wang11affiliation: Department of Astronomy, School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China , H.H.
Wang11affiliation: Department of Astronomy, School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China , L.C.-C.
Lin22affiliation: Department of Physics, UNIST, Ulsan 44919, Korea , C.-P.
Hu33affiliation: Department of Physics, National Changhua University of
Education, Changhua 50007, Taiwan , K.L. Li44affiliation: Department of
Physics, National Cheng Kung University, 70101 Tainan, Taiwan , and A.K.H.
Kong55affiliation: Institute of Astronomy, National Tsing Hua University,
Hsinchu 30013, Taiwan<EMAIL_ADDRESS>
###### Abstract

We report a study of the X-ray emission from the white dwarf/M-type star
binary system AR Scorpii using archival data taken in 2016-2020. It has been
known that the X-ray emission is dominated by the optically thin thermal
plasma emission, and its flux level varies significantly over the orbital
phase. The X-ray emission also contains a component that modulates with the
beat frequency between the white dwarf’s spin frequency and orbital frequency.
In this new analysis, the 2020 data taken by NICER shows that the X-ray
emission is modulating with the spin frequency as well as the beat frequency,
indicating that part of the X-ray emission is coming from the white dwarf’s
magnetosphere. It is found that the signal of the spin frequency appears only
at a specific orbital phase, while the beat signal appears over the orbital
phase. We interpret the X-ray emission modulating with the spin frequency and
the beat frequency as a result of the synchrotron emission from electrons with
a smaller and larger pitch angle, respectively. In a long-term evolution, the
beat pulse profile averaged over the orbital phase changed from a single-peak
structure in 2016/2018 to a double-peak structure in 2020. The observed X-ray
flux levels measured in 2016/2017 are higher than those measured in 2018/2020.
The plasma temperature and amplitude of the orbital waveform might vary with
time too. These results indicate that the X-ray emission from AR Scorpii
evolves on a timescale of years. This long-term evolution would be explained
by a super-orbital modulation related to, for example, a precession of the
white dwarf, or a fluctuation of the system related to activity of the
companion star.

##  1\. Introduction

AR Scorpii is a binary system composed of a white dwarf (hereafter WD) and a
M-type main-sequence star, and its distance from the Earth is estimated to be
$d\sim 116$pc (Marsh et al., 2016; Bailer-Jones et al., 2018; Peterson et al.,
2019). The emission from radio to X-ray bands is modulating with the orbital
period of $P_{o}\sim 3.56$ hours (Marsh et al., 2016; Stanway et al., 2018;
Takata et al., 2018). Marsh et al. (2016) report that in addition to the
orbital modulation, the radio/optical emission from AR Scorpii is modulating
in a beat frequency, $\nu_{b}=\nu_{s}-\nu_{o}\sim 8.46$ mHz, where
$\nu_{s}\sim 8.54$ mHz is the spin frequency of the WD and $\nu_{o}\sim 0.078$
mHz is the orbital frequency. AR Scorpii is the first WD that shows a
modulation of the radio emission related to the WD’s spin. The pulsed emission
up to soft X-ray bands has been confirmed (Takata et al., 2018), and no
confirmation of the GeV emission has been reported (Kaplan et al., 2019; Singh
et al., 2020). A spin down of the WD has been measured by the optical
observations (Marsh et al., 2016; Potter & Buckley, 2018a; Stiller et al.,
2018; Gaibor et al., 2020). With a longer base line of the optical
observations, Gaibor et al. (2020) measure a spin down rate as
$\dot{\nu}_{s}\sim(-4.82\pm 0.18)\times 10^{-17}{\rm~{}Hz~{}s^{-1}}$,
indicating the spin down power of $L_{sd}\sim 2\times
10^{33}{\rm~{}erg~{}s^{-1}}$.

Marsh et al. (2016) report that the optical emission is minimum around the
inferior conjunction of the companion star’s orbit, where the companion star
is between the WD and Earth, while it achieves the maximum level slightly
earlier than the superior conjunction, where the companion star is behind the
WD. Littlefield et al. (2017) find that the orbital modulation has been fairly
stable with 78 days of observations, but suggest with a longer base line
optical data that the orbital maximum shifts toward earlier phases relative to
the superior conjunction with time. Katz (2017) interprets this observed
orbital waveform and the long-term evolution with a precession model of the
WD’s spin axis. Peterson et al. (2019), on the other hand, report that the
orbital waveform of AR Scorpii has remained constant across about one century,
and constrain the first time derivative of the orbital frequency as
$\dot{\nu}_{o}\leq 3.8\times 10^{-20}{\rm~{}Hz~{}s^{-1}}$.

Marcote et al. (2017) carry out a radio observation with a high angular
resolution and find that the radio emission region of AR Scorpii is compact
with a scale of $\leq 4R_{\odot}$. They also report no evidence for a radio
outflow or collimated jets. Then, Stanway et al. (2018) confirm that the
orbital modulation in radio bands shows a similar trend to that of the optical
emission. This indicates that inner-surface of the companion star is probably
heated by an interaction with the WD magnetosphere or by an irradiation of
synchrotron radiation from the WD’s magnetosphere (Garnavich et al., 2019),
and the emission from radio to optical bands is mainly produced at the region
near or on the inner-surface.

The pulse profile folded with the beat frequency shows a double-peak structure
in the radio, optical and UV bands. In optical/UV bands, the pulse shape can
be described by a main peak plus a small secondary peak, for which the phase
separation is $\sim 0.5$ (Marsh et al., 2016). The pulsed fraction reaches
$>95$% in the UV band. Marsh et al. (2016) also find a signal of the spin
frequency of the WD in the periodograms, suggesting that the optical pulsed
emission originates from a region co-rotating with the WD. The optical
emission from AR Scorpii is also unique with a strong linear polarization.
Buckley et al. (2017) show that the degree of the linear polarization reaches
up to $\sim 40\%$ at the pulse peak, and the position angle swings through
$180^{\circ}$. These properties of the linear polarization indicate an origin
of synchrotron emission from the relativistic electrons spiraling along the
magnetic field of the WD (Potter & Buckley, 2018b; Takata & Cheng, 2019; du
Plessis et al., 2019).

In radio bands, Marsh et al. (2016) show the pulse profile can be described by
a double-peak structure with a pulsed fraction of $\sim 10\%$, which is
smaller than those in the optical/UV bands. Stanway et al. (2018) report that
the beat signal decreases with the radio frequency, and no beat signal is
found at 1.5GHz. The radio emission does not show an evident signal of the WD
spin in the periodograms, which indicates that the pulsed radio emission
originates from the region in the magnetosphere of the companion star rather
than in the magnetosphere of the WD. Stanway et al. (2018) also measure a weak
linear polarization with a degree $<1\%$ but a large circular polarization,
which reaches to $\sim 30$% at the orbital phase $\phi_{o}\sim 0.8$, where
$\phi_{o}=0$ is defined at the inferior conjunction of the companion star’s
orbit. They suggest that a non-relativistic cyclotron emission dominates the
radio emission at lower frequencies.

The broadband spectrum from the radio to UV bands can be described by a
thermal optical emission from the companion star plus non-thermal component
with a total luminosity of $L\sim 2\times 10^{32}{\rm~{}erg~{}s^{-1}}$, which
is about 10% of the spin down power. Garnavich et al. (2019) measure the slope
of the optical spectrum and find that the power-law index $s$,
$F_{\nu}\propto\nu^{-s}$, shows a significant variation in $0.8<s<2.0$
depending on the orbital phase and the optical brightness. They suggest that
the variation in the power-law index is related to the synchrotron cooling and
frequency of the replenishment of the energy injection. This non-thermal
component has an energy peak at around 0.01eV and probably extends down to
radio bands with a slope of $F_{\nu}\sim\nu^{1/3}$ (Stanway et al., 2018),
which is consistent with the slope of low energy tail of the synchrotron
emission (Geng et al., 2016). In the radio bands, the spectrum becomes
shallower below $\sim 1$GHz (Marsh et al., 2016; Stanway et al., 2018). This
also indicates that an additional component (probably cyclotron emission from
non-relativistic electrons) dominates at low radio frequency bands.

In the higher energy band, we (Takata et al., 2018) study the X-ray emission
from AR Scorpii with an XMM-Newton observation carried out in 2016. The
observed spectrum shows an iron-K line that is well described by an optically
thin thermal plasma emission with several different temperatures; there is no
signature of the emission from an accreting matter in the observed spectrum.
The observed X-ray flux varies significantly over the orbital phase (Figure
1), and it achieves maximum and minimum at around the superior conjunction and
inferior conjunction of the companion star, respectively. The hardness ratio
of the observed X-ray emission, on the other hand, is almost constant over the
orbital phase. These emission properties are compatible to the thermal process
of the plasmas heated to $0.1-10$ keV on or near the inner surface of the
companion star, and the orbital variation of the observed flux is caused by a
variation of the visible emission region over the orbital phase rather than an
absorption by the accretion column. The X-ray luminosity is of the order of
$L_{X}\sim 4\times 10^{30}~{}\rm{erg~{}s^{-1}}$, and it is $\sim 1$% of the
optical/UV luminosity of the non-thermal component. Hence, a small fraction of
the released energy is used for heating up the plasma to a temperature of
$\sim$keV, and majority is for creating a population of the non-thermal
electrons.

Takata et al. (2018) find a modulation with the beat frequency in the X-ray
data below $2$ keV. The pulse profile averaged over the orbital phase shows a
single-peak (or probably double-peak with a small secondary peak). On the
other hand, the orbital resolved analysis shows that the pulse structure
evolves over the orbital phase, and the pulse profile at around the inferior
conjunction of the companion star can be described by a double-peak structure
with a prominent secondary peak. The averaged pulsed fraction is about $14\%$,
suggesting that the unpulsed component is the main component in the X-ray
bands. By comparing the pulse profile measured by the EPIC camera with the
profile by the OM camera onboard XMM-Newton, we find that the X-ray peak and
optical/UV peak are in phase. This suggests that the non-thermal component in
optical/UV bands extends to the soft X-ray bands. The phase-resolved spectrum
of the pulsed component in the 0.15-2 keV band shows evidence of the power-law
component with a photon index $\Gamma=2.3\pm 0.5$ and an unabsorbed flux
$F_{\rm{0.15-2}}\sim 4\times 10^{-13}{\rm erg~{}cm^{-2}~{}s^{-1}}$.

After the discovery of the radio/optical pulsed emission from AR Scorpii in
2016, more X-ray data have been collected. It enables us to study a long-term
evolution of the X-ray emission from AR Scorpii. In this paper, we study the
X-ray emission properties of AR Scorpii using archival data taken by XMM-
Newton, Chandra, NuSTAR and NICER during 2016-2020 (Table 1). In section 2, we
summarize the data reduction for each observation. We present the results of
the timing analysis in section 3 and of spectral analysis in section 4. In
section 5, we discuss the X-ray emission process in the WD magnetosphere based
on the results of the study.

##  2\. Data reduction

Table 1Journal of the observation data used in this study. | Obs. ID | Date | Exposure |
---|---|---|---|---
|  |  | (ks) |
XMM | 0783940101 | 2016 Sep. 10 | 38.7 | X16
-Newton | 0795720101 | 2018 Feb. 19 | 30.0 | X18
Chandra | 19711 | 2017 June 23 | 26.1 | Ch17
NuSTAR | 30301027002 | 2018 Feb. 18 | 134.7 | Nu18
NICER | 3589010101 | 2020 June 14 | 5.6 |
| 3589010102 | 2020 June 14 | 19.3 |
| 3589010103 | 2020 June 16 | 15.4 | Ni20
| 3589010104 | 2020 June 16 | 9.5 |
| 3589010105 | 2020 June 18 | 6.9 |

###  2.1. $XMM$-Newton

We analyze the archival XMM-Newton data (MOS1/2 and PN) taken at 2016
September 10 (Obs. ID:0783940101, PI: Steeghs) and at 2018 February 19. (Ob.
ID: 0795720101, PI: Pavlov) with total exposure times of 39 ks and 30 ks,
respectively (Table 1). The observations were operated under the small window
mode for the MOS1/MOS2 and the large window mode for the PN in 2016 (hereafter
X16), while with full frame mode for all EPIC cameras in 2018 (hereafter X18).
Event list for each detector and observation is produced in the standard way
using XMM-Newton Science Analysis Software (`XMMSAS`, version 18.0.0). A point
source is significantly detected by the `XMMSAS` task `detect_chain`. We
generate the source spectra within a radius of $20^{\prime\prime}$ circle
centered at source position (R.A.,
Dec.)=($16^{\rm{h}}21^{\rm{m}}47^{\rm{s}}$.28,$-22^{\circ}53^{\prime}10^{\prime\prime}.4$)
and the background spectra from a source free region. When we generate the
spectra, we remove epoch with background flaring by limiting the count rate
above 10 keV to $<0.35~{}{\rm count~{}s^{-1}}$ for the MOS1/MOs2 data and
$<0.4~{}{\rm count~{}s^{-1}}$ for the PN data. We create response files with
the `XMMSAS` tasks `rmfgen` and `arfgen` with updated calibration files. We
group the spectral bins to achieve a signal-to-noise ratio of $S/N\geq 3$ in
each bin using the `XMMSAS` task `specgrp`. We use the `XMMSAS` task `barycen`
to obtain barycenter corrected arrival times (ephemeris DE405) and create the
background subtracted light curve with the task `epiclccorr` to examine the
orbital variation. We generate the orbital light curves in the 0.3-10keV band,
and compare it with the results of the Chandra and NICER data (Figure 1).

The X18 data recorded a large flare-like background event covering a
significant part of the observation. After removing such a background event,
the exposure of the data used in the spectral analysis is about $10$ ks for
the PN data and $20$ ks for the MOS1/MOS2 data, respectively. Figure 2 shows
the orbital phase distribution of the exposure time for the data used in the
spectral analysis; the maximum exposure time is normalized to unity. We can
see that a large exposure loss because of the background events causes a non-
uniform orbital coverage of the X18 data. For the X16 data, the exposure loss
due to the flare-like background events is less serious compared to the case
of the X18 data.

###  2.2. Chandra

A Chandra observation for AR Scorpii (Obs. ID: 19711, PI:Pavlov) was operated
under ACIS Timed Exposure (1/8 sub-array) mode and an exposure time of 24.8 ks
(hereafter Ch17). We re-process the data with standard tasks in the `CIAO`
Version 4.11. We use the `CIAO` task `specextract` to generate the source
spectrum and associated response files from a circular region centered at
aforementioned position with a radius $3.6^{\prime\prime}$, and create a
background spectrum from a vicinity source free region. It is known that pile-
up of a Chandra data is significant when the count rate is $\geq
0.1~{}\rm{count~{}s^{-1}}$ (Davis, 2001). In section 4.1.3, we examine the
pile-up fraction as $\sim 5$% for this observation. For a timing analysis, we
apply the `CIAO` task `axbary` to obtain barycen-corrected arrival time of
each event, and use `dmextract` task to create a background subtracted light
curve in the 0.3-10keV band to examine the orbital variation. We check the
background light curve and find no indication of flaring. For searching a beat
frequency in the data, we use the `Xselect` to extract events from the source
region and obtain a photon count of $\sim$4200.

###  2.3. NuSTAR

To examine the hard X-ray emission from AR Scorpii, we analyze the archival
NuSTAR data taken at 2018 February (Obs. ID: 30301027002, PI: Pavlov) with an
exposure $\sim$134 ks (hereafter Nu18). Event list and source/background
spectra are produced in the standard way using `nupipline` and `nuproducts`
tasks of the `Heasoft`. Barycenter correction is applied in the `nuproducts`
task. A point source is clearly seen in the generated image. We create the
source spectrum with the data extracted from a circular region with a radius
of $50^{\prime\prime}$ centered at aforementioned source position, and we
obtain a background spectrum from a vicinity source free region. We use
`grppha` task of the `FTOOLS` to group the spectral bins such that new
grouping contains a minimum of 30 counts in each bin. We create a background
subtracted light curve in the 3-78 keV band. We extract $\sim 10,000$ events
from the source region, and search for the beat frequency in the data.

###  2.4. NICER

NICER observed AR Scorpii at 2020 June 14-16 (Obs. IDs: 589010101-358901015
PI: Mori) with a total exposure $\sim 57$ ks (hereafter Ni20). We merge the
`mpu7_ufa` files for the five observations using the task `nimpumerge` of the
`Heasoft`, and screen the data using the tasks `nicerclean` and `nimaketime`.
We filter the raw data with the standard criteria for analysis of the
beat/spin modulation, and we apply an additional filtering based on the
magnetic cut-off rigidity to generate spectrum and orbital light curve. We
accept the time with $>1.5$GeV/c in current analysis and find that a more
tight filtering does not affect much to the result of the analyses. We use the
latest calibration files (`xti20200722`) and recalibrate the energy scale
using the file `nixtiflightpi20170601v004.fits`. We apply a barycentric time
correction to the created event list using `barycor` task of the `Heasoft`.
After standard screening, we also remove flare-like background events using
the `Xselect`.

We generate the spectra of the source region and background modeling with
`nibackgeng3C50` task implemented in the `HEASoft`; we create a filter file
with an option `coltypes="base,3c50"` in the task
`niprefilter2`111https://heasarc.gsfc.nasa.gov/lheasoft/ftools/headas/niprefilter2.html.
We use the `grppha` task to group the generated spectrum such that new
grouping contains a minimum of 40 counts in each bin. We note that the best-
fit parameters for the spectrum obtained with the original event file is
consistent with those obtained with the `nibackgeng3C50` product within the
errors. In this study, therefore, we will present the fitting result with the
`nibackgeng3C50` product.

Figure 2 shows that the 2020 NICER observation does not cover uniformly the
orbital phase. We generate an exposure corrected light curve in the 0.3-10 keV
band using the `Xselect`. Using the source and background spectra obtained by
the `nibackgeng3C50` task, we estimate that the background emission
contributes to $\sim 55$% of total emission in the 0.3-10 keV band. By
removing 55% of the total counts from the generated light curve, we create a
background subtracted light curve, which is compared with other observations
(Figure 1).

##  3\. Timing analysis

###  3.1. Properties of the orbital modulation

Figure 1.— Orbital modulation measured by each data set. Background
contribution is subtracted and exposure loss for each bin is corrected (Figure
2). For Ni20 data, we assume the background level as 55% of the total emission
in the 0.3-10 keV band (section 2.4). The blue solid curve in each panel
represents the best-fit function with one Gaussian component plus a constant
component.

Figure 2.— Distribution of the exposure fraction (filled circles) over the
orbital phase. Histograms shows the orbital variation after correction of the
exposure loss. For X18 data, the large exposure loss in the orbital phase
$\phi_{o}=0.4-1$ is caused by large flare-like background events. For Ni20
data, a non-continuous monitoring of the source produces a large fluctuation
of the exposure fraction over the orbital phase.

We create folded light curve with the orbital ephemeris reported by Marsh et
al. (2016): $\nu_{o}=0.077921380$ mHz and $T_{0}(\rm{MJD})=57264.09615$ when
the companion star is located at the inferior conjunction. Figure 1 shows the
orbital light curves in the 0.3-10 keV band for the X16/18, Ch17 and Ni20 data
and in the 3-78 keV band for the Nu18 data; the background level is subtracted
for each light curve. We can see that all light curves show a significant
orbital modulation, and the observed flux levels achieve the maximum and
minimum at around the superior conjunction ($\phi_{o}=0.5$) and inferior
conjunction ($\phi_{o}=0$), respectively. As discussed in Takata et al.
(2018), this pattern of the orbital modulation indicates that X-ray emission
from AR Scorpii mainly originates from the region located near or on the inner
surface of the companion star.

To examine a long-term evolution of the orbital waveform, we fit the obtained
light curve with one Gaussian component plus a constant component (Figure 1).
We define amplitude of the orbital variation as
$(f_{max}-f_{min})/(f_{max}+f_{min})$, where $f_{max}$ and $f_{min}$ are the
maximum and minimum counts obtained from the best-fit function, respectively
(İnam et al., 2004; Rea et al., 2007). We obtain amplitudes of $\sim 0.22\pm
0.01$ for the X16 data, $\sim 0.16\pm 0.03$ for the C1h7 data, $\sim 0.21\pm
0.01$ for the X18 data and $\sim 0.18\pm 0.02$ for the Nu18 data. For the Ni20
data, we obtain an amplitude of $\sim 0.32\pm 0.01$ by assuming $55$% of the
background level, and we find that the amplitude is larger than those of other
data. This would provide evidence of the long-term evolution of the orbital
waveform in the X-ray bands, although additional observations with an imaging
capability is necessary to confirm it. By comparing the amplitude of X18 data
with that of Nu18 data, we find less dependency of the orbital waveform on the
photon energy. This constant amplitude among different energy bands also
indicates that the orbital modulation is caused by a variation of the visible
emission region over the orbital phase.

###  3.2. Searching for pulsation

Figure 3.—  $Z_{2}^{2}$ periodogram for the XMM-Newton 2018 data. $\nu_{b}$
and $\nu_{s}=\nu_{b}+0.0792{\rm mHz}$ represents the beat frequency and spin
frequency of the WD, respectively.

Figure 4.—  $Z_{2}^{2}$ periodogram for the NICER 2020 data. $\nu_{b}$ and
$\nu_{s}$ represents the beat frequency and spin frequency of the WD,
respectively.

Figure 5.—  Pulse profiles in the 0.3-2 keV band folded with the beat
frequency for the three data sets. Each event is folded with the timing
solution described in section 3.3, and the background contribution is
subtracted for each light curve. The 55% of the background level in the total
emission is assumed for Ni20 data. “ON” and “OFF” in each panel indicate the
on-pulse phase and off-pulse phase to generate the spectrum of the pulsed
component, respectively (section 4.3).

We search the beat frequency $\nu_{b}\sim 8.6$mHz (Marsh et al., 2016) using
the $Z_{2}^{2}$-test (Buccheri et al., 1983). For the X16 and X18 data sets,
the extracted photon events from the source region yield $\sim 45,000$ and
$\sim 20,000$ counts, respectively. We use $\sim 4,200$ counts for Ch17, $\sim
10,000$ counts for Nu18 FPMA/FPMB and $\sim 135,000$ photons for Ni20 to
search the beat frequency. In addition to the X16 data reported by Takata et
al. (2018), a significant pulsed signal is found in X18 and Ni20 data sets.
Since the pulsed fraction is $\sim 14\%$ (Takata et al., 2018), the counts of
the pulsed emission in the Ch17 data will be only $\sim 600$ or less, and it
is unlikely to find a beat signal in the data. Although we combine all
observations, which include the Nu20 data, to search a pulsation in the $>2$
keV band, we do not find any significant signal at the beat frequency. This
non-detection of the beat frequency with the number of the collected photon
counts is expected, because only $\sim 2$% of the total emission from AR
Scorpii is pulsating in the $>2$ keV band (section 4.3).

Figures 3 and 4 present the $Z_{2}^{2}$-periodograms for the X18 and Ni20
data, respectively [see Takata et al. (2018) for the $Z_{2}^{2}$ periodogram
of the X16 data]. We can see a strong signal at the beat frequency
$\nu_{b}\sim 8.46$mHz and the 1st harmonic in the Ni20 data. In addition to
the beat frequency, we can also confirm a significant signal at the spin
frequency of the WD, $\nu_{s}\sim 8.53$mHz and the 1st/2nd harmonics. The
signal of the spin frequency with $Z_{2}^{2}\sim 125$ of the Ni20 data is much
stronger than $Z_{2}^{2}\sim 50$ of the X16 data (Takata et al., 2018) and
$Z_{2}^{2}\sim 20$ of the X18 data. The Ni20 data therefore provides a
concrete evidence that the pulsed X-ray emission originates from the WD
magnetosphere and its emission region is co-rotating with the WD.

###  3.3. Long-term evolution of pulse profile

Figure 6.—  Dynamic pulse profiles in the 0.3-2 keV band for X16 (left) and
Ni20 (right) data. The data is folded with t he beat phase and the pulse
profile at each orbital phase is created by the data of successive 0.1 orbital
phase. Total count (area of the profile) of each pulse profile is normalized
to unity.

Figure 7.— Orbital resolved pulse profiles in the 0.3-2 keV band for X16
(left) and Ni20 (right) data. Bottom right two panels show the pulse profiles
in the phase segments $\phi_{o}=0.1-0.6$ and $0.6-1.1$, respectively.

Figure 8.— Left: Dynamic pulse profile folded with the spin frequency for Ni20
data set. Right: Pulse profiles in the 0.3-2 keV band folded with the spin
frequency for Ni20 data set.

To compare between the pulse profiles measured by XMM-Newton and NICER
observations, we restrict the photon energy in the 0.3-2 keV band, and create
folded light curves with the beat ephemeris provided by Gaibor et al. (2020):
$\nu_{b}=8.460297018$mHz, $\dot{\nu}_{b}=-4.82\times 10^{-17}{\rm
Hz~{}s^{-1}}$ and $T_{0}(\rm{MJD})=57941.16888790$. For X16/18 data, we
determine the background level from a nearby source free region. For Ni50
data, we generate the background subtracted light curve by assuming that the
background level is 55% of the total emission.

Figure 5 shows the pulse profiles for the three data sets, and we find that
the pulse profile of the Ni20 data is significantly different from those of
the X16/18 data. The pulse profiles of the X16/X18 data show a single-peak
structure. For the N20 data, on the other hand, the pulse profile shows a
double-peak structure with a phase separation of $\delta\phi_{b}\sim 0.4-0.5$.
The difference in the pulse profiles between the X16/X18 and Ni20 data may
indicate a temporal evolution of the pulsed X-ray emission. Takata et al.
(2018), on the other hand, find that the pulse profile evolves over the
orbital phase, and it can be represented by either a single peak structure or
a double-peak structure. The pulse profile obtained with all data of each
observation therefore may depend on how the observation covers the orbital
phase. In fact, we can see a large fluctuation of the exposure time over the
orbital phase for the Ni20 data (Figure2).

We examine the orbital evolution of the pulse profiles to investigate the
cause of the difference between the pulse profiles of the X16/18 and Ni20
data. Figure 6 is a dynamic pulse profile folded with the beat frequency over
the orbital phase using the X16 data (left panel) and the Ni20 data (right
panel); the pulse profile at each orbital phase is created with the data of
successive 0.1 orbital phase, and the total count (area of the profile) is
normalized to unity. We can see that the dynamic pulse profiles in two data
sets are similar to each other, except for feature of the secondary peak
appeared at the beat phase $\phi_{b}\sim 0$ (equivalently $\phi_{b}\sim 1$)
and during the orbital phase $\phi_{o}\sim 0.5-1$.

Figure 7 shows the orbital resolved pulse profiles for the X16 (left) and Ni20
(right) data. In both data sets, the pulse profile at the orbital phase
$\phi_{o}\sim 0.3$ can be described by a single-peak structure appeared at the
beat phase $\phi_{b}\sim 0.5$. Then the peak position migrates to earlier beat
phase during the orbital phase $\phi_{o}\sim 0.3-0.6$, as the dynamic pulse
profiles indicate. After the orbital phase $\phi_{o}>0.6$, we can see a double
peak structure in both the X16 and N20 data. For the X16 data, the original
peak at the beat phase $\phi_{b}\sim 0.5$ always dominates the secondary peak
at $\phi_{b}\sim 0$. For Ni20 data, on the other hand, the secondary peak has
a similar or larger intensity compared to the original peak. In the bottom row
of the figure, we present the pulse profiles integrated over the orbital phase
$\phi_{o}=0.1-0.6$ and $\phi_{o}=0.6-1.1$, respectively. We find in the figure
that the pulse profiles in $\phi_{o}=0.1-0.6$ for the X16 and Ni20 both can be
described by a single-peak structure. For $\phi_{o}=0.6-1.1$, on the other
hand, the pulse profile of the Ni20 data shows a double-peak structure with a
secondary peak stronger than the original peak. For the Ni20 data, we can see
that the pulse profile averaged over the orbital phase using the normalized
pulse profiles presented in Figure 6 is also described by the double-peak
structure. These results therefore suggest that the difference between the
pulse profiles of the X16/X18 data and Ni20 data is an intrinsic feature, and
the pulsed X-ray emission experiences a long-term evolution.

We define the pulsed fraction with the equation,
$(f_{max}-f_{min})/(f_{max}+f_{min})$, with $f_{max}$ and $f_{min}$ being the
maximum and minimum counts in the light curve. In 0.3-2.0 keV bands, we obtain
$0.13\pm 0.02$ for X16, $0.11\pm 0.03$ for X18, and $0.20\pm 0.02$ for Ni20
data, respectively. We find that the pulsed fraction of the Ni20 data is
higher than those of the X16/X18 data. The pulsed fraction may show a long-
term evolution, as well as the structure of the pulse profile. Although we
could not find the pulsed emission above 2 keV bands, we estimate a pulsed
fraction less than 2% on the basis of the spectral analysis (section 4.3).

###  3.4. Mixture between spin and beat modulations

Figure 9.— Orbital resolved $Z_{2}^{2}$ periodogram for Ni20 data set.
$\nu_{b}$ and $\nu_{s}$ represents the beat frequency and spin frequency of
WD, respectively. The dotted lines in the left panel show the position of the
average between two frequencies.

In the dynamic pulse profile in Figure 6, we can see that a migration of the
pulse peak position at the orbital phase $\phi_{o}\sim 0.2-0.5$. This
migration suggests an existence of a signal with another period. In fact, the
$Z^{2}_{2}$ periodogram for the Ni20 data (Figure 4) clearly indicates that a
component of the X-ray emission from AR Scorpii modulates with not only the
beat frequency but also the WD’s spin frequency. Figure 8 presents dynamic
pulse profile folded with the spin frequency (left panel) and integrated pulse
profile (right panel) of the Ni20 data. In the dynamic pulse profile, we can
see a relatively constant peak position in the spin phase during the orbital
phase $\phi_{o}\sim 0.2-0.5$. This suggests that the signal of the WD’s spin
frequency appears in the X-ray data only at a specific orbital phase interval.
As the right panel in Figure 8 shows, the spin pulse profile integrated over
the orbital phase can be described by a single peak structure, which is
different from the double-peak structure of the beat pulse profile.

We divide the Ni20 data into two orbital intervals, $\phi_{o}=0.1-0.6$ and
$\phi_{o}=0.6-1.1$, and perform the $Z^{2}_{2}$-test for each interval (Figure
9). The periodogram of the orbital phase $\phi_{o}=0.1-0.6$ shows the signal
of the spin frequency that is stronger than the signal of the beat frequency.
During $\phi_{o}=0.1-0.6$, moreover, we can see a strong signal at the average
of two frequencies (dotted vertical lines), and the first harmonics of the
averaged frequency has a power stronger than those of the $2\nu_{b}$ and
$2\nu_{s}$ in the periodogram. This supports the fact that the short-term
modulation during $\phi_{o}=0.1-0.6$ is caused by the mixture between the beat
frequency and spin frequency. For $\phi_{o}=0.6-1.1$ (right panel in Figure
9), on the other hand, no signal of the spin frequency is found in the
$Z_{2}^{2}$ periodogram, and the short-term modulation is dominated by the
beat frequency. We find that the fist harmonic of the beat frequency during
the orbital phase $\phi_{o}=0.6-1.1$ shows a $Z_{2}^{2}$ power similar to or
stronger than the power of the fundamental frequency. This is consistent with
the double-peak structure with a phase separation $\sim 0.5$ appeared in the
pulse profile folded with the beat frequency.

We note that the dynamic pulse profile of the X16 data (left panel in Figure
6) also shows a migration of the peak position and therefore indicates an
appearance of the spin frequency at the orbital phase interval $\phi_{o}\sim
0.2-0.5$. The similar behavior of the X16 and N20 data suggests that the
appearance of signal of the WD’s spin frequency only at a specific orbital
phase interval is regularly repeated by every orbit. In section 5, we discuss
the mixture between the spin and beat frequencies in the short-term modulation
and its dependency on the orbital phase.

##  4\. Spectral analysis

The spectrum of the X-ray emission from AR Scorpii can be described by an
optically thin thermal plasma emission with several different temperatures
(Takata et al., 2018). In this study, we fit the generated spectra with a
multi-temperature VMEKAL model (Mewe et al., 1985, 1986) using `Xspec` version
12.10.1f. We apply `TBABS` model for modeling photoelectric absorption.

Table 2Fittings for 2016 XMM-Newton data (0.2-12 keV bands) with three-temperature VMEKAL models.  | Takata et al. (2018) aaExcept for iron, the abundances are fixed to the solar abundances [Anders and Grevesse (1989)] | Model 1 aaExcept for iron, the abundances are fixed to the solar abundances [Anders and Grevesse (1989)] | Model 2 bbThe abundances of O, Ne, Mg, S and Si are linked to that of Fe.
---|---|---|---
$N_{\rm H}$ ($10^{20}{\rm cm^{-2}}$) | $3.4^{+0.8}_{-0.8}$ | $1.2^{+0.3}_{-0.2}$ | $2.8^{+0.5}_{-0.5}$
$F_{e}$ | $0.67^{+0.29}_{-0.21}$ | $0.45^{+0.07}_{-0.07}$ | $0.28^{+0.06}_{-0.05}$
$kT_{1}$ (keV) | $8.0^{+2.8}_{-1.6}$ | $7.6^{+6.0}_{-1.9}$ | $6.3^{+1.3}_{-1.0}$
$kT_{2}$ (keV) | $1.7^{+0.4}_{-0.3}$ | $2.4^{+0.9}_{-0.7}$ | $1.5^{+0.2}_{-0.2}$
$kT_{3}$ (keV) | $0.6^{+0.1}_{-0.1}$ | $0.63^{+0.04}_{-0.05}$ | $0.32^{+0.04}_{-0.03}$
$\chi_{red}^{2}$ (dof) | 1.03 (404) | 1.20 (406) | 1.08 (406)

###  4.1. Analysis for individual data set

|
---|---
|
Figure 10.— X-ray spectrum of AR Scorpii. Top left: 2016 XMM-Newton
observation. The black, red and green correspond to PN,MOS1 and MOS2,
respectively. Top right: 2017 Chandra observation. Bottom left: 2018 XMM-
Newton/NuSTAR observations. For XMM-Newton data, the black, red and green show
the spectra of PN, MOS1 and MOS2, respectively. For NuSTAR data, the blue and
cyan correspond to FPMA and FPMB data, respectively. Bottom right: 2020 NICER
observation.

####  4.1.1 Reanalysis for X16 data

Although Takata et al. (2018) fit the spectrum of the X16 data with a three-
temperature VMEKAL model, we reanalyze the data with the updated calibration
files of XMM-Newton observation. Table 2 compares the best-fit parameters
reported in Takata et al. (2018) with those in the new analysis. In Takata et
al. (2018) (the second column), a reasonable fit is obtained by thawing iron
(Fe) abundance and by fixing other elements to the solar abundance. In our new
analysis, however, we find that the reduced $\chi^{2}$ of the fitting with the
model (Model 1 in third column) becomes worse comparing to the previous
result, and the best-fit parameters show a smaller hydrogen column density
$(N_{H}$).

To find a better fitting model, we thaw each element step by step. We find
that abundances of another five elements (O, Ne, Mg, S and Si) smaller than
the solar abundance significantly improve the fitting results with a $F$
statistic value of 8.2, which means that the probability of this improvement
being caused by chance is $2.3\times 10^{-7}$. Since the predicted abundances
of those components are consistent with that of the iron within the errors,
Table 2 (Model 2) shows the best-fit parameters by linking the abundances
between the iron and the five elements. We find that the hydrogen column
density and two higher temperatures are consistent with previous result, but
the abundance of the Fe and the lowest temperature decrease from the previous
one.

####  4.1.2 NICER data

Since NICER has no imaging capability, we check whether the spectra obtained
with the process described in section 2.4 is indeed consistent with the X-ray
emission from AR Scorpii. We restrict the photon energy in the 0.3-2 keV band
for the spectral analysis of the NICER data because of strong background
contamination in lower/higher energy bands. We use the Model 2 of Table 2, in
which the abundances of O, Ne, Mg, S, Si and Fe elements are linked. We find
that the background subtracted spectrum is reasonably fitted by a two-
temperature VMEKAL model with the best-fit hydrogen column density
$N_{H}=4.6^{+0.8}_{-0.7}\times 10^{20}{\rm cm^{-2}}$, temperatures
$kT_{1}=2.6^{+0.4}_{-0.4}$ keV and $kT_{2}=0.67^{+0.09}_{-0.06}$ keV
($\chi_{red}^{2}=1.15$ for 163 dof: third temperature component is not
necessary). The Fe abundance of $0.27^{+0.11}_{-0.08}$ (solar abundance) is
consistent with the result of the X16 data in Table 2. A simple black-body
radiation model or a power-law model cannot provide a reasonable fit. Hence,
the generated spectrum likely represents the X-ray emission from AR Scorpii.

The pulse profiles folded with the beat frequency of the Ni20 data shows a
single-peak structure in $\phi={0.1-0.6}$ and a double-peak structure in
$\phi={0.6-1.1}$ (Figure 7). We therefore create the orbitally resolved
spectra for $\phi_{o}={0.1-0.6}$ and $\phi_{o}={0.6-1.1}$, and we check any
orbital evolution of the best-fit parameters with the two-temperature VMEKAL
model. We do not find a significant difference in the best-fit parameters of
two orbital phase intervals. This is not unexpected because the pulsed
emission is only $\sim$10% of the total emission, and the spectrum is
dominated by the un-pulsed component.

####  4.1.3 Pile-up fraction of Chandra data

First, we use the `Web PIMMS` Version 4.11 to estimate pile-up fraction of the
Chandra data. Since the `Web PIMMS` implements the APEC222http://atomdb.org/
model for the plasma emission, we fit the generated spectrum with a two-
temperature APEC model with the `Xspec`. Since the data cannot constrain the
hydrogen column density, we fix it at $N_{H}=3.0\times 10^{20}{\rm cm^{-2}}$.
A reasonable fit is obtained with the temperatures $kT_{1}=5.1^{+1.2}_{-0.5}$
keV, $kT_{2}=1.0^{+0.4}_{-0.3}$ keV and the iron abundance $0.6^{+0.2}_{-0.3}$
($\chi_{red}^{2}=1.0$ for 107 dof). With the frame time of the observation,
$0.8$ s, the `Web PIMMS` estimates the pile-up fraction as $\sim 5$ percent.
Second, we use the `CIAO`’s fitting package, `Sherpa`, to estimate the pile-up
fraction. We fit the spectrum with a two-temperature VMEKAL components plus
pile-up model (Davis, 2001). The fitting result shows the plasma temperatures
$kT_{1}=4.1^{+0.5}_{-0.5}$ keV, $kT_{2}=0.8^{+0.3}_{-0.4}$ keV and abundance
Fe$=0.67^{+0.14}_{-0.13}$. The estimated pile-up fraction is $\sim 5$ percent
(so called $\alpha$-parameter and $f$-parameter are $\sim 0.8$ and $0.89$,
respectively), which is consistent with the result of the `Web PIMMS`.

####  4.1.4 Joint fitting for X18 and Nu18 data

Table 3Join fitting for the spectra (0.2-78 keV bands) of 2018 XMM-Newton and NuSTAR observations.  | Model 2 | Model 3 aaTwo VMEKAL + Power-law model. | Model 4 bbTwo VMEKAL + Power-law plus exponential cut-off model.
---|---|---|---
$N_{\rm H}$ ($10^{20}{\rm cm^{-2}}$) | $4.0^{+2.1}_{-1.4}$ | $3.4^{+1.6}_{-1.1}$ | $2.6^{+0.6}_{-0.6}$
$F_{e}$ | $0.22^{+0.05}_{-0.05}$ | $0.28^{+0.08}_{-0.08}$ | $0.31^{+0.07}_{-0.06}$
$kT_{1}$ (keV) | $6.8^{+0.8}_{-0.6}$ | $4.1^{+0.4}_{-0.4}$ | $2.7^{+0.3}_{-0.3}$
$kT_{2}$ (keV) | $1.5^{+0.2}_{-0.2}$ | $0.82^{+0.07}_{-0.07}$ | $0.66^{+0.14}_{-0.09}$
$kT_{3}$ (keV) | $0.34^{+0.32}_{-0.06}$ | - | -
$\alpha$ccPhoton index of the power-law component. | - | $1.6^{+0.1}_{-0.3}$ | $-0.20^{+0.5}_{-1.0}$
$E_{c}$ddCut-off energy of the power-law plus cut-off component. (keV) | - | - | $4.2^{+1.9}_{-1.1}$
$\chi_{red}^{2}$ (dof) | 1.07 (546) | 1.14 (546) | 1.07 (545)

The observations for the X18 and Nu18 data were carried out on 2018 February
18 and 19, respectively. This joint observation enables us to investigate
X-ray spectrum of AR Scorpii over a wide energy range in the 0.2-78keV bands.
In particular, we can investigate the behavior of the X-ray spectrum above
10keV bands. We perform a joint fit for the two spectra by introducing a
`constant` factor to account for cross-calibration mismatch between the XMM-
Newton and NuSTAR observations. During the fitting process, we fix the
`constant` factor for the X18 data at unity, and obtain a constant factor of
$\sim 1.1$ for the Nu18 data.

First we fit the spectrum with a three-temperature VMEKAL model (Model 2 in
Table 3), in which we link the abundances of six elements (O, Ne, Mg, S, Si
and Fe) and obtain the best-fit parameters as $N_{H}=4.0^{+2.1}_{-1.4}\times
10^{20}{\rm cm^{-2}}$, $kT_{1}=6.8^{+0.8}_{-0.6}$ keV,
$kT_{2}=1.5^{+0.2}_{-0.2}$ keV, $kT_{3}=0.34^{+0.32}_{-0.06}$ and
$F_{e}=0.22^{+0.05}_{-0.05}$ with $\chi_{red}^{2}=1.07$ for 546 dof., which
are consistent with the results of X16 data within the errors (Table 2).

To examine the spectral behavior above 10keV bands, we replace one VMEKAL
component to a power-law model, and fit the spectrum with
`constant*tbabs*(vmekal+vmekal+powerlaw)` (Model 3 in Table 3). We obtain the
best-fitting photon index of $\alpha=1.6^{+0.1}_{-0.3}$, but we find that the
$\chi_{red}^{2}\sim 1.14$ of the fitting becomes worse in comparison to the
three-temperature VMEKAL model. This would suggest that the observed spectrum
does not extend in the hard X-ray band with a single power-law form.

Finally we consider possibility of a spectral cut-off feature by replacing the
power-law component to a power-law plus an exponential cut-off model. We fit
the data with `constant*tbabs*(vmekal+vmekal+cutoffpl)` (Model 4 in Table 3),
and obtain an acceptable fit with $\chi_{red}^{2}\sim 1.07$ for 545 dof, which
is comparable to that of three VMEKAL model. The best-fit parameters, however,
show a very hard photon index of $\alpha\sim-0.20^{+0.5}_{-1.0}$ with a cut-
off energy $E_{c}\sim 4.2^{+1.9}_{-1.1}$ keV, and indicate that the spectrum
of the third component in higher energy band has a convex shape. This spectral
property is consistent with a thermal emission rather than a non-thermal
emission, and therefore the optically thin thermal plasma emissions with three
different temperatures describes better the X-ray emissions from AR Scorpii
over the 0.2-78keV bands.

###  4.2. Investigation for long-term evolution

Table 4Simultaneous fittings for the spectra (0.2-78 keV bands) with three-temperature VMEKAL models and by linking three temperatures. | X16 | Ch17 | X18/Nu18 | Ni20
---|---|---|---|---
$N_{\rm H}$ ($10^{20}{\rm cm^{-2}}$) | $3.8^{+0.4}_{-0.4}$
$F_{e}$ | $0.25^{+0.03}_{-0.03}$
$kT_{1}$ (keV) | $6.4^{+0.6}_{-0.4}$
$kT_{2}$ (keV) | $1.4^{+0.19}_{-0.05}$
$kT_{3}$ (keV) | $0.30^{+0.02}_{-0.02}$
$N_{1}$aaNormalization of the VMEKAL component in units of $10^{-14}/(4\pi d^{2})\int n_{e}n_{H}dV$ , where $d$ (cm) is the distance to the source. $(10^{-4})$ | $12.7^{+0.5}_{-1.0}$ | $13.1^{+1.0}_{-1.0}$ | $11.0^{+0.6}_{-1.1}$ |  $8.4^{+0.8}_{-0.9}$
$N_{2}$aaNormalization of the VMEKAL component in units of $10^{-14}/(4\pi d^{2})\int n_{e}n_{H}dV$ , where $d$ (cm) is the distance to the source. $(10^{-4})$ | $5.0^{+1.3}_{-0.7}$ | $6.9^{+1.6}_{-1.5}$ | $5.3^{+1.3}_{-0.8}$ | $4.6^{+0.9}_{-0.8}$
$N_{3}$aaNormalization of the VMEKAL component in units of $10^{-14}/(4\pi d^{2})\int n_{e}n_{H}dV$ , where $d$ (cm) is the distance to the source. $(10^{-4})$ | $2.5^{+0.4}_{-0.4}$ | $1.6^{+1.4}_{-1.4}$ | $1.6^{+0.4}_{-0.4}$ | $0.8^{+0.2}_{-0.2}$
$F_{0.3-10}$bbUnabsorbed flux in the 0.3-10 keV band in units of $10^{-12}~{}{\rm erg~{}s^{-1}cm^{-2}}$ . | $2.88^{+0.03}_{-0.03}$ | $3.05^{+0.09}_{-0.09}$ | $2.54^{+0.04}_{-0.04}$ | $1.97^{+0.03}_{-0.03}$
$\chi_{red}^{2}$ (dof) | 1.12 (1232)
Table 5Simultaneous fitting for the spectra (0.2-78 keV bands) with three-
temperature VMEKAL models and by unlinking lower two temperatures.

$N_{\rm H}$ ($10^{20}{\rm cm^{-2}}$) & $3.6^{+0.4}_{-0.3}$
$F_{e}$ $0.26^{+0.03}_{-0.03}$
$kT_{1}$ (keV) $6.8^{+0.6}_{-0.5}$

$kT_{2}$ (keV) $1.6^{+0.1}_{-0.2}$ $1.6^{+0.6}_{-0.3}$
$1.6^{+0.2}_{-0.2}$$1.6^{+0.3}_{-0.2}$

$kT_{3}$ (keV) $0.29^{+0.02}_{-0.02}$ $0.3^{+0.7}_{-0.2}$
$0.32^{+0.07}_{-0.04}$ $0.61^{+0..08}_{-0.11}$

$N_{1}$ $(10^{-4})$ $11.8^{+1.0}_{-1.0}$$11.9^{+1.6}_{-2.5}$
$10.1^{+1.1}_{-1.1}$ $7.6^{+1.6}_{-2.0}$

$N_{2}$ $(10^{-4})$ $6.0^{+1.1}_{-1.1}$ $8.0^{+2.6}_{-2.2}$
$6.2^{+1.3}_{-1.3}$ $4.9^{+1.9}_{-1.6}$

$N_{3}$ $(10^{-4})$ $2.2^{+0.4}_{-0.3}$ $1.7^{+100.0}_{-1.4}$
$1.4^{+0.4}_{-0.3}$ $1.0^{+0.3}_{-0.2}$

$F_{0.3-10}$ $2.86^{+0.03}_{-0.03}$ $3.05^{+0.09}_{-0.09}$
$2.51^{+0.04}_{-0.04}$ $1.93^{+0.04}_{-0.04}$

$\chi_{red}^{2}$ (dof) 1.09 (1226)

Table 6Fittings for 2016 and 2018 XMM-Newton data (0.2-12keV bands) with three-temperature MEKAL models. | 2016 | 2018
---|---|---
| 0-0.35/0.85-1 | 0.35-0.85 | 0-0.35/0.85-1 | 0.35-0.85
$N_{\rm H}$ ($10^{20}{\rm cm^{-2}}$) | $2.9^{+0.5}_{-0.4}$
$F_{e}$bbfootnotemark:  | $0.27^{+0.05}_{-0.04}$
$kT_{1}$ (keV) | $7.0^{+1.5}_{-1.2}$
$kT_{2}$ (keV) | $1.6^{+0.1}_{-0.2}$
$kT_{3}$ (keV) | $0.32^{+0.04}_{-0.02}$
$N_{1}$ $(10^{-4})$ | $10.3^{+1.5}_{-1.2}$ | $12.3^{+1.7}_{-1.4}$ | $9.9^{+1.4}_{-1.2}$ | $10.4^{+1.7}_{-1.4}$
$N_{2}$aafootnotemark:  $(10^{-4})$ | $5.6^{+1.6}_{-1.9}$ | $6.5^{+1.8}_{-2.1}$ | $5.1^{+1.6}_{-1.7}$ | $6.7^{+1.9}_{-3.0}$
$N_{3}$aafootnotemark:  $(10^{-4})$ | $1.4^{+0.4}_{-0.3}$ | $2.2^{+0.5}_{-0.4}$ | $1.2^{+0.4}_{-0.4}$ | $1.1^{+0.5}_{-0.5}$
$F_{0.3-10}$ccfootnotemark:  | $2.53^{+0.04}_{-0.04}$ | $3.03^{+0.04}_{-0.04}$ | $2.38^{+0.05}_{-0.05}$ | $2.62^{+0.06}_{-0.06}$
$\chi_{red}^{2}$ (dof) | 0.99 (1229)

We fit all spectra of X16, Ch17, X18/Nu18 and Ni20, simultaneously, with a
three-temperature VMEKAL model by linking the hydrogen column densities in
different data sets. We apply the pile-up model implemented in the `Xspec`,
and fit all spectra with the model
`pileup*constant*tbabs*(vmekal+vmekal+vmekal)`, for which we turn on the pile-
up model only for the Chandra data. Since the pile-up parameters $\alpha$ and
$f$ are not constrained by the current data, we fix the those parameters to be
the best-fit values $0.8$ and $0.89$, respectively, obtained in section 4.1.3.
The `constant` factors for the X16, Ch17, X18 and Ni20 data are fixed to
unity. We link the abundances of the six elements (O, Ne, Mg, S, Si and Fe)
and fix other elements to the solar abundances, as discussed in section 4.1.1.
A long-term evolution of the spectral properties is discussed with following
three cases.

First, by assuming a constant plasma temperature and flux with time, we link
the three temperatures and normalization factors for all data sets. We find
that the fitting model is unacceptable with a reduced chi-square of
$\chi_{red}^{2}\sim 2.3$ for 1232 dof. Then we unlink the normalization
factors, but keeping the linking for the temperatures. We obtain a reasonable
fit with $\chi_{red}^{2}\sim 1.12$ for 1216 dof. Table 4 and Figure 10
summarize the best-fit parameters and spectra, respectively. We find that the
best-fit parameters are consistent with the results of the fitting for the X16
data only (Model 2 in Table 2). In Table 4, we can see that the 0.3-10 keV
unabsorbed flux deviates from a constant flux with time. One possibility is an
intrinsic temporal evolution of the X-ray emission of AR Scorpii as indicated
by the evolution of the pulse profile discussed in section 3. For example, we
can see in Table 4 that the normalization factors ($N_{1}$, $N_{2}$ and
$N_{3}$) of the Ni20 data are smaller than the values of the X16 data,
indicating the volume of the emission region evolves with time.

We also unlink the abundances for the different data sets to see how the best-
fit parameters evolve. We obtain the best-fit abundances as $F_{e}=0.21\pm
0.04$ for X16, $0.33^{+0.22}_{-0.15}$ for Ch17, $0.22\pm 0.04$ for X18/Nu18
and $0.54^{\pm 0.18}_{-0.13}$ for Ni20, respectively. The abundance of Ni20
data is higher than those of the other data sets. This may be due to the fact
that the spectrum of Ni20 is limited in the 0.3-2.0 keV bands, which does not
cover the 6.8 keV line emission from He-like $Fe$, and hence the abundance of
Fe may be overestimated. We can see that other best-fit parameters are
consistent with the results in Table 4 within the errors.

Second, to investigate another possibility that the change of the X-ray flux
is associated with a change of the temperature of the emitting plasmas, we re-
fit the data by unlinking the lower two temperature components ($kT_{2}$ and
$kT_{3}$). For the highest temperature, since the Ch17 and Ni20 data sets
cannot constrain the value, we keep the linking. An improvement of the fitting
is found with a $F$ static value 5.5, which means that the probability of this
improvement being caused by chance is $\sim 10^{-5}$. Table 4.2 summarizes the
best-fit parameters of this model. The X16, Ch17 and X18/Nu18 data are
consistent with a constant temperature, although the error range for $kT_{2}$
of Ch17 is large. We find however that the lowest temperature $kT_{3}\sim
0.61$ keV of the Ni20 data is higher than the values $kT_{3}\sim 0.3$ keV of
the other data sets, suggesting the plasma temperature shows a temporal
evolution. We find that this tendency of the temperatures does not change,
even though we unlink the abundances for the different data sets. In such a
case, we obtain the abundances of $F_{e}=0.25\pm 0.05$ for X16,
$0.50^{+0.28}_{-0.23}$ for Ch17, $0.23\pm 0.05$ for X18/Nu18, and
$0.45^{+0.18}_{-0.12}$ for Ni20, respectively.

Finally, we can see in Figure 2 that the different observations were carried
out with a different exposure distribution over the orbital phase. Since the
observed X-ray emission shows an orbital modulation with a large amplitude, it
is possible that the different coverage over the orbit produces an apparent
long-term evolution of the flux. To investigate such a possibility, we check
orbital phase-resolved spectra for the X16 and X18 data. As we can see in
Figure 2, the exposure distribution at the orbital phase $\phi_{o}\sim
0.35-0.85$ is relatively uniform for both X16 and X18 data sets. Hence, we
divide the orbital phase into two intervals, $0-0.35/0.85-1$ and $0.35-0.85$,
and obtain a spectrum for each orbital interval.

Table 6 summarizes the best-fit parameters with a three-temperature VMEKAL
model; we link the temperatures of the two data sets and the abundances for 6
elements (O, Ne, Mg, S, Si and Fe). In each observation, we find that the flux
in $\phi_{o}=0.35-0.85$ is higher than that in $\phi_{o}=0-0.35/0.85-1$, and
the flux averaged over the two phase intervals is consistent with the flux
obtained from the whole data of each observation (Table 4). In each orbital
phase interval, the flux of the X16 data is higher than that of the X18 data.
These results suggest that the difference in the exposure distribution over
the orbital phase is not the main reason of the difference between the fluxes
of the X16 and X18 data sets. Although more concrete evidence is required with
more data, the non-pulsed X-ray emission, as well as the pulsed emission, of
AR Scorpii likely experiences a long-term evolution on a timescale of years.

###  4.3. Pulsed spectrum

Figure 11.— Spectrum of the pulsed component. The red, green and black
point/model show the X16, X18 and Ni20 data sets, respectively.

Takata et al. (2018) examine the spectrum of the pulsed component with the X16
data and find evidence of a power-law component with a photon index
$\Gamma=2.3\pm 0.5$. In the new analysis, we include the X18 and Ni20 data
set, but exclude MOS1/MOS2 data of the X16/X18 data because of a small amount
of the photon counts. We define the on-pulse phase and off-pulse phase based
on the pulse profiles folded with the beat frequency, and we generate the
spectrum of the pulsed component by subtracting the off-pulse spectrum from
the on-pulse spectrum. As indicated in Figure 5, we define the on-pulse phase
to cover the main peak for the X16/18 data, and to cover main/secondary peaks
for the Ni20 data. We generate on-pulse and off-pulse spectra using the
`XMMSAS` for the X16/X18 data and using the `Xselect` for the N20 data. We use
the `grppha` task to group the spectral bins so that the count in each bin
after subtracting the off-pulse emission remains $>$40 photons for the X16/X18
data and $>100$ photons for the Ni20 data.

During the fitting, we fix the hydrogen column density at $N_{H}=3.5\times
10^{20}~{}{\rm cm^{-2}}$ inferred from the phase averaged spectrum. The power-
law model provides a reasonable fit with a photon index of
$\Gamma=2.80^{+0.28}_{-0.25}$ ($\chi_{red}^{2}=0.8$) and an unabsorbed flux of
$F_{0.2-2~{}keV}=2.37^{+0.29}_{-0.28}\times 10^{-13}~{}{\rm
erg~{}cm^{-2}s^{-1}}$. Such a large photon index is consistent with a slope
expected from the ratio of the fluxes in optical and X-ray bands. An optically
thin plasma emission model (e.g. two-temperature MEKAL model) also provides a
comparable goodness of the fit with temperatures $kT_{1}=1.8^{+1.5}_{-1.3}$
keV and $kT_{2}=0.11^{+0.03}_{-0.08}$ keV, but we cannot constrain the
abundance. The black-body radiation model predicts an emission size of $\sim
0.1$ km with a source distance $d\sim 0.1$ kpc, which is unreasonably small.

Garnavich et al. (2019) find that the photon index of the optical emission
varies significantly with the orbital phase and with the emission magnitude.
With the current X-ray data quality, however, it is difficult to investigate
such a dependency of the properties of the pulsed emission on the beat phase
or on the orbital phase.

As discussed in section 3, we cannot find the beat signal in the $>2$ keV
band. If the spectrum of the pulsed emission could extend beyond 2 keV band
with a power-law distribution, the expected flux in the 2-10 keV band is
$F_{2-10~{}keV}=3.2^{+1.2}_{-0.8}\times 10^{-14}~{}{\rm erg~{}cm^{-2}s^{-1}}$,
which is $\sim 2$% of the total emission in the 2-10 keV band. Using all data
sets, we have a total of $\sim$55,000 events in the $>$2 keV band. Since the
background emission dominates in the $>$2 keV band of the Ni20 data, the
events related to the pulsed emission will be much less than $\sim$1,100. With
current collected photon counts, therefore, the significance of the pulsation
in the $>2$ keV band will be too small to be detected (i.e., $<5\sigma$), even
if the power-law component extends to that energy band.

##  5\. Sumary and Discussion

We have studied the X-ray emission from AR Scorpii with the archival data
taken in 2016-2020. In the short-term modulation, the periodogram using the
Ni20 data clearly indicates a modulation of the spin frequency of WD,
indicating that part of the X-ray emission is coming from the WD’s
magnetosphere. The short term modulation during the orbital phase
$\phi_{o}\sim 0.1-0.6$ is caused by the mixture between the beat frequency and
the spin frequency, while the intensity peak during $\phi_{o}\sim 0.6-1.1$ is
regularly repeated with the beat frequency. The X16 data also indicates such a
modulation behavior over the orbital phase. Hence, the appearance of the WD’s
spin frequency in the X-ray emission only at the specific orbital phase is
repeated by every orbit. In the long-term variation, we find that the observed
flux level in 2016/2017 was higher than those measured in 2018/2020, and the
amplitude of the orbital waveform of the 2020 NICER data might be larger than
those of other data sets taken in 2016-2018. Moreover, the beat pulse profile
averaged over the orbital phase changed from a single-peak structure in
2016/2018 to a double-peak structure in 2020. These results suggest that the
X-ray emission from AR Scorpii experiences a long-term evolution on a
timescale of years. The spectrum of the pulsed emission can be fit by a power-
law function with a photon index $\Gamma\sim 0.25-0.31$.

A mixture of the beat frequency and the spin frequency in the short-term
modulation is also observed in the optical emission from AR Scorpii (Marsh et
al., 2016; Potter & Buckley, 2018a, b). According to the figure 3 in Potter &
Buckley (2018b) (the periodogram of the optical data), the beat signal is much
stronger than the spin signal in the optical band. In the X-ray bands, the
signal of the WD’s spin frequency during the orbital phase $\phi_{o}\sim
0.1-0.6$ is similar strength to or may be stronger than that of the beat
frequency, while the signal of the spin frequency disappears during the
orbital phase $\phi_{o}\sim 0.6-1,1$, as we discussed in section 3. This shows
that a component of the X-ray emission modulating with the WD spin frequency
appears and disappears in the observation over the orbital phase. We may
interpret the signal of the beat frequency and the spin frequency in the X-ray
emission as a result of the synchrotron emission from two kind of the electron
populations characterized by a larger and a smaller pitch angles,
respectively. A brief discussion of this interpretation is as follows.

We assume a non-thermal emission process for the observed X-ray emission
modulating with the spin and/or beat frequencies. The non-thermal process of
AR Scorpii has been discussed that an interaction between the magnetosphere of
the WD and companion star may form a bow shock or cause a magnetic dissipation
process, which could accelerate the electrons from the companion star (Geng et
al., 2016; Buckley et al., 2017; Takata et al., 2017; Singh et al., 2020).
Garnavich et al. (2019) and Lyutikov et al. (2020), on the other hand, propose
that it is related to the magnetic reconnection event in the magnetosphere of
the WD because of the interaction with the flow from the companion star.
Bednarek (2018) proposes a hadronic model, in which the electrons and protons
are accelerated in a strongly magnetized turbulent region around the M-type
star, and the primary electrons and/or secondary electron/positron pairs
produce a non-thermal emission from AR Scorpii. We assume that the accelerated
electrons at the interaction region are trapped in the WD magnetosphere and,
produce the non-thermal emission through the synchrotron radiation.

The modulation of the beat frequency has been interpreted that the energy
and/or particle injection for the synchrotron radiation is enhanced when the
magnetic pole points toward the companion star (Geng et al., 2016; Potter &
Buckley, 2018b). Takata et al. (2017) suggest on the other hand that the
synchrotron emission around the magnetic mirror points, where the pitch angle
of the electrons becomes $90$ degree, produces a pulsed emission. In their
model, if the WD’s magnetic pole is inclined from the spin axis, the electrons
injected into the WD’s magnetosphere at different spin/beat phases migrate to
the first magnetic mirror points with different travel times. With the effect
of the travel times, the emission from different magnetic mirror points can
enhance the observed intensity at certain spin/beat phase, and results in the
formation of the pulse. This model expects that the pulse profile evolves over
the orbital phase, because the position of the WD relative to the companion
star measured from the Earth depends on the orbital phase.

In the observed X-ray emisison, the modulation with the WD’s spin frequency
appears only at a specific orbital phase interval. Potter & Buckley (2018b)
discuss synchrotron radiation from a region fixed in the WD’s magnetosphere to
explain the optical pulsed emission from AR Scorpii. They demonstrate that by
choosing an appropriate pitch angle of the synchrotron radiation, the pulsed
emission modulating with the spin frequency appears brighter during the
orbital phase $\phi_{o}\sim 0.2-0.45$ only, because of a beaming effect of the
synchrotron radiation. Our speculative interpretation therefore is that the
pulsed X-ray emission of AR Scorpii are produced by two kinds of the electron
populations; electrons having a larger pitch angle and one having a smaller
pitch angle. The emission that modulates with the beat frequency is produced
by the population with a larger pitch angle, while the emission that modulates
with the spin frequency is produced by the population with a smaller pitch
angle.

The electrons that are injected from the interaction region will initially
move toward the WD surface along the magnetic filed line. Because of the
conservation of the first adiabatic invariance, the pitch angle increases as
the electrons approach to the WD surface, because the strength of the magnetic
field increases. If the electrons with a pitch angle $\theta_{0}$ are injected
into the WD magnetosphere at a radius distance $a$ from the WD, the pitch
angle will evolve with the radial distance, $r$, measured from the WD as

$\sin\theta(r)=\left(\frac{a}{r}\right)^{3/2}\sin\theta_{0},$ (1)

where we assume a magnetic dipole field and ignore a change of the Lorentz
factor due to the synchrotron radiation loss. The electrons injected with a
larger pitch angle will lose their energy at or near the magnetic mirror
point. Since the emission from the electrons with a larger pitch angle covers
a wider area in the sky, it can be observable over the orbital phase. We
expect that such an emission modulates with the beat frequency.

For the electrons injected with a smaller pitch angle, on the other hand, they
can move in the deeper WD’s magnetosphere while keeping a small pitch angle.
If the magnetic field is strong enough, the electrons may lose their energy
through the synchrotron radiation before reaching the magnetic mirror points.
A simple quantitative estimation is the following. As an electron moves along
the dipole magnetic field line toward the WD, the synchrotron loss timescale
($\tau_{syn}$) decreases since the pitch angle and magnetic field strength
increase. When the synchrotron loss timescale becomes of the order of the
dynamic timescale, $\tau_{c}=r/c$, the electrons significantly lose the
energy. By equating $\tau_{syn}=\tau_{c}$, and using the equation (1), we
estimate the Lorentz factor of the electron (i) that are injected with a pitch
angle $\theta_{0}$ at $r=a$ and (ii) that emit the typical synchrotron photon
energy $E_{syn}$ as

$\gamma=\frac{2}{3}\left[\frac{16\pi}{27}r_{e}^{9/16}r_{g}^{-1/8}R_{WD}^{3/8}a^{3/16}\left(\frac{E_{syn}}{hc}\right)\right]^{16/23}\sin^{2/23}\theta_{0},$
(2)

where $r_{e}=e^{2}/m_{e}c^{2}$ and $r_{g}=m_{e}c^{2}R_{WD}^{3}/(e\mu_{WD})$
with $\mu_{WD}$ and $R_{WD}$ being the dipole moment and radius of the WD,
respectively. The radial distance, $r_{s}$, at the emission region and the
pitch angle at $r_{s}$ become

$r_{s}=\left[\frac{16\pi}{27}r_{e}^{2}r_{g}^{-3}R_{WD}^{9}a^{9/2}\left(\frac{E_{syn}}{hc}\right)\right]^{2/23}\sin^{6/23}\theta_{0},$
(3)

and

$\sin\theta(r_{s})=\left[\frac{16\pi}{27}r_{e}^{2}r_{g}^{-3}R_{WD}^{9}a^{-7}\left(\frac{E_{syn}}{hc}\right)\right]^{-3/23}\sin^{14/23}\theta_{0},$
(4)

respectively.

We may apply $\mu_{WD}=5\times 10^{34}~{}\rm{G~{}cm^{3}}$ and $R_{WD}=10^{9}$
cm for the WD parameters of AR Scorpii (Buckley et al., 2017), and $a=8\times
10^{8}$ cm of the separation between centers of two stars as the injection
distance (since the companion star, whose radius is $\sim 2\times 10^{10}$ cm,
occupies a large portion of the space between centers of the two stars, this
distance of the injection may be a crude approximation). Then, if the initial
pitch angle is very small, for example $\sin\theta_{0}=0.01$, we obtain
$\gamma\sim 9000$, $\sin\theta(r_{s})\sim 0.06$ and $r_{s}=2\times 10^{10}$cm,
respectively, for the photon energy $E_{syn}=1$ keV. Such an electron
therefore loses the energy while keeping the small pitch angle. The emission
from the small pitch angle approximately directs along the magnetic field
line, and it is observable by the observer whose line of sight is parallel to
the direction of the magnetic field at the radial distance $r_{s}$. We expect
that the X-ray emission from such a small pitch angle is observed only at a
specific $spin~{}phase$ interval (like the pulsar emission) and hence the
observed emission modulates with the spin frequency.

The detectability of the emission from a small pitch angle for AR Scorpii will
depend on the orbital phase because of the beaming effect. We assume that new
high-energy electrons on a WD’s magnetic field line are injected when the
magnetic field line interacts with the companion star magnetosphere/outflow.
Since the dynamical timescale ($\sim$ one second) and the synchrotron
timescale at the emission point are much shorter than the WD’s spin period, we
anticipate that, in terms of the spin phase, the ejected electrons immediately
lose most of their energy after the injection. At a specific orbital phase
interval, therefore, if there is no electron injection when the magnetic field
at the distance $r_{s}$ directs toward the observer, we cannot observe the
X-ray emission modulating with the WD’s spin. This beaming effect could
explain why the signal of the WD’s spin frequency in the observed X-ray
emission appears and disappears over the orbital phase.

Finally, another unique property is that AR Scorpii shows a long-term
evolution of the X-ray emission. In particular, the X-ray light curve folded
with the beat frequency clearly indicates an evolution of the pulsed component
on a time-scale of years. This evolution could be caused either by an
fluctuation or secular evolution of the state of the system. The change of the
pulse profile and orbital modulation have been discovered in the X-ray
emission from an intermediate polar (IP), XY Arietis (Norton & Mukai, 2007),
which is the accreting magnetic white dwarf, and the observed X-ray emission
mainly originates from the white dwarf surface hit by the accretion flow. XY
Arietis experienced an outburst because of a substantial increase of the
accretion rate in 1996 (Hellier et al., 1997), and the increase of the mass
transfer was probably caused by a result of a thermal-viscous disk instability
(Hameury & Lasota, 2017) or an enhancement mass transfer rate from the
companion star. Hellier et al. (1997) report that the spin pulse profile
evolved from a double-peaked profile to a single-peaked profile during the
outburst. Norton & Mukai (2007) find that the spin pulse profile and the
orbital waveform of XY Arietis also show a long-term evolution in the
quiescent state. They interpret the long-term evolution in quiescent state
with a precession model of an accretion disk.

AR Scorpii is a different type of WD’s binary system from the IP in the sense
that there is no accretion disk in the system and the pulsed X-ray emission is
likely produced by the synchrotron radiation of the relativistic electrons in
the WD’s magnetosphere. It has been expected that the modulation peak in the
beat phase occurs when the magnetic poles direct toward the companion star.
The double peak structure of the Ni20 data (and optical data) with a phase
separation $\sim 0.5$ implies the emission from the northern and southern
hemispheres. The single-peak structure observed in the X16 and X18 data may
indicate that the X-ray emission from one hemisphere in 2016/2018 was
intrinsically dim because of an anisotropic energy injection from the
interaction region. If this interpretation is true, the emerging secondary
peak in the Ni20 data might be due to an increase of the energy injection to
the hemisphere.

Another possible interpretation can be that the X-ray emission from one
hemisphere was geometrically invisible in 2016/2018 observations, and it
became visible in 2020. This interpretation could be possible if the spin axis
of the WD is precessing. Schwarzenberg-Czerny (1992) estimate the precession
period would be $\sim 10^{4}-10^{5}$ times the WD’s spin period, which depends
on the oblateness of the star. Katz (2017) applies the precession model to
explain the orbital waveform in the optical bands, and predicts a precession
period of 20-200 years for AR Scorpii. Tovmassian et al. (2007, 2012) report
an evidence of the WD’s precession in the X-ray emission from two IPs. If this
precession scenario is true for AR Scorpii, the transition between a single-
peak and a double-peak structures in the pulse profile folded with the beat
frequency will be regularly repeated with the precession frequency. Since the
pulse peak in the X-ray bands is aligned with the peak in the optical/UV
bands, more long-term multi-wavelength observations will be required to
understand the origin of the long-term evolution of the emission from AR
Scorpii.

We thank to referee for his/her useful comments and suggestions. J.T., W.X.F
and W.H.H. are supported by the National Science Foundation of China (NSFC)
under 11573010, 11661161010, U1631103 and U1838102. L.C.C.L. is supported by
NRFK through the grant number of 2016R1A5A1013277. K.L.L is supported by the
Ministry of Science and Technology of the Republic of China (Taiwan) through
grants 108-2112-M-007-025-MY3 and 109-2636-M-006-017, and he is a Yushan
(Young) Scholar of the Ministry of Education of the Republic of China
(Taiwan). AKHK is supported by the Ministry of Science and Technology of the
Republic of China (Taiwan) through grants 108-2628-M-007-005-RSP and
109-2628-M-007-005-RSP.

## References

  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58
  * Bednarek (2018) Bednarek, W. 2018, MNRAS, 476, L10
  * Buccheri et al. (1983) Buccheri, R., Bennett, K., Bignami, G. F., et al. 1983, A&A, 128, 245
  * Buckley et al. (2017) Buckley, D. A. H., Meintjes, P. J., Potter, S. B., Marsh, T. R., & Gänsicke, B. T. 2017, Nature Astronomy, 1, 0029
  * Davis (2001) Davis, J. E. 2001, ApJ, 562, 575
  * du Plessis et al. (2019) du Plessis, L., Wadiasingh, Z., Venter, C., & Harding, A. K. 2019, ApJ, 887, 44
  * Gaibor et al. (2020) Gaibor, Y., Garnavich, P. M., Littlefield, C., Potter, S. B., & Buckley, D. A. H. 2020, MNRAS, 496, 4849
  * Garnavich et al. (2019) Garnavich, P., Littlefield, C., Kafka, S., et al. 2019, ApJ, 872, 67
  * Geng et al. (2016) Geng, J.-J., Zhang, B., & Huang, Y.-F. 2016, ApJ, 831, L10
  * Hameury & Lasota (2017) Hameury, J. M., & Lasota, J. P. 2017, A&A, 602, A102
  * Hellier et al. (1997) Hellier, C., Mukai, K., & Beardmore, A. P. 1997, MNRAS, 292, 397
  * İnam et al. (2004) İnam, S. Ç., Baykal, A., Matthew Scott, D., Finger, M., & Swank, J. 2004, MNRAS, 349, 173
  * Kaplan et al. (2019) Kaplan, Q., Meintjes, P. J., Singh, K. K., et al. 2019, arXiv e-prints, arXiv:1908.00283
  * Katz (2017) Katz, J. I. 2017, ApJ, 835, 150
  * Littlefield et al. (2017) Littlefield, C., Garnavich, P., Kennedy, M., et al. 2017, ApJ, 845, L7
  * Lyutikov et al. (2020) Lyutikov, M., Barkov, M., Route, M., et al. 2020, arXiv e-prints, arXiv:2004.11474
  * Marcote et al. (2017) Marcote, B., Marsh, T. R., Stanway, E. R., Paragi, Z., & Blanchard, J. M. 2017, A&A, 601, L7
  * Marsh et al. (2016) Marsh, T. R., Gänsicke, B. T., Hümmerich, S., et al. 2016, Nature, 537, 374
  * Mewe et al. (1985) Mewe, R., Gronenschild, E. H. B. M., & van den Oord, G. H. J. 1985, A&AS, 62, 197
  * Mewe et al. (1986) Mewe, R., Lemen, J. R., & van den Oord, G. H. J. 1986, A&AS, 65, 511
  * Norton & Mukai (2007) Norton, A. J., & Mukai, K. 2007, A&A, 472, 225
  * Peterson et al. (2019) Peterson, E., Littlefield, C., & Garnavich, P. 2019, AJ, 158, 131
  * Potter & Buckley (2018a) Potter, S. B., & Buckley, D. A. H. 2018a, MNRAS, 478, L78
  * Potter & Buckley (2018b) —. 2018b, MNRAS, 481, 2384
  * Rea et al. (2007) Rea, N., Nichelli, E., Israel, G. L., et al. 2007, MNRAS, 381, 293
  * Schwarzenberg-Czerny (1992) Schwarzenberg-Czerny, A. 1992, A&A, 260, 268
  * Singh et al. (2020) Singh, K. K., Meintjes, P. J., Kaplan, Q., Ramamonjisoa, F. A., & Sahayanathan, S. 2020, Astroparticle Physics, 123, 102488
  * Stanway et al. (2018) Stanway, E. R., Marsh, T. R., Chote, P., et al. 2018, A&A, 611, A66
  * Stiller et al. (2018) Stiller, R. A., Littlefield, C., Garnavich, P., et al. 2018, AJ, 156, 150
  * Takata & Cheng (2019) Takata, J., & Cheng, K. S. 2019, ApJ, 875, 119
  * Takata et al. (2018) Takata, J., Hu, C. P., Lin, L. C. C., et al. 2018, ApJ, 853, 106
  * Takata et al. (2017) Takata, J., Yang, H., & Cheng, K. S. 2017, ApJ, 851, 143
  * Tovmassian et al. (2012) Tovmassian, G., Zharikov, S., & Neustroev, V. 2012, Mem. Soc. Astron. Italiana, 83, 713
  * Tovmassian et al. (2007) Tovmassian, G. H., Zharikov, S. V., & Neustroev, V. V. 2007, ApJ, 655, 466


#  Deep Q-Learning-based Distribution Network Reconfiguration for Reliability
Improvement

Mukesh Gautam, _Student Member, IEEE_ , Narayan Bhusal, _Student Member, IEEE_
,
and Mohammed Benidris _Senior Member, IEEE_ ,
Department of Electrical & Biomedical Engineering,
University of Nevada, Reno, Reno, NV 89557
(emails: {mukesh.gautam<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>
###### Abstract

Distribution network reconfiguration (DNR) has proved to be an economical and
effective way to improve the reliability of distribution systems. As optimal
network configuration depends on system operating states (e.g., loads at each
node), existing analytical and population-based approaches need to repeat the
entire analysis and computation to find the optimal network configuration with
a change in system operating states. Contrary to this, if properly trained,
deep reinforcement learning (DRL)-based DNR can determine optimal or near-
optimal configuration quickly even with changes in system states. In this
paper, a Deep Q Learning-based framework is proposed for the optimal DNR to
improve reliability of the system. An optimization problem is formulated with
an objective function that minimizes the average curtailed power. Constraints
of the optimization problem are radial topology constraint and all nodes
traversing constraint. The distribution network is modeled as a graph and the
optimal network configuration is determined by searching for an optimal
spanning tree. The optimal spanning tree is the spanning tree with the minimum
value of the average curtailed power. The effectiveness of the proposed
framework is demonstrated through several case studies on 33-node and 69-node
distribution test systems.

###### Index Terms:

Deep Q Network, distribution system reliability, network reconfiguration,
reinforcement learning, spanning tree.

##  I Introduction

The goal of electric utilities is to reliably and economically supply
electricity to their customers through utilization of available resources.
Since power interruptions are mostly due to failure of the distribution system
components [1], enhancing the reliability of distribution systems is
inevitable to provide uninterruptible electric power supply to customers.
Reliability of distribution systems can be enhanced in two ways: (a) optimal
utilization of available resources using smart grid technologies and (b)
installation of redundant resources. The option of installing redundant
resources is not economical and waste of resources. Therefore, sophisticated
smart grid technologies should be developed to optimally utilize distribution
system resources in an optimum manner. In this context, distribution network
reconfiguration (DNR) is one of such smart grid technologies to provide
economic and reliable supply of electricity. DNR can optimize existing
resources by modifying the configuration of distribution networks through
changing status of sectionalizing and tie-switches.

Several analytical and population-based intelligent search approaches have
been proposed in the literature to solve the DNR problem. A genetic algorithm-
based DNR has been proposed in [2] for power quality and reliability
improvement. Similarly, genetic algorithm has been used for DNR to improve the
reliability and optimal placement of distributed generators in [3]. In [4],
mixed-integer quadratic programming has been used for reliability constrained
power loss minimization. In [5], neighbor search algorithm has been used for
DNR to improve the reliability and reduce the loss by taking into account the
uncertainties of data. A review on DNR approaches that improve the reliability
and reduce the power loss is provided in [6]. Analytical and population-based
intelligent search methods used for DNR to improve reliability has the
following shortfalls. Accuracy and effectiveness of analytical-based methods
for DNR depend upon the accuracy of models used, where accurate models impose
scalability challenges. Also, mathematical models are usually derived based on
several approximations and they require complete system information.
Population based methods, on the other hand, are computationally expensive due
to the large search space, especially when system sizes increase.

In addition to population-based intelligent search techniques and analytical
methods, learning-driven approaches are gaining significant attention for
optimal DNR. In [7], batch-constrained reinforcement learning has been used
for the dynamic DNR with the objective of minimizing network operational
costs. In [8], network power loss and number of switching actions in
distribution network are minimized using deep learning. Reinforcement learning
has also been used to simultaneously reduce network power loss and improve
voltage profiles [9], where loop-based encoding is leveraged with Noisy Net
deep Q learning to improve training effectiveness and computational
efficiency. In [10], deep Q learning has been implemented to minimize
switching actions while performing the DNR. The work presented in [10] has
been tested for the computational cost and the scalability as compared to
brute-force search algorithm and the genetic algorithm. Although there are
several similarities between DNR for different objectives, optimum DNR for
reliability improvement is a challenging task since it requires estimating
energy not supplied for each possible configuration. Therefore, developing
intelligent learning-based approaches for DNR to improve the reliability is
pivotal.

This paper proposes a deep reinforcement learning (DRL)-based framework for
DNR to improve the reliability of distribution systems. In the proposed
optimization framework, the average curtailed power is used as an objective
function. In addition to operation and technical constraints, all-node-
traversing and radiality constraints are considered. In the training phase of
the proposed framework, Q values are predicted using forward propagation of a
deep neural network (DNN). Actions are selected using Epsilon-Greedy
algorithm. When actions are passed through the training environment, the DRL
agent gets rewarded (or penalized) based on its performance. Target Q values
are calculated based on the reward. The mean squared error (MSE), which is the
most commonly employed loss function for regression, is computed using the
predicted and target Q values. Errors are then back-propagated to update the
weights of DNN. The trained DRL agent is then used to find the best network
configuration. The proposed framework is validated through case studies on
several distribution test systems, and the results show that the proposed
framework can effectively find a network configuration with high reliability
level.

The main contributions of this paper are summarized as follows:

  * •

The proposed work represents an advancement over existing DNR approaches that
either use approximate mathematical models or heuristic techniques. In
contrary to existing analytical and population-based methods, learning-driven
methods provide a flexibility for optimal DNR for changing system states
(change of loads, addition of distributed energy sources, etc.) since they do
not need to repeat the entire calculation processes when there is a change in
system states. If trained properly, DRL can find the optimal or near-optimal
network configuration in almost real-time.

  * •

This papers provides a framework for optimal network reconfiguration with an
objective of reliability enhancement. Formulating an objective function for
reliability enhancement through DNR has been a challenging task because the
the reliability level needs to be evaluated for each system configuration.

The rest of the paper is arranged as follows. Section II explains problem
formulation for proposed DNR. Section III describes the proposed framework and
solution approach. Section IV validates the proposed work through case
studies. Section V provides concluding remarks.

##  II Problem Formulation

This section explains the formulation of the objective function and the
constraints of the DNR problem.

###  II-A Objective Function

Reliability is one of the major factors that indicates the performance of the
system. Reliability of distribution systems can be quantified using several
reliability indices including system average interruption frequency index
(SAIFI), system average interruption duration index (SAIDI), customer average
interruption duration index (CAIDI), and average curtailed power. The average
curtailed power is taken as the objective function for the problem under
consideration since it can capture the severity of the outages and is directly
affected by the topology or configuration of a distribution network.
Therefore, in this work, the objective of the DNR is to minimize the total
annual energy not supplied by the system. Mathematically, the average
curtailed power can be expressed as follows.

$\mbox{Average Curtailed Power}=\sum_{k\in\Omega_{k}}P_{d,k}U_{k}\mbox{,}$ (1)

where $P_{d,k}$ is the power demand at node $k$; $\Omega_{k}$ is the set of
nodes with power demand; and $U_{k}$ is the average annual power
unavailability duration at node $k$, which can be defined as follows.

$U_{k}=\sum_{l\in\Omega_{lk}}\lambda_{l}r_{l}\mbox{,}$ (2)

where $\lambda_{l}$ is the failure rate of branch (or edge) $l$; $r_{l}$ is
the outage duration (or repair time) of branch (or edge) $l$; and
$\Omega_{lk}$ is the set of branches (or edges) between substation node and
node $k$.

###  II-B Constraints

During the process of selection of the optimal configuration of the
distribution network, two constraints (radiality constraint and all-node-
traversing constraint) should always be satisfied.

  * •

Radiality constraint: Radiality constraint is always maintained in a
distribution system in order to design the protection coordination schemes.
Each candidate configuration should have a radial topology since most of the
practical distribution systems do not have loop structure.

  * •

All-node-traversing constraint: A distribution system operator should always
configure the network in such a way that all loads are supplied with power in
non-contingent scenarios. Therefore, for each candidate configuration, all-
node-traversing constraint should always be satisfied.

If we search for a spanning tree, both of the aforementioned constraints are
satisfied. In order to search for the spanning tree, we have to represent the
distribution network as an undirected graph
$\mathscr{G}=(\mathscr{N},\mathscr{E})$, where $\mathscr{N}$ is a set of nodes
(or vertices) and $\mathscr{E}$ is a set of edges (or branches). For the
graph, a node-branch incidence matrix can be constructed after satisfying all-
node-traversing constraint. If the node-branch incidence matrix is full
ranked, then the radiality constraint is satisfied.

##  III Proposed Framework

This work leverages recently advanced reinforcement learning techniques for
DNR to improve reliability of distribution systems. This section provides a
brief overview of Deep Q learning, reward function, and training attributes of
the Deep Q learning.

###  III-A Deep Q Learning

A reinforcement learning (RL) is a branch of machine learning in which an
_agent_ learns to take suitable _actions_ to maximize cumulative _reward_ it
gets from an uncertain _environment_. In general, an RL system consists of
four main integrants: policy, reward, value functions, and environment model.
An agent decides the action to be taken based on the policy. The policy maps
states to actions. When the agent takes an action, it gets rewarded (or
penalized). Value function calculates the expected value of cumulative reward
that a agent gets when it follows a certain policy. There are different
algorithms for RL. The choice of an algorithm depends on many factors such as
continuous/discrete nature of states, continuous/discrete action-space, etc.
For the DNR problem under consideration, the action-space is discrete in
nature, which makes Q-Learning a suitable candidate for the problem. However,
a basic Q-Learning needs large sized look-up tables where state-action values
are stored. To avoid the use of large sized look-up tables, a deep neural
network (DNN) is used as an action-value function approximator. The addition
of DNN in the basic Q-Learning make the framework a Deep Q Network (DQN). The
update rule for action-value function in Q-learning is defined as follows
[11].

$\displaystyle Q(S_{t},A_{t})\leftarrow Q(S_{t},A_{t})+\alpha\times[R_{t+1}$
(3)
$\displaystyle+\gamma\times\max\limits_{a}Q(S_{t+1},A_{t+1})-Q(S_{t},A_{t})]$

where $A_{t}$ and $S_{t}$ are the action and state of an agent at
$t^{\text{th}}$ iteration; $Q(S_{t},A_{t})$ is the action-value function at
$t^{\text{th}}$ iteration; $Q(S_{t+1},A_{t+1})$ is the action-value function
at $(t+1)^{\text{th}}$ iteration; $\alpha$ is the learning rate; and $\gamma$
is the reward discount factor.

Instead of iteratively updating the action-value function, the DNN is trained
and the parameters of the action-value function are optimized to minimize the
mean-squared error (MSE) loss function (i.e., regression loss function), which
is expressed as follows [12].

$L(\theta)=\mathbb{E}[(Q(S_{t},A_{t}|\theta)-y_{t})^{2}]$ (4)

where $\mathbb{E}$ denotes expectation operator; $\theta$ denotes the
parameter of action-value function $Q(S_{t},A_{t})$; and $y_{t}$ denotes the
target action-value function, which is defined a follows.

$y_{t}=R(S_{t},A_{t})-\gamma\times\max\limits_{a}Q(S_{t+1},A_{t+1}|\theta^{\prime})$
(5)

In (5), $R(S_{t},A_{t})$ denotes reward function at $t^{\text{th}}$ iteration
and $\theta^{\prime}$ denotes the parameter of action-value function
$Q(S_{t+1},A_{t+1})$.

For the DNR problem under consideration, the status of sectionalizing and tie
switches denotes the state; and the indices of opened branches (or edges)
denote the actions.

###  III-B Reward Function

A reward function based on average curtailed power is designed to evaluate the
actions taken by the agent. For each action (here, set of edges to be opened),
firstly, all-node-traversing constraint is checked. If this constraint is not
satisfied, a high negative reward (or penalty) is given to the agent. After
this, a node-branch incidence matrix is computed for the network
configuration. If the node-branch incidence matrix is full-ranked, this means
that the network is radial. Therefore, if radiality constraint is not
satisfied, a high negative reward (or penalty) is given to the agent. After
both the constraints are satisfied, the average curtailed power of the network
configuration is calculated using annual power unavailability duration
($U_{k}$) and power demand ($P_{d,k}$) at each node $k$.

###  III-C Training Attributes

The training of DQN is performed for a certain number of episodes ($n_{ep}$).
The initial state of the system is the state with tie-switches open. The
weights of DNN are initialized with some random values. In each episode, the
predicted Q values corresponding to each edge of the system is computed based
on forward propagation of DNN. For the selection of actions, the Epsilon-
Greedy (exploration-exploitation) algorithm [13] is used. The value of
exploration rate i.e., epsilon ($\varepsilon$) is initialized at 1. The
epsilon is updated after each episode as follows.

$\varepsilon_{new}=\varepsilon_{old}-\frac{\varepsilon_{old}-\varepsilon_{min}}{n_{ep}}$
(6)

where $\varepsilon_{min}$ is the minimum exploration rate. The target Q value
of the DQN is computed using (5). The MSE losses for each episode are computed
based on (4) using the actual and target Q-values. These MSE losses are back-
propagated to update the weights of the DNN.

Figure 1: Flowchart showing the training phase of the proposed framework

The flowchart shown in Fig. 1 describes the training phase of the proposed
framework.

##  IV Case studies and Discussions

The proposed framework is implemented on the $33$-node and $69$-node
distribution systems. The case study parameters for both systems are as
follows. The failure rate of each branch is assumed proportional to its branch
impedance. The highest failure rate (in this case 0.4 failures/year or f/yr)
is assigned to the branch with the largest impedance; the lowest failure rate
(in this case 0.1 f/yr) is assigned to the branch with the smallest impedance;
and linear interpolation is used to determine the failure rates of remaining
branches. Regarding the outage duration (or repair rate) of each branch, its
value is assumed to be constant (6 hr is used for all case studies). Normally-
open switches (or tie switches) are assumed to be fully reliable. The annual
power unavailability duration ($U_{k}$) of each node is determined based on
the values of failure and repair rates for each branch.

###  IV-A Case I: 33-node Distribution System

The $33$-node distribution test system is $100$ kVA, $12.66$ kV radial
distribution system with $33$ nodes, $32$ branches and $5$ tie-lines.
Therefore, the total number of branches in this system is $37$. All branches
(including tie-lines) are numbered from $1$ to $37$ following our previous
work [14]. The total load of the system is 3.71 MW. The detailed data of the
system is provided in [15].

Figure 2: Rewards of training episodes for the 33-node system Figure 3: MSE of
losses of training episodes for the 33-node system Figure 4: Average curtailed
power of training episodes for the 33-node system

The training of the DQN for the 33-node system is performed for 10,000
episodes. The initial state of the system is the state with all tie-switches
(i.e., edges 33, 34, 35, 36, and 37) are open. Initially, the rewards are very
low but then they increase as the number of episodes increases. Fig. 2 shows
the actual rewards and running mean (100-episode window) of actual rewards as
the episode progresses. It can be seen from Fig. 2 that as the number of
episodes increases, the running mean of the reward increases and saturates
after nearly 8,000 episodes.

Fig. 3 shows the actual values and running mean (100-episode window) of MSE
losses as the episode progresses. Fig. 3 also shows that initially MSE loss is
very high but its running mean remains low after a few episodes. Fig. 4 shows
the plot of running mean (100-episode window) of the average curtailed power
as the episode progresses. From Fig. 4 we can see that the average curtailed
power decreases as the number of episodes increases and settles down after
nearly 8,000 episodes. The final value of the average curtailed power obtained
during the training of the proposed framework is 23.96 MWh/yr. For this value
of the average curtailed power, the open edges (branches) are 7, 14, 26, 33,
and 34.

###  IV-B Case II: 69-node Distribution System

The 69-node distribution test system is a $12.66$ kV radial distribution
system with $69$ nodes, $68$ branches and $5$ tie-lines. Therefore, the total
number of branches in this system is $73$. All branches (including tie-lines)
are numbered from $1$ to $73$ following our previous work [14]. The total load
of the system is 3.80 MW. Detailed data of the system are given i [16].

Figure 5: Rewards of training episodes for the 69-node system Figure 6: MSE of
losses of training episodes for the 69-node system Figure 7: Average curtailed
power of training episodes for the 69-node system

The training of the DQN for the 69-node system is performed for 10,000
episodes. The initial state of the system is the state with all tie-switches
(i.e., edges 69, 70, 71, 72, and 73) are open. Initially, rewards are very low
but then they increase as the number of episodes increases. Fig. 5 shows the
actual rewards and running mean (100-episode window) of actual rewards as the
episode progresses. We can see from the Fig. 5 that as the number of episode
increases, the running mean of the rewards increases and saturates after
nearly 8,000 episodes. Fig. 6 shows the actual values and running mean (100
episode window) of MSE losses as the episode progresses. Fig. 6 also shows
that initially MSE loss is very high but its running mean remains low after a
few episodes. Fig. 7 shows the plot of running mean (100 episode window) of
the average curtailed power as the episode progresses. From Fig. 7 we can see
that the the average curtailed power decreases as the number of episodes
increases and settles down after nearly 8,000 episodes. The final value of the
average curtailed power obtained during the training of the proposed framework
is 28.48 MWh/yr. For this value of the average curtailed power, the open edges
(branches) are 14, 18, 21, 58, and 69.

##  V Conclusion

This paper has proposed a DRL-based framework to determine the configuration
of a distribution network with the least value of the average curtailed power.
During the training of the proposed framework, the weights of DNN were
initialized with random values and the system state was initialized with the
configuration with tie-switches open. The target Q values were computed based
on reward the DRL agent gets from the environment. The predicted and target Q
values were used to update the weights of DNN. Case studies were performed on
33-node and 69-node distribution test systems. The results exhibit the
effectiveness of the proposed framework to improve the reliability level of
distribution systems.

## Acknowledgement

This work was supported by the U.S. National Science Foundation (NSF) under
Grant NSF 1847578.

## References

  * [1] A. A. Chowdhury and D. Koval, _Power Distribution System Reliability: Practical Methods and Applications_. J. Wiley Inc., 2009.
  * [2] N. Gupta, A. Swarnkar, and K. Niazi, “Distribution network reconfiguration for power quality and reliability improvement using genetic algorithms,” _International Journal of Electrical Power & Energy Systems_, vol. 54, pp. 664–671, 2014.
  * [3] Y. Tian, M. Benidris, S. Sulaeman, S. Elsaiah, and J. Mitra, “Optimal feeder reconfiguration and distributed generation placement for reliability improvement,” in _2016 International Conference on Probabilistic Methods Applied to Power Systems (PMAPS)_ , 2016, pp. 1–7.
  * [4] J. Jose and A. Kowli, “Reliability constrained distribution feeder reconfiguration for power loss minimization,” in _2016 National Power Systems Conference (NPSC)_ , 2016, pp. 1–6.
  * [5] P. Zhang, W. Li, and S. Wang, “Reliability-oriented distribution network reconfiguration considering uncertainties of data by interval analysis,” _International Journal of Electrical Power & Energy Systems_, vol. 34, no. 1, pp. 138–144, 2012.
  * [6] B. Sultana, M. Mustafa, U. Sultana, and A. R. Bhatti, “Review on reliability improvement and power loss reduction in distribution system via network reconfiguration,” _Renewable and sustainable energy reviews_ , vol. 66, pp. 297–310, 2016.
  * [7] Y. Gao, W. Wang, J. Shi, and N. Yu, “Batch-constrained reinforcement learning for dynamic distribution network reconfiguration,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 6, pp. 5357–5369, 2020.
  * [8] X. Ji, Z. Yin, Y. Zhang, B. Xu _et al._ , “Real-time autonomous dynamic reconfiguration based on deep learning algorithm for distribution network,” _Electric Power Systems Research_ , vol. 195, p. 107132, 2021.
  * [9] B. Wang, H. Zhu, H. Xu, Y. Bao, and H. Di, “Distribution network reconfiguration based on noisynet deep q-learning network,” _IEEE Access_ , vol. 9, pp. 90 358–90 365, 2021.
  * [10] S. H. Oh, Y. T. Yoon, and S. W. Kim, “Online reconfiguration scheme of self-sufficient distribution network based on a reinforcement learning approach,” _Applied Energy_ , vol. 280, p. 115900, 2020.
  * [11] A. Zai and B. Brown, _Deep reinforcement learning in action_. Manning Publications, 2020.
  * [12] D. Cao, W. Hu, J. Zhao, G. Zhang, B. Zhang, Z. Liu, Z. Chen, and F. Blaabjerg, “Reinforcement learning and its applications in modern power and energy systems: A review,” _Journal of Modern Power Systems and Clean Energy_ , vol. 8, no. 6, pp. 1029–1042, 2020.
  * [13] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [14] M. Gautam, N. Bhusal, M. Benidris, and S. J. Louis, “A spanning tree-based genetic algorithm for distribution network reconfiguration,” in _2020 IEEE Industry Applications Society Annual Meeting_. IEEE, 2020, pp. 1–6.
  * [15] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” _IEEE Power Engineering Review_ , vol. 9, no. 4, pp. 101–102, 1989.
  * [16] J. Savier and D. Das, “Impact of network reconfiguration on loss allocation of radial distribution systems,” _IEEE Transactions on Power Delivery_ , vol. 22, no. 4, pp. 2473–2480, 2007.