A Sparkler in the Fireworks Galaxy: Discovery of an Ultraluminous X-Ray Transient with a Strong Oxygen Line in NGC 6946

Chandra observed NGC 6946 nine times between 2001 and 2017 (in most cases, to follow supernova explosions, which occur quite frequently in this galaxy). The observation log is shown in Table 1. All observations were made with the S3 chip of the Advanced CCD Imaging Spectrometer (ACIS) array, except for the ACIS observation of 2012 (ObsID 13435), in which the target was placed on the S2 chip. After downloading the data from the public archive, we reprocessed and analyzed them with the Chandra Interactive Analysis of Observations (ciao) software version 4.10 (Fruscione et al. 2006) with calibration database version 4.7.9. Specifically, we rebuilt level 2 event files with the task chandra_repro and filtered out background flares with the task deflare. We created images in multiple energy bands for each epoch with dmcopy. We searched for point sources in each epoch with with wavdetect. Our target was undetected at every epoch until its first appearance in 2012 May. It was also later detected in the subsequent Chandra observations of 2016 September and 2017 June.

Table 1.  Chandra, XMM-Newton, and Swift Observations of NGC 6946 and Luminosity of CXOU J203451.1+601043

ObsID	Good Time Interval	Observation Date	Observed X-Ray Fluxa	X-Ray Luminositya
	(ks)		(erg cm−2 s−1)	(erg s−1)
Chandra/ACIS
1043	58.3	2001 Sep 7	<1.2 × 10−15	<2.1 × 1037
4404	28.7	2002 Nov 25	<1.4 × 10−15	<2.5 × 1037
4631	28.4	2004 Oct 22	<1.1 × 10−15	<2.0 × 1037
4632	25.2	2004 Nov 6
4633	26.6	2004 Dec 3
13435	20.4	2012 May 21	${4.0}_{-0.6}^{+0.6}\times {10}^{-14}$	${1.1}_{-0.2}^{+0.2}\times {10}^{39}$
17878	40.0	2016 Sep 28	${2.7}_{-0.2}^{+0.2}\times {10}^{-13}$	${4.8}_{-0.4}^{+0.4}\times {10}^{39}$
19887	18.5	2016 Sep 28
19040	9.8	2017 Jun 11	${0.6}_{-0.1}^{+0.1}\times {10}^{-14}$	${1.7}_{-0.3}^{+0.3}\times {10}^{39}$
XMM-Newton/EPIC
0093641501	0.6	2003 Apr 18	$\lt 1\times {10}^{-14}$	$\lt 2\times {10}^{38}$
0093641601	2.2	2003 May 17
0093641701	1.2	2003 Jun 18
0200670101	3.9	2004 Jun 9	<3.1 × 10−15	<3 × 1037
0200670201	12.7	2004 Jun 11
0200670301	11.3	2004 Jun 13
0200670401	8.8	2004 Jun 25
0401360101	18.7	2006 May 23	<2 × 10−15	<1.5 × 1037
0401360201	4.7	2006 Jun 2
0401360301	4.9	2006 Jun 18
0500730101	26.0	2007 Nov 8	<1.2 × 10−15	<2 × 1037
0500730201	31.7	2007 Nov 2
0691570101	109.3	2012 Oct 21	${1.6}_{-0.1}^{+0.1}\times {10}^{-13}$	${1.6}_{-0.1}^{+0.1}\times {10}^{39}$
0794581201	43.1	2017 Jun 1	${0.6}_{-0.1}^{+0.1}\times {10}^{-13}$	${2.8}_{-0.1}^{+2.0}\times {10}^{39}$
Swift/XRT
31113001–31113004	10	2008 Feb 4–2008 Feb 14	$\lt 2\times {10}^{-14}$	< 4 × 1038
49820001–49820003	7	2013 May 31–2013 Jun 4	${1.4}_{-0.4}^{+0.4}\times {10}^{-13}$	${2.4}_{-0.5}^{+0.5}\times {10}^{39}$
10130001–10130029	44	2017 May 13–2017 Sep 17	${1.1}_{-0.2}^{+0.2}\times {10}^{-13}$	${3.6}_{-0.4}^{+0.4}\times {10}^{39}$
94059001–94059044	43	2018 Apr 1–2018 Dec 27	${2.2}_{-0.3}^{+0.3}\times {10}^{-13}$	${2.8}_{-0.4}^{+0.4}\times {10}^{39}$
94059045–94059072	17	2019 Jan 6–2019 Apr 6	${1.7}_{-0.3}^{+0.3}\times {10}^{-13}$	${3.2}_{-0.4}^{+0.4}\times {10}^{39}$

Note.

^aFor the Chandra observations, we estimated observed fluxes and intrinsic luminosities (or their respective upper limits) in the 0.3–10 keV band with the ciao task srcflux, assuming a power-law model with photon index ${\rm{\Gamma }}=2.5$ and total column density ${N}_{{\rm{H}}}=4\times {10}^{21}$ cm⁻² (twice the Galactic line-of-sight value). This model was chosen because it approximates the best-fitting power-law model for the 2016 Chandra spectrum (Table 2). For XMM-Newton observations, we also used a simple power-law model to convert from count rates to fluxes and luminosities in this table. When we had enough counts for a detailed fit (2012 and 2017 XMM-Newton observations), we used the best-fitting values (Table 2); for the nondetections, we assumed Γ = 2.5 and N_H = 4 × 10²¹ cm⁻². For Swift, we fixed N_H = 4 × 10²¹ cm⁻², then used pimms to estimate the power-law photon index that best approximates the observed (1.5–10)/(0.3–1.5) hardness ratio of each observation. We used those indices to infer fluxes and luminosities of the respective observations; when not enough counts were available, we again assumed Γ = 2.5. See Table 2 for a comparison of flux and luminosity conversions using power laws versus more complex spectral models.

Download table as:  ASCIITypeset image

In order to check whether the source is spatially extended, we computed the point-spread functions (PSFs) at the off-axis locations of the source in 2012 and 2016 using the Chandra Ray Tracer (ChaRT⁵ ) and simulated PSF files with the marx software.⁶ We compared the profiles of the source with those of the respective simulated PSFs with the ciao tool srcextent. We confirm that the transient source is consistent with being pointlike. This is further confirmed by the 2017 observation (ObsID 19040), the only one in which the target is almost on-axis: its observed FWHM of ≈2'' is consistent with the expected width of an on-axis PSF.

For each observation (both those with a detection and those with no detection), we extracted the source events from a circular aperture of either 5'' or 2'' radius (for ObsID 19040) to match the size of the PSF at the respective locations. Background events were extracted from nearby source-free regions at least three times the size of the source region. We used the ciao task scrflux to estimate an upper limit to the source flux in the epochs when it was not detected and an approximate flux in the two epochs (2012 and 2017) when it was detected but with a small number of counts. For the 2016 data set, we combined the spectra from the two exposures taken on 2016 September 28 and obtained enough counts for meaningful spectral fitting; we extracted source and background spectra with specextract, which also generates the appropriate auxiliary response files (ARFs) and response matrix files (RMFs) for the subsequent spectral analysis. Spectra were grouped to a minimum of 15 counts bin^–1 for χ² fitting. We also repeated our spectral analysis on the same spectra grouped to 1 count bin^–1 using the Cash statistics (Cash 1979). Moreover, to check for short-term variability in each observation in which the source was detected, we extracted background-subtracted light curves using the ciao tool dmextract.

We did subsequent data analysis with NASA's High Energy Astrophysics Software (HEASOFT): ds9 (Joye & Mandel 2003) version 8.0 for imaging and photometry, ftools/Xronos (Blackburn 1995) version 6.25 for timing analysis, and xspec (Arnaud 1996) version 12.9.1 for spectral modeling. The reported errors are 90% confidence intervals for the fitting parameters. For the 2016 Chandra spectra, all of the best-fitting parameters obtained from χ² and Cash statistics fitting (cstat in xspec) agree well within their error ranges; thus, in Section 3.3 and Table 2, we report only the values from χ² fitting for simplicity.

Table 2.  Best-fitting Parameters of the EPIC Spectra from 2012 and 2017 and the ACIS Spectrum from 2016

Model Parameters	Values
	2012	2016	2017
tbabs × tbabs × (po + Gaussian)
NH,Gal (1022 cm−2)	[0.20]	[0.20]	[0.20]
NH,int (1022 cm−2)	${0.14}_{-0.03}^{+0.03}$	${0.20}_{-0.12}^{+0.13}$	${0.23}_{-0.06}^{+0.07}$
Γ (keV)	${2.17}_{-0.06}^{+0.07}$	${2.63}_{-0.20}^{+0.22}$	${3.50}_{-0.26}^{+0.31}$
${N}_{\mathrm{po}}$ (10−5 ph keV−1 cm−2 s−1 at 1 keV)	${5.8}_{-0.4}^{+0.4}$	${14.2}_{-2.9}^{+3.9}$	${6.3}_{-1.1}^{+1.4}$
${E}_{\mathrm{line}}$ (keV)	${0.61}_{-0.02}^{+0.02}$	⋯	${0.66}_{-0.01}^{+0.01}$
${\sigma }_{\mathrm{line}}$ (keV)	[0]	⋯	$\lt 0.030$
${N}_{\mathrm{line}}$ (10−5 ph cm−2 s−1)	${1.2}_{-0.7}^{+1.2}$	⋯	${2.7}_{-1.4}^{+2.3}$
${\chi }^{2}/$dof	346.8/230 (1.51)	68.0/61 (1.12)	89.3/64 (1.40)
${f}_{0.3-10}$ (10−13 erg cm−2 s−1)a	${1.68}_{-0.06}^{+0.06}$	${2.49}_{-0.20}^{+0.22}$	${0.62}_{-0.04}^{+0.04}$
${L}_{0.3-10}$ (1039 erg cm−2 s−1)b	${2.23}_{-0.14}^{+0.16}$	${4.86}_{-0.97}^{+1.50}$	${3.10}_{-0.84}^{+1.45}$
tbabs × tbabs × (bbodyrad + diskbb + Gaussian)
${N}_{{\rm{H}},\mathrm{Gal}}$ (1022 cm−2)	[0.20]	[0.20]	[0.20]
${N}_{{\rm{H}},\mathrm{int}}$ (1022 cm−2)	${0.08}_{-0.08}^{+0.12}$	$\lt 0.15$	${0.26}_{-0.18}^{+0.23}$
${{kT}}_{\mathrm{bb}}$ (keV)	${0.15}_{-0.02}^{+0.03}$	${0.27}_{-0.09}^{+0.08}$	${0.15}_{-0.02}^{+0.03}$
${N}_{\mathrm{bb}}$ (km2)c	${11.3}_{-9.2}^{+41.5}$	${1.7}_{-1.0}^{+8.7}$	${89}_{-79}^{+908}$
${{kT}}_{\mathrm{in}}$ (keV)	${1.25}_{-0.07}^{+0.07}$	${1.14}_{-0.20}^{+0.40}$	${0.72}_{-0.11}^{+0.14}$
${N}_{\mathrm{dbb}}$ (10−3 km2)d	${3.7}_{-0.9}^{+1.1}$	${6.3}_{-4.8}^{+10.0}$	${13.3}_{-8.1}^{+19.7}$
${E}_{\mathrm{line}}$ (keV)	${0.61}_{-0.04}^{+0.04}$	⋯	${0.66}_{-0.02}^{+0.02}$
${\sigma }_{\mathrm{line}}$ (keV)	[0]	⋯	$\lt 0.033$
${N}_{\mathrm{line}}$ (10−5 ph cm−2 s−1)	${0.48}_{-0.26}^{+1.44}$	⋯	${3.2}_{-2.1}^{+13.5}$
${\chi }^{2}/$dof	$235.5/228$ (1.03)	$68.3/59$ (1.16)	$72.5/62$ (1.17)
${f}_{0.3-10}$ (10−13 erg cm−2 s−1)a	${1.57}_{-0.05}^{+0.05}$	${2.28}_{-0.10}^{+0.11}$	${0.61}_{-0.04}^{+0.04}$
${L}_{0.3-10}$ (1039 erg cm−2 s−1)b	${1.75}_{-0.30}^{+0.76}$	${2.22}_{-0.17}^{+0.76}$	${2.52}_{-1.44}^{+6.62}$
tbabs × tbabs × (bbodyrad + comptt + Gaussian)
${N}_{{\rm{H}},\mathrm{Gal}}$ (1022 cm−2)	[0.20]	[0.20]	[0.20]
${N}_{{\rm{H}},\mathrm{int}}$ (1022 cm−2)	${0.05}_{-0.01}^{+0.01}$	${0.23}_{-0.10}^{+0.33}$	${0.28}_{-0.02}^{+0.03}$
${{kT}}_{\mathrm{bb}}$ (keV)	${0.16}_{-0.01}^{+0.01}$	${0.07}_{-0.07}^{+0.06}$	${0.13}_{-0.01}^{+0.06}$
${N}_{\mathrm{bb}}$ (km2)c	${7.3}_{-1.0}^{+13.7}$	(unconstrained)	${117}_{-12}^{+258}$
kT0 (keV)=${{kT}}_{\mathrm{bb}}$	$[{0.16}_{-0.01}^{+0.01}]$	$[{0.07}_{-0.07}^{+0.06}]$	$[{0.13}_{-0.01}^{+0.06}]$
kTe (keV)	${1.01}_{-0.10}^{+0.15}$	$\gt 0.98$	${0.65}_{* }^{+0.89}$
τ	${13.7}_{-0.3}^{+0.4}$	${5.0}_{-5.0}^{+3.4}$	${13.6}_{-0.8}^{+7.3}$
${N}_{{\rm{c}}}$ (10−5)	${9.7}_{-0.2}^{+0.2}$	${63}_{* }^{+606}$	${11.9}_{-1.1}^{+2.4}$
${E}_{\mathrm{line}}$ (keV)	${0.61}_{-0.04}^{+0.03}$	−–	${0.66}_{-0.02}^{+0.01}$
${\sigma }_{\mathrm{line}}$ (keV)	[0]	⋯	$\lt 0.025$
${N}_{\mathrm{line}}$ (10−5 ph cm−2 s−1)	${0.36}_{-0.27}^{+0.27}$	⋯	${3.7}_{-1.4}^{+1.4}$
${\chi }^{2}/$dof	$234.6/227$ (1.03)	$64.1/58$ (1.11)	$72.6/61$ (1.19)
${f}_{0.3-10}$ (10−13 erg cm−2 s−1)a	${1.56}_{-0.05}^{+0.05}$	${2.59}_{-0.32}^{+1.85}$	${0.61}_{-0.04}^{+0.04}$
${L}_{0.3-10}$ (1039 erg cm−2 s−1)b	${1.63}_{-0.04}^{+0.05}$	($\gtrsim 3.8$)e	${2.83}_{-0.08}^{+2.03}$
tbabs × tbabs × (diskpbb + Gaussian)
NH,Gal (1022 cm−2)	[0.20]	[0.20]	[0.20]
NH,int (1022 cm−2)	< 0.02	${0.03}_{-0.03}^{+0.06}$	${0.03}_{-0.02}^{+0.18}$
kTin (keV)	${1.49}_{-0.11}^{+0.13}$	${1.40}_{-0.31}^{+0.29}$	${0.68}_{-0.09}^{+0.11}$
p	${0.60}_{-0.02}^{+0.01}$	${0.50}_{* }^{+0.04}$	${0.50}_{* }^{+0.04}$
Ndpbb (10−3 km2)d	${0.92}_{-0.31}^{+0.44}$	${0.89}_{-0.53}^{+2.17}$	${6.7}_{-3.4}^{+8.1}$
${E}_{\mathrm{line}}$ (keV)	${0.61}_{-0.04}^{+0.04}$	⋯	${0.67}_{-0.02}^{+0.01}$
σline (keV)	[0]	⋯	$\lt 0.029$
Nline (10−5 ph cm−2 s−1)	${0.28}_{-0.19}^{+0.23}$	⋯	${0.80}_{-0.36}^{+0.45}$
χ2/dof	252.7 /229 (1.10)	65.9/60 (1.10)	102.5/63 (1.63)
f0.3−10 (10−13erg cm−2 s−1)a	${1.58}_{-0.05}^{+0.06}$	${2.41}_{-0.20}^{+0.19}$	${0.59}_{-0.04}^{+0.04}$
L0.3−10 (1039 erg cm−2 s−1)b	${1.50}_{-0.05}^{+0.06}$	${1.71}_{-0.14}^{+0.13}$	${0.93}_{-0.06}^{+0.06}$

Notes.

^aObserved fluxes in the 0.3–10 keV band. ^bFor all spectral models, the de-absorbed luminosities L_0.3−10 (0.3–10 keV band) were defined as 4πd² times the de-absorbed fluxes. ^cN_bb = (R_bb/D₁₀)², where R_bb is the source radius in km and D₁₀ is the distance to the source in units of 10 kpc (here D₁₀ = 770). ^dN_dbb = (R_in/D₁₀)² cos θ, where R_in is the apparent inner disk radius in km, D₁₀ is the distance to the source in units of 10 kpc, and θ is our viewing angle (θ = 0 is face-on). N_dpbb is defined exactly as N_dbb. ^eUpper limit unconstrained because of the degeneracy between absorption column density and seed blackbody luminosity.

Download table as:  ASCIITypeset images: 1 2

There are 14 XMM-Newton observations of NGC 6946 with the European Photon Imaging Camera (EPIC) in the full-frame mode. We downloaded the data from the XMM-Newton Science Archive and reduced them with the Science Analysis Software (sas) version 17.0. The EPIC-pn and EPIC-MOS events were processed using the XMM-Newton pipeline and associated calibration files. We removed time intervals of high background from the event files with evselect using a "RATE ≤ 0.4" threshold for the pn and "RATE ≤ 0.35" for MOS1 and MOS2.

Between 2003 and 2007, CXOU J203451.1+601043 was not detected in any of the observations. We grouped the nondetection observations by year to increase the signal-to-noise ratio. For each year, we built a stacked pn and MOS image and measured the total counts inside a circle with 20'' radius at the position of the source and the background counts from surrounding regions. We then applied the Bayesian method of Kraft et al. (1991) for Poisson-distributed counts to obtain the 90% upper limit to the net counts and count rates.⁷ For the conversion from count rates to fluxes, we used the EPIC exposure maps to account for vignetting, and we assumed that for an equal effective exposure time between the three instruments, the pn contributes to about 61% of the count rate and MOS1+MOS2 to about 39%. We obtained this relative ratio between the EPIC instruments with the Portable, Interactive Multi-Mission Simulator (pimms) version 4.9 for a range of plausible spectral models; the relative contribution of pn and MOS changes at most by 2% or 3% from model to model, which is negligible for the purpose of our analysis. We converted the upper limits on the pn+MOS count rates to flux upper limits with pimms using a power-law model with photon index Γ = 2.5 and total absorbing column density N_H = 4 × 10²¹ cm⁻² (Table 1).

In contrast, CXOU J203451.1+601043 was detected in 2012 and 2017 with luminosities exceeding 10³⁹ ergs s⁻¹. For the 2012 and 2017 detections, source events were extracted from a circular aperture of 20'' radius. We generated background-subtracted light curves binned to 0.1 s with the sas task evselect followed by epiclccorr. Background events were extracted from point source–free, circular regions on the same chip approximately three times larger than the source region. We created spectral files and associated instrumental responses with the script multiespecget, which extracts the spectra of all three detectors and combines them into a single EPIC spectrum. We chose to combine the three EPIC spectra to increase the signal-to-noise ratio of possible line features in the soft X-ray band. For spectral extraction, we used the standard filtering conditions "(FLAG==0) && (PATTERN ≤ 4)" for the pn and "(#XMMEA_EM && (PATTERN ≤ 12)" for MOSs. Spectra were grouped to a minimum of 25 counts bin^–1 for χ² fitting.

As for the Chandra data, subsequent imaging, timing, and spectral analysis were carried out with HEASOFT packages (ds9, ftools, and xspec) and astropy (Astropy Collaboration et al. 2018) specifically for period searches (Section 3.2).

The X-ray Telescope (XRT) on board Swift has monitored NGC 6946 frequently over the years, but typical exposure times are very short (∼1 ks). For this work, we used all of the XRT observations from 2008 onward, stacked into data sets for individual years (Table 3). The 2008 observations are the last ones from any X-ray observatory in which the transient ULX is not detected; it is detected in all subsequent Swift observations. We used standard HEASOFT packages (version 6.25) for spectral extraction; we created a combined spectrum and exposure map with xselect and an ancillary response function with xrtmkarf, and the ready-made response file comes from the XRT Calibration Database.⁸

Table 3.  Alternative Set of Models for the 2017 EPIC Spectra

Model Parameters	Values
tbabs × tbabs × (bbodyrad + diskbb + vapec)
NH,Gal (1022 cm−2)	[0.20]
NH,int (1022 cm−2)	${0.13}_{-0.13}^{+0.19}$
${{kT}}_{\mathrm{bb}}$ (keV)	${0.14}_{-0.04}^{+0.06}$
${N}_{\mathrm{bb}}$ (km2)a	${22}_{-21}^{+358}$
${{kT}}_{\mathrm{in}}$ (keV)	${0.72}_{-0.12}^{+0.14}$
${N}_{\mathrm{dbb}}$ (10−3 km2)b	${7.5}_{-4.3}^{+9.5}$
kTvapec (keV)	${0.41}_{-0.10}^{+0.09}$
O abundance (solar units)	>10
Ne = Mg = Si abundances (solar units)	>10
Other abundances (solar units)	[1]
Nvapec (10−6)c	${0.07}_{-0.03}^{+4.3}$
χ2/dof	$68.6/61$ (1.13)
${f}_{0.3-10}$ (10−13 erg cm−2 s−1)d	${0.62}_{-0.04}^{+0.04}$
${L}_{0.3-10}$ (1039 erg cm−2 s−1)e	${1.32}_{-0.57}^{+2.51}$
tbabs × tbabs × (bbodyrad + comptt + vapec)
${N}_{{\rm{H}},\mathrm{Gal}}$ (1022 cm−2)	[0.20]
${N}_{{\rm{H}},\mathrm{int}}$ (1022 cm−2)	${0.12}_{-0.02}^{+0.16}$
${{kT}}_{\mathrm{bb}}$ (keV)	${0.15}_{-0.02}^{+0.01}$
${N}_{\mathrm{bb}}$ (km2)c	${17.5}_{-2.8}^{+39.9}$
kT0 (keV)=${{kT}}_{\mathrm{bb}}$	$[{0.15}_{-0.02}^{+0.01}]$
kTe (keV)	${0.59}_{* }^{+0.46}$
τ	${26.4}_{-3.1}^{+4.1}$
${N}_{{\rm{c}}}$ (10−5)	${4.3}_{-0.5}^{+0.7}$
${{kT}}_{\mathrm{vapec}}$ (keV)	${0.43}_{-0.04}^{+0.05}$
O abundance (solar units)	${385}_{-140}^{+155}$
Ne = Mg = Si abundances (solar units)	${230}_{-95}^{+110}$
Other abundances (solar units)	[1]
${N}_{\mathrm{vapec}}$ (10−6)c	${0.18}_{-0.05}^{+0.65}$
${\chi }^{2}/$dof	$68.4/60$ (1.14)
${f}_{0.3-10}$ (10−13 erg cm−2 s−1)d	${0.61}_{-0.04}^{+0.04}$
${L}_{0.3-10}$ (1039 erg cm−2 s−1)e	${1.23}_{-0.07}^{+2.71}$
tbabs × tbabs × (diskpbb + vapec)
${N}_{{\rm{H}},\mathrm{Gal}}$ (1022 cm−2)	[0.20]
${N}_{{\rm{H}},\mathrm{int}}$ (1022 cm−2)	${0.05}_{-0.05}^{+0.06}$
${{kT}}_{\mathrm{in}}$ (keV)	${0.93}_{-0.16}^{+0.23}$
p	${0.50}_{* }^{+0.12}$
${N}_{\mathrm{dpbb}}$ (10−3 km2)b	${1.2}_{-0.6}^{+6.7}$
${{kT}}_{\mathrm{vapec}}$ (keV)	${0.40}_{-0.04}^{+0.07}$
O abundance (solar units)	>5.5
Ne = Mg = Si abundances (solar units)	>6.3
Other abundances (solar units)	[1]
${N}_{\mathrm{vapec}}$ (10−6)c	${2.8}_{-2.7}^{+4.5}$
${\chi }^{2}/$ dof	$71.7/62$ (1.16)
${f}_{0.3-10}$ (10−13 erg cm−2 s−1)d	${0.62}_{-0.04}^{+0.04}$
${L}_{0.3-10}$ (1039 erg cm−2 s−1)e	${0.95}_{-0.20}^{+0.23}$

Notes.

^a ${N}_{\mathrm{bb}}={({R}_{\mathrm{bb}}/{D}_{10})}^{2}$, where ${R}_{\mathrm{bb}}$ is the source radius in km and D₁₀ is the distance to the source in units of 10 kpc (here ${D}_{10}=770$). ^b ${N}_{\mathrm{dbb}}={({R}_{\mathrm{in}}/{D}_{10})}^{2}\cos \theta$, where R_in is the apparent inner disk radius in km, D₁₀ is the distance to the source in units of 10 kpc, and θ is our viewing angle ($\theta =0$ is face-on). N_dpbb is defined exactly as N_dbb. ^c ${N}_{\mathrm{vapec}}=\displaystyle \frac{{10}^{-14}}{4\pi \,{d}^{2}}\int {n}_{e}{n}_{{\rm{H}}}{dV}$, where d is the angular diameter distance to the source (cm), and n_e and n_H are the electron and hydrogen densities (cm⁻³). ^dObserved fluxes in the 0.3–10 keV band. ^eFor all spectral models, the de-absorbed luminosities ${L}_{0.3-10}$ (0.3–10 keV band) were defined as 4πd² times the de-absorbed fluxes

Download table as:  ASCIITypeset image

We also searched for NuSTAR observations that covered the position of the transient and found two: one from 2017 May (ObsID 90302004002; 66.8 ks) and the other from 2017 June (ObsID 90302004004; 47.8 ks). We examined the focal plane module A and focal plane module B data, processed by the pipeline task nupipeline; however, we did not detect any source at the ULX position, either in the individual images or in a stack of the two images.

The field containing our target source was observed (Table 4) with the Advanced Camera for Surveys (ACS) Wide Field Channel (WFC) on 2004 July 29 with the F814W filter (exposure time of 120 s). It was then reobserved with ACS-WFC on 2016 October 26 in the F606W and F814W filters (exposure times of 2430 and 2570 s, respectively). It was later imaged with the Wide Field Camera 3 (WFC3) Ultraviolet and VISible light camera (UVIS) on 2018 January 5 in the F555W band for 710 s and in F814W for 780 s.

We retrieved calibrated, geometrically corrected images (.drc files) from the Mikulski Archive for Space Telescopes. We used ds9 to inspect the images and perform aperture photometry of the candidate counterparts, with a source extraction radius of 015 and a local annular background region. We then converted these small-aperture measurements to infinite-aperture values with the help of the online tables of encircled energy fractions for ACS-WFC and WFC3-UVIS. Finally, we applied the corresponding zero-points to convert our photometric measurements into magnitudes in the Vega system.

The first detection of this transient as a bright new X-ray source was in the Chandra observations of 2012 May; its luminosity was 2 orders of magnitude higher than typical previous nondetection limits (Figure 1 and Table 3). It has remained very luminous in all subsequent Chandra, XMM-Newton, and Swift observations, including the most recent ones (Figure 2). The transient resides in one of the spiral arms (Figure 3); the projected galactocentric radius is ≈90'', which is ≈3.4 kpc at the assumed distance of 7.7 Mpc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At first, we improved the astrometry of the stacked Chandra/ACIS images with the help of a small sample of sources that have a counterpart in the Gaia and 2MASS catalogs (available in ds9). We estimate that our transient source is located at R.A. (J2000) $=\,{20}^{{\rm{h}}}{34}^{{\rm{m}}}51\buildrel{\rm{s}}\over{.} 1$, decl. (J2000) = 60°10'436, but the uncertainty remained large, ≈05. The reason for this large positional error is caused by the off-axis location of the X-ray source (hence, a distorted PSF) in three of the four observations in which it is detected; the only Chandra observation in which the ULX is almost on-axis is also short (ObsID 19040; 10 ks), and there are no direct X-ray/Gaia or X-ray/2MASS associations for that short exposure alone.

Therefore, we used a different method in our search for an optical counterpart. In the on-axis observation 19040, there are two bright X-ray sources within 2' of the transient ULX, with a well-identified, pointlike optical counterpart in the HST images. One is SN 2017eaw (Wiggins 2017), located at R.A. (J2000) = 20^h34^m4424, decl. (J2000) = 60°11'359; the other is the ULX inside the MF16 nebula (Roberts & Colbert 2003).⁹ We determined the centroids of the X-ray emission from the three sources with wavdetect applied to the ACIS image from ObsID 19040. The uncertainty in the centroid position is ≲01 for each of the three sources. The transient ULX is located 514 east and 526 south of SN 2017eaw. It is also located 718 west and 473 south of the MF16 ULX.

Using those relative offsets in the HST images, we constrained the error circle of the transient ULX. The SN 2017eaw (but not MF16) is in the field of view of the WFC3 image from 2018; its relative offset gives the transient ULX location in those images. In the 2004 and 2016 ACS images, both MF16 and the position of SN 2017eaw are in the field of view¹⁰ ; this gives us two reference offsets for the relative position of the transient ULX. Both offsets point to the same location with a difference of <01 between them. In summary, we constrained the relative position of the transient ULX on the HST images with an error radius of ≲02.

Finally, we also refined the absolute astrometry of the HST images based on Gaia and Sloan Digital Sky Survey associations. This was a much simpler and straightforward task and reduced the uncertainty on the absolute astrometry of the HST images to ≲01. The most accurate position for the transient ULX is then R.A. (J2000) = 20^h34^m5112, decl. (J2000) = 60°10'433 (±02).

First, we studied the long-term X-ray variability of CXOU J203451.1+601043. We determined the net count rate or 90% upper limit in each Chandra/ACIS observation in the 0.3–7 keV band. In some cases, we stacked observations taken a few weeks apart to reach a deeper detection limit. We converted 0.3–7 keV count rates to 0.3–10 keV fluxes assuming the same model for all Chandra observations: a power law with photon index Γ = 2.5 and intrinsic neutral absorption column density N_H = 2 × 10²¹ cm⁻², in addition to the line-of-sight Galactic N_H = 2 × 10²¹ cm⁻². We chose these model parameters because they are a good approximation to those found from a detailed fit to the 2016 Chandra data (Table 2 and Section 3.3). If we use a "standard" photon index Γ = 1.7 and only line-of-sight absorption, the inferred luminosities or upper limits will be ≈75% of those reported in Table 3.

We did a similar analysis for the XMM-Newton/EPIC observations, stacking the exposures from individual years to increase the signal-to-noise ratio. Count rates were extracted from the 0.3–10 keV band. For the 2012 and 2017 observations, we had enough counts to fit multicomponent models to the data (as we shall discuss in Section 3.3). We adopted the best-fitting Comptonization model to convert from net count rates to the fluxes and luminosities listed in Table 3. For the XMM-Newton observations in which the source was not detected, we determined upper count-rate limits from the combined EPIC images, and we used our fiducial power-law model (Γ = 2.5 and total column density N_H = 4 × 10²¹ cm⁻²) to convert to flux and luminosity limits.

For the stacked Swift/XRT observations from 2008 (total of ≈10 ks between February 4 and 14) and 2013 (total of ≈7 ks between May 31 and June 4), we used our fiducial power-law model (Γ = 2.5, N_H = 4 × 10²¹ cm⁻²) to determine fluxes and luminosities or their upper limits. For the other three sets of stacked Swift observations (44 ks in 2017, 43 ks in 2018, and 17 ks so far in 2019), we had enough counts to determine the hardness ratios between the 1.5 and 10 keV band and the 0.3–1.5 keV band. We fixed the total column density to N_H = 4 × 10²¹ cm⁻² and used pimms to estimate the power-law photon indices that most closely reproduce the observed hardness ratios; the values are Γ = 3.2 ± 0.5 in 2017, Γ = 2.1 ± 0.4 in 2018, and Γ = 2.6 ± 0.5 in 2019. We then used those photon indices to convert from net count rates to fluxes and luminosities in the respective observations.

The resulting long-term light curve is shown in Figure 4. In all observations between 2001 and 2008, CXOU J203451.1+601043 was always undetected, and it was always detected from 2012 to 2019. The upper limit to the nondetections is typically ≈few × 10³⁷ erg s⁻¹ for individual years and <10³⁷ erg s⁻¹ if all of the observations with nondetections are stacked up. We also know that the source was undetected in the ROSAT High Resolution Imager in 1994 May (Schlegel et al. 2000), with an upper limit to the de-absorbed 0.3–10 keV luminosity of ≈10³⁸ erg s⁻¹ (after converting from the model and distance used in that paper to those used for this work). Since 2012, the source has been hovering at luminosities ≈1–4 × 10³⁹ erg s⁻¹ (depending on the choice of spectral model).

Zoom In Zoom Out Reset image size

We also determined and examined the background-subtracted light curves from individual Chandra and XMM-Newton observations with detections; for this, we used the ftools tasks lcurve and lcstats. We found that the 2016 Chandra observations, binned to 500 s, show statistically significant intra-observation variability by a factor of 2 (Figure 5, top panels), with a χ² probability of constant rate <1%. The intra-observational variability during the 2012 and 2017 XMM-Newton observations is more marginal (Figure 5, bottom panels). We searched for periods or quasi-periodic oscillations in both the Chandra and XMM-Newton light curves using powspec, efsearch, and efold, but none of the signal peaks is significant above the noise level. In particular, we searched for periods of ∼1 s in the 2012 and 2017 EPIC-pn light curves (binned to 0.1 s), by analogy with typical periods found in ULX pulsars (e.g., Bachetti et al. 2014; Rodríguez Castillo et al. 2019; Sathyaprakash et al. 2019). For this, we used the LombScargle routine in astropy version 3.2.1 (Astropy Collaboration et al. 2018). For the 2012 light curve, we found a probability >40% that any of the peaks in the Lomb–Scargle periodogram are due to random fluctuations of photon counts; for the 2017 data set, that probability is >50%. Thus, we cannot detect any significant period in this source. The relatively low number of counts and signal-to-noise ratio of these observations are not sufficient for any more detailed variability analysis on such short timescales.

Zoom In Zoom Out Reset image size

We have enough counts for detailed modeling of the spectra from XMM-Newton in 2012 and 2017 and Chandra in 2016. As a first step, we tried a simple power-law model (tbabs × tbabs × pow) with a fixed Galactic absorption component N_H,Gal = 2 × 10²¹ cm⁻² (Kalberla et al. 2005) and a free intrinsic component. The model provides a good fit for the 2016 spectrum (partly because of the more limited band coverage of ACIS) but leaves significant systematic residuals in the 2012 and 2017 spectra. One finding that is already obvious from a power-law fit is that the 2017 spectrum is significantly softer (Γ ≈ 3.5) than the other two. Another finding is a strong residual feature consistent with an unresolved emission line at E ≈ 0.66 keV in the 2017 and (at a weaker level) 2012 data. The presence of this feature is significant for every choice of continuum model (power-law and every other continuum component tested in the rest of this section); we will discuss this line in more detail later.

Next, we tried another simple one-component model: tbabs × tbabs × diskbb. The disk–blackbody model does not improve on the power-law model for any of the three spectra; in simple terms, it is too curved for the observed data points. We improved on the disk–blackbody model in two alternative ways: by adding a second (softer) thermal component (tbabs × tbabs × (diskbb + bbodyrad)) and by using a p-free disk model (tbabs × tbabs × diskpbb) with p < 0.75.

The double thermal model provides a good fit (χ²_ν < 1.2) at all three epochs (Table 2). The lower-temperature blackbody has a characteristic temperature kT_bb ≈ 0.15–0.3 keV, typical of the soft excess observed in many ULXs (Gladstone et al. 2009; Kajava & Poutanen 2009; Kaaret et al. 2017; Zhou et al. 2019), and a characteristic radius of ∼10³ km for all three epochs (again, a common feature for this class of systems). The best-fitting inner disk color temperatures in 2012 (kT_in ≈ 1.25 keV) and 2016 (kT_in ≈ 1.14 keV) are consistent with the temperature expected from the inner disk of a stellar-mass BH at the Eddington limit. The physical inner disk radius R_in is assumed to be the innermost stable circular orbit and is defined as R_in ≈ 1.19r_in, where r_in is the best-fitting radius in the diskbb model in xspec (Kubota et al. 1998); the standard scaling factor 1.19 accounts for the hardening factor (≈1.7) and the fact that the peak effective temperature occurs slightly outward from the inner disk radius, at R ≈ (49/36) R_in (Kubota et al. 1998). In our case, we find ${R}_{\mathrm{in}}\sqrt{\cos \theta }\approx 55$ km in 2012 and ${R}_{\mathrm{in}}\sqrt{\cos \theta }\approx 70$ km in 2016; again, these values are typical of stellar-mass BHs at high luminosities. Instead, the best-fitting temperature kT_in ≈ 0.7 keV, determined for the 2017 spectrum, is too low for the observed luminosity. For comparison, in other BH X-ray binaries in which the accretion disk emission is still the dominant component of the X-ray spectrum at L_X ≈ L_Edd ≈ 10³⁹ erg s⁻¹, the peak color temperature kT_in ≈ 1.2–1.5 keV (Kubota & Makishima 2004; Soria 2007; Sutton et al. 2013). This suggests that although the curvature of a disk model provides a formally good fit to the data, it does not provide the correct physical interpretation of the emission at that epoch. The de-absorbed 0.3–10 keV luminosity, defined as 4πd² times the de-absorbed flux f¹¹ , is L_0.3−10 ≈ 2 × 10³⁹ erg s⁻¹ in both 2012 and 2016 (Table 2).

The p-free disk model is a simple approximation to slim disk models (Abramowicz et al. 1988; Watarai et al. 2001), suitable for stellar-mass BHs near or slightly above their Eddington limit. The definition of the parameter p is that the effective temperature on the disk scales as T(R) ∝ R^−p; the standard disk case is p = 0.75, while near-Eddington BHs tend to have p ≲ 0.6 (Sutton et al. 2017). In our case, the p-free disk provides statistically equivalent fits to the double thermal model for the 2012 and 2016 spectra (${\chi }^{{2}_{\nu }}\approx 1.1$ at both epochs). The best-fitting temperature kT_in ≈ 1.5 keV in 2012 and kT_in ≈ 1.4 keV in 2016, consistent with typical values expected for slim disks near the Eddington limit. For a p-free disk, the physical inner radius is usually defined as R_in ≈ 3.19r_in (Vierdayanti et al. 2008); in our case, we obtain ${R}_{\mathrm{in}}\sqrt{\cos \theta }\approx 75$ km in both 2012 and 2016. De-absorbed luminosities are ≈1.5–2 × 10³⁹ erg s⁻¹ in both epochs. Instead, a p-free disk is not a good fit for the 2017 spectrum (χ²_ν ≈ 1.6); specifically, it does not model well the characteristic downturn seen in the data above 2 keV.

For our fourth attempt to model the data, we used a Comptonization model: tbabs × tbabs × (bbodyrad + comptt). Here the bbodyrad component plays the role of seed thermal emission. Replacing bbodyrad with a diskbb seed component does not change our results, because we are only looking at the Wien section of the curve. This model provides a good fit for all three epochs. For 2012 and 2016, it is statistically equivalent to the double thermal and p-free models; for 2017, it is as good as the double thermal model. The seed photon temperature is kT₀ ≈ 0.15 keV in 2012 and 2017, while it is unconstrained (≲0.1 keV) in 2016. The electron temperature in the Comptonizing cloud, which determines the location of the high-energy downturn, is ∼1 keV in 2012 and 2017, while it is unconstrained in 2016, where the presence of a high-energy downturn is not statistically significant (mostly because of the more limited band coverage in Chandra). Low temperatures and high optical depths are among the defining properties of ULXs (Gladstone et al. 2009) when their spectra are fitted with Comptonization models. The most common interpretation for this kind of spectra is that the harder photons from the inner disk region are down-scattered in a thick disk outflow, seen at high inclination angles. The de-absorbed luminosity for the Comptonization model in 2012 and 2017 is ≈1.5–3 × 10³⁹ erg s⁻¹, similar to the luminosity inferred from the other models. Instead, the luminosity is unconstrained in the 2016 Chandra spectrum because the temperature of the seed thermal component is too low and its normalization is unconstrained; if we neglect the contribution of the seed photons, the luminosity is ≈4 × 10³⁹ erg s⁻¹ in 2016. Alternatively, we fixed the temperature of the blackbody and seed thermal components at kT₀ = kT_bb ≡ 0.10 keV to reduce the degeneracy between temperature, normalization, and absorption column. The best-fitting parameters are essentially unchanged; the de-absorbed luminosity has a best-fitting value of ≈5.4 × 10³⁹ erg s⁻¹ (of which ≈3.5 × 10³⁹ erg s⁻¹ is from the Comptonized component), with a 90% lower limit of ≈2.6 × 10³⁹ erg s⁻¹ and a badly constrained upper limit of ≈3.3 × 10⁴⁰ erg s⁻¹.

In summary, all three spectra are consistent with a mildly super-Eddington stellar-mass BH. The EPIC spectrum from 2012 is more consistent with a broadened disk regime (Gladstone et al. 2009; Sutton et al. 2013) but can also be the result of Comptonization. The 2016 ACIS spectrum is softer than that in 2012 but cannot be reliably classified because of the limited band coverage. The 2017 EPIC spectrum is significantly softer than the other two, is best modeled with a Comptonization model, and belongs to the soft ultraluminous regime. The de-absorbed luminosity was ≈1.5–2 × 10³⁹ erg s⁻¹ in 2012, ≈2–4 × 10³⁹ erg s⁻¹ in 2016, and ≈2–3 × 10³⁹ erg s⁻¹ in 2017. The softness of the 2017 spectrum compared with the other two is obvious when we plot unfolded spectra (Figure 6), based on Comptonization models for consistency.

Zoom In Zoom Out Reset image size

For all continuum spectral models described above, we found a significant emission line residual at E = 0.66 ± 0.01 keV in the 2017 spectrum (Table 2); a similar residual but with much lower significance is also seen in the 2012 spectrum. We modeled the residual with a Gaussian emission line. We found an equivalent width of ${95}_{-60}^{+55}$ eV in 2017 (almost independent of the choice of continuum model) and an intrinsic line luminosity of $\approx 2.4$ $\times {10}^{38}$ erg s⁻¹, which is ≈10% of the total intrinsic luminosity in the EPIC band. The FWHM of the line (σ_line in Table 2) is consistent with zero (that is, the intrinsic velocity broadening of the line is much less than the instrumental broadening) and constrained to be ≲30 eV (≈14,000 km s⁻¹) to the 90% confidence level in the 2017 spectrum, with only small differences between the different continuum models.

To assess the significance of the line component, first of all, we note that the 90% lower limit of the line normalization parameter (photon flux) is >0 in both the 2012 and 2017 spectral fits (Table 2). The improvement in the fit statistic when we include the Gaussian emission line component is ${\rm{\Delta }}{\chi }^{2}\gt 11.4$ for all of the continuum spectral models. We obtained a more rigorous statistical constraint with the likelihood ratio test lrt in xspec: we ran 10,000 simulations for various pairs of models with and without the line component. The significance of the line in 2012 is only ≈70%, while in 2017, it is >99.8%, regardless of the continuum model. In fact, Figure 7 shows that the 0.66 keV emission line in the 2017 spectrum is evident even simply by eye. We also built and inspected EPIC-pn and MOS images in the 0.60–0.70 keV band (Figure 8) for the 2017 data set to make sure that the emission is dominated by the pointlike ULX and not contaminated by, for example, diffuse hot gas in a spiral arm. We compared the narrow-band images from 2017 with those from earlier epochs, and from Chandra when the ULX was not detected, to exclude the possibility of a preexisting supernova remnant at that location. As a result of these tests, we are confident that the ≳10³⁸ erg s⁻¹ line emission seen does come from the ULX.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The most likely identification is the O viii Lyα line (rest-frame energy of 0.654 keV). Such a line was strongly detected in other soft ULXs and is interpreted as a signature of fast outflows (Pinto et al. 2016, 2017; Kosec et al. 2018). Although it is not a surprise to detect this line in a ULX, it is unusual to see it so strong and significant even in a moderate-quality spectrum at CCD resolution (see also Section 4.2 for a comparison with other ULXs).

Another line that has been found associated with ULX outflows is the O vii triplet at 0.56 keV (Pinto et al. 2016, 2017; Kosec et al. 2018). In our spectra, we do not have enough counts to detect or place meaningful constraints on the presence of this line. The main reason for the lower signal-to-noise ratio is that photoelectric absorption is stronger at 0.56 than at 0.65 keV, and, conversely, the pn and MOS effective areas are lower. Nonetheless, we tried adding a narrow line with energy fixed at 0.56 keV and estimated the 90% upper limit to its normalization for each of the four spectral models described before. The result of this test is that a 0.56 keV with the same photon flux as the 0.66 keV line (i.e., ∼10⁻⁵ photons cm⁻² s⁻¹, or EW ≲ 75 eV) is still consistent with our data, even though this additional line would not stand out "by eye" in the spectra. The lack of strong constraints on the strength of the O viii triplet means that we cannot exclude a (minor) contribution to the 0.66 keV line from the O viii Heβ line at 0.67 keV.

The 2017 EPIC spectrum also shows hints of other emission line features in the soft band (Figure 7, bottom right panel), for example, the Mg xi triplet at 1.35 keV. The presence of other possible line features motivated us to try and replace the single Gaussian component in our 2017 spectral fits with a thermal plasma component, vapec in xspec, suitable for collisionally ionized gas. Keeping the metal abundances of all elements at the solar value does not provide an improvement compared with the simpler model fits without a Gaussian line. Similarly, leaving all abundances free but locked together does not improve the fit. Instead, we found that we can obtain good fits (Table 3) by leaving the abundance of some elements (O, Ne, Mg, and Si) free and keeping the other abundances at the solar level. The very high relative abundances of the α elements required to fit our spectrum may be a clue for the type of donor star, as we shall discuss later (Section 4.3). The de-absorbed luminosity in the vapec component is ≈3 × 10³⁸ erg s⁻¹.

For the 2017 spectrum, we then replaced the solar-metallicity intrinsic absorption component (tbabs) with variable abundance absorption models (we tried tbvarabs and vphabs). We did not find any significant improvement or change in the emission line properties when the oxygen abundance in the intrinsic absorber is a free parameter. In fact, the oxygen abundance in the intrinsic absorber is constrained to be ≲2 times solar at the 90% confidence level. This suggests that the gas responsible for the intrinsic absorption is not the same oxygen-rich gas responsible for the strong 0.66 keV emission line; the intrinsic absorption component may be located in the halo and disk of NGC 6946, rather than within the binary system.

Finally, we modeled the stacked spectra from the Swift/XRT observations in 2017, 2018, and 2019. Because of the low number of counts, the detailed models described above are degenerate and cannot be adequately constrained; however, we can still get useful information on the hardness evolution. To do so, we fixed the intrinsic column density N_H = 2 × 10²¹ cm⁻², as derived from most of the Chandra and XMM-Newton spectral fits, and fitted the power-law slope and normalization. We found that the ULX spectrum has hardened in recent years, from Γ = 3.2 ± 0.5 in 2017 (consistent with the contemporaneous XMM-Newton observations), to a more moderate Γ = 2.1 ± 0.4 in 2018, to Γ = 2.6 ± 0.5 in 2019. At the same time, the X-ray luminosity has remained approximately constant at ≈3 × 10³⁹ erg s⁻¹.

We have already described (Section 3.1) how we used the relative offsets to MF16 and SN 20017eaw to pinpoint the location of the X-ray transient in the HST images. We detect only one faint optical source inside the error circle (Figure 9), in the 2016 ACS images, in both the F606W and F814W bands. The same source is too faint to be detectable in the 2004 and 2017 observations (shorter exposure times). Thus, we cannot constrain the source variability. Aside from the positional coincidence, there is no direct evidence that this optical source is the actual donor star for the transient X-ray source.

Zoom In Zoom Out Reset image size

We measured an apparent brightness m_F606W = (26.35 ± 0.15) and m_F814W = (26.05 ± 0.15) mag (corrected to infinite aperture). After correcting for line-of-sight extinction (A_F606W = 0.85, A_F814W = 0.52 mag; Schlafly & Finkbeiner 2011), we obtain absolute magnitudes M_F606W = (−3.95 ± 0.15) and M_F814W = (−3.90 ± 0.15) mag (Table 4). This is consistent with an early B star (main-sequence or subgiant), while it rules out O stars and supergiants. It is a type of optical counterpart very common in ULXs (Gladstone et al. 2013; Motch et al. 2014). However, in Section 4.3 we will propose a white dwarf scenario for the donor star, in which case this faint optical source is unrelated to the ULX.

There are at least two types of transient X-ray sources that reach the ultraluminous regime. One type, well exemplified by CXOU J203451.1+601043, turns on as a ULX after years of nondetection (typically at least 2 orders of magnitude fainter) and then remains in the super-Eddington regime for many years. Another example of this behavior is the transient ULX in M83 (Soria et al. 2012, 2015). In our galaxy, GRS 1915+105 turned on in 1992 (Castro-Tirado et al. 1992) and has been near the Eddington limit ever since, at least until mid-2018, when it started to decline toward the low/hard state (Negoro et al. 2018; Miller. et al. 2019; Rodriguez et al. 2019). The second type of ultraluminous transient is well represented by another source recently discovered in NGC 6946, labeled ULX-4 in Earnshaw et al. (2019a); its outburst lasted only ∼10 days. Another example of a short-duration ULX transient is the first ULX discovered in M31, which went back to quiescence after a few months (Middleton et al. 2013) . In the Milky Way, V404 Cyg provides the most notable example of a short-duration outburst reaching the Eddington luminosity and then declining to quiescence, between 2015 June and August (Sivakoff et al. 2015; Kimura et al. 2016; Muñoz-Darias et al. 2016; Motta et al. 2017).

The most commonly accepted theoretical interpretation of X-ray outbursts in Galactic X-ray binaries is based on the thermal-viscous disk instability (King & Ritter 1998; Dubus et al. 2001; Lasota 2001), which depends on the existence of an outer disk region where hydrogen is mostly neutral. The natural extension of this model for long-duration super-Eddington transients (such as CXOU J203451.1+601043) requires the simultaneous presence of two ingredients: a sufficiently high irradiation of the outer disk to keep it in the ionized state and a sufficiently high accretion rate to keep the source at or above the Eddington limit for years. If the donor is a low-mass star, which cannot keep a persistently high mass transfer rate through the L1 Lagrangian point, a long-duration phase of super-Eddington accretion may still be achieved if the disk is very large¹² and able to store enough mass prior to an outburst (the mass stored in the disk is proportional to R_out³, and the peak luminosity scales as R_out²; King & Ritter 1998). For example, in the case of GRS 1915+105, the estimated disk size is ∼10¹² cm, and the mass in the disk is ∼10²⁸ g, which suffices to keep the X-ray source at or above 10³⁹ erg s⁻¹ for decades (Done et al. 2004). As for the other requirement (i.e., a strong irradiation of the outer disk), broadband spectral modeling of observational data suggests (Sutton et al. 2014) reprocessing fractions of ∼10⁻³ for sub-Eddington stellar-mass X-ray binaries and some ULXs and ∼10⁻² for other ULXs, specifically those with soft X-ray spectra. The reason for the enhanced irradiation factor in some ULXs may be that their strong disk wind intercepts and scatters a fraction of photons emitted in the polar funnel and redirects them onto the outer disk (Sutton et al. 2014).

This is not the only viable explanation for transient behavior in ULXs. If a ULX is powered by an NS rather than a BH, transitions between the propeller and the accretor state will cause transient behavior (Dall'Osso et al. 2015; Tsygankov et al. 2016; Earnshaw et al. 2018). The drop in luminosity happens when the magnetospheric radius of the NS (the radius at which the magnetic pressure "stops" the inflowing matter) becomes larger than the corotation radius of the disk, thus creating a centrifugal barrier (Illarionov & Sunyaev 1975; Stella et al. 1986). The location of the magnetospheric radius depends on the NS spin, mass accretion rate, and strength of the NS magnetic field. The existence of a luminosity gap in the long-term light curve of a transient ULX, between ∼a few 10³⁷ and ∼10³⁹ erg s⁻¹, would be a clue in favor of an accretor/propeller switch model, rather than a thermal-viscous disk instability (where we expect the system to evolve smoothly between higher and lower luminosities without discrete jumps). For CXOU J203451.1+601043, all we can say at the moment is that the system was never detected at any other luminosity <10³⁹ erg s⁻¹; thus, the accretor/propeller model is still viable. Future observations of CXOU J203451.1+601043 will be needed to determine whether, when, and how it will decline.

We showed (Section 3.3) that the EPIC spectrum from 2012 is consistent with either a broadened disk regime or Comptonization, while the 2017 spectrum is significantly softer, is not consistent with disk models, and suggests down-scattering of the direct X-ray emission in a cooler medium, such as the disk outflow. A high-energy downturn already at a photon energy of ≈3 keV (or, equivalently, a steep photon index Γ ≈ 3.5 when fitted in the 0.3–8 keV band) makes the 2017 state of this ULX even softer than NGC 5408 X-1, i.e., the source usually taken as a standard for the soft ultraluminous regime (Sutton et al. 2013; Middleton et al. 2014). The 2017 spectrum puts it in the same class as NGC 55 X-1 (Stobbart et al. 2004; Pinto et al. 2017) and NGC 247 X-1 (Feng et al. 2016), i.e., the two sources that are in a transitional state between the classical ULX and supersoft regimes. The very soft spectral appearance found in 2017 is also reminiscent of the eclipsing ULX CXOM51 J132940.0+471237 in M51 (Urquhart & Soria 2016). The high scattering optical depth (τ ≈ 13) fitted to the comptt model is consistent with the relation between optical depth and coronal temperature found by Pintore et al. (2014), again at the very soft end of the sequence.

The spectral and timing properties of the soft ultraluminous regime are generally attributed to our viewing angle passing through the thick disk outflow (Kawashima et al. 2012; Sutton et al. 2013; Pintore et al. 2014; Middleton et al. 2015a; Narayan et al. 2017; Pinto et al. 2017). This reduces or suppresses our detection of hard X-ray photons (directly emitted from the inner disk region) and increases the contribution of down-scattered soft X-ray photons. Spectral softening may be caused by an increase of the radiatively driven wind mass-loss rate and, therefore, its optical depth, which is a function of, for example, the accretion rate (Poutanen et al. 2007). In the case of CXOU J203451.1+601043, the de-absorbed luminosity fitted to the observed spectral data points varies only by a factor of 2 between the Chandra, XMM-Newton, and Swift observations; however, those luminosity estimates correct only for the effect of cold absorption, not for the down-scattering of hard X-ray photons into the soft X-ray or far-UV band. It is plausible that the intrinsic luminosity in 2017 would be much higher if we could see the direct emission unaffected by the down-scattering outflow (for example, if we had a pole-on view). Precession of the viewing angle may also change the scattering optical depth seen by distant observers for constant wind properties.

If our interpretation of CXOU J203451.1+601043 as an extreme example of the soft ultraluminous regime is correct, we expect two other properties associated with this regime. One is the large amount of short-term variability (Middleton et al. 2015a). Moderate intra-observational variability is observed in our X-ray timing analysis; however, the observed count rate is too low to constrain the rms fractional variability at high frequencies. Thus, we cannot make any firm conclusions on whether this source is more variable than other ULXs.

The second property that should be associated with a thick down-scattering outflow is the presence of spectral residuals in the soft X-ray band, caused by blends of emission and absorption lines (Middleton et al. 2015b; Pinto et al. 2016, 2017). Indeed, we have shown the presence of at least one strong O viii emission line at E = (0.66 ± 0.01) in the 2017 EPIC spectrum, with a luminosity of ≈2 × 10³⁸ erg s⁻¹ and an equivalent width of ≈100 eV. The line was significantly stronger than that in the 2012 spectrum. We have also shown hints of other likely residuals around 1.35 (Mg xi) and 1.7 (Mg xii) keV. The luminosity of the O viii emission line in 2017 was an order of magnitude higher than that of analogous lines detected in other ULXs, such as NGC 1313 X-1 and NGC 5408 X-1 (Middleton et al. 2015b; Pinto et al. 2016), and also an order of magnitude higher than that of the O viii line in NGC 55 X-1 (Pinto et al. 2017), despite the similarity in the X-ray continuum.

Why is the oxygen line so strong? We speculate that the donor star is oxygen-rich, well above solar abundance. One scenario is that the donor star is an oxygen-rich Wolf-Rayet (WO subclass; Crowther 2007; McClelland & Eldridge 2016). The optical luminosity of a WO star is low enough (Sander et al. 2019) to be consistent with our upper limit of ≈−4 mag for the optical counterpart. Our spectral modeling suggests that the source is likely viewed at a high inclination angle. Consequently, if the donor star was a Wolf-Rayet, we would expect to see strong sinusoidal variability or even eclipsing behavior in the X-ray flux, by analogy with other Wolf-Rayet X-ray binaries (Qiu & Soria 2019; Qiu et al. 2019). Instead, no such variability is detected in the light curves of this ULX (Figure 5). A more plausible scenario is that the donor star is a CO or (preferably) an O–Ne–Mg white dwarf; the latter are the most massive subclass of white dwarfs (Truran & Livio 1986; Shara & Prialnik 1994), formed from B-type stars with initial masses just below the limit (≈8 M_⊙) for supernova explosions. In order to form a luminous X-ray binary, the white dwarf must be filling its Roche lobe in an ultracompact system (e.g., van Haaften et al. 2012). Strong, relativistically broadened O viii Lyα lines attributed to Compton reflection have been seen from some (sub-Eddington) ultracompact X-ray binaries (UCXBs) in the Milky Way (Madej et al. 2010, 2014; Madej & Jonker 2011). The origin of the emission line in CXOU J203451.1+601043 may be different (an outflow rather than Compton reflection on the inner disk surface), but we mention this analogy simply as an indication of an oxygen-rich accretion flow. It is also theoretically possible that UCXBs reach super-Eddington luminosities: ULXs in globular clusters, such as the one in the RZ 2109 cluster of NGC 4472 (Maccarone et al. 2007), have been interpreted as UCXBs. In particular, a strong [O iii] λ5007 emission line was detected in the optical spectra of the RZ 2109 source (Steele et al. 2014) and interpreted as evidence of a hydrogen-poor, oxygen-rich donor and an outflow powered by the accreting compact object. In the soft X-ray band, several XMM-Newton and Chandra spectra of the UCXB in RX 2109, presented by Dage et al. (2018), also show emission residuals at 0.6–0.7 keV consistent with O viii Lyα emission.

A possible issue with the UCXB scenario is that the candidate ULX UCXBs suggested in the literature so far are all in globular clusters, where dynamical formation is strongly enhanced; instead, CXOU J203451.1+601043 is in the field, in or near a spiral arm, almost certainly not in a globular cluster, because its optical counterpart is fainter than M_I ≈ −4 mag. The typical optical luminosity distribution of old globular clusters spans $-11\lesssim {M}_{I}(\mathrm{mag})\lesssim -5$ (e.g., Secker 1992; Barmby et al. 2000; Jordán et al. 2007). On the other hand, several Galactic UCXBs are also found in the field outside globular clusters or the bulge (Cartwright et al. 2013). The transient nature of CXOU J203451.1+601043 is also a puzzle: the disk in an ultraluminous UCXB is too hot (and therefore fully ionized) to undergo thermal-viscous instabilities. Instead, the transient behavior could be caused by mass transfer instabilities from the donor star or an accretor/propeller switch if the compact object is a magnetized NS.