# NGC 7793 P9: An Ultraluminous X-Ray Source Evolved from a Canonical Black
Hole X-Ray Binary

Chin-Ping Hu  _Department of Physics, The University of Hong Kong, Pokfulam
Road, Hong<EMAIL_ADDRESS>Albert K. H. Kong  _Institute of Astronomy,
National Tsing Hua University, Hsinchu 30013, Taiwan_ C.-Y. Ng  _Department of
Physics, The University of Hong Kong, Pokfulam Road, Hong<EMAIL_ADDRESS>K. L. Li  _Department of Physics and Astronomy, Michigan State University,
East Lansing, MI 48824, USA_

(Received 2018 March 31; Revised 2018 July 19; Accepted 2018 July 22)

###### Abstract

Transient ultraluminous X-ray sources (ULXs) provide an important link
bridging transient low-mass X-ray binaries and ULXs. Here we report the first
discovery of both a canonical sub-Eddington outburst and an ultraluminous
super-Eddington outburst from an unusual transient ULX, NGC 7793 P9 with a
variability factor higher than $10^{3}$. Its X-ray spectrum switches between
the typical high/soft state and the steep power-law state during the canonical
outburst. The inner radius of the accretion disk and the disk
temperature–luminosity correlation suggest that P9 harbors a stellar-mass
black hole (BH). At the beginning of the ultraluminous outburst, we observe a
cool outer disk with a hard Comptonized spectrum, implying a transition to the
ULX regime. When the luminosity increases to $L\gtrsim 3\times 10^{39}$ erg
s-1, P9 shows a significantly curved spectrum that can be described by either
a slim disk or a strongly curved Comptonized corona. The phenomenological
model suggests that the hot disk observed near the peak of the ultraluminous
outburst is coincidentally consistent with the extension of the thermal track.
Utilizing more physical Comptonized disk models, we suggest that the corona
cools down and the apparent disk-like spectrum is dominated by soft
Comptonization. The significant variability above 1 keV supports this two-
component scenario. The spectral evolution can also be interpreted with the
supercritical accretion model. All these indicate that a canonical black hole
X-ray binary can show properties of a ULX.

accretion, accretion disks – black hole physics – X-rays: binaries – X-rays:
individual (NGC 7793 P9)

††facilities: _CXO_ (ACIS), _XMM_ (MOS, PN), _NuSTAR_ , _Swift_ (XRT, UVOT)

##  1 Introduction

Accreting black holes (BHs) are bright X-ray sources powered by the
gravitational potential energy of accreting materials. The luminosities of
Galactic BHs are strongly limited to the Eddington luminosity
($L_{\rm{Edd}}$), at which the radiation pressure balances the inward
gravitational force of the accreting materials. Thanks to the revolution in
astronomical instrumentation over the past several decades, the discovery of
ultraluminous X-ray sources (ULXs) challenged the Eddington limit (see reviews
by Frank et al., 2002; Feng & Soria, 2011; Kaaret et al., 2017). ULXs are
extragalactic off-nuclear X-ray point sources with luminosities exceeding the
Eddington limit of a stellar-mass BH. Most ULXs are believed to be powered by
stellar-mass BHs with mild beaming and/or moderate super-Eddington accretion
rates. They show ultraluminous states distinct from the canonical spectral
states of BH X-ray binaries (BHXBs, Feng & Soria, 2011; Gladstone et al.,
2009; Kaaret et al., 2017). Unlike the spectral behavior of a BHXB in the
steep power-law (SPL) state that shows a hot disk and an optically thin
corona, a genuine ULX usually shows a much cooler disk with a larger inner
radius plus an optically thick corona and has a luminosity of $L\gtrsim
3\times 10^{39}$ erg s-1. ULX states can be further divided into the hard
ultraluminous (HUL), soft ultraluminous (SUL), and supersoft ultraluminous
spectral states (Sutton et al., 2013; Kaaret et al., 2017). In addition, those
low-luminosity ULXs with $L\lesssim 3\times 10^{39}$ erg s-1usually show a
curved spectrum that can be described with a broadened disk (BD) model,
implying that the disk component is no longer a geometrically thin accretion
disk but the luminosity remains sub-Eddington (Sutton et al., 2013).

The best studied ULXs are persistent sources with low variability factors
(Feng & Kaaret, 2006; Kaaret & Feng, 2009). A stable high mass accretion rate
is necessary to interpret the long lifetime of the ultraluminous state. As
long as the companion star of a ULX is a massive star with $M\gtrsim
10M_{\odot}$, super-Eddington mass transfer occurs when the companion evolves
to the Hertzsprung gap without a common envelope (King et al., 2001; Rappaport
et al., 2005; Pavlovskii et al., 2017). Most of the currently known ULXs,
including the ultraluminous pulsar P13 in NGC 7793, have early-type companions
with high optical luminosities (Bachetti et al., 2014; Motch et al., 2014; Lau
et al., 2016). It has also been proposed that some ULXs may exhibit a low-mass
companion and be powered by a long-lasting super-Eddington outburst (King,
2002). In contrast, BHXBs with low-mass companions are usually transient X-ray
sources. Therefore, transient ULXs is key to testing this scenario and showing
that a low-mass X-ray binary (LMXB) can evolve to a ULX with typical UL
spectral behaviors. The Galactic LMXBs GRS 1915+105 and XTE J1550$-$564 showed
optically thick Comptonized spectra during the peak of the outburst with a
luminosity of $\sim 10^{39}$ erg s-1. They could fill in the link bridging
BHXBs and ULXs (Kubota & Done, 2004; Soria, 2007; Vierdayanti et al., 2010).
Several transient ULXs were reported as evolved from canonical BHXBs although
they showed either canonical BHXB spectral states or ultraluminous outburst
only (see, e.g., Middleton et al., 2012; Kaaret & Feng, 2013). Therefore, it
remains unclear whether a canonical BHXB showing sub-Eddington outbursts can
evolve into a ULX with a much higher luminosity and typical UL spectral
behaviors. In this paper, we report the discovery of a transient ULX, NGC 7793
P9 (hereafter P9), that showed both a canonical BHXB outburst and an
ultraluminous outburst that fill in the gap.

NGC 7793 is an SA(s)d type galaxy in the Sculptor Group with a distance of
3.6–3.9 Mpc (Karachentsev et al., 2003; Crowther, 2013; Radburn-Smith et al.,
2011; Tully et al., 2016). It has a bulge with filamentary spiral structures
and several Hii regions (Read & Pietsch, 1999). P9 is located at the northeast
edge of one filament extending from the host galaxy and has no obvious
association with any Hii regions or supernova remnants (Read & Pietsch, 1999).
P9 was first characterized as a highly variable point source in a _ROSAT_
survey (Read & Pietsch, 1999). The X-ray spectrum of P9 was well fit with
either a PL or a thermal bremsstrahlung model. The absorption column density
($N_{\textrm{H}}$) was rather large. A deep _Chandra_ observation (ObsID
http://cda.harvard.edu/chaser/viewerContents.do?obsid=3954 (catalog 3954))
revealed a faint source, CXOU J235808.7$-$323403, at the position of P9, but
its spectrum was not well characterized (Pannuti et al., 2011). Therefore, the
luminosity was not well constrained. By simply assuming a PL with a photon
index $\Gamma=1.5$ and $N_{\textrm{H}}=1.15\times 10^{20}$ cm-2, the
luminosity of P9 was estimated to be $4.5\times 10^{37}$ erg s-1. No optical
or radio counterparts have been found (Read & Pietsch, 1999; Pannuti et al.,
2011).

We report a detailed analysis of the transient ULX P9 based on the _Neil
Gehrels Swift_ (hereafter _Swift_) monitoring data and several dedicated
observations made with _Chandra_ , _XMM-Newton_ , and _NuSTAR_. We introduce
the observation and data analysis in Section 2. The detailed spectral
analysis, including the variability of $\Gamma$, the stacked _Swift_
spectroscopy, and tests of more physical models during the ultraluminous
outburst are presented in Section 3. We further discuss the implications of
the spectral variability, especially the evolution of the disk component and
the accretion models, in the super-Eddington accretion regime in Section 4.
Finally, we summarize this work in Section 5.

##  2 Observations and Data Reduction

Figure 1: The _Swift_ XRT light curve of NGC 7793 P9. The black diamonds are
the data points determined with _Swift_ XRT, where the gray crosses denote the
upper limits. The _Chandra_ (blue squares) and _XMM-Newton_ (red filled
circles) data are also plotted on the same figure. Table 1: Observation Log of
NGC 7793 P9.

Epoch | Observatory | ObsID | Instrument | MJD | Exposurea (ks)
---|---|---|---|---|---
C1 | _Chandra_ | http://cda.harvard.edu/chaser/viewerContents.do?obsid=3954 (catalog 3954) | ACIS-S | 52889 | 49
C2 | _Chandra_ | http://cda.harvard.edu/chaser/viewerContents.do?obsid=14231 (catalog 14231) | ACIS-S | 55787 | 59
C3 | _Chandra_ | http://cda.harvard.edu/chaser/viewerContents.do?obsid=13439 (catalog 13439) | ACIS-S | 55920 | 58
C4 | _Chandra_ | http://cda.harvard.edu/chaser/viewerContents.do?obsid=14378 (catalog 14378) | ACIS-S | 55926 | 25
X1 | _XMM-Newton_ | http://nxsa.esac.esa.int/nxsa-web/#obsid=0693760101 (catalog 0693760101) | PN | 56061 | 23
|  |  | MOS |  | 36
X2 | _XMM-Newton_ | http://nxsa.esac.esa.int/nxsa-web/#obsid=0693760401 (catalog 0693760401) | PN | 56621 | 46
|  |  | MOS |  | 48
X3 | _XMM-Newton_ | http://nxsa.esac.esa.int/nxsa-web/#obsid=0748390901 (catalog 0748390901) | PN | 57002 | 47
|  |  | MOS |  | 49
X4b | _XMM-Newton_ | http://nxsa.esac.esa.int/nxsa-web/#obsid=0781800101 (catalog 0781800101) | PN | 57528 | 24
|  |  | MOS |  | 48
| _NuSTAR_ | 80201010002 | FPMA+FPMB | 57528 | 114

aafootnotetext: Effective exposure time after filtering out the time intervals
with flaring background.

bbfootnotetext: A joint observation with _NuSTAR_.

###  2.1 _Swift_ Observations

The host galaxy NGC 7793 has been monitored with _Swift_ since 2010. Four
series of observations have been carried out in late 2010, from mid 2011 to
late 2013, in late 2014, and from 2016 April to 2017 June. We extracted the
XRT light curve with the online XRT
pipeline111http://www.swift.ac.uk/user_objects/ (Evans et al., 2007, 2009). P9
is clearly detected in 88 of the total 122 observations. We converted the
0.3–10 keV count rates to the luminosities with
WebPIMMS222https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl by
assuming $N_{\rm{H}}=10^{21}$ cm-2, a PL spectrum with $\Gamma=2$, and a
distance of 3.9 Mpc. We will show that $\Gamma$ varied between $\sim$1.5 and
$\sim$2.5 in a later section. The value $\Gamma=2$ is a conventional average
value to estimate the luminosity variability in Figure 1. We searched for
periodic signals between $P=10$–$1000$ days to check if P9 contained orbital
or superorbital modulations like other ultraluminous pulsars but we did not
find significant periodicity (Kong et al., 2016; Walton et al., 2016; Hu et
al., 2017). The luminosity of P9 varied between $\sim 10^{38}$ and $\sim
10^{39}$ erg s-1, well below the ULX regime before 2013 (see Figure 1).
Interestingly, P9 showed an ultraluminous outburst with a peak luminosity of
$\sim 5\times 10^{39}$erg s-1 after 2014, and the luminosity returned to
$\lesssim 10^{39}$ erg s-1 in 2017. We further downloaded all of the _Swift_
UVOT images taken between 2010 and 2017, and stacked them to obtain deep UV
images in four _Swift_ bands (_u_ , _w1_ , _m2_ , and _w2_). No UV counterpart
was seen within the error circle of the _Chandra_ position.

###  2.2 _Chandra_ Data

NGC 7793 was observed with _Chandra_ in 2003 (C1: ObsID
http://cda.harvard.edu/chaser/viewerContents.do?obsid=3954 (catalog 3954)) and
2011 (C2: ObsID http://cda.harvard.edu/chaser/viewerContents.do?obsid=14231
(catalog 14231), C3: ObsID
http://cda.harvard.edu/chaser/viewerContents.do?obsid=13439 (catalog 13439),
C4: ObsID http://cda.harvard.edu/chaser/viewerContents.do?obsid=14378 (catalog
14378)). The observation log of P9 is listed in Table 1. These four
observations were made with the Advanced CCD Imaging Spectrometer array
(ACIS-S) in the timed exposure mode. We reprocessed the data using the
pipeline chandra_repro in the _Chandra_ Interactive Analysis of Observations
(CIAO) version 4.9 with the most recent calibration database (CALDB) version
4.7.3 (Fruscione et al., 2006). We extracted the source events from a circular
aperture of 2″ radius. The background events were extracted from a nearby
source-free region. We estimated the source luminosity with the srcflux
command and plotted the results in Figure 1. The luminosities of the 2011
observations are $\sim 5\times 10^{38}$ erg s-1, consistent with the _Swift_
light curve.

We extracted the light curves of P9 to check for any variability in each
observation. We used the task dmextract to create the background-subtracted
light curves, and the results are plotted in Figure 2. The light curves show
limited fluctuations, and no significant trends were observed.

###  2.3 _XMM-Newton_ Data

_XMM-Newton_ has observed NGC 7793 nine times with the EPIC CCD in the full-
frame mode, and four of them are publicly available. We reduced the data using
the _XMM-Newton_ Science Analysis Software (SAS version 16.0.0) and used
epproc and emproc to reprocess the PN and MOS data with the latest calibration
files, respectively. We filtered out the time intervals with a flaring
particle background. The effective exposure times are listed in Table 1. The
2012 observation (X1: ObsID http://nxsa.esac.esa.int/nxsa-
web/#obsid=0693760101 (catalog 0693760101)) shows that the luminosity of P9
was $\sim 2\times 10^{38}$ erg s-1. In 2013, P9 was not detected with _XMM-
Newton_ (X2: ObsID http://nxsa.esac.esa.int/nxsa-web/#obsid=0693760401
(catalog 0693760401)). We estimated a 3$\sigma$ upper limit as $\sim 8\times
10^{36}$ erg s-1. On the other hand, P9 was very bright with luminosities
exceeding $10^{39}$ erg s-1 in 2014 (X3: ObsID http://nxsa.esac.esa.int/nxsa-
web/#obsid=0748390901 (catalog 0748390901)) and 2016 (X4: ObsID
http://nxsa.esac.esa.int/nxsa-web/#obsid=0781800101 (catalog 0781800101)). The
nondetection of P9 in 2013 suggests that P9 was in quiescence, and the entire
time span can be divided into two epochs: a less luminous outburst from 2011
to 2012 and an ultraluminous outburst from 2014.

We then extracted the background-subtracted light curves of P9 from individual
observations. We extracted the source events from a circular aperture of 15″
radius to minimize the possible contamination from CXOU J235810.4$-$323357
(Pannuti et al., 2011). We used epiclccorr to subtract the background light
curve and apply the relative and absolute corrections. Figure 2 shows the
light curves of the _XMM-Newton_ observations. No significant variability was
seen in the X1 and X3 data sets. However, the light curve of X4 shows a
decreasing trend. We calculated the fractional variability ($F_{\rm{var}}$)
for X3 and X4 in which P9 was bright enough. $F_{\rm{var}}$ can be defined as

$F_{\rm{var}}=\frac{1}{\left<X\right>}\sqrt{S^{2}-\left<\sigma_{err}^{2}\right>}\rm{,}$
(1)

where $S^{2}$ is the total variance of the light curve,
$\left<\sigma_{err}^{2}\right>$ is the mean square of the uncertainty of each
data point, and $\left<X\right>$ is the mean count rate (Edelson et al.,
2002). For individual data sets, we calculated $F_{\rm{var}}$ for the full
energy range (0.3–10 keV), the soft X-ray band (0.3–1 keV), and the hard X-ray
band (1–10 keV). We did not observe significant variability in all the three
bands of X3, and the 2$\sigma$ upper limit is estimated as
$F_{\rm{var}}<0.06$. On the other hand, X4 shows a clear variability with
$F_{\rm{var}}=0.12\pm 0.01$ in 0.3–10 keV, and the variability is dominated by
the X-ray emission in 1–10 keV with $F_{\rm{var}}=0.15\pm 0.02$.
$F_{\rm{var}}$ in the 0.3–1 keV band was less obvious with a value of $0.05\pm
0.02$. We corrected the photon arrival time to the solar barycenter and
searched for periodicity in the range of 0.1–1000 s with both the $H$-test (de
Jager et al., 1989) and the Fourier analysis, but no periodicity was found.

###  2.4 _NuSTAR_ Data

NGC 7793 was observed with both _NuSTAR_ and _XMM-Newton_ on 2016 May 20
(ObsID 80201010002). We reduced the data using the _NuSTAR_ analysis software
NUSTARDAS (version 1.9.1) with the calibration database of version 20170614.
We used the script nupipeline to apply data calibration and screening, and
then extracted the images, light curves, and spectra using nuproducts. The
X-ray emission of P9 was detected with both the FPMA and FPMB below $\sim$20
keV. Figure 3 shows the stacked images obtained by all of the _XMM-Newton_
EPIC detectors in the 0.3–10 keV band, and _NuSTAR_ FPMA+FPMB in the 3–10 keV,
10–20 keV, and 20–79 keV ranges. The emission of P9 is dominated by soft
X-rays below 10 keV. For the 10–20 keV band, P9 was marginally detected with a
detection significance of 8.8$\sigma$ using the wavelet source detection
algorithm wavdetect.

We extracted the light curve of each _NuSTAR_ orbit to check if P9 showed
long-term variability. We extracted the source events from a circular aperture
of 60″ radius and background events from an annulus with an inner radius of
80″ and an outer radius of 120″. The result is plotted in Figure 2 together
with the scaled X4 light curve obtained with _XMM-Newton_. The _NuSTAR_ flux
dropped significantly after $\sim$MJD 57529.5.

Figure 2:  (Left) _Chandra_ 0.5–7 keV light curves of NGC 7793 P9. (Right)
_XMM-Newton_ (0.3–10 keV) and _NuSTAR_ (3–79 keV) light curves of NGC 7793 P9.
The X4 light curve on the third panel was scaled to fit the _NuSTAR_ count
rate. The bottom panel showed the soft (0.1–1 keV) and hard (1–10 keV) X-ray
light curves of X4 observation.

##  3 X-Ray Spectroscopy

We started the spectral analysis with the _Chandra_ , _XMM-Newton_ , and
_NuSTAR_ data sets, which have much better signal-to-noise ratios than the
_Swift_ data sets. We used the absorption model
tbnew333http://pulsar.sternwarte.uni-erlangen.de/wilms/research/tbabs/ (Wilms
et al., 2000), fixed the Galactic absorption column density to
$N_{\rm{H}}(\rm{gal})=1.23\times 10^{20}$ cm-2, and set an additional
absorption for the external galaxy (Kalberla et al., 2005). We followed the
empirical classification strategy of ULX spectral states proposed by Sutton et
al. (2013). We first employed simple phenomenological models, the multicolor
disk blackbody (MCD; Mitsuda et al., 1984) and the power-law (PL), to fit the
spectra. Then, we tried a combination of them (MCD+PL) to check if adding an
additional component improves the fit. For this case, we further check the
empirical classification of the UL spectra by computing the photon index
($\Gamma$) and the fraction of the PL flux in 0.3–1 keV. After that, we fit
the spectrum with more physical models like the slim disk and the Comptonized
models to further characterize the spectral behaviors of those data sets with
the UL spectral features. The best-fit parameters for each _Chandra_ and _XMM-
Newton_ observations with a simple MCD+PL model were summarized in Table 2.

###  3.1 Phenomenological Models

####  3.1.1 _Chandra_ Spectroscopy

The _Chandra_ 0.5–7 keV spectra were generated using the specextract command,
which creates source and background spectra, corresponding response files
(rmf), and ancillary files (arf). The spectra were then grouped to a minimum
of 15 photons per bin for C1 and 25 photons per bin for C2–C4 for spectral
fitting. The C1 observation made in 2003 can be well described with a simple
PL but $N_{\rm{H}}=1.3^{+1.7}_{-1.3}\times 10^{21}$ cm-2 is not well
constrained. The uncertainties of the spectral parameters are obtained within
a 90% confidence interval. The photon index is $\Gamma=2.4^{+0.5}_{-0.4}$,
which is quite soft but with a large uncertainty. To doubly confirm the result
obtained with the $\chi^{2}$ statistic, we tried the Cash statistic without
grouping the photons and yielded a consistent result (Cash, 1979). The C2
spectrum is well fit by the MCD with an inner-disk temperature
$kT_{\rm{in}}=0.96_{-0.06}^{+0.07}$ keV. The normalization constant indicates
an apparent inner radius of $R_{in}\sqrt{\cos{\theta}}\approx 50$ km where
$\theta$ is the inclination angle. We tried the MCD+PL model but obtained no
significant improvement. The C3 and C4 observations were taken four months
after the C2 observation and had a time separation of five days. An absorbed
PL model with $\Gamma\approx 2.3$ and $N_{\rm{H}}\sim 10^{21}$ cm-2 can well
describe both spectra, and adding an MCD does not improve the fit statistic
significantly. All of the _Chandra_ spectra and the corresponding best-fit
models are shown in Figure 4. It is clear that the C2 spectrum, which was best
described with an MCD model, has a stronger curvature than other three data
sets that can be described by a PL model. We checked the spectral fits with
Cash statistics and the results are fully consistent with those obtained with
$\chi^{2}$ statistics.

####  3.1.2 _XMM-Newton_ Spectroscopy

To extract high-quality _XMM-Newton_ spectra, we filtered out the background
flaring epochs and chose the events with good grades by restricting FLAG$=$0.
We set PATTERN$\leq$4 for the PN data and PATTERN$\leq$12 for the MOS data.
The source and background regions were normalized using the backscale command,
and the corresponding response files and ancillary files were created using
rmfgen and arfgen. The spectra were then grouped to a minimum of 15 photons
per bin for X1 and 25 photons per bin for the X3 and X4 observations using
specgroup.

Figure 3: X-ray images of NGC 7793 P9 obtained with (a) _XMM-Newton_ in 0.3–10
keV, and _NuSTAR_ in (b) 3–10 keV,(c) 10–20 keV, and (d) 10–79 keV. The _XMM-
Newton_ image is stacked from MOS1, MOS2, and PN images, and the _NuSTAR_
images are stacked from the FPMA and FPMB images. The green cross denotes the
_Chandra_ position of P9.

We fit the PN, MOS1, and MOS2 spectra in 0.3–10 keV simultaneously with all
the parameters tied between instruments. To further correct for any
differences in calibrations between the detectors, we added an overall
normalization constant for the MOS1 and MOS2 spectra. The discrepancies in the
following analysis are well below $\lesssim$5%. The X1 observation was taken
in the tail of the canonical outburst, roughly five months after the _Chandra_
C3 and C4 observations. The spectra can be well described with an MCD model.
The best-fit $kT_{\rm{in}}$ is $0.79\pm 0.06$ keV and the apparent inner
radius of the accretion disk is $R_{in}\sqrt{\cos{\theta}}=53_{-8}^{+9}$ km,
consistent with the C2 spectrum taken in 2011. Similar to all the _Chandra_
observations, the fitting cannot be significantly improved by two-component or
more complex models. The spectral change from C2 to X1 suggests that P9 had a
canonical outburst in 2011–2012, and its spectrum switched between an MCD with
an inner radius of $\sim 50$ km and a PL with $\Gamma\sim 2.3$.

Table 2: Best-fit Spectral Parameters of NGC 7793 P9 obtained from the
_Chandra_ , _XMM-Newton_ , and _NuSTAR_ Data Sets.

Obs. | Model | $L_{\textrm{0.3--10 keV}}$ | $L_{\rm{MCD}}$ | $N_{\rm{H}}$ | $R_{\rm{in}}\sqrt{\cos\theta}$ | $kT_{\rm{in}}$ | $\Gamma$ | $\chi^{2}/\rm{dof}$ | State
---|---|---|---|---|---|---|---|---|---
|  | ($10^{39}$ erg s-1) | ($10^{39}$ erg s-1) | ($10^{22}$ cm-2) | (km) | (keV) |  |  |
C1 | PL | $0.10_{-0.02}^{+0.03}$ |  $\cdots$ | $0.13_{-0.13}^{+0.17}$ |  $\cdots$ |  $\cdots$ | $2.4_{-0.4}^{+0.5}$ | $6.8/14$ | SPL
C2 | MCD | $0.5\pm 0.1$ | $0.5\pm 0.1$ | $0.07\pm 0.04$ | $49^{+7}_{-6}$ | $0.96^{+0.07}_{-0.06}$ |  $\cdots$ | $62.7/57$ | Thermal
C3 | PL | $0.79\pm 0.06$ |  $\cdots$ | $0.09\pm 0.04$ |  $\cdots$ |  $\cdots$ | $2.2\pm 0.1$ | $62.2/67$ | SPL
C4 | PL | $0.83_{-0.08}^{+0.10}$ |  $\cdots$ | $0.14_{-0.07}^{+0.08}$ |  $\cdots$ |  $\cdots$ | $2.3\pm 0.02$ | $27.1/31$ | SPL
X1 | MCD | $0.26_{-0.07}^{+0.08}$ | $0.26_{-0.07}^{+0.08}$ | $0.024_{-0.014}^{+0.018}$ | $53^{+8}_{-7}$ | $0.78\pm 0.02$ |  $\cdots$ | $77.8/93$ | Thermal
X3 | MCD+PL | $1.1\pm 0.2$ | $0.36^{+0.07}_{-0.08}$ | $0.07\pm 0.02$ | $370^{+90}_{-60}$ | $0.34\pm 0.04$ | $1.8\pm 0.2$ | $171.9/189$ | UL (HUL)
X4a | MCD+PL | $3.4_{-0.4}^{+0.6}$ | $1.4^{+0.5}_{-0.4}$ | $0.11_{-0.02}^{+0.03}$ | $34^{+7}_{-6}$ | $1.5\pm 0.2$ | $2.0_{-0.1}^{+0.2}$ | $308.3/311$ | UL (BD?)
| MCD+PLb | $4_{-2}^{+5}$ | $0.6_{-0.4}^{+0.8}$ | $0.06\pm 0.04$ | $100_{-60}^{+90}$ | $0.8_{-0.2}^{+0.5}$ | $1.4_{-0.8}^{+0.4}$ | $302.1/309$ |

aafootnotetext: The _NuSTAR_ data set was also included in the fitting.

bbfootnotetext: A high-energy cutoff was included in this PL component with a
cutoff energy of $4_{-4}^{+1}$ keV and the folding energy of $5\pm 5$ keV.

The X3 observation was made in the early phase of the ultraluminous outburst.
We found that neither the MCD model nor the PL model provides acceptable
fitting. On the other hand, a two-component MCD+PL model can significantly
improve the fitting with a $\chi^{2}/dof=171.9/189$ (see Figure 4). The best-
fit $N_{\rm{H}}$ is $(7\pm 2)\times 10^{20}$ cm-2. The temperature and the
radius of the inner disk are $kT_{\rm{in}}=0.34\pm 0.04$ keV and
$R_{\rm{in}}\sqrt{\cos\theta}=370^{+90}_{-60}$ km, which imply a much cooler
disk with an innermost radius much farther away from the central compact
object compared to the canonical thermal and SPL states of a BHXB. The PL
component has $\Gamma=1.8(2)$, clearly harder than the value of the canonical
SPL spectrum.

The X4 observation was made near the peak of the ultraluminous outburst. We
included the _NuSTAR_ 3–20 keV data in the fitting to increase the statistic
and constrain the spectral behavior in the hard X-ray band. Since the _XMM-
Newton_ X4 observation was made before this drop, we extracted the _NuSTAR_
spectra during the X4 observation time interval. We fit all the five spectra
(PN, MOS1, MOS2, FPMA, FPMB) simultaneously. Similar to the X3 observation, a
two-component MCD+PL model provides a good statistic with
$\chi^{2}/\rm{dof}=308.3/311$. The disk temperature is $kT_{\rm{in}}=1.5\pm
0.2$ keV and $\Gamma=2.0\pm 0.2$. In the soft X-ray band, the apparent
spectral curvature caused by a highly optically thick Comptonized component
may be classified as a disk-like spectrum (Sutton et al., 2013). Therefore, we
used an MCD plus PL with a high-energy cutoff to fit the data. The fitting
statistic is equally good with a $\chi^{2}/\rm{dof}=302.1/309$. The
$N_{\rm{H}}$ value is $(7\pm 4)\times 10^{20}$ cm-2, consistent with that
obtained in the previous section. The cutoff energy is $5\pm 5$ keV,
indicating a strong spectral curvature caused by an optically thick corona but
not well constrained.

####  3.1.3 Determining the Extragalactic Absorption

Most of the best-fit models yield a similar $N_{\rm{H}}$ with a value of
($\sim$5–10)$\times 10^{20}$ cm-2. Therefore, we assume that $N_{\rm{H}}$ did
not significantly vary during the entire observation time span. We fit all the
_Chandra_ , _XMM-Newton_ , and _NuSTAR_ data simultaneously with a tied
$N_{\rm{H}}$ and their corresponding best-fit models to further constrain the
$N_{\rm{H}}$. The resulting $\chi^{2}/\rm{dof}$ is $767.3/789$, and
$N_{\rm{H}}$ is determined as $(7\pm 1)\times 10^{20}$ cm-2. We also tried to
use the abundances with the absorption model phabs with the corresponding
abundance and obtained $N_{\rm{H}}=(6\pm 1)\times 10^{20}$ cm-2, slightly
lower than that obtained with tbnew (Anders & Grevesse, 1989). Other
parameters are fully consistent with those listed in Table 2. The $N_{\rm{H}}$
is consistent with that of P13, indicating that P9 is unlikely a heavily
absorbed source compared to other point sources in NGC 7793 (Israel et al.,
2017b).

###  3.2 More Physical Models

From the spectral shapes in Figure 4 and fitting parameters in Table 2, we
noticed that the X3 and X4 spectra likely represent the ultraluminous states.
X3 can be well described with a cool disk plus a hard PL without a significant
high-energy cutoff below 10 keV. The spectral shape is similar to BHXB spectra
in the low/hard state but the luminosity is much higher. This is the typical
spectral feature of hard ULXs. The low $kT_{\rm{in}}<0.5$ keV and hard
$\Gamma<2$ implies an HUL spectrum according to the empirical classification
(Sutton et al., 2013). The X4 spectrum is disk-like but much broader than a
standard MCD. This spectrum can be classified as a BD. However, X4 shows a
significant difference from X3 only above $\sim$1 keV, implying that the major
difference could be contributed by the PL tail. This heavily Comptonized
spectra are sometimes misclassified as a BD. Hence, we carefully examine the
X3 and X4 spectra with more physical Comptonized models and slim disk models
to explore the underlying physics.

####  3.2.1 Slim Disk

We first fit X4 with a $p$-free accretion disk (DiskPBB), in which the radial
advection is taken into account (Mineshige et al., 1994; Hirano et al., 1995).
This model is an approximation of the slim disk. The radial temperature
profile in this model has a PL dependence $T(r)\propto r^{-p}$, where $r$ is
the radius and $p$ is a free parameter. A standard MCD corresponds to
$p=0.75$, and a slim disk can be approximated with $p<0.75$. The best-fit
$\chi^{2}/\rm{dof}$ is $308.3/312$, slightly better than the statistic of the
MCD+PL model. The inner-disk temperature is $kT_{\rm{in}}=2.2\pm 0.2$ keV,
much higher than that in the high/soft state of the normal outburst (see Table
3). The low $p=0.57\pm 0.02$ is consistent with other ULXs with a BD spectrum
(Kaaret et al., 2017). The apparent inner radius,
$R_{\rm{in}}\sqrt{\cos\theta}=16\pm 3$ km, is extremely small. We further
tried a table model for a slim disk (Kawaguchi, 2003), in which included
corrections for local BB contribution, Comptonized spectrum, or relativistic
effects. We found that the spectrum can be described by the model ID 4, which
includes a slim disk and a Comptonized local spectrum. The resulting
$\chi^{2}/\rm{dof}$ is $304.7/311$, which is equally good as the $p$-free disk
model. The BH mass is $4^{+27}_{-2}M_{\odot}$, $N_{\rm{H}}$ is $(7\pm 2)\times
10^{20}$ cm-2, the mass accretion rate is $\dot{m}=6^{+5}_{-1}$
$L_{\rm{Edd}}/c^{2}$, but the viscosity is not well constrained. $N_{\rm{H}}$
is fully consistent with previous measurements. We also tried to use the
$p$-free disk model to fit X3 spectra; however, the statistic is unacceptably
high since the slim disk is unable to describe the high-energy tail. We
therefore add a PL component and found that the disk component has a small
value of $p$ but the uncertainty is large.

####  3.2.2 Comptonized Model: SIMPL

The disk corona usually contributes a high-energy tail in the X-ray spectrum
of a BHXB. We described the tail with the convolution model SIMPL. This model
empirically calculates the Comptonized spectrum from the scattering of a
fraction of seed photons, i.e., the MCD component (Steiner et al., 2009). The
seed photons are redistributed to both higher and lower energies. The
normalization factor is the scattered fraction ($f_{\rm{sc}}$) We fit the X3
data set, which results in an acceptable statistic of
$\chi^{2}/dof=172.9/189$. The disk is partially scattered with a fraction of
$f_{\rm{sc}}=0.33$, and has an apparent inner-disk temperature of
$kT_{\rm{in}}=0.33\pm 0.04$ keV and a radius of
$R_{\rm{in}}\sqrt{\cos\theta}=500_{-100}^{+160}$ km. This result still shows a
truncated disk similar to that obtained from phenomenological models.

We fit the X4 data with the same model, but the best-fit statistic
($\chi^{2}/\rm{dof}=366.2/315$) is much worse than that yielded by MCD+PL. We
obtained $f_{\rm{sc}}=1$, implying that all the disk photons are scattered to
the Comptonized component. The seed MCD component has $kT_{\rm{in}}=0.92\pm
0.6$ keV and $R_{\rm{in}}\sqrt{\cos\theta}=109_{-12}^{+14}$ km, a bit closer
to the BH than that obtained from the X3 data set.

####  3.2.3 Comptonized Model: NTHCOMP

In addition to the convolution model, we used an additive Comptonized model,
NTHCOMP, to describe the Comptonized continuum (Zdziarski et al., 1996; Życki
et al., 1999). The high-energy rollover, which is sharper than an exponential
cutoff, is parameterized by the electron temperature $kT_{e}$. We set the
maximum seed photon temperature $kT_{bb}$ to the inner-disk temperature
$kT_{\rm{in}}$ of the MCD component. The statistic of X3 is insensitive to the
electron temperature $kT_{e}$ because no significant high-energy cutoff is
observed below 10 keV. Hence, we fixed it to 100 keV. The fit statistic is
acceptable and the parameter suggests a cool disk with a large inner radius
(see Table 3). The true disk radius can be re-estimated using the equation
provided by Kubota & Makishima (2004) for the presence of disk Comptonization.
Assuming an inclination angle of 60∘, $R_{\rm{in}}$ can be estimated as
$770_{-140}^{+300}$ km, slightly larger than that directly obtained from the
MCD component.

The Comptonized component of X4 can also be described with this model.
However, the normalization of the disk component is consistent with zero,
implying that MCD is negligible. Therefore, we used NTHCOMP to fit the spectra
and set $kT_{bb}$ as a free parameter. The statistic is as good as that
obtained from phenomenological models. The best-fit
$kT_{e}=1.9_{-0.3}^{+0.5}$, roughly corresponds to an exponential cutoff at
4–7 keV. This is consistent with that obtained from the MCD+PL model with a
high-energy cutoff. The inner-disk radius can be re-estimated as
$460_{-140}^{+280}$ km, consistent with that obtained from the X3 data set.

Figure 4: (Left) The _Chandra_ spectra of NGC 7793 P9. The black, red, blue,
and green points represent the C1, C2, C3, and C4 spectra, respectively. The
solid lines are the corresponding best-fit models. (Right) The _XMM-Newton_
and _NuSTAR_ spectra of NGC 7793 P9. The solid lines are the corresponding
best-fit MCD, MCD+PL, and DiskPBB models for the X1, X3, and X4 data sets,
respectively.

####  3.2.4 Disk–Corona Coupling: OPTXAGNF

The above two models account for the partial Comptonization of a standard MCD.
However, the size of the corona may play an important role if the disk is only
partially covered. The innermost accretion disk could be fully or partially
Comptonized, while the outer disk remains standard and contributes a
significant amount of soft X-ray emission. Done & Kubota (2006) proposed the
DKBBFTH model to deal with the spectrum from the Comptonized corona over the
inner disk. This model has been applied to the spectra of several ULXs to
recover the properties of the inner disk (Gladstone et al., 2009). An upgraded
model, OPTXAGNF, was proposed to fit the broadband spectra of narrow-line
Seyfert 1 galaxies with the same aspect (Done et al., 2012). Given a similar
accretion property, this model can be used to describe the UL and SPL states
of BHXBs. The disk component is a canonical MCD where the radius is larger
than the coronal radius. The disk inside the corona is not fully thermalized.
Instead, the emission was Comptonized to form both soft excess and a hard PL
tail. There are 11 parameters in OPTXAGNF: the BH mass ($M_{\rm{BH}}$),
distance, luminosity in units of Eddington ratio ($\log L/L_{\rm{Edd}}$), BH
spin ($a$), coronal radius ($R_{\rm{cor}}$ in units of
$r_{g}=GM_{\rm{BH}}/c^{2}$), outer radius of the disk ($R_{\rm{out}}$),
electron temperature ($kT_{e}$) and optical depth ($\tau$) for the soft
Comptonization component, photon index of the hard Comptonization component
($\Gamma$), fraction ($f_{\rm{pl}}$) of the power below $R_{\rm{cor}}$ that is
redistributed to the hard PL, and the redshift. We fit the X3 and X4 data sets
simultaneously and tied the BH mass with each one, and assumed
$N_{\rm{H}}=7\times 10^{20}$ cm-2 and $a=0$. We found that the $R_{\rm{cor}}$
for these two data sets hit the hard limit. We therefore set
$R_{\rm{cor}}=100$ $r_{g}$. A BH mass of $M_{\rm{BH}}=22\pm 3$ $M_{\odot}$ can
result in a good statistic (see Table 3). The luminosities are $1.1\pm 0.2$
$L_{\rm{Edd}}$ and $1.3\pm 0.2$ $L_{\rm{Edd}}$ for the X3 and X4 data sets,
respectively. The Comptonized component has a high $\tau$ for both data sets,
implying an optically thick corona. We also tried to set the BH spin to the
maximum value of $a=0.998$. In this fitting, $R_{\rm{cor}}$ is well
constrained but $\tau$ is not. The fitting is insensitive to $kT_{e}$. Hence,
we fixed it to 1 keV to account for the soft excess. We found that the BH mass
is $120\pm 20$ $M_{\odot}$, six times larger than that of the $a=0$ case. The
corona has a radius of $\sim$15 $r_{g}$ for both the X3 and X4 spectra,
roughly the same size as $\sim$100 $r_{g}$ of a $\sim$20 $M_{\odot}$ BH.
Similar to the best-fit result for a non-spinning BH, a highly Comptonized,
optically thick corona is required for both data sets.

###  3.3 Spectral Variability of _Swift_ Observations

To obtain an overall hardness–intensity relationship, we fit the individual
_Swift_ , _Chandra_ , and _XMM-Newton_ data sets with a simple PL model and
fixed $N_{\rm{H}}$ to $7\times 10^{20}$ cm-2 (see Figure 5). From the overall
$\Gamma$–$L$ relationship obtained from individual _Swift_ observations, we
found that $\Gamma$ varied over a wide range when the luminosity is lower than
$10^{39}$ erg s-1, likely caused by the transition between the SPL and the
high/soft states as seen in canonical BHXBs. In addition, $\Gamma$ shows a
break and a decreasing trend when the luminosity is higher than $10^{39}$ erg
s-1. We then stacked the _Swift_ observations with similar luminosities and
spectral behaviors to increase the signal-to-noise ratio (see Table 4 for
definition and result). We fit the six data sets individually with
$N_{\rm{H}}$ fixed to $7\times 10^{20}$ cm-2 using the the Cash statistic. We
also binned the spectra and fit them with the minimum $\chi^{2}$ algorithm,
and the results were fully consistent with those obtained from Cash statistic.
We fit the spectra with MCD+PL to get hints about the spectral behavior. We
did not attempt to use more physical models to fit the _Swift_ spectra due to
a limited number of photons.

Table 3: Best-fit Spectral Parameters of the X3 and X4 Data Sets with More
Physical Models.

Model | Parameters | X3 | X4
---|---|---|---
|  $N_{\rm{H}}$ ($10^{22}$ cm-2) | $0.12\pm 0.05$ | $0.09\pm 0.02$
|  $T_{\rm{in}}$ (keV) | $0.5_{-0.1}^{+0.2}$ | $2.2\pm 0.2$
DiskPBB+PL |  $R_{\rm{in}}\sqrt{\cos\theta}$ (km) | $120_{-80}^{+190}$ | $16\pm 3$
| $p$ | $0.47_{-0.06}^{+0.20}$ | $0.57\pm 0.02$
| $\Gamma$ | $1.5_{-0.5}^{+0.3}$ |  $\cdots$
| $\chi^{2}/\rm{dof}$ | $168.2/188$ | $308.3/312$
|  $N_{\rm{H}}$ ($10^{22}$ cm-2) | $0.05\pm 0.02$ | $0.007\pm 0.007$
|  $T_{\rm{in}}$ (keV) | $0.33\pm 0.04$ | $0.92_{-0.06}^{+0.07}$
SIMPL(MCD) |  $R_{\rm{in}}\sqrt{\cos\theta}$ (km) | $500_{-100}^{+170}$ | $110_{-12}^{+14}$
| $\Gamma$ | $1.8\pm 0.2$ | $2.9_{-0.2}^{+0.3}$
| $f_{\rm{sc}}$ | $0.33_{-0.07}^{+1.0}$ | $1.0^{a}$
| $\chi^{2}/\rm{dof}$ | $172.9/189$ | $367.4/314$
|  $N_{\rm{H}}$ ($10^{22}$ cm-2) | $0.05\pm 0.02$ | $0.04_{-0.01}^{+0.02}$
|  $T_{\rm{in}}$ (keV) | $0.33\pm 0.04$ | $0.4\pm 0.1^{b}$
|  $R_{\rm{in}}\sqrt{\cos\theta}$ (km) | $400_{-70}^{+100}$ |  $\cdots$
MCD+NTHCOMP |  $kT_{e}^{d}$ (keV) | $100^{a}$ | $1.9_{-0.3}^{+0.5}$
| $\Gamma$ | $1.8\pm 0.2$ | $1.9\pm 0.1$
| $\chi^{2}/\rm{dof}$ | $173.2/189$ | $310.7/313$
| BH Mass ($M_{\odot}$) | $22\pm 3$
|  $\log L$ ($L_{\rm{Edd}}$) | $0.02\pm 0.07$ | $0.10\pm 0.08$
|  $R_{\rm{cor}}$ ($r_{g}$) | $100^{a}$ | $100^{a}$
OPTXAGNF | $\tau$ | $100^{a}$ | $15_{-2}^{+7}$
($a=0$) | $\Gamma$ | $1.60_{-0.05}^{+0.04}$ | $1.9\pm 0.3$
|  $kT_{e}^{e}$ (keV) | $0.5\pm 0.1$ | $1.2\pm 0.2$
| $f_{\rm{pl}}$ | $0.97\pm 0.01$ | $0.4\pm 0.2$
| $\chi^{2}/\rm{dof}$ | $500.6/502$
| BH Mass ($M_{\odot}$) | $120\pm 20$
|  $\log L$ ($L_{\rm{Edd}}$) | $-0.7_{-0.1}^{+0.2}$ | $-0.6_{-0.1}^{+0.2}$
|  $R_{\rm{cor}}$ ($r_{g}$) | $15_{-2}^{+3}$ | $15_{-6}^{+12}$
OPTXAGNF | $\tau$ | $24_{-7}^{+76\,c}$ | $100^{a}$
($a=0.998$) | $\Gamma$ | $1.61\pm 0.06$ | $2.0\pm 0.2$
|  $kT_{e}^{e}$ (keV) | $1.0^{a}$ | $1.0^{a}$
| $f_{\rm{pl}}$ | $1.0^{a}$ | $0.71_{-0.08}^{+0.01}$
| $\chi^{2}/\rm{dof}$ | $510.6/503$

aafootnotetext: This parameter is frozen.

bbfootnotetext: This temperature is $kT_{bb}$ in the NTHCOMP model because the
MCD component is not included.

ccfootnotetext: The upper bound of $\tau$ is not well constrained.

ddfootnotetext: The electron temperature in the NTHCOMP model denotes the
high-energy rollover. An exponential cutoff could be roughly estimated as
$E_{c}=$2–3 $kT_{e}$.

eefootnotetext: The electron temperature in the OPTXAGNF model denotes the
soft excess caused by Comptonization. The electron temperature for the hard
Comptonization component is fixed to 100 keV.

Figure 5: Photon indices obtained by all the _Swift_ individual observations
(gray), _Chandra_ (blue circle), _XMM-Newton_ (red square), and stacked
_Swift_ observations (magenta triangle).

The S1 spectrum ($L<10^{39}$ erg s-1) is obviously disk-like. It can be well
described by the MCD model with $kT_{\rm{in}}=1.3_{-0.2}^{+0.3}$ keV and
$R_{\rm{in}}\sqrt{\cos\theta}=33^{+11}_{-9}$ km. We attempted to fit the S1
spectrum with the DiskPBB model and found $p=0.7_{-0.1}^{+0.2}$, fully
consistent with that of an MCD. The S2–S4 data sets are better described by a
simple PL compared to the MCD. We obtained an anticorrelation between $\Gamma$
and flux from these three data sets (see Figure 5). We tried to fit the
spectra with the MCD+PL model, and the statistics were slightly improved. The
best-fit MCD component suggested that the disk is cool with $kT_{\rm{in}}\sim
0.5$ keV and the inner radius is larger than the innermost stable circular
orbit. We tried the DiskPBB model and yielded equally good statistics. The
best-fit $p$-values are much smaller than those obtained from X4 and close to
the hard limit of the DiskPBB model.

On the other hand, the S5 and S6 spectra are disk-like and can be well
described with an MCD. Adding a PL component to S5 slightly improved the
statistic but the change in the MCD component is limited and the PL component
was not well constrained. The MCD fit of S5 has
$kT_{\rm{in}}=1.4_{-0.3}^{+0.4}$ keV and
$R_{\rm{in}}\sqrt{\cos\theta}=46^{+14}_{-13}$ km, similar to the MCD component
in X4. This spectrum is equally well described with a $p$-free slim disk. The
best-fit $p=0.65_{-0.07}^{+0.07}$ is consistent with that obtained in the X4
data set but much closer to that of a standard MCD. For the S6 spectrum,
adding a PL component or fitting $p$ could not improve the fit statistic. The
parameters of the PL components are not constrained, and the $p$-value is
fully consistent with 0.75. We plot all of the stacked _Swift_ spectra with
the best-fit model in Figure 6 and rebin the spectral channels to have at
least 5$\sigma$ significance. It is clear that the spectral shape changes from
PL-like to disk-like as flux increased when $L\gtrsim 10^{39}$ erg s-1. S3 and
S4 are similar to X3. They show a hint of a hard PL tail although the
significance is low due to the limited energy range and low photon statistic.
S5 and S6 are disk-like, but they do not have an obvious deviation from the S3
and S4 spectra below $\sim 1$ keV. This behavior is similar to that of X4,
implying that the apparent disk-like spectrum could be a misleading.

##  4 Discussion

Table 4: Best-fit Spectral Parameters of Stacked _Swift_ Spectra During the
Ultraluminous Outburst.

Bins | Luminosity Range | Model | $L_{\textrm{0.3--10 keV}}^{a}$ | $L_{\textrm{MCD}}$ | $R_{\rm{in}}\sqrt{\cos\theta}$ | $kT_{\rm{in}}$ |  $\Gamma$ or $p$ | $C^{2}/\rm{dof}$
---|---|---|---|---|---|---|---|---
| ($10^{39}$ erg s-1) |  | ($10^{39}$ erg s-1) | ($10^{39}$ erg s-1) | (km) | (keV) |  |
S1 |  $L<1$, $\Gamma<1.8$ | MCD | $0.8\pm 0.1$ | $0.8\pm 0.1$ | $33^{+11}_{-9}$ | $1.3^{+0.3}_{-0.2}$ |  $\cdots$ | $143.9/167$
S2 |  $L<1$, $\Gamma>1.8$ | MCD+PL | $0.86\pm 0.8$ | $0.3\pm 0.1$ | $150_{-100}^{+120}$ | $0.5_{-0.1}^{+0.3}$ | $1.9(5)$ | $224.8/254$
|  | DiskPBB |  |  $\cdots$ | $<48$ | $1.1_{-0.2}^{+0.4}$ | $0.53_{-0.03}^{+0.04\,b}$ | $223.2/255$
S3 |  $L=1$–1.3 | MCD+PL | $1.2\pm 0.8$ | $0.4\pm 0.2$ | $300_{-150}^{+180}$ | $0.39_{-0.08}^{+0.12}$ | $1.4_{-0.9}^{+0.5}$ | $149.1/184$
|  | DiskPBB |  |  $\cdots$ | $<25$ | $1.3_{-0.4}^{+1.5}$ | $0.53_{-0.04}^{+0.05\,b}$ | $154.3/184$
S4 |  $L=1.3$–2 | MCD+PL | $1.5\pm 0.2$ | $0.6_{-0.3}^{+0.5}$ | $150_{-80}^{+90}$ | $0.6_{-0.1}^{+0.2}$ | $1.5_{-0.9}^{+0.4}$ | $258.3/267$
|  | DiskPBB |  |  $\cdots$ | $23_{-21}^{+18}$ | $1.4_{-0.3}^{+0.8}$ | $0.56_{-0.03}^{+0.04\,b}$ | $264.8/267$
S5 |  $L=2$–3 | MCD+PL | $2.2\pm 0.2$ | $1.7_{-0.6}^{+0.3}$ | $46^{+14}_{-13}$ | $1.4^{+0.2}_{-0.3}$ | $2.3^{+6.2}_{-1.2}$ | $280.0/324$
|  | DiskPBB |  |  $\cdots$ | $34^{+20}_{-17}$ | $1.5^{+0.4}_{-0.2}$ | $0.65^{+0.07\,b}_{-0.05}$ | $279.8/325$
S6 | $L>3$ | MCD | $3.6_{-0.3}^{+0.4}$ | $3.6_{-0.3}^{+0.4}$ | $56^{+11}_{-10}$ | $1.5^{+0.2}_{-0.1}$ |  $\cdots$ | $253.5/312$

aafootnotetext: Derived from the best-fit model.

bbfootnotetext: This is the $p$-value of the $p$-free disk.

Figure 6: _Swift_ stacked spectra and the best-fit MCD+PL model P9.

###  4.1 Implications from the MCD+PL Model

####  4.1.1 Spectral Evolution of P9

During the canonical outburst, the X-ray spectrum of P9 switched between an
MCD and a PL. The MCD spectra were observed in mid 2011 and mid 2012. The disk
had an apparent inner radius of $R_{\rm{in}}\sqrt{\cos\theta}\approx 50$ km
and an inner-disk temperature $kT_{\rm{in}}=0.7$–0.9 keV. This radius
corresponds to a true inner radius of 90 km after applying corrections for the
color–temperature and inner boundary condition of an MCD and assuming an
inclination angle of 60∘ (Shimura & Takahara, 1995; Kubota et al., 1998).
Assuming a non-spinning BH, the mass can be estimated to be $\sim 10$
$M_{\odot}$. For an extremely spinning BH, the apparent value corresponds to a
BH mass of $\sim 60$ $M_{\odot}$. In addition, the disk could be weakly
Comptonized and result in a smaller $R_{\rm{in}}$ (Kubota & Makishima, 2004).
The above mass estimation represents the lower limit of the BH mass in P9.
This indicates that P9 has a regular or massive stellar-mass BH. Its X-ray
spectrum switched between the SPL and high/soft states of a canonical outburst
before 2013 (Tetarenko et al., 2016).

When the luminosity exceeds $10^{39}$ erg s-1, $\Gamma$ shows a clear
decreasing trend with the luminosity (see Figure 5), and the high-quality
_XMM-Newton_ spectra in this ultraluminous state cannot be well fit by single-
component models. They are significantly different from those in the canonical
outburst. When the luminosity is slightly higher than $L\approx 10^{39}$ erg
s-1 (X3, S3, and S4), P9 has a cool disk and a hard PL tail. Moreover, the
spectral curvature increases with flux when $L\gtrsim 2\times 10^{39}$ erg s-1
(X4, S5, S6). These data sets can be described by disk-like spectral models.

As described in Section 1, the BD usually dominates the spectra of ULXs with
$L\lesssim 3\times 10^{39}$ erg s-1, and those of ULXs with $L\gtrsim 3\times
10^{39}$ erg s-1 showed two-component spectra like HUL and SUL (Gladstone et
al., 2009; Sutton et al., 2013). However, some ULXs like NGC 1313 X-2 and
Holmberg IX X-1 showed an inverse transition, i.e., an HUL spectrum in the
low-luminosity state, and a disk-like BD spectrum near their peak luminosities
(Pintore & Zampieri, 2012; Middleton et al., 2015; Luangtip et al., 2016). The
BD feature of these sources is likely to be misclassified due to the strong
spectral curvature caused by the highly optically thick Comptonized spectrum
(Sutton et al., 2013). Similar to these two sources, P9 has an apparent
inverse transition between these two states. However, the X4 light curve shows
a significant variability dominated by the emission above 1 keV
($F_{\rm{var}}=0.15(2)$ for 1–10 keV and $F_{\rm{var}}=0.05(2)$ for 0.3–1
keV). It is difficult to interpret with a single slim disk model since the
variability should be independent of the energy for a single emission
component. Hence, we suggest that the apparent BD spectral shape observed in
X4 consists of two components like NGC 1313 X-2 and Holmberg IX X-1.

####  4.1.2 Variability of the Disk Component

P9 evolved from the canonical BHXB regime to the ULX regime when its
luminosity crossed the Eddington limit of a $\sim 10$ $M_{\odot}$ BH. To
obtain the disk evolution, we calculated the 0.3–10 keV luminosity of the disk
to check its correlation with the inner-disk temperature. Assuming a fixed
inner radius, the luminosity of a standard disk shows a positive correlation
with the disk temperature:

$L_{\rm{MCD}}=\frac{\pi\sigma
G^{2}M^{2}T_{\rm{in}}^{4}}{6c_{0}^{4}f_{\rm{col}}^{4}c^{4}},$ (2)

where $\sigma$ is the Stefan–Boltzmann constant, $M$ is the mass of the
central compact object, $c_{0}\approx 0.1067$, and $f_{\rm{col}}^{4}$ is the
color–temperature correction, which is usually larger than unity (Gierliński
et al., 1999; Gierliński & Done, 2004). Many ULXs (especially supersoft ones)
show an anticorrelation between the disk luminosity and the temperature (see,
e.g., Kajava & Poutanen, 2009; Urquhart & Soria, 2016). On the other hand,
some ULXs show a positive correlation close to $L_{\rm{disk}}\propto T^{4}$
(Miller et al., 2003). Since P9 has both canonical and ultraluminous
outbursts, it provides a unique chance to test the accretion disk model.

We collected the 0.3–10 keV luminosity of the MCD component and disk
temperature from both the high/soft state (C2, X1, and S1), HUL state (X3, S3,
S4, and possibly S2), and disk-like state in the ultraluminous outburst (X4,
S5, and S6). We show the disk luminosity–temperature relation of P9 together
with that of a BHXB, XTE J1550$-$564, in Figure 7. The disk configuration of
P9 during the high/soft state in the canonical outburst (X1 and C2) was fully
consistent with the thermal track of XTE J1550$-$564\. We attempted to add an
NTHCOMP component to these two data sets but the Componized component
contributed negligible emission. The disk component in the disk-like spectral
states of P9 during the ultraluminous outburst well lies on the extension of
the thermal track. This is likely a coincidence since they can also be well
described by Comptonized models. The MCD component in the HUL state clustering
at the lower-left corner of this figure is similar to the ultraluminous branch
of XTE J1550$-$564\. However, P9 shows a much harder ($\Gamma\sim$1.6–1.8) PL
tail compared to XTE J1550$-$564 ($\Gamma\sim$2.5–2.8). The presence of this
state and the apparent disk-like state makes P9 as a unique source among
canonical BHXBs.

Figure 7: Comparison of the $kT_{\rm{in}}$–$L_{\rm{MCD}}$ relation between P9
and XTE J1550$-$564 (Soria, 2007). The red dashed line is the best-fit
$L_{\rm{MCD}}\propto T_{\rm{in}}^{4}$ relation for XTE J1550$-$564\. The red
squares are obtained from the high/soft state of XTE J1550$-$564 and the green
ones are from the SPL state and ultraluminous branch. The black data points
denote the thermal and BD states of P9, and the blue ones denote the HUL
spectral state. The black circle with the dashed line is the best-fit thermal
component of X4 spectra with a two-component model containing an MCD and a PL
with a high-energy cutoff.

###  4.2 Implications from the Comptonized Disk

The simple MCD+PL model is helpful for classifying the spectral states. The
result suggests that the innermost region of the accretion disk was covered by
the Comptonized corona. We further used physical models to deal with the
coupling of both the disk and the corona. We noticed that all of the spectra
with luminosities above $L\gtrsim 10^{39}$ erg s-1, including X3, X4, and
S3–S6, show clear spectral difference above $\sim$1 keV (see Figures 4 and 6).
This implies that the mass accretion rate during the ultraluminous outburst
has no significant difference, and the observed variability over (1–5)$\times
10^{39}$ erg s-1 is dominated by the change of the Comptonized component. If
we use NTHCOMP to describe the Comptonized corona, then a cool inner-disk
temperature, a large re-estimated radius, a hard PL with $\Gamma\sim 1.8$, and
a large optical depth are necessary for both X3 and X4 spectra. In contrast,
X4 shows a higher normalization and a clear cutoff at 4–6 keV. This implies
that the plasma temperature of the corona is hot at the beginning of the
ultraluminous outburst, and a huge amount of disk photons are scattered to
hard X-ray. As the corona temperature cools down, fewer photons can be
scattered to hard X-rays. The spectrum shows a cutoff in the _XMM-Newton_
energy range and an apparently higher luminosity in 0.3–10 keV.

We also described the UL spectra with the OPTXAGNF model. The BH mass can be
estimated as $\sim 20$ $M_{\odot}$ for a non-spinning BH or $\sim 120$
$M_{\odot}$ for an extremely fast spinning BH. This result is a bit higher
than that obtained from the simple estimation based on the apparent inner-disk
radius in the high/soft state. In the case of $M_{\textrm{BH}}\sim
20M_{\odot}$, most of the disk emission below $R_{\rm{cor}}$ is scattered into
the hard Comptonized component ($f_{\rm{pl}}\sim 1$) and produces a hard PL
with $\Gamma=1.6$ in X3 spectrum. When the X-ray luminosity increases to near
the peak of the ultraluminous outburst, a significant power is emitted in the
low-temperature soft Comptonization component. This result is consistent with
the cooling of the corona, which is implied from the NTHCOMP model. The
luminosity is close to or slightly higher than the Eddington luminosity, which
suggests that P9 belongs to the major ULX population: the high-luminosity tail
of BHXBs. Of course, we could not exclude the explanation that X4 shows a slim
disk since the spectrum can be well fit by both the DiskPBB model and the slim
disk model. If this is true, the corona just fades away when the luminosity is
close to the peak. However, the short-term variability is only observed above
1 keV in the X4 observation. Therefore, the apparent BD spectrum likely
contains two components.

A BH with a mass of $\sim 120M_{\odot}$ can also produce the observed spectra
with a smaller coronal radius (in units of $r_{g}$). The luminosity is roughly
0.2 $L_{\rm{Edd}}$, consistent with that of canonical outbursts in BHXBs. If
this is true, P9 is consistent with a massive stellar-mass BH that was
recently found via gravitational wave observations (Abbott et al., 2016,
2017). However, the hard PL tail observed in X3 is not expected in canonical
outbursts. Hence, this scenario is less favored.

###  4.3 Supercritical Accretion

The UL spectra can also be physically interpreted with the supercritical
accretion model. The inner disk is geometrically thick and an optically thick
outflowing wind is blown. A funnel-shaped cavity forms along the rotational
axis of the BH (King et al., 2001; Ohsuga & Mineshige, 2007; Middleton et al.,
2015). In this case, the soft emission component is contributed by the wind,
and the hard emission component is from the inner accretion disk or the
accretion flow, although detailed modeling is still under development
(Poutanen et al., 2007; Kajava et al., 2012). The beaming factor was
proportional to the mass accretion rate. Hence, the inclination angle is
expected to be low since the hard component dominates the observed flux during
the HUL state. The mass accretion rate is maximum during the HUL state (X3),
where the beaming factor is large and we observed the innermost accretion
flow. When the mass accretion rate decreased, the inner-disk temperature
slightly cooled down and the beaming factor decreased. The resulting spectrum
is best described by a BD if the emission from the wind becomes less dominant.
As the disk luminosity drops to sub-Eddington levels, we observe a thin disk
in the tail of the outburst. We attempted to fit the hard component with
DiskPBB in the X3 data set, but the disk temperature is extremely high
($kT\sim 4$ keV) and the disk radius is not well constrained. Therefore, this
scenario is less favored although not entirely excluded. Of course, a corona
can possibly form around the supercritical accretion disk and cause a much
more complex X-ray spectrum.

Finally, P9 shows the HUL spectrum when $L\approx 10^{39}$ erg s-1, much lower
than several genuine ULXs showing the same spectral behavior (e.g.,
1–1.2$\times 10^{40}$erg s-1 for Holmberg IX X-1, Sutton et al., 2013). This
could imply that those genuine ULXs have much higher BH masses if the spectral
state classification only depends on the Eddington ratio. However, the
luminosity of a ULX can be strongly affected by the beaming factor and viewing
angle. An important counterexample is NGC 5907 ULX-1, which shows an HUL
spectral behavior at $L\gtrsim 10^{40}$ erg s-1. Recently, it was found to be
powered by a neutron star (Israel et al., 2017a). Therefore, the spectral
state and the mass of the central compact object cannot be directly
associated. The empirical luminosity range of the spectral states may stop
being valid as the ULX sample grows. Comprehensively analyzing the timing and
spectral behaviors of ULXs with next-generation X-ray telescopes and improving
the supercritical accretion spectral model will be helpful for describing the
accretion disk structure of ULXs.

##  5 Summary

We carried out a systematic analysis of the spectral variability of the newly
discovered transient ULX NGC 7793 P9 with _Swift_ XRT monitoring and dedicated
_Chandra_ , _XMM-Newton_ , and _NuSTAR_ observations. P9 is an LMXB and showed
a canonical outburst behavior before 2013. The ultraluminous outburst after
2014 suggests that a BHXB can evolve to a ULX, supporting the scenario that
some ULXs are powered by long-lasting super-Eddington outbursts. P9 is an
important example bridging the BHXB and the ULX and shows the most fruitful
spectral states among all the ULXs. The spectral analysis implies that the
compact object is a stellar-mass BH. The simple MCD+PL model suggests that P9
shows an ultraluminous branch, similar to XTE J1550$-$564\. However, the
spectrum is dominated by a much harder PL compared to XTE J1550$-$564\. This
is the feature of an HUL spectrum of a ULX. We suggest that the inner
accretion disk is covered by a corona that scatters most of the inner-disk
emission to a hard PL. In addition, P9 evolved to an apparent BD spectral
state when $L\gtrsim 3\times 10^{39}$ erg s-1 although it is likely
misclassified due to the strong curvature of the Comptonized component. We
suggested that the corona cools down but the disk component remains unchanged.
The soft Comptonization dominated the observed spectrum. The spectral
evolution can also be interpreted with the supercritical accretion model,
although a detailed spectrum of this model is still under development. The
theoretical models can be further tested by tracking the spectral evolution of
more samples of transient ULXs.

We thank Prof. Chris Done for useful suggestions in modeling the X-ray spectra
and classifying the spectral states. We thank the anonymous referee for the
comments that improved this paper. This research is in part based on the data
obtained from the Chandra Data Archive, and has made use of software provided
by the Chandra X-ray Center (CXC) in the application packages CIAO, ChIPS, and
Sherpa. This research has used the observations obtained with _XMM-Newton_ and
ESA science mission with instruments and contributions directly funded by the
ESA member states and NASA. This research has also made use of data obtained
with _NuSTAR_ , a project led by Caltech, funded by NASA and managed by
NASA/JPL, and has utilized the NUSTARDAS software package, jointly developed
by the ASDC (Italy) and Caltech (USA). This work made use of data supplied by
the UK Swift Science Data Centre at the University of Leicester. C.P.H. and
C.Y.N. are supported by a GRF grant from the Hong Kong Government under HKU
17300215P. A.K.H.K. is supported by the Ministry of Science and Technology of
the Republic of China (Taiwan) through grants 103-2628-M-007-003-MY3 and
105-2112-M-007-033-MY2.

## References

  * (1)
  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review Letters, 116, 061102
  * Abbott et al. (2017) —. 2017, Physical Review Letters, 118, 221101
  * Anders & Grevesse (1989) Anders, E., & Grevesse, N. 1989, GeCoA, 53, 197
  * Bachetti et al. (2014) Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nature, 514, 202
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939
  * Crowther (2013) Crowther, P. A. 2013, MNRAS, 428, 1927
  * de Jager et al. (1989) de Jager, O. C., Raubenheimer, B. C., & Swanepoel, J. W. H. 1989, A&A, 221, 180
  * Done et al. (2012) Done, C., Davis, S. W., Jin, C., Blaes, O., & Ward, M. 2012, MNRAS, 420, 1848
  * Done & Kubota (2006) Done, C., & Kubota, A. 2006, MNRAS, 371, 1216
  * Edelson et al. (2002) Edelson, R., Turner, T. J., Pounds, K., et al. 2002, ApJ, 568, 610
  * Evans et al. (2007) Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379
  * Evans et al. (2009) —. 2009, MNRAS, 397, 1177
  * Feng & Kaaret (2006) Feng, H., & Kaaret, P. 2006, ApJ, 650, L75
  * Feng & Soria (2011) Feng, H., & Soria, R. 2011, New A Rev., 55, 166
  * Frank et al. (2002) Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (3rd ed; Cambridge: Cambridge Univ. Press), 398
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, Proc. SPIE, 6270, 62701V
  * Gierliński & Done (2004) Gierliński, M., & Done, C. 2004, MNRAS, 347, 885
  * Gierliński et al. (1999) Gierliński, M., Zdziarski, A. A., Poutanen, J., et al. 1999, MNRAS, 309, 496
  * Gladstone et al. (2009) Gladstone, J. C., Roberts, T. P., & Done, C. 2009, MNRAS, 397, 1836
  * Hirano et al. (1995) Hirano, A., Kitamoto, S., Yamada, T. T., Mineshige, S., & Fukue, J. 1995, ApJ, 446, 350
  * Hu et al. (2017) Hu, C.-P., Li, K. L., Kong, A. K. H., Ng, C.-Y., & Chun-Che Lin, L. 2017, ApJ, 835, L9
  * Israel et al. (2017a) Israel, G. L., Belfiore, A., Stella, L., et al. 2017a, Science, 355, 817
  * Israel et al. (2017b) Israel, G. L., Papitto, A., Esposito, P., et al. 2017b, MNRAS, 466, L48
  * Kaaret & Feng (2009) Kaaret, P., & Feng, H. 2009, ApJ, 702, 1679
  * Kaaret & Feng (2013) —. 2013, ApJ, 770, 20
  * Kaaret et al. (2017) Kaaret, P., Feng, H., & Roberts, T. P. 2017, ARA&A, 55, 303
  * Kajava & Poutanen (2009) Kajava, J. J. E., & Poutanen, J. 2009, MNRAS, 398, 1450
  * Kajava et al. (2012) Kajava, J. J. E., Poutanen, J., Farrell, S. A., Grisé, F., & Kaaret, P. 2012, MNRAS, 422, 990
  * Kalberla et al. (2005) Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775
  * Karachentsev et al. (2003) Karachentsev, I. D., Grebel, E. K., Sharina, M. E., et al. 2003, A&A, 404, 93
  * Kawaguchi (2003) Kawaguchi, T. 2003, ApJ, 593, 69
  * King (2002) King, A. R. 2002, MNRAS, 335, L13
  * King et al. (2001) King, A. R., Davies, M. B., Ward, M. J., Fabbiano, G., & Elvis, M. 2001, ApJ, 552, L109
  * Kong et al. (2016) Kong, A. K. H., Hu, C.-P., Lin, L. C.-C., et al. 2016, MNRAS, 461, 4395
  * Kubota & Done (2004) Kubota, A., & Done, C. 2004, MNRAS, 353, 980
  * Kubota & Makishima (2004) Kubota, A., & Makishima, K. 2004, ApJ, 601, 428
  * Kubota et al. (1998) Kubota, A., Tanaka, Y., Makishima, K., et al. 1998, PASJ, 50, 667
  * Lau et al. (2016) Lau, R. M., Kasliwal, M. M., Bond, H. E., et al. 2016, ApJ, 830, 142
  * Luangtip et al. (2016) Luangtip, W., Roberts, T. P., & Done, C. 2016, MNRAS, 460, 4417
  * Middleton et al. (2015) Middleton, M. J., Heil, L., Pintore, F., Walton, D. J., & Roberts, T. P. 2015, MNRAS, 447, 3243
  * Middleton et al. (2012) Middleton, M. J., Sutton, A. D., Roberts, T. P., Jackson, F. E., & Done, C. 2012, MNRAS, 420, 2969
  * Miller et al. (2003) Miller, J. M., Fabbiano, G., Miller, M. C., & Fabian, A. C. 2003, ApJ, 585, L37
  * Mineshige et al. (1994) Mineshige, S., Hirano, A., Kitamoto, S., Yamada, T. T., & Fukue, J. 1994, ApJ, 426, 308
  * Mitsuda et al. (1984) Mitsuda, K., Inoue, H., Koyama, K., et al. 1984, PASJ, 36, 741
  * Motch et al. (2014) Motch, C., Pakull, M. W., Soria, R., Grisé, F., & Pietrzyński, G. 2014, Nature, 514, 198
  * Ohsuga & Mineshige (2007) Ohsuga, K., & Mineshige, S. 2007, ApJ, 670, 1283
  * Pannuti et al. (2011) Pannuti, T. G., Schlegel, E. M., Filipović, M. D., et al. 2011, AJ, 142, 20
  * Pavlovskii et al. (2017) Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092
  * Pintore & Zampieri (2012) Pintore, F., & Zampieri, L. 2012, MNRAS, 420, 1107
  * Poutanen et al. (2007) Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G., & Abolmasov, P. 2007, MNRAS, 377, 1187
  * Radburn-Smith et al. (2011) Radburn-Smith, D. J., de Jong, R. S., Seth, A. C., et al. 2011, ApJS, 195, 18
  * Rappaport et al. (2005) Rappaport, S. A., Podsiadlowski, P., & Pfahl, E. 2005, MNRAS, 356, 401
  * Read & Pietsch (1999) Read, A. M., & Pietsch, W. 1999, A&A, 341, 8
  * Shimura & Takahara (1995) Shimura, T., & Takahara, F. 1995, ApJ, 445, 780
  * Soria (2007) Soria, R. 2007, Ap&SS, 311, 213
  * Steiner et al. (2009) Steiner, J. F., Narayan, R., McClintock, J. E., & Ebisawa, K. 2009, PASP, 121, 1279
  * Sutton et al. (2013) Sutton, A. D., Roberts, T. P., & Middleton, M. J. 2013, MNRAS, 435, 1758
  * Tetarenko et al. (2016) Tetarenko, B. E., Sivakoff, G. R., Heinke, C. O., & Gladstone, J. C. 2016, ApJS, 222, 15
  * Tully et al. (2016) Tully, R. B., Courtois, H. M., & Sorce, J. G. 2016, AJ, 152, 50
  * Urquhart & Soria (2016) Urquhart, R., & Soria, R. 2016, MNRAS, 456, 1859
  * Vierdayanti et al. (2010) Vierdayanti, K., Mineshige, S., & Ueda, Y. 2010, PASJ, 62, 239
  * Walton et al. (2016) Walton, D. J., Fürst, F., Bachetti, M., et al. 2016, ApJ, 827, L13
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914
  * Zdziarski et al. (1996) Zdziarski, A. A., Johnson, W. N., & Magdziarz, P. 1996, MNRAS, 283, 193
  * Życki et al. (1999) Życki, P. T., Done, C., & Smith, D. A. 1999, MNRAS, 309, 561
  * (67)


# Wheeler-DeWitt quantization of point-particles

Ward Struyve111Instituut voor Theoretische Fysica, KU Leuven,
Belgium∗222Centrum voor Logica en Filosofie van de Wetenschappen, KU Leuven,
Belgium

###### Abstract

We present the Hamiltonian formulation of a relativistic point-particle
coupled to Einstein gravity and its canonical quantization à la Wheeler-
DeWitt. In the resulting quantum theory, the wave functional is a function of
the particle coordinates and the 3-metric. It satisfies a particular
Hamiltonian and diffeomorphism constraint, together with a Klein-Gordon-type
equation. As usual in the Wheeler-DeWitt theory, the wave function is time-
independent. This is also reflected in the Klein-Gordon-type equation, where
the time derivative is absent. Before considering gravity, we consider the
coupling of a particle with electromagnetism, which is treated similarly, but
simpler.

##  1 Introduction

Canonical quantization of Einstein gravity leads to the Wheeler-DeWitt theory
where the state is described by a wave functional $\Psi(h_{ij})$, which is a
function of a 3-metric $h_{ij}({\bf x})$ [1]. A common way to introduce matter
is through a scalar field, so that the wave functional also is a function of a
scalar field $\phi({\bf x})$, i.e., $\Psi=\Psi(h_{ij},\phi)$. A remarkable
feature is that the wave functional turns out to be time-independent. The wave
equations are merely constraints on the wave functional: the Hamiltonian and
the diffeomorphism constraint. The time-independence of $\Psi$ leads to the
question of how time evolution can be accounted for. There are various
attempts to answer this question, but there is far from a consensus on the
right approach [2, 3, 1, 4].

In this paper, we consider the coupling of gravity to point-particles instead
of a scalar field. This allows for the study of quantum gravity for a fixed
number of particles without having to consider quantum field theory. This
system has been studied in the context of the problem of time by Rovelli [5,
6] and Pavšič [7, 8].

Classically, the history of a point-particle is described by its world-line
$X^{\mu}(\lambda)$, which depends on an arbitrary parameter $\lambda$, while
the space-time metric $g_{\mu\nu}(x)$ is a function of space-time $x=(t,{\bf
x})$.333Units are assumed such that $\hbar=c=1$. The particle and field are
hence parameterized differently, with $X^{0}(\lambda)=t$. To pass to the
Hamiltonian picture, a common temporal parameter for the field and the
particle needs to be chosen. This is already familiar from the case of the
free particle in Minkowski space-time. In that case, the disadvantage of
choosing the parameter $t$ is that one needs to take a square root to express
the Hamiltonian in terms of the momenta, resulting in the Hamiltonian
$H=\pm\sqrt{{\bf p}^{2}+m^{2}}$. Choosing the parameter $\lambda$ instead
leads to a zero Hamiltonian, but with the constraint $p_{\mu}p^{\mu}+m^{2}=0$.
Quantization of the latter leads to the Klein-Gordon equation, while the
former leads to the positive or negative frequency part of the Klein-Gordon
equation.

The situation in the case of gravity is similar. Rovelli [5, 6] considered the
parameter $t$ and ended up taking a square root. Pavšič seems to use $t$ for
the field, while using $\lambda$ for the particle. This does not amount to a
standard Hamiltonian formulation and hence makes this quantization of the
theory questionable. In this paper, the Hamiltonian formulation is presented
using $\lambda$ as the common temporal parameter, together with the
corresponding quantum theory. This avoids needing to take the square root and
as such generalizes Rovelli’s treatment. The wave function is again time-
independent; It does not depend on $X^{0}$, just as in Rovelli’s treatment,
but contrary to the findings of Pavšič. So $X^{0}$ does not appear as a
possible candidate for the evolution parameter (contrary to for example the
treatment of dust [9]). It will appear however that Pavšič’s theory (when
suitably interpreted) can be derived from ours in a mixed Schrödinger-
Heisenberg picture.

To start with, the treatment of point-particles coupled to an electromagnetic
field is considered, in Minkowski space-time, which shares many features with
the gravitational case.

##  2 Electrodynamics

To formulate the quantum theory of a point-particle interacting with an
electromagnetic field, we start from the classical theory, pass to the
Hamiltonian picture, and apply the usual canonical quantization methods [10,
11, 12, 13, 14].

Consider a classical point-particle with worldline $X^{\mu}(\lambda)$
interacting with an electromagnetic field with vector potential
$A^{\mu}=(A_{0},{\bf A})$. Derivatives with respect to $\lambda$ and $t$ will
be denoted by respectively by primes and dots throughout, e.g.
$X^{\prime\mu}=dX^{\mu}/d\lambda$. It is also assumed throughout that
$X^{\prime 0}>0$. The action is

${\mathcal{S}}={\mathcal{S}}_{M}+{\mathcal{S}}_{F},$ (1)

with

${\mathcal{S}}_{M}=-m\int
d\lambda\sqrt{{X^{\prime}}^{\mu}(\lambda)X^{\prime}_{\mu}(\lambda)}-e\int
d\lambda{X^{\prime}}^{\mu}(\lambda)A_{\mu}(X(\lambda))$ (2)

and

${\mathcal{S}}_{F}=-\frac{1}{4}\int d^{4}xF_{\mu\nu}(x)F^{\mu\nu}(x),$ (3)

where $F^{\mu\nu}=\partial^{\mu}A^{\nu}-\partial^{\nu}A^{\mu}$ is the
electromagnetic field tensor. The temporal parameter in the matter action and
the field action are different. In the matter action, the parameter $\lambda$
acts as temporal parameter, while $t$ does it for the electromagnetic field.
In order to pass to a Hamiltonian formulation, a common temporal coordinate is
needed. We will first use $t$ and then $\lambda$.

To write the matter action in terms of $t$, the identity $1=\int
dt\delta(t-X^{0}(\lambda))$ can be inserted in the action and integrate over
$\lambda$. Defining ${\widetilde{X}}^{\mu}(t)=X^{\mu}(\lambda(t))$, where
$X^{0}(\lambda(t))\equiv 1$, leads to

${\mathcal{S}}_{M}=-m\int dt\sqrt{1-\dot{\widetilde{\bf X}}(t)^{2}}+e\int
dt\dot{\widetilde{\bf X}}(t)\cdot{\bf A}(t,{\widetilde{\bf X}}(t))-e\int
dtA_{0}(t,{\widetilde{\bf X}}(t)).$ (4)

Writing ${\mathcal{S}}=\int dtL({\widetilde{\bf X}},\dot{\widetilde{\bf
X}},A,\dot{A})$, with $L$ the Lagrangian, the canonical momenta are

${\widetilde{\bf P}}=\frac{\partial L}{\partial\dot{\widetilde{\bf
X}}}=m\frac{{\widetilde{\bf X}}}{\sqrt{1-\dot{\widetilde{\bf X}}^{2}}}+e{\bf
A}({\widetilde{\bf X}}),$ (5) $\Pi_{0}({\bf x})=\frac{\delta
L}{\delta\dot{A}_{0}({\bf x})}=0,\qquad{\boldsymbol{\Pi}}({\bf
x})=\frac{\delta L}{\delta\dot{\bf A}({\bf x})}=\dot{\bf A}({\bf
x})+{\boldsymbol{\nabla}}A_{0}({\bf x}).$ (6)

The canonical Hamiltonian is

$\displaystyle H$ $\displaystyle={\bf P}\cdot\dot{\widetilde{\bf X}}+\int
d^{3}x(\pi_{0}\dot{A}_{0}+{\boldsymbol{\Pi}}\cdot\dot{\bf A})-L$
$\displaystyle=\frac{m}{\sqrt{1-\dot{\widetilde{\bf
X}}^{2}}}+eA_{0}({\widetilde{\bf X}})+H_{F},$ (7)

where $H_{F}$ is the usual Hamiltonian for the free electromagnetic field:

$H_{F}=\int
d^{3}x\left[\frac{1}{2}{\boldsymbol{\Pi}}^{2}+\frac{1}{2}({\boldsymbol{\nabla}}\times{\bf
A})^{2}-A_{0}{\boldsymbol{\nabla}}\cdot{\boldsymbol{\Pi}}\right].$ (8)

As is well known, there are two constraints: the primary constraint

$\Pi_{0}=0,$ (9)

and the secondary constraint (which is the Gauss law)

${\boldsymbol{\nabla}}\cdot{\boldsymbol{\Pi}}-e\delta({\bf x}-{\widetilde{\bf
X}})=0.$ (10)

The matter part of the Hamiltonian (7) is not yet expressed in terms of the
momenta. Using

$\frac{m^{2}}{1-\dot{\widetilde{\bf X}}(t)^{2}}=m^{2}+\left({\widetilde{\bf
P}}-e{\bf A}({\widetilde{\bf X}})\right)^{2},$ (11)

we obtain

$H=\pm\sqrt{m^{2}+\left({\widetilde{\bf P}}-e{\bf A}({\widetilde{\bf
X}})\right)^{2}}+eA_{0}({\widetilde{\bf X}})+H_{F},$ (12)

which involves the choice of a square root.

Quantization in the Schrödinger picture (and dropping the tildes) leads to the
following equations for $\Psi({\bf X},A,t)$:

$\textrm{i}\partial_{t}\Psi=\pm\sqrt{m^{2}-\left({\boldsymbol{\nabla}}-\textrm{i}e{\bf
A}({\bf X})\right)^{2}}\Psi+eA_{0}({\bf X})\Psi+{\widehat{H}}_{F}\Psi,$ (13)
$\frac{\delta\Psi}{\delta
A_{0}}=0,\qquad\textrm{i}{\boldsymbol{\nabla}}\cdot\frac{\delta\Psi}{\delta{\bf
A}({\bf x})}+e\delta({\bf x}-{\bf X})\Psi=0,$ (14)

where

${\widehat{H}}_{F}=\frac{1}{2}\int d^{3}x\left(-\frac{\delta^{2}}{\delta{\bf
A}({\bf x})^{2}}+[{\boldsymbol{\nabla}}\times{\bf A}({\bf x})]^{2}\right).$
(15)

The equations (14) are the constraints (9) and (10) that are imposed as
operator constraints. The square root could be eliminated by considering the
square of (13) to obtain a Klein-Gordon-like equation. If $\lambda$ rather
than $t$ is taken as a common temporal parameter then this is indeed the
equation that will be obtained, as we will now show.

To use $\lambda$ as the temporal parameter, $1=\int
d\lambda\delta(\lambda-\lambda(t))$ is inserted in the field action. Defining
${\widetilde{A}}(\lambda,{\bf x})=A(X^{0}(\lambda),{\bf x})$, so that
${\widetilde{A}}^{\prime}=\dot{A}X^{\prime 0}$, the field action can be
written as

${\mathcal{S}}_{F}=\int dtL_{F}(A,\dot{A})=\int d\lambda
L^{*}_{F}({\widetilde{A}},{\widetilde{A}}^{\prime}),$ (16)

with

$L^{*}_{F}({\widetilde{A}},{\widetilde{A}}^{\prime})=X^{\prime
0}L_{F}\left({\widetilde{A}},{\widetilde{A}}^{\prime}/X^{\prime 0}\right).$
(17)

With ${\mathcal{S}}=\int d\lambda
L^{*}(X,X^{\prime},{\widetilde{A}},{\widetilde{A}}^{\prime})$, the conjugate
momenta are

${\widetilde{\Pi}}_{\mu}({\bf x})=\frac{\delta
L^{*}}{\delta{\widetilde{A}}^{\prime\mu}({\bf x})}=\frac{\delta
L}{\delta\dot{A}^{\mu}({\bf x})}=\Pi_{\mu}({\bf x}),$ (18)
$P_{\mu}=\frac{\partial L^{*}}{\partial
X^{\prime\mu}}=-m\frac{X^{\prime}_{\mu}}{\sqrt{{X^{\prime}}^{\nu}X^{\prime}_{\nu}}}-e{\widetilde{A}}_{\mu}({\bf
X})-\delta^{0}_{\mu}H_{F}({\widetilde{A}},{\widetilde{\Pi}}).$ (19)

So the expressions for the momenta for the field are just the same as before,
cf. (6). The particle momentum $P_{0}$ gets a contribution from field
Lagrangian which is just minus the field Hamiltonian (8).

The canonical Hamiltonian is

$H^{*}=P_{\mu}X^{\prime\mu}+\int d^{3}x{\widetilde{\Pi}}_{\mu}({\bf
x}){\widetilde{A}}^{\prime\mu}({\bf x})-L^{*}=0$ (20)

and is zero because of the reparameterization invariance of the action. There
are two primary constraints: $\chi_{1}=\Pi_{0}=0$ (as before) and

$\chi_{2}=\left[P_{\mu}+e{\widetilde{A}}_{\mu}({\bf
X})+\delta^{0}_{\mu}H_{F}({\widetilde{A}},{\widetilde{\Pi}})\right]\left[P^{\mu}+e{\widetilde{A}}^{\mu}({\bf
X})+\delta^{\mu}_{0}H_{F}({\widetilde{A}},{\widetilde{\Pi}})\right]-m^{2}=0.$
(21)

There is a secondary constraint which follows from the requirement that the
Poisson bracket of the two primary constraints with the total Hamiltonian

$H^{*}_{T}=H^{*}+\int d^{3}x\lambda_{1}({\bf x})\chi_{1}({\bf
x})+\lambda_{2}\chi_{2}$ (22)

vanishes, which results in $[\chi_{1}({\bf x}),\chi_{2}]_{P}=0$ or

$\frac{\delta
H_{F}({\widetilde{A}},{\widetilde{\Pi}})}{\delta{\widetilde{A}}_{0}({\bf
x})}+e\delta({\bf x}-{\bf X})=0,$ (23)

which amounts to the Gauss constraint (10).

Quantization (again dropping the tildes) leads to the following equation for
$\Psi(X,A)$:

$\left(\partial_{\mu}+\textrm{i}eA_{\mu}({\bf
X})+\textrm{i}\delta^{0}_{\mu}{\widehat{H}}_{F}\right)\left(\partial^{\mu}+\textrm{i}eA^{\mu}({\bf
X})+\textrm{i}\delta^{\mu}_{0}{\widehat{H}}_{F}\right)\Psi+m^{2}\Psi=0,$ (24)

with ${\widehat{H}}_{F}$ as before in (15), together with the constraints
(14). This Klein-Gordon-type equation is just the square of (13), provided
$X^{0}$ is identified with $t$.

So far we have dealt with just a single particle. How to extend the theory to
many particles? The theory with the square root Hamiltonian (13) is directly
generalized to many particles (13), with now $\psi=\psi({\bf X}_{1},\dots,{\bf
X}_{n},A,t)$:

$\textrm{i}\partial_{t}\Psi=\sum_{k}\left[\sqrt{m^{2}-\left({\boldsymbol{\nabla}}_{k}-\textrm{i}e{\bf
A}({\bf X}_{k})\right)^{2}}+eA_{0}({\bf
X}_{k})\right]\Psi+{\widehat{H}}_{F}\Psi,$ (25) $\frac{\delta\Psi}{\delta
A_{0}}=0,\qquad\textrm{i}{\boldsymbol{\nabla}}\cdot\frac{\delta\Psi}{\delta{\bf
A}({\bf x})}+e\sum_{k}\delta({\bf x}-{\bf X}_{k})\Psi=0.$ (26)

We have chosen the positive square root Hamiltonian for each particle. The
Klein-Gordon-type equation (24) is not that straightforwardly generalizable to
many particles because it is not of the Schrödinger form like (13). One option
is to run through the quantization procedure again for many particles, but
this seems rather complicated. Another option is to recast the Klein-Gordon-
type equation (24) in a Schrödinger form, which could be done using the Kemmer
formulation [15] (which actually concerns a Dirac-like equation for spin-0),
and which at least in the case of external field is directly extendable to
many particles [16]. Yet another option is to use the multi-time picture where
the wave function depends on the time-component of each of the particles [17].
We will not pursue this further here.

##  3 Gravity

The analysis of gravity proceeds completely analogously as that of
electromagnetism. The action for a classical point-particle coupled to gravity
is

${\mathcal{S}}={\mathcal{S}}_{M}+{\mathcal{S}}_{G},$ (27)

with

${\mathcal{S}}_{M}=-m\int
d\lambda\sqrt{g_{\mu\nu}(X(\lambda)){X^{\prime}}^{\mu}(\lambda)X^{\prime\nu}(\lambda)}$
(28)

and

${\mathcal{S}}_{G}=-\frac{1}{\kappa}\int d^{4}x\sqrt{-g}R,$ (29)

with $\kappa=16\pi G$, is the Einstein-Hilbert action.

Before considering the Hamiltonian picture in terms of the temporal parameter
$\lambda$, we recall the usual Hamiltonian formulation for the Einstein-
Hilbert action in terms of $t$ (following the conventions of [18]). It is
supposed that space-time can be foliated in terms of space-like hypersurfaces
such that the space-time manifold is diffeomorphic to
${\mathbb{R}}\times\Sigma$, with $\Sigma$ a 3-surface. Coordinates
$x^{\mu}=(t,{\bf x})$ are chosen such that the time coordinate $t$ labels the
leaves of the foliation and ${\bf x}$ are the coordinates on $\Sigma$. In
terms of these coordinates, the metric and its inverse are written as

$g_{\mu\nu}=\begin{pmatrix}N^{2}-N_{i}N^{i}&-N_{i}\\\
-N_{i}&-h_{ij}\end{pmatrix}\,,\qquad
g^{\mu\nu}=\begin{pmatrix}\frac{1}{N^{2}}&\frac{-N^{i}}{N^{2}}\\\
\frac{-N^{i}}{N^{2}}&\frac{N^{i}N^{j}}{N^{2}}-h^{ij}\end{pmatrix}\,,$ (30)

where $N$ and $N^{i}$ are respectively the lapse and the shift vector, and
$h_{ij}$ is the Riemannian metric on $\Sigma$, with $h$ its determinant.
Spatial indices of $h_{ij}$ and $N^{i}$ and the corresponding momenta are
raised and lowered with this spatial metric.

With $L_{G}$ the gravitational Lagrangian, the canonically conjugate momenta
are

$\pi^{ij}=\frac{\delta
L_{G}}{\delta{\dot{h}}_{ij}}=-\frac{1}{\kappa}\sqrt{h}(K^{ij}-Kh^{ij})$ (31)

(which is a tensor density of weight $-1$), with

$K_{ij}=\frac{1}{2N}({D_{i}N_{j}+D_{j}N_{i}-\dot{h}}_{ij}),\qquad
K=K_{ij}h^{ij},$ (32)

the extrinsic curvature, and

$\pi=\frac{\delta L_{G}}{\delta N}=0,\qquad\pi_{i}=\frac{\delta L_{G}}{\delta
N^{i}}=0.$ (33)

The Hamiltonian is

$H_{G}=\int_{\Sigma}d^{3}x\left(N{\mathcal{H}}+N^{i}{\mathcal{H}}_{i}\right),$
(34)

with

${\mathcal{H}}=\kappa
G_{ijkl}\pi^{ij}\pi^{kl}+{\mathcal{V}},\qquad{\mathcal{H}}_{i}=-2h_{ik}D_{j}\pi^{jk},$
(35)

where $G_{ijkl}=(h_{ij}h_{jl}+h_{il}h_{jk}-h_{ij}h_{kl})/2\sqrt{h}$ is the
DeWitt metric, $D_{i}$ is the covariant derivative corresponding to $h_{ij}$,
and ${\mathcal{V}}=-\sqrt{h}R^{(3)}/\kappa$ is the gravitational potential
density. Apart from the primary constraints (33), there are also the secondary
constraints

${\mathcal{H}}=0,\qquad{\mathcal{H}}_{i}=0.$ (36)

In order to use the temporal variable $\lambda$, we introduce
${\widetilde{g}}_{\mu\nu}(\lambda,{\bf x})=g_{\mu\nu}(X^{0}(\lambda),{\bf x})$
and similarly for other variables. The total action is

${\mathcal{S}}=\int d\lambda
L^{*}(X,X^{\prime},{\widetilde{g}},{\widetilde{g}}^{\prime}).$ (37)

The momenta for the metric are the same as before, i.e.,

${\widetilde{\pi}}^{ij}=\frac{\delta
L^{*}}{\delta{\widetilde{h}}^{\prime}_{ij}}=\frac{\delta
L_{G}}{\delta\dot{h}_{ij}}=\pi^{ij},\quad{\widetilde{\pi}}=\frac{\delta
L^{*}}{\delta{\widetilde{N}}}=0,\quad{\widetilde{\pi}}_{i}=\frac{\delta
L^{*}}{\delta{\widetilde{N}}^{i}}=0,$ (38)

and the momentum for the particle becomes (indices of $X^{\mu}$ and $P_{\mu}$
are raised and lowered with the metric ${\widetilde{g}}_{\mu\nu}({\bf X})$)

$P_{\mu}=\frac{\partial L^{*}}{\partial
X^{\prime\mu}}=-m\frac{X^{\prime}_{\mu}}{\sqrt{{X^{\prime}}^{\nu}X^{\prime}_{\nu}}}-\delta^{0}_{\mu}{\widetilde{H}}_{G},$
(39)

where

${\widetilde{H}}_{G}=H_{G}({\widetilde{N}},{\widetilde{N}}^{i},{\widetilde{g}}_{ij},{\widetilde{\pi}}^{ij}).$
(40)

The canonical Hamiltonian is zero, i.e.,

$H^{*}=P_{\mu}X^{\prime\mu}+\int_{\Sigma}d^{3}x{\widetilde{\pi}}^{ij}({\bf
x}){\widetilde{g}}^{\prime}_{ij}({\bf x})-L^{*}=0.$ (41)

There are three primary constraints:

${\widetilde{\pi}}=0,\qquad{\widetilde{\pi}}_{i}=0,$ (42)
$\chi={\widetilde{g}}^{\mu\nu}({\bf
X})\left[P_{\mu}+\delta^{0}_{\mu}{\widetilde{H}}_{G}\right]\left[P_{\nu}+\delta^{0}_{\nu}{\widetilde{H}}_{G}\right]-m^{2}=0.$
(43)

There are two secondary constraints corresponding to
$\delta\chi/\delta{\widetilde{N}}=0$ and
$\delta\chi/\delta{\widetilde{N}}^{i}=0$, resulting in

${\widetilde{\mathcal{H}}}({\bf x})=\delta({\bf x}-{\bf
X})\left({\widetilde{N}}({\bf X})p^{0}+\frac{1}{{\widetilde{N}}({\bf
X})}{\widetilde{H}}_{G}\right)$ (44)

and

${\widetilde{\mathcal{H}}}_{i}({\bf x})=\delta({\bf x}-{\bf X})p_{i},$ (45)

where (as in (40)) the tilde on ${\widetilde{\mathcal{H}}}$ and
${\widetilde{\mathcal{H}}}_{i}$ denotes the functions (35) evaluated for
fields with tildes. There are no further constraints.

By multiplying (44) with ${\widetilde{N}}({\bf x})$ and (45) with
${\widetilde{N}}^{i}({\bf x})$, integrating over all space, and using
$p^{0}=(p_{0}-{\widetilde{N}}^{i}({\bf X})p_{i})/{\widetilde{N}}({\bf
X})^{2}$, it follows that

$p_{0}=0.$ (46)

Hence, another interesting consequence is that

${\widetilde{N}}({\bf x}){\widetilde{\mathcal{H}}}({\bf
x})+{\widetilde{N}}^{i}({\bf x}){\widetilde{\mathcal{H}}}_{i}({\bf
x})=\delta({\bf x}-{\bf X}){\widetilde{H}}_{G}.$ (47)

Summarizing, the Hamiltonian dynamics is completely determined by the
constraints (42)-(45). Using (46), they can be simplified to (dropping the
tildes):

$\pi=0,\qquad\pi_{i}=0,$ (48) $\left[H_{G}-N^{i}({\bf
X})p_{i}\right]^{2}-N({\bf X})^{2}\left[h^{ij}({\bf
X})p_{i}p_{j}+m^{2}\right]=0,$ (49) ${\mathcal{H}}({\bf x})=\delta({\bf
x}-{\bf X})\frac{1}{N({\bf X})}\left[-N^{i}({\bf X})p_{i}+H_{G}\right],$ (50)
${\mathcal{H}}_{i}({\bf x})=\delta({\bf x}-{\bf X})p_{i}.$ (51)

Quantization in the Schrödinger picture leads to the following equations for
$\Psi({\bf X},h_{ij})$:444We have chosen the Laplace-Beltrami operator
ordering for the particle but not for gravity [19].

$\left[\textrm{i}N^{i}({\bf
X})\nabla_{i}+{\widehat{H}}_{G}\right]^{2}\Psi-N({\bf
X})^{2}\left(-\nabla^{2}+m^{2}\right)\Psi=0,$ (52)
${\widehat{\mathcal{H}}}({\bf x})\Psi=\delta({\bf x}-{\bf X})\frac{1}{N({\bf
X})}\left[\textrm{i}N^{i}({\bf X})\nabla_{i}+{\widehat{H}}_{G}\right]\Psi,$
(53) ${\widehat{\mathcal{H}}}_{i}({\bf x})\Psi=-\textrm{i}\delta({\bf x}-{\bf
X})\nabla_{i}\Psi,$ (54)

where $\nabla_{i}$ is the covariant derivative with respect to the metric
$h_{ij}({\bf X})$, with $\nabla_{i}\psi=\partial_{i}\psi$,
$\nabla^{2}=\nabla_{i}\nabla^{i}$ is the Laplacian, and

${\widehat{H}}_{G}=\int_{\Sigma}d^{3}x\left(N{\widehat{\mathcal{H}}}+N^{i}{\widehat{\mathcal{H}}}_{i}\right),$
(55) ${\widehat{\mathcal{H}}}=-\kappa G_{ijkl}\frac{\delta^{2}}{\delta
h_{ij}\delta
h_{kl}}+{\mathcal{V}}(h,\phi),\qquad{\widehat{\mathcal{H}}}_{i}=2\textrm{i}h_{ik}D_{j}\frac{\delta}{\delta
h_{jk}}.$ (56)

The wave function $\Psi$ does not depend on $N$, $N^{i}$ and $X^{0}$, because
of the operator constraints following from (46) and (48). The latter means
that the wave functional does not depend on time as is familiar in Wheeler-
DeWitt quantization.

The wave functional does not depend on $N$ and $N^{i}$, but they still appear
in the wave equations. Choosing $N=1$ and $N^{i}=0$, results in

$\left({\widehat{H}}_{G}^{2}+\nabla^{2}-m^{2}\right)\Psi=0,$ (57)
${\widehat{\mathcal{H}}}({\bf x})\Psi=\delta({\bf x}-{\bf
X}){\widehat{H}}_{G}\Psi,$ (58) ${\widehat{\mathcal{H}}}_{i}({\bf
x})\Psi=-\textrm{i}\delta({\bf x}-{\bf X})\nabla_{i}\Psi.$ (59)

Applying the quantization recipe using $t$ as temporal coordinate leads to the
following quantum theory [5, 6]:

${\widehat{\mathcal{H}}}({\bf x})\Psi=\pm\delta({\bf x}-{\bf
X})\sqrt{-\nabla^{2}+m^{2}}\Psi,\qquad{\widehat{\mathcal{H}}}_{i}({\bf
x})\Psi=-\textrm{i}\delta({\bf x}-{\bf X})\nabla_{i}\Psi,$ (60)

with $\Psi({\bf X},h_{ij})$. This theory also follows from (48)-(51) if the
square root is taken in (49). As in the case of electromagnetism, the
extension to many particles is straightforward, leading to

${\widehat{\mathcal{H}}}({\bf x})\Psi=\sum_{k}\delta({\bf x}-{\bf
X}_{k})\sqrt{-\nabla^{2}_{k}+m^{2}}\Psi,\qquad{\widehat{\mathcal{H}}}_{i}({\bf
x})\Psi=-\textrm{i}\sum_{k}\delta({\bf x}-{\bf X}_{k})\nabla_{ki}\Psi,$ (61)

with $k$ the particle label, $\Psi({\bf X}_{1},\dots,{\bf X}_{n},h_{ij})$. The
extension of the Klein-Gordon form (57)-(59) to many particles requires a bit
more effort (though apparently less than in the case of electromagnetism due
to absence of the time derivative). The natural generalization seems to be

$\left({\widehat{O}}^{2}_{k}+\nabla^{2}_{k}-m^{2}\right)\Psi=0,$ (62)
${\widehat{\mathcal{H}}}({\bf x})\Psi=\sum_{k}\delta({\bf x}-{\bf
X}_{k}){\widehat{O}}_{k}\Psi,$ (63) ${\widehat{\mathcal{H}}}_{i}({\bf
x})\Psi=-\textrm{i}\sum_{k}\delta({\bf x}-{\bf X}_{k})\nabla_{ki}\Psi.$ (64)

The operators ${\widehat{O}}_{k}$ can be determined almost completely by self-
consistency. Namely, for configurations $({\bf X}_{1},\dots,{\bf X}_{N})$ with
all the ${\bf X}_{k}$ different, we can integrate (63) over a test function
$f({\bf x})$ such that $f({\bf X}_{k})=1$ and $f({\bf X}_{l})=0$ for $l\neq
k$, to obtain ${\widehat{O}}_{k}\Psi=\int d^{3}xf({\bf
x}){\widehat{\mathcal{H}}}({\bf x})\Psi$. If not all the ${\bf X}_{k}$ are
different, the action of the ${\widehat{O}}_{k}$ could be defined similarly by
requiring symmetry. For example, for a configuration with ${\bf X}_{k}$ =
${\bf X}_{l}$, we can require that
${\widehat{O}}_{k}\Psi={\widehat{O}}_{l}\Psi$. This generalization is
consistent with (61), in the sense that (61) follows by taking particular
square roots in (62).

The Wheeler-DeWitt quantization of the relativistic particle was considered
before by Pavšič [7]. However, in the quantization procedure no common
temporal parameter was used. In the definition of the particle momentum the
parameter $\lambda$ was used, while in the definition of the field momenta the
parameter $t$ was used. As such the quantization deviates from the usual
recipe and the resulting quantum theory does not agree with the one presented
here. However, the equations obtained by Pavšič follow from ours when
interpreted in a partial Heisenberg picture (and as such is reminiscent of the
Dirac-Fock-Podolsky formulation of a particle interacting with an
electromagnetic field [20]). That is, defining

$\Phi(X,h_{ij})=\textrm{e}^{\textrm{i}{\widehat{H}}_{G}X^{0}}\Psi({\bf
X},h_{ij}),\qquad
h_{ij}(X)=\textrm{e}^{\textrm{i}{\widehat{H}}_{G}X^{0}}h_{ij}({\bf
X})\textrm{e}^{-\textrm{i}{\widehat{H}}_{G}X^{0}},$ (65)
${\widehat{\mathcal{H}}}({\bf
x},X^{0})=\textrm{e}^{\textrm{i}{\widehat{H}}_{G}X^{0}}{\widehat{\mathcal{H}}}({\bf
x})\textrm{e}^{-\textrm{i}{\widehat{H}}_{G}X^{0}},\qquad{\widehat{\mathcal{H}}}_{i}({\bf
x},X^{0})=\textrm{e}^{\textrm{i}{\widehat{H}}_{G}X^{0}}{\widehat{\mathcal{H}}}_{i}({\bf
x})\textrm{e}^{-\textrm{i}{\widehat{H}}_{G}X^{0}},$ (66)

the equations (57)-(59) reduce to

$\left(\partial^{2}_{0}-{\widetilde{\nabla}}^{2}+m^{2}\right)\Phi=0,$ (67)
${\widehat{\mathcal{H}}}({\bf x},X^{0})\Phi=-\textrm{i}\delta({\bf x}-{\bf
X})\partial_{0}\Phi,$ (68) ${\widehat{\mathcal{H}}}_{i}({\bf
x},X^{0})\Phi=-\textrm{i}\delta({\bf x}-{\bf X}){\widetilde{\nabla}}_{i}\Phi,$
(69)

where ${\widetilde{\nabla}}$ and ${\widetilde{\nabla}}^{2}$ concern the
covariant derivative with respect to the metric $h_{ij}(X)$. (Pavšič also
chooses a different operator ordering, but that difference is not essential.)
For Pavšič the appearance of the time derivative in (68) is interesting the
light of the problem of time. However, the proper quantum theory is time-
independent.

##  4 Acknowledgments

This work is supported by the Research Foundation Flanders (Fonds
Wetenschappelijk Onderzoek, FWO), Grant No. G066918N. It is a pleasure to
thank Christian Maes and Kasper Meerts for useful discussions.

## References

  * [1] C. Kiefer, Quantum Gravity, International Series of Monographs on Physics 124, Clarendon Press, Oxford (2004).
  * [2] K.V. Kuchar̆, “Time and interpretations of quantum gravity”, in Proceedings of the 4th Canadian Conference on General Relativity and Relativistic Astrophysics, eds. G. Kunstatter, D. Vincent and J. Williams, World Scientific, Singapore (1992), reprinted in Int. J. Mod. Phys. D 20, 3-86 (2011).
  * [3] C. Rovelli, Quantum Gravity, Cambridge University Press, Cambridge (2004).
  * [4] N. Pinto-Neto and W. Struyve, “Bohmian quantum gravity and cosmology”, in Applied Bohmian Mechanics: From Nanoscale Systems to Cosmology, 2nd edition, eds. X. Oriols and J. Mompart, 607-656 (2019) and arXiv:1801.03353 [gr-qc].
  * [5] C. Rovelli, “What is observable in classical and quantum gravity?”, Class. Quantum Grav. 8, 297-316 (1991).
  * [6] C. Rovelli, “Quantum reference systems”, Class. Quantum Grav. 8, 317-331 (1991).
  * [7] M. Pavšič, “Klein–Gordon–Wheeler–DeWitt–Schrödinger equation”, Phys. Lett. B 703, 614-619 (2011) and arXiv:1106.2017 [gr-qc].
  * [8] M. Pavšič, “A novel approach to quantum gravity in the presence of matter without the problem of time”, Int. J. Mod. Phys. A 35, 2050015 (2020) and arXiv:1906.03987 [physics.gen-ph].
  * [9] J.D. Brown and K.V. Kuchař, “Dust as a standard of space and time in canonical quantum gravity”, Phys. Rev. D 51, 5600 (1995) and arXiv:gr-qc/9409001.
  * [10] P.A.M. Dirac, Lectures On quantum Mechanics, Belfer Graduate School of Science, Yeshiva University, New York (1964), reprinted by Dover Publications, New York (2001).
  * [11] A. Hanson, T. Regge and C. Teitelboim, Constrained Hamiltonian Systems, Accademia Nazionale Dei Lincei, Roma (1976).
  * [12] K. Sundermeyer, Constrained Dynamics, Lecture Notes in Physics 169, Springer-Verlag, Berlin (1982).
  * [13] D.M. Gitman and I.V. Tyutin, Quantization of Fields with Constraints, Springer-Verlag, Berlin Heidelberg (1990).
  * [14] M. Henneaux and C. Teitelboim, Quantization of Gauge Systems, Princeton University Press, New Jersey (1991).
  * [15] N. Kemmer, “The particle aspect of meson theory”, Proc. R. Soc. A 173, 91-116 (1939).
  * [16] W. Struyve, W. De Baere, J. De Neve and S. De Weirdt, “On the uniqueness of paths for spin-0 and spin-1 quantum mechanics”, Phys. Lett. A 322, 84-95 (2004) and arXiv:quant-ph/0311098.
  * [17] S. Petrat and R. Tumulka, “Multi-time equations, classical and quantum”, Proc. R. Soc. A 470, 20130632 (2014) and arXiv:1309.1103 [quant-ph].
  * [18] S. Carlip, General Relativity, Oxford University Press, Oxford (2019).
  * [19] T. Christodoulakis and J. Zanelli, “Operator Ordering in Quantum Mechanics and Quantum Gravity”, Nuovo Cimento B 93, 1-21 (1986).
  * [20] P.A.M. Dirac, V.A. Fock and B. Podolsky, “On quantum electrodynamics”, Physikalische Zeitschrift der Sowjetunion 2, 468–479 (1932); reprinted in ‘Selected Papers on Quantum Electrodynamics”, ed. J. Schwinger, Dover, New York (1958).