# The Black Hole Candidate Swift J1728.9$-$3613 and the Supernova Remnant
G351.9$-$0.9

Mayura Balakrishnan Department of Astronomy, The University of Michigan, 1085
S. University Ave., Ann Arbor, MI, 48109, USA  Paul A. Draghis Department of
Astronomy, The University of Michigan, 1085 S. University Ave., Ann Arbor, MI,
48109, USA  Jon M. Miller  Department of Astronomy, The University of
Michigan, 1085 S. University Ave., Ann Arbor, MI, 48109, USA  Joe Bright
Department of Astronomy, University of California at Berkeley, 501 Campbell
Hall #3411, Berekeley, CA, 94720, USA  Department of Physics, University of
Oxford, Keble Rd, Oxford, UK  Robert Fender  Department of Physics, University
of Oxford, Keble Rd, Oxford, UK  Mason Ng  MIT Kavli Institute for
Astrophysics and Space Research, 70 Vassar St, Cambridge, MA, 02139, USA
Edward Cackett  Department of Physics & Astronomy, Wayne State University, 666
W Hancock St, Detroit, MI, 48201, USA  Andrew Fabian  Institute of Astronomy,
University of Cambridge, Madingley Rd, Cambridge, UK  Kip Kuntz  NASA Goddard
Space Flight Center, Greenbelt, MD, 20771, USA  The Henry A. Rowland
Department of Physics and Astronomy, Johns Hopkins University, 3701 San Martin
Drive, Baltimore, MD, 20218, USA  James C.A. Miller-Jones  International
Centre for Radio Astronomy Research, Curtin University, GPO Box U1987, Perth,
WA, 6845, Australia  Daniel Proga  Department of Physics, University of Nevada
at Las Vegas, NV, 89154, USA  Paul S. Ray  U.S. Naval Research Laboratory,
4555 Overlook Ave., SW. Washington, 20375, DC, USA  John Raymond  Center of
Astrophysics at Harvard, 60 Garden St.,Cambridge, MA, 02138, USA  Mark
Reynolds  Department of Astronomy, The University of Michigan, 1085 S.
University Ave., Ann Arbor, MI, 48109, USA  Department of Astronomy, The Ohio
State University, 140 West 18th Avenue, Columbus, OH, 43210, USA  Abderahmen
Zoghbi  Deparment of Astronomy, University of Maryland, College Park, MD,
20742, USA  HEASARC, Code 6601, NASA/GSFC, Greenbelt, MD 20771  CRESST II,
NASA Goddard Space Flight Center, Greenbelt, MD, 20771, USA

###### Abstract

A number of neutron stars have been observed within the remnants of the core-
collapse supernova explosions that created them. In contrast, black holes are
not yet clearly associated with supernova remnants. Indeed, some observations
suggest that black holes are “born in the dark”, i.e. without a supernova
explosion. Herein, we present a multi-wavelength analysis of the X-ray
transient Swift J1728.9$-$3613, based on observations made with Chandra, ESO-
VISTA, MeerKAT, NICER, NuSTAR, Swift, and XMM-Newton. Three independent
diagnostics indicate that the system likely harbors a black hole primary.
Infrared imaging signals a massive companion star that is broadly consistent
with an A or B spectral type. Most importantly, the X-ray binary lies within
the central region of the catalogued supernova remnant G351.9$-$0.9. Our deep
MeerKAT image at 1.28 GHz signals that the remnant is in the Sedov phase; this
fact and the non-detection of the soft X-ray emission expected from such a
remnant argue that it lies at a distance that could coincide with the black
hole. Utilizing a formal measurement of the distance to Swift J1728.9$-$3613
($d=8.4\pm 0.8$ kpc), a lower limit on the distance to G351.9$-$0.9 ($d\geq
7.5$ kpc), and the number and distribution of black holes and supernova
remnants within the Milky Way, extensive simulations suggest that the
probability of a chance superposition is $<1.7\%$ (99.7% credible interval).
The discovery of a black hole within a supernova remnant would support
numerical simulations that produce black holes and remnants, and thus provide
clear observational evidence of distinct black hole formation channels. We
discuss the robustness of our analysis and some challenges to this
interpretation.

X-rays: stellar-mass black holes — X-rays: binaries — Supernova Remnants –
X-ray transients — Radio observations

††journal: ApJ††software: CIAO (v4.13; Fruscione et al. (2006)), SAS (v1.3;
Gabriel et al. (2004)) v20.0.0 with a development version of ESAS based on
Snowden et al. (2008) that will be released with SAS v21 , HEAsoft (v6.28;
Nasa High Energy Astrophysics Science Archive Research Center (2014)
(Heasarc)), nibackgen3C50 (Remillard et al., 2022), nicer-background
(v0.4.t1.200e; Zoghbi et al. in prep; https://github.com/zoghbi-a/nicer-
background) stingray (Huppenkothen et al., 2019, 2019), XSPEC (Arnaud, 1996),
oxkat (Heywood, 2020), CASA (McMullin et al., 2007b), CARTA (Comrie et al.,
2021), PRESTO (Ransom et al., 2002), Scipy (Virtanen et al., 2020)

##  1 Introduction

While neutron stars have been detected within the remnants of the core-
collapse explosion that created them (Shternin et al., 2011; Heinz et al.,
2013), clear evidence of a black hole within the supernova remnant of a
progenitor star has been elusive. It has been argued that Cygnus X-1, arguably
the most famous black hole X-ray binary, formed without a supernova explosion
(Mirabel & Rodrigues, 2003). The detection of a black hole in a supernova
remnant would broaden our understanding of massive stellar evolution, by
demonstrating an independent formation channel that also allows the creation
of a remnant. An understanding of black hole formation will help us better
grasp the demographics of X-ray binaries and gravitational wave source
populations, as well as the final stages of the lives of massive stars.

As many as three sources may be examples of black holes within supernova
remnants. SS 433 is an X-ray binary associated with the supernova remnant W50.
The combination of dust along the line of sight to SS 433 and a complex local
environment give rise to very high obscuration. This makes it extremely hard
to constrain the mass of the compact object using the motions of its massive
companion star. Indeed, spectroscopy and timing of the compact object itself
may be impossible. Several estimates based on optical photometry and
spectroscopy yield masses above the neutron star range (Blundell et al., 2008;
Kubota et al., 2010; Cherepashchuk et al., 2018). An association between the
black hole candidate MAXI J1535-571 and the remnant G323.7$-$1.0 has recently
been claimed (Maxted et al., 2020). However, differences in the column of
obscuring gas along the line of sight suggest that the black hole and its low
mass stellar companion may lie behind the remnant (Saji et al., 2018; Miller
et al., 2018). If they are associated, a specialized evolutionary channel may
be required to allow a low-mass companion to fill its Roche lobe without
significant binary evolution (Maxted et al., 2020). Finally, the asymmetric
supernova remnant W49B is relatively young, with an age of only $\tau\simeq
5000$ years, yet no cooling neutron star or pulsar has been detected within
it. This absence, a rare morphology, and unusual metal abundances within W49B
leave open the possibility that the explosion formed a black hole (Lopez et
al., 2013).

Core-collapse supernovae explosions leave behind a compact object, along with
diffuse emission in the form of a supernova remnant. The shock produced by the
explosion can heat material to temperatures that make it visible in X-rays.
Gamma-rays can be produced by energetic electrons through inverse Compton
scattering, or by energetic ions through collisions in dense gas that produce
pions (Bykov et al., 2018). Neutron stars and pulsars emit radiation on their
own, while the detection of stellar mass black holes is contingent upon the
black hole accreting material from a companion star and the source going into
outburst. In this paper, we report that the transient X-ray binary Swift
J1728.9$-$3613 likely harbors a black hole primary, consistent with results
from Saha et al. (2022), and may be physically associated with the supernova
remnant G351.9$-$0.9.

On January 25th 2019, the Swift Burst Alert Telescope (BAT) triggered and
located a new transient in outburst that was named Swift J1728.9-3613
(Barthelmy et al., 2019). Subsequent monitoring with MAXI suggested that the
source is an accreting X-ray pulsar or a hard-state black hole X-ray binary,
based on the time scale of the flux rise and the spectral hardness of the
source (Negoro et al., 2019). Observations with NICER revealed a highly
absorbed spectrum, in addition to band-limited variability below 1 Hz and a
narrow quasi-periodic oscillation (QPO) at 5.5 Hz, consistent with type-B QPOs
that are often seen in the soft-intermediate state of black holes (Enoto et
al., 2019). The MeerKAT radio telescope also observed Swift J1728.9$-$3613,
detecting the source six days after discovery (Bright et al., 2019). After a
relatively short bright phase, this outburst of Swift J1728.9$-$3613 decayed
and the source returned to quiescence by mid-2019. No additional outbursts
from the source have been reported in any wavelength.

An extensive analysis of the spectral and timing properties of Swift J1728
using Swift and NICER observations concludes that the system hosts a black
hole (Saha et al., 2022). The authors in Saha et al. (2022) report further on
the detected type-B QPOs and note a partial hysteresis in the RMS-intensity
diagram that supports the conclusion that Swift J1728 is a black hole X-ray
binary.

The position of Swift J1728.9$-$3613 (Barthelmy et al., 2019) is spatially
coincident with supernova remnant G351.9$-$0.9, previously discovered in radio
observations (Whiteoak & Green, 1996). G351.9$-$0.9 is largely unremarkable:
it lacks the spiral cocoon structure of W50 or the barrel structure of W49B.
Rather, G351.9$-$0.9 has a defined edge with heightened emission, and fainter,
smoothly varying emission throughout its interior. Like other remnants that
are highly obscured by foreground dust and gas, it is not detected in optical
wavelengths. The remnant is also not detected in the Fermi-Lat SNR catalog
(Acero et al., 2016).

In Section 2 of this paper, we describe our observations and initial data
reduction. In Section 3, we analyze the nature of the X-ray binary Swift
J1728.9$-$3613\. Section 4 describes the nature of the spatially coincident
supernova remnant, G351.9$-$0.9. In Section 5, we investigate the probability
of coincidental overlap between Swift J1728.9$-$3613 and G351.9$-$0.9. In
Section 6, we discuss the strengths and weaknesses of our analysis, and the
implications of an association between Swift J1728.9$-$3613 and G351.9$-$0.9.

##  2 Observations and Preliminary Data Reduction

We made observations of Swift J1728.9 in outburst using Chandra, MeerKAT,
NICER, NuSTAR, and Swift. After Swift J1728.9 entered quiescence, we also made
X-ray observations of G351.9$-$0.9 using Swift and XMM-Newton, ensuring that
the flux from a bright point source did not obscure or distort any diffuse
emission from the remnant. Radio observations with MeerKAT were made while
Swift J1728.9$-$3613 was in outburst, but this did not impact imaging of
G351.9$-$0.9. NICER observations were used to construct a color-color diagram,
and all NICER observations were searched for specific phenomena; for brevity,
these are not listed within the table. The key observations in our campaign
are listed in Table 7.

###  2.1 Chandra

We observed Swift J1728.9$-$3613 with Chandra on six occasions. In this
analysis, we have focused on just two of these; a full spectroscopic analysis
of the remaining observations will be the subject of a subsequent paper. The
first five Chandra observations were made using the High Energy Transmission
Grating Spectrometer (HETGS), dispersed onto the ACIS-S array of CCDs. The
highest flux was observed in the first observation, OBSID 21213, obtained on
2019 February 02 (MJD 58526.4). We focused on this spectrum to determine the
column density and distance to Swift J1728.9 (see Section 3.2). The final
Chandra observation of Swift J1728.9$-$3613 was obtained when the source was
faint, allowing us to use a standard imaging mode. In OBSID 22289 (obtained on
2019 July 28, or MJD 58692.21), the HETG was not deployed, but the ACIS-S3
chip was operated in a 1/8 sub-array to mitigate photon pile-up (e.g. Davis,
2001), reducing the nominal frame time from 3.2 s to 0.4 s. Both of these
observations were reduced using the standard tools and procedures available in
the CIAO suite, version 4.13, and the associated calibration files. Prior to
analysis, all observations were reprocessed using the CIAO tool chandrarepro.
Further details regarding our preparation and analysis of these observations
follow below.

###  2.2 NuSTAR

Swift J1728.9-3613 was observed with NuSTAR on 2019 February 03 (MJD
58517.62), near the peak of the outburst. The net exposure of the observation
was 16.7 ks. A count rate of $\sim 280~{}{\rm counts}~{}{\rm s}^{-1}$ was
recorded in each of the two NuSTAR Focal Plane Module (FPM) detectors: FPMA
and FPMB. These detectors deliver both imaging and spectroscopy; however,
owing to the brightness of the source, nearly the entire detector plane is
illuminated and any weak diffuse emission is overwhelmed. Using the standard
routines in HEASOFT version 6.28 through the NuSTARDAS pipeline v2.0.0 and
CALDB version 20210524, we extracted time-averaged spectra from the entire
3–79 keV pass band from the two NuSTAR detectors, using circular extraction
regions of radius 120” centered on the position of the source. We extracted
the background rate from regions of the same size as the source regions,
placed at the edges of the detectors. The full details of the observation and
our analysis of the spectra is the subject of companion work by Draghis et al.
(2023). In the sections that follow, we utilize the key results of that
analysis.

###  2.3 XMM-Newton

After Swift J1728.9$-$3613 entered quiescence, XMM-Newton observed
G351.9$-$0.9 on three occasions (see Table 7). The first observation, OBSID
0860140101, was an 18 ks snapshot that started on 2020 September 04 (MJD
59116.17). The next two observations were much deeper, each lasting 83 ks;
OBSID 0901540101 started on 2022 February 28 (MJD 59638.07), and was
immediately followed by OBSID 0910540201 starting on 2022 March 04 (MJD
59642.06).

The central CCDs in the EPIC-MOS cameras are able to cover most of
G351.9$-$0.9, potentially enabling an image with minimal disruptions owing to
chip gaps and other considerations. In all three observations, then, the PN
and MOS cameras were run in “full frame” mode, providing the largest possible
field of view. The “thin” EPIC optical blocking filter was utilized in the
initial snapshot, whereas the ”medium” filter was used in the deeper
exposures. Owing to the large column density along this line of sight, our
analysis was restricted to the EPIC cameras. All of the data were reduced
using XMMSAS version 19.1.0 and the corresponding “CCF” calibration files. The
specifics of how the data were prepared and analyzed for different purposes
are detailed in the sections that follow.

###  2.4 NICER

The NICER telescope’s primary instrument, the X-ray Timing Instrument (XTI),
is made up of 52 live detectors with an operating range of 0.2–12 keV. NICER
observed Swift J1728.9$-$3613 at a high cadence throughout its 2019 outburst.
For brevity, we list in Table 7 only the 48 observations featured in Figure 1,
although a spectroscopic analysis was performed on the initial 80
observations. In order to conduct a search for Type-1 X-ray bursts in the
light curve, and to search for coherent pulsations, we constructed and
analyzed light curves from all of the NICER observations of Swift
J1728.9$-$3613 available in the public archive. All of the observations were
processed using the standard pipeline tool, nicerl2 in HEASOFT version 6.28,
and the files in CALDB version 20210728. We utilized three background
estimation tools that are publicly available. Additional details follow in the
sections below.

###  2.5 Neil Gehrels Swift Observatory

Swift observed Swift J1728.9-3613 six times: once in windowed timing (WT) mode
and five times in photon counting (PC) mode (see Table 7). We processed all of
the observations using tools available in HEASOFT version 6.28 and the
associated calibration files. Initial event cleaning was carried out using the
task XRTPIPELINE. Using the standard tools, we constructed light curves for
the purpose of searching for X-ray bursts. In addition, Swift observed
G351.9$-$0.9 on 2022 September 2022 (MJD 59844.56) for 1.5 ks. This
observation was obtained in PC mode and was also processed using the standard
tools.

###  2.6 MeerKAT

We observed the field of Swift J1728.9$-$3613 with the MeerKAT radio telescope
(Mauch et al., 2020) as part of the ThunderKAT large survey project weekly XRB
monitoring program (Fender et al., 2016) beginning on 2019 February 12 (MJD
58533). We obtained a total of 13 observations amounting to an on-source time
of 255 minutes. We initially reduced and imaged each observation independently
using the oxkat (Heywood, 2020) reduction pipeline, a set of semi-automated
scripts that perform phase reference calibration, amplitude and phase self-
calibration, and wide-field imaging (McMullin et al., 2007a; Offringa et al.,
2014; Kurtzer et al., 2017; Kenyon et al., 2018; Makhathini, 2018). From each
of these images we measured the flux density from G351.9$-$0.9.

In addition to imaging each observation individually, we also created a deep
image of the field using the entire 255 minutes of on-source observing time.
We concatenated the individual self-calibrated measurement sets in CASA and
imaged the combined data using WSClean (Offringa et al., 2014). We cleaned the
field using $8\sigma$ automasking, then re-imaged the field using a custom
mask based on the initial clean image. We employed multiscale cleaning to
attempt to account for extended structure in the field and set a stopping
threshold of 20$\mu\rm{Jy}$. This is significantly larger than the expected
thermal noise from a 255 minute MeerKAT observation, however extended
structures in the field make reaching the thermal noise impossible.

##  3 X-ray Binary Swift J1728.9–3613

###  3.1 Source Position

After its discovery with the Swift/BAT (Barthelmy et al., 2019), follow-up
observations with the XRT recorded a source position of RA,DEC $=$ 17h 28m
58.64s, -36∘ 14’ 37.7” (J2000), with an error of 1.7” (Kennea, 2019). This
places the source within the region of the sky occupied by the catalogued
supernova remnant G351.9$-$0.9 (Whiteoak & Green, 1996). In order to determine
the best Chandra position of the source, we ran the CIAO tool wavdetect on
Chandra OBSID 22289. Swift J1728.9$-$3613 is detected at (RA,DEC) = 17h 28m
58.76s, -36∘ 14’ 36.20”. The statistical error on this position is negligible,
but a systematic error of 0.3” is reasonable based on the limits of Chandra’s
attitude control system (e.g. Miller et al., 2017). The nominal positions
differ by 2” but are formally consistent within errors. We adopt the Chandra
measurement as the position of Swift J1728.9 in the work that follows.

###  3.2 Column Density and Distance

The high column density along the line of sight to Swift J1728.9$-$3613 makes
optical measures of its distance impossible. If the source is once again
observed in radio during a prolonged outburst, it may be possible to determine
its proper motion and distance using very long baseline interferometry. In the
meantime, the same column density that complicates other distance measures can
be turned into a distance constraint in X-rays. In short, the distance to a
given source can be determined by integrating the gas density distribution
along the line of sight until the accumulated column density matches that
measured in X-rays. Using a combination of measured galactic $n_{H}$ radial
density profiles (Kalberla & Kerp, 2009; Marasco et al., 2017), we obtained a
continuous distribution of neutral hydrogen in the galactic midplane as a
function of galactocentric radius. For a given line of sight, we calculate the
galactocentric radius and get the equivalent hydrogen number density, $n_{H}$
(cm-3), and integrate along that line of sight to obtain a column density
$N_{H}$ (cm-2).

To determine the column density to Swift J1728.9$-$3613, we made joint fits to
the combined MEG and combined HEG spectra obtained in Chandra OBSID 21213
using XSPEC version 12.11.1 (Arnaud, 1996). Prior to fitting, the spectra were
binned according to the “optimal binning” algorithm detailed in Kaastra &
Bleeker (2016), using the tool ftgrouppha. Fits to the MEG spectrum were
restricted to the 1–5 keV range, and fits to the HEG spectrum were restricted
to the 1.2–8 keV range, owing to poor sensitivity below and above these
bounds. The spectra were fit minimizing the Cash statistic (Cash, 1979),
allowing a normalizing constant to float between them. The tbabs (Wilms et
al., 2000) model was used to describe line-of-sight absorption, after adopting
the required cross sections and abundances (“Verner” abundances and “Wilms”
cross-sections, within XSPEC). The continuum is well described with only a
disk blackbody (Mitsuda et al., 1984). This model measures $N_{H}=4.22\pm
0.02\times 10^{22}~{}{\rm cm}^{-2}$, $kT=1.27\pm 0.01$ keV, and a flux
normalization of $K=242\pm 5$, and gives a fit statistic of $C=3270.4$ for
$\nu=2365$ degrees of freedom ($1\sigma$ errors). The constant indicates that
the MEG normalization is 8% higher than the HEG in OBSID 21213.

This value of the column density to Swift J1728.9$-$3613 is broadly consistent
with that measured using NuSTAR, $N_{H}=3.7^{+0.2}_{-0.1}\times 10^{22}~{}{\rm
cm}^{-2}$ (Draghis et al. 2023). In that work, they considered a series of
more sophisticated models, including thermal emission from the accretion disk,
non-thermal emission from the corona, and reflection from the innermost region
of the disk. These measurements are not in formal agreement at the $1\sigma$
level, but the degree of agreement is excellent considering the different pass
bands and point spread functions of each telescope, and the different models
required to fit the spectra.

Integrating along the line of sight, the column density measured with Chandra
nominally translates to a distance of $d=9.08^{+0.11}_{-0.09}$ kpc, and the
column density measured with NuSTAR indicates a distance of
$d=7.9^{+0.34}_{-0.27}$ kpc. Although it is likely the case that Chandra is
better suited to this measurement, we conservatively regard the difference as
a systematic error, and adopt a distance constraint of 7.6 kpc $\leq d\leq$
9.2 kpc, or $d=8.4\pm 0.8$ kpc.

In order to verify that this technique delivers reasonable distance estimates,
we collected published column densities from a number of X-ray binaries
located close to the plane of the Milky Way that have independent distance
constraints. For each source, we calculated the distance to the source based
on its column density, and compared that value to the independent constraint.
In general, the two are in good agreement; at least within the plane, it
appears that our technique is robust. Please see Appendix A for additional
details.

###  3.3 X-ray Colors

Figure 1: Recreations of Figure 8 from Done & Gierliński (2003), where black
holes are represented by black half-circles, and the filled circles refer to
neutron stars. The red, cyan, and blue filled circles correspond to atolls, Z
sources, and Circinus X-1, respectively. The orange hatched region represents
the area on the plot inaccessible by neutron stars due to the physical
boundary at their surfaces. We overplot values from our analysis: the violet
half-circles are measured from NICER spectra, while the green diamond point
represents the NuSTAR observation. Each panel corresponds to NICER hard and
soft color values derived using different background estimation routines.
Please see the text for details.

Although accreting neutron stars and black holes both have accretion disks and
non-thermal coronae, the build-up of accreted material in a boundary layer on
the neutron star surface is potentially a distinguishing feature. Done &
Gierliński (2003) analyzed a large set of RXTE observations of X-ray binaries
to show that the absence of boundary layer emission may be used to infer that
the compact object in a given X-ray binary is a black hole, and that these
differences are evident in X-ray colors. In that work, no neutron star systems
show a hard color (defined as the ratio of 9.7 - 16 keV flux to 6.4 - 9.7 keV
flux) between 0.25 and 1, and a soft color (defined as the ratio of 4 - 6.4
keV flux to 3 - 4 keV flux) that is lower than $\sim$ 1.4. For reference, we
have recreated most of Figure 8 from Done & Gierliński (2003) in Figure 1,
with black semi-circles indicating black hole systems, filled circles
indicating neutron star systems, and the orange hatched region signifying
X-ray color ranges inaccessible to neutron stars. We do not replicate the
tracks shown in the original figure that arise from models with different
combinations of disk temperature and photon index, as these are not relevant
to this work. NICER obtained the densest monitoring of Swift J1728.9$-$3613,
and its effective area is much higher than that of Swift close to 10 keV. In
order to apply the diagnostics developed by Done & Gierliński (2003), we
therefore analyzed the first 50 NICER observations of Swift J1728.9$-$3613\.
These represent the brightest observations, from which the most robust colors
can be measured. For each observation, we filtered for events in 1–10 keV
energy range and removed instrumental flaring intervals. Source spectra were
extracted from cleaned event files, and the background spectra were generated
using nibackgen3C50 (Remillard et al., 2021). The standard NICER response
matrix and ancillary files were used from CALDB version 20210728. Prior to
fitting each spectrum in XSPEC, the data were binned using the “optimal
binning” scheme developed by Kaastra & Bleeker (2016) and implemented in the
tool ftgrouppha. Each spectrum was fit over the 1–10 keV range with a simple
model consisting of a disk blackbody and a power-law, modified by line of
sight absorption: tbabs $\times$ (diskbb + powerlaw). The hardest band in the
Done & Gierliński (2003) diagnostic, 9.7–16.0 keV, extends above the NICER
pass band. Since this scheme relies on fluxes, not count rates, it is possible
to extrapolate models into this regime. However, caution is required, because
the NICER background is complicated by variable space weather, its position
aboard the International Space Station, and the fact that NICER is not a
focusing instrument. Background flux is inevitably highest at the edges of an
instrumental pass band. We therefore repeated the above procedure, using two
alternative derivations of the NICER background, the “space weather” method
(also supported by the NICER team) and nicer-background v0.4.t1.200e (Zoghbi
et al. in prep; https://github.com/zoghbi-a/nicer-background). We grouped the
spectra using oggrpha
(https://github.com/zoghbi-a/aztools/blob/master/scripts/ogrppha.py), which
allows us to take into account the high background by binning the spectra to a
minimum S/N of 3 and oversampling the detector resolution by a factor of 3.
Again following the work of Done & Gierliński (2003), we used the counts in
each band as a proxy for the errors, allowing the fact that we were using 3 -
10 keV spectra to be extrapolated up to 16 keV to be conveyed in the errors on
the hard color. To illustrate, the error on the soft color was calculated as
follows:

$\frac{\sigma_{soft}}{\text{color}_{\rm
soft}}=\sqrt{\left(\frac{\sqrt{\text{counts}_{\rm 3-4keV}}}{\text{counts}_{\rm
3-4keV}}\right)^{2}+\left(\frac{\sqrt{\text{counts}_{\rm
4-6.4keV}}}{\text{counts}_{\rm 4-6.4keV}}\right)^{2}}$

Figure 1 shows the location of these NICER observations of Swift
J1728.9$-$3613 on the color–color diagram with violet half-circles. The NuSTAR
observation is indicated by the green diamond. The three panels in the figure
correspond to the three different background methods that were utilized. For
clarity, only those observations yielding colors with fractional errors less
than 10% are plotted (29 of 50 observations). The strong majority of the
points fall within the zone that Done & Gierliński (2003) identify as only
possible for black holes. Including the observations with less certain colors
does not change this result. Similarly, the results are robust against the
implementation of different binning schemes within the simple spectral fits.

The observation of Swift J1728.9$-$3613 made with NuSTAR is particularly
important to our analysis, and the NuSTAR pass band is more similar to the
RXTE spectra considered by Done & Gierliński (2003). Adopting the simple
tbabs*(diskbb + powerlaw) model, fitting the full 3–79 keV pass band, and
deriving errors based on the counts in each band, we derive a hard color ratio
of 0.652 $\pm$ 0.004 and soft color ratio of 1.258 $\pm$ 0.006. This point is
shown in Figure 1 as a green diamond that also lies within the zone that is
only accessible to black hole X-ray binaries. Draghis et al. (2023) fit the
NuSTAR spectrum of Swift J1728.9 with a set or disk reflection models; the
hard and soft colors derived from these models differ from the above values by
less than 1% and still fall within the black hole zone.

###  3.4 Spin Parameter

The dimensionless angular momentum or “spin parameter” of a compact object,
$a=cJ/GM^{2}$, can act as a direct probe of its nature. Neutron stars have an
observed spin limit of approximately $a\simeq 0.2$ (Miller & Miller, 2015;
Reynolds, 2021), while the theoretical upper limit on neutron star spin due to
break up is $a\leq 0.7$ (Lo & Lin, 2011). Previous timing and spectral
analysis (Saha et al., 2022) and our measurements of X-ray colors observed in
Swift J1728.9$-$3613 are suggestive of a black hole primary, but a spin
parameter that demands a black hole primary is potentially more definitive.

X-ray measures of black hole spin all rely on the premise that the accretion
disk obeys the innermost stable circular orbit (ISCO) set by the spin
parameter, when the Eddington fraction is moderate. The most broadly
applicable method of determining black hole spin is based on measuring
relativistic distortions to the reflection spectrum of the inner accretion
disk. The extremity of these distortions is set by the ISCO, and therefore by
the spin of the black hole. This is now a well-established spectral modeling
methodology that has resulted in numerous precise spin measurements,
especially when the broad pass band and high energy sensitivity of NuSTAR is
utilized (see, e.g. Draghis et al., 2022).

In a companion work, Draghis et al. (2023) modeled the relativistic reflection
features present in the NuSTAR observation (OBSID 90501303002). Fits to the
3–79 keV NuSTAR spectra with the “relxill” family of disk reflection models
(e.g. Dauser et al., 2013; García et al., 2013) yield a black hole spin
parameter of $a=0.86\pm 0.02$ ($1\sigma$ statistical errors). This value
nominally excludes the theoretical maximum for neutron stars at the $5\sigma$
level of confidence. The highest observed neutron star spin parameters are
excluded at much higher levels of confidence.

This spin measurement is a lower limit; determined by the innermost radius
consistent with the observed spectrum. Inner regions of the disk that are not
illuminated would not be measured, and therefore the spin of the black hole
could be higher.

###  3.5 Search for Pulsations

Although Swift J1728.9$-$3613 appears to harbor a black hole, the detection of
coherent pulsations would unambiguously signal a neutron star primary. For
this reason, it is important to understand the sensitivity of the NICER data
to coherent pulsations similar to those observed in other X-ray binaries.
These dense monitoring obtained with NICER, and its pass band, afford the best
chance of detecting pulsations.

Figure 2 shows a hardness–intensity diagram, constructed using NICER
observations obtained in between MJD 58513–58548, 58559–58613, 58633–58650,
and 58651–58768. Specifically, the total count rate is plotted versus the
ratio of counts in the 4–10 keV to 2–4 keV band. In this HID, the data can be
grouped into four spectral states. The time series was broken up into 128
second segments, and for each segment, we constructed 1 s-binned light curves;
if the light curve had a duty cycle of at least 90%, it was admitted into the
calculation of the averaged spectrum.

For each admitted 128 s segment of data from these observations, the time
series was then binned to 0.5 ms (i.e., Nyquist frequency of 1000 Hz), and we
calculated the Leahy-normalized power spectrum. The total power spectrum is
then averaged across all of the admitted segments (Bartlett, 1948). In the
end, we had 152, 17, 182, and 195 segments corresponding to the four spectral
states identified for the search. No coherent pulsations were detected above
the noise.

Figure 2: Hardness intensity diagram from NICER observations of Swift
J1728.9$-$3613\. The different boxes show how the observations were split into
four spectral states in order to do the pulsation search. No pulsations are
found, with upper limits ruling out pulse amplitudes in the range typically
detected in neutron star X-ray binaries.

We then calculated the 3-sigma upper limits on the fractional pulsation
amplitude using the formalism laid out in Vaughan et al. (1994), modifying it
for the case of a stacked power spectrum. The maximum possible pulse
amplitudes (3-$\sigma$ upper limits) for Swift J1728.9$-$3613 at 10, 100, and
400 Hz are 0.0051%, 0.0078%, 0.36%, and 0.053%, respectively. This is lower
than typical pulse amplitudes for millisecond pulsars, which range from 1-10 %
(Patruno & Watts, 2021).

Acceleration searches are excellent for detecting pulsations from pulsars in
binary systems, because the Doppler motion due to the orbit modulates the NS
pulsation frequency, and the signal that was spread across multiple Fourier
bins can be combined (Ransom et al., 2002). We initially ‘bunched up’ adjacent
GTIs; that is, if the interval between the end of one GTI and the start of the
next GTI is sufficiently short, the two GTIs were combined. With a gap of 10
s, we had 301 GTIs in the end, and we only performed acceleration searches on
GTIs that were longer than 200 s (175 such GTIs). The acceleration searches
were done with accelsearch in PRESTO over 0.01–100 Hz (Ransom et al., 2002).
We also set $z_{\rm max}$ to 100, where we surmised that any Doppler-modulated
signal would have been smeared across a maximum of 100 Fourier bins (standard
choice in the field). The range of accelerations searched is given by
$a=zc/(f_{0}T^{2})$ where $a$ is the orbital acceleration (${\rm m\,s^{-2}}$),
$z$ is the number of Fourier bins the signal would smear across, $c$ is the
speed of light (m/s), $f_{0}$ is the pulsar frequency (Hz), and $T$ is the
span of the segment being searched. Upon visual inspection of the candidates
and their diagnostic plots (e.g., subintegration plots), there were no
significant candidates that could be identified as a coherent signal.

###  3.6 A Search for Type I X-ray Bursts

Type I bursts are thermonuclear runaway events that result from the build-up
of accreted material on the surface of a neutron star. Therefore, the
detection of one or more bursts in a given X-ray binary unambiguously signals
that it harbors a neutron star. No bursts have been reported in Swift
J1728.9$-$3613, but it is important to make a systematic search and to
understand whether or not prior observations were capable of detecting bursts
with the properties and recurrence times observed in other sources. It is also
important to note that not all neutron star binaries display type-I bursts.

The frequency, strength, and timescales of Type I bursts depend on the
composition of the burning material and its metallicity. There are variations
from source to source because of differences in the core temperature of
neutron stars and the accretion rate (Galloway et al., 2008). In particular,
Type I bursts are typically seen at Eddington fractions of $\lambda_{\rm
Edd}=0.003-0.3$ when looking at the 2.5 – 25 keV flux band, (Galloway et al.,
2008).

Using the distance estimate to Swift J1728.9$-$3613 determined in this paper
($d=8.4\pm 0.8\rm\ kpc$) and the spectra extracted from all of the
observations made using Chandra, Swift, NICER, and NuSTAR, we selected
observations that fall within this window (extrapolating the pass band, where
necessary). $L_{edd}$ was calculated with an assumed mass of 1.4 M⊙. NICER
observations dominate this set, owing to the dense monitoring it achieved
across the outburst of Swift J1728.9$-$3613\. As noted above, NICER has a
complex background; we therefore generated background light curves in the
0.4–15 keV band from telescope environment data and compared them to the
source light curves. This was repeated for all three background formulations.
All potential bursts were due to increases in the detector background, and not
due to flaring or bursting in Swift J1728.9$-$3613\. We note that we repeated
this procedure for all NICER observations of Swift J1728.9$-$3613, without
respect to Eddington fraction, and we again find no evidence of bursts from
the source. In their timing analysis, Saha et al. (2022) also report no
bursts.

###  3.7 The Binary Companion Star

The Chandra and NuSTAR spectra of Swift J1728.9$-$ 3613 in outburst are
broadly consistent with a column density of $N_{H}=4\times\rm 10^{22}\
cm^{-2}$. The color excess is given by: $E(B-V)=\frac{N_{H}}{5.8\times
10^{21}\rm atoms\ cm^{-2}mag^{-1}}$ (Cox, 2000). For typical conditions in the
ISM of the Milky Way, the optical extinction, AV, is related to the color
excess by: $\frac{A_{V}}{E(B-V)}\sim 3.1$. Our column density measurements
then correspond to A${}_{V}\sim 21$ magnitudes of extinction. Even an A0
supergiant with MV = $-$6.5 would then have an apparent magnitude $\sim$
29-30. This is at the very limit of what can be detected in the deepest
observations with Hubble; it is well beyond what can be detected with the
Swift UVOT telescope and ground-based telescopes. We conclude that the optical
counterpart of Swift J1728.9$-$3613 has not been detected, and that optical
observations cannot reveal the nature of the companion star and binary system.

###  3.8 Infrared Observations

Figure 3: Finder chart showing the IR emission from ESO-VISTA of Swift J1728.
The red circle represents a region of radius equal to 1.3 arc seconds,
centered at the known position of the source. The image was composed by
displaying the $K_{s}$ band in red, the $H$ band in green, and the $J$ band in
blue. Figure 4: _Left_ : The IR light curve of Swift J1728. Red points
represent $K_{s}$ magnitudes, green points represent $H$ magnitudes, and blue
points represent $J$ magnitudes. The horizontal lines show the median
magnitude prior to the outburst, in the three bands. The vertical dashed line
shows the date of the Swift detection of the outburst. _Right_ : Apparent
magnitudes for main sequence stars assumed to be at a distance of $8.4\pm
0.8~{}\rm kpc$ (dashed line with shaded regions). Values are computed assuming
$N_{H}=4\times 10^{22}\;cm^{-2}$ and $R_{V}=A_{V}/E(B-V)=3.1$.

We examined archival infrared observations in $J,H,$ and $K_{s}$ bands from
the VIRCAM camera on the ESO-VISTA telescope to search for emission at the
Chandra X-ray position of Swift J1728.9$-$3613\. Figure 3 was composed using a
$K_{s}$ frame to illustrate the color red, an $H$ frame to show color green,
and a $J$ frame to show color blue. The red circle has a radius of 1.3” and is
centered at the Chandra position. A candidate source is clearly
distinguishable from the background. We therefore downloaded a total of 31
$K_{s}$ frames, 7 $J$ frames, and 3 $H$ frames (all fully reduced) from the
ESO Science Archive and ran Gaussian-weighted aperture photometry at the
position of Swift J1728.9$-$3613 and of roughly 4000 other point sources with
known $J,H,$ and $K_{s}$ magnitudes obtained from the Two Micron All Sky
Survey (2MASS) catalog (Skrutskie et al., 2006). We required the catalog
sources to have an apparent J magnitude between $11.5\leq m_{J}\leq 14.5$, in
order to avoid saturated sources and still have reliable magnitude uncertainty
measurements. The aperture photometry was performed via the NOAO IRAF software
suite, using circular extraction regions with 5-pixel radius. The Gaussian
weighting was set to be equal to the FWHM of the individual images. We
extracted the background from annuli centered at the position of the sources,
with an inner radius and width of 10 pixels, and calculated the background
rate as the mode of the values of the pixels in the annulus. In order to
account for possible small calibration errors in the WCS registration of the
frames, we allowed a centering algorithm to shift the position of the source
and background regions by at most one pixel. We excluded sources where a
region shift of more than one pixel was needed. Last, we performed
differential photometry in order to calculate the apparent magnitude of the
source at the position of Swift J1728.9$-$3613\. The left panel in Figure 4
shows the apparent magnitude in the $K_{s},H,$ and $J$ filters in red, green,
and blue (respectively) as a function of time, over the past 11 years. Non-
detections and measurements with an uncertainty larger than 0.7 magnitudes
were not plotted or used in subsequent analysis. The dashed vertical black
line corresponds to the time of the Swift discovery of the source. The
horizontal lines correspond to the median $K_{s},H,$ and $J$ magnitudes over
the measurements, prior to the 2019 outburst. The median apparent magnitudes
correspond to $m_{K_{s}}\sim 18.2$, $m_{H}\sim 18.7$, and $m_{J}\sim 20.9$.
Shortly after the 2019 outburst, the $J$ apparent magnitude of the source is
observed to be 2 magnitudes brighter than the median before the outburst;
thereafter, the $J$ flux decays in time. About 200 days after the Swift
detection, the $K_{s}$ apparent magnitude is $\sim 0.7$ magnitudes brighter
than prior to the outburst. This suggests a significant IR brightening either
coincidental to the outburst or shortly after that could be associated with
jet production in harder spectral states or disk brightening, followed by a
decay to the value observed 200 days after the Swift discovery.

Figure 5: Stacked $1.28$ GHz MeerKAT image of G351.9$-$0.9, with white
contours at 25 $\mu$Jy $\times$ (-3, 3, 6, 9, 12, 15) $\sigma$. The radio
color bar indicates the flux density scale in Janskys. The Chandra position of
Swift J1728.9-3613 is indicated with a filled blue circle.

The right panel in Figure 4 shows the expected $K_{s},H$, and $J$ apparent
magnitudes of main sequence stars of spectral types ranging from O9 to M0, in
red, green, and blue respectively. The dashed lines show expected apparent
magnitudes for different stellar types (values taken from Cox 2000),
calculated for a distance of $d=8.4\pm 0.8$ kpc and a hydrogen column density
of $N_{H}=4\times 10^{22}$ cm-2 while accounting for the gas-phase, the grain-
phase, and molecular interstellar medium and a ratio of total to selective
extinction $R_{V}=A_{V}/E(B-V)=3.1$. While the pre-outburst emission does not
perfectly correspond to the expected magnitudes of a catalog main sequence
star, the values of the apparent magnitudes in the three bands and the colors
are broadly consistent with a late B or A class star.

##  4 Supernova Remnant G351.9–0.9

###  4.1 Radio Imaging

The deep, 255-minute MeerKAT image of the field containing Swift
J1728.9$-$3613 yielded an excellent view of the supernova remnant G351.9$-$0.9
shown in Figure 5. The Chandra position of Swift J1728.9$-$3613 clearly places
the source within the central region of the remnant. G351.9$-$0.9 has an
elliptical shape, with a major axis that is approximately 12’ long, and a
minor axis that is approximately 9’ long. It is center-filled, but retains a
clear edge structure with enhanced emission along its perimeter. The
appearance of G351.9$-$0.9 in the deep MeerKAT image resembles that of the
Cygnus Loop, which is in the Sedov phase (Fesen et al., 2018). Both remnants
have higher brightness patches where the supernova blast wave is interacting
with denser clouds and the shock has become radiative, increasing the radio
emission due to the compression (Raymond et al., 2020). This implies the
remnant is still in the Sedov phase, and therefore less than $\tau=30,000$
years old (Cioffi et al., 1988). A pressure-driven post-Sedov supernova
remnant would appear more disjoined and more strongly limb-brightened, while a
pre-Sedov remnant would be smaller and brighter given the distance constraint
provided by the column density. At a distance of 8.4 kpc (coinciding with
Swift J1728), the supernova remnant would have a diameter of $\sim$ 20-30 pc,
which is reasonable for a remnant in the Sedov phase.

To constrain the total flux density contribution from the remnant, we used
image analysis tools from CARTA (Comrie et al., 2021). In our final clean
image we measure a sum of $64\,\rm{Jy~{}beam}^{-1}$ (corresponding to a flux
density of $1.8\,\rm{Jy}$) from $\sim 3.4\times 10^{5}$ pixels from a manually
drawn region encompassing the structure. Within the same region we measure a
sum of $5.3\,\rm{Jy~{}beam}^{-1}$ (corresponding to a flux density of
$0.2\,\rm{Jy}$) in the residual image, implying $\sim 1.6\,\rm{Jy}$ within the
sky model of the extended emission.

Figure 6: LEFT: The combined MOS-2 image of the field containing G351.9$-$0.9
in the 2.0-7.5 keV band. RIGHT: The combined MOS-2 image of the same field in
the 0.3-1.0 keV band. In both cases, strong point sources were excised, soft
proton flaring was removed on a chip-by-chip basis, an estimate of the
quiescent particle background has been subtracted, and the images are
exposure-corrected and then smoothed with a Gaussian ($\sigma=$10 pixels). The
contours shown in green trace the outline of the remnant in radio bands, as
detected with MeerKAT. In both panels, the Chandra position of the black hole
candidate Swift J1728.9-3613 is indicated with a filled yellow circle. There
is no evidence of X-ray emission from the remnant, potentially consistent with
a spectrally soft, young remnant that is obscured by the relatively high
foreground column density.

There are a number of point sources embedded within the diffuse radio
structure, including Swift J1728.9$-$3613\. However, we did not remove them
before calculating the total flux density as they each contribute $<<1\%$ of
the total. In addition to the diffuse structure, there are a number of bright
extended structures detected outside of the remnant by MeerKAT; these cause
the image noise to be highly non-Gaussian and non-uniform. It is therefore
difficult to determine if the remaining flux density in the residual image
corresponds to the elevated noise level, or to genuine emission, or indeed if
the clean model contains exclusively emission from the nebula itself. We
therefore place a tentative range of $1.6$-$1.8\,\rm{Jy}$ on the total diffuse
emission. Our estimates are consistent with other reported measurements of
G351.9-0.9: 10.7 Jy at 300 MHz (uGMRT; Veena et al., 2019) and 2.0 Jy at 843
MHz (MOST; Whiteoak & Green, 1996). These measurements translate to $\sim$ 0.9
- 1.6 Jy at 1.28 GHz for the 1 $\sigma$ range of flux densities and spectral
indices reported.

Additionally, we note that, due to MeerKAT being a radio interferometer, our
images are affected by the ‘zero-spacing’ problem. This results in flux
density not being recovered beyond a certain angular scale determined by the
shortest baseline length. Angular scales $\lesssim 10\arcmin$ should be well
recovered by Meerkat (Knowles et al., 2022) due to the short minimum baseline
and abundance of short spacings. This is the approximate size of G351.9$-$0.9,
so while our flux density measurement is formally a lower limit it should
provide a reasonable estimate for the total emission from the remnant at
$1.28\,\rm{GHz}$.

###  4.2 X-ray Imaging

The EPIC cameras have a very large field of view, approximately 30’. This
feature, and the sensitivity afforded by the collecting area of these cameras,
make them excellent for surveys of point sources and bright extended sources.
However, imaging of faint extended sources is more difficult. Instrumental
flaring due to soft proton flux, the quiescent particle background of the
instruments, and the diffuse astrophysical foreground caused by solar wind
charge exchange all pose significant difficulties. The ESAS software suite
(Snowden et al., 2004, 2008) was developed to handle imaging of faint diffuse
structures; it is particularly effective when a source is large enough that
local backgrounds are difficult or impossible to obtain. This suite combines
standard SAS routines and new routines to produce images and spectra that are
as free from these complications as possible. We found that some of the tools
within the ESAS suite are incompatible with recent SAS releases, and therefore
utilized a development version of the ESAS suite that has been updated to work
with the latest libraries and compilers.

G351.9$-$0.9 is slightly larger than the central CCD in each MOS camera,
requiring the use of adjacent CCDs. The loss of two of these CCDs in the MOS-1
camera prevents robust exposure correction and background estimation at the
expected source boundary. Similarly, while the PN camera remains the best
instrument for point sources, the growing number of bad pixels and columns –
especially close to the chip boundaries that run through the center of
G351.9$-$0.9 given its position within the field – it is likely not suited to
studies of faint extended emission. In contrast, all of the CCDs in the MOS-2
camera continue to operate, and the detector is relatively free of bad pixels
and columns. Particularly since XMM-Newton does not dither, this is important
to robust imaging of faint emission. For these reasons, our analysis of the
field containing G351.9$-$0.9 is restricted to the MOS-2 camera, although
mosaics including the MOS-1 and PN cameras were constructed and consulted. The
deeper mosaics using the more disparate data did not reveal any further
insight.

We ran the standard SAS task emchain on the raw MOS-2 event list from each
observation. The ESAS tool espfilt was then run on the resulting event lists,
to identify periods of soft proton flaring in each CCD and times when
individual CCDs were operating in “anomalous” states with elevated
backgrounds. One output of espfilt is a GTI file that can be applied to remove
the flagged intervals. The ESAS tool cheese was then run in order to detect
and excise point sources from the data. Given the high column density expected
along the line of sight to G351.9$-$0.9, we ran this tool in soft and hard
bands; bright point sources were primarily detected in the soft band,
indicating that they lie in the foreground rather than within or behind
G351.9$-$0.9. The ESAS tools mosspectra and mosback were run on each event
list to create spectra and images of the quiescent particle background. These
tools rely on calibration files that are essentially archival observations
that were taken with the filter wheel in the “closed” position. It must be
noted that the quiescent particle background is not evenly distributed across
the face of the detector; rather, it is diffuse but concentrated close to the
optical axis, and can easily be mistaken for diffuse emission from the
targeted source. This is particularly true for the MOS cameras, since the
quiescent particle background spectra include strong Al K$\alpha$ and Si
K$\alpha$ lines at 1.49 keV and 1.75 keV, respectively. Even using the ESAS
tools, these lines cannot be entirely removed from resultant images and
spectra, and they are unfortunately similar to lines expected in clusters of
galaxies and supernova remnants. Again owing to the high column density along
the line of sight to G351.9$-$0.9, we did not correct for solar-wind charge
exchange emission, since this contributes most at very low energy and the
remnant should only be evident at higher energy.

Finally, the ESAS tool combimage was used to combine the products created from
each observation, to create summed, filtered, fluxed images, corrected for
exposure, soft proton flaring, the quiescent particle background, and point
sources. These images are shown in Figure 6. The left panel shows the combined
MOS-2 image in the 2.0–7.0 keV band, with a representative set of contours
from the radio image of G351.9$-$0.9 marking its location. In this band, there
is clearly no excess X-ray flux in the region of the remnant. For comparison,
the right panel shows the image in the 0.3–1.0 keV band; again, there is no
indication of enhanced flux in the region of the remnant. We therefore
conclude that the remnant was not detected in X-rays, even in the deep
combined exposure.

In order to estimate an upper limit on the X-ray flux from G351.9$-$0.9, we
extracted spectra from each observation, using standard tools and an
elliptical extraction region containing the remnant. The spectra were then
added using the FTOOLS addascaspec and addrmf, and binned using the “optimal
binning” algorithm of Kaastra & Bleeker (2016). The summed spectrum was fit
over the 0.3-10.0 keV band, minimizing a $\chi^{2}$ statistic. We find that a
model consisting of two thermal plasmas (via “mekal” in XSPEC) with
temperatures of $kT=0.88~{}keV$ and $kT=0.45$ keV, two Gaussian lines centered
at $E=1.49$ keV and $E=1.75$ keV, and a broken power-law ($\Gamma_{1}=0.46$,
$E_{br}=5.7$ keV, $\Gamma_{2}=-3.5$) can describe the flux. This model does
not achieve an acceptable fit ($\chi^{2}=584$ for $\nu=131$ degrees of
freedom), but it does capture the contributions from solar-wind charge
exchange, the quiescent particle background, and residual soft proton flaring.
This model measures a flux of $F=2.2\times 10^{-11}~{}{\rm erg}~{}{\rm
cm}^{-2}~{}{\rm s}^{-1}$ (0.3-10.0 keV), and we therefore set this as an upper
limit on the X-ray flux of G351.9$-$0.9. These results were checked using the
Swift snapshot observation of G351.9$-$0.9. The shape of the astrophysical
background is similar, but the detector backgrounds differ. The background
observed in the much deeper XMM-Newton exposures offers much tighter
constraints despite its complexity.

The interaction between the neutron star X-ray binary Cir X-1 and its
surrounding SNR makes its association clear (Heinz et al., 2013). As a final
check, then, we examined the Chandra OBSID 22289 for evidence of interaction
between Swift J1728.9$-$3613 and its local environment. Using the CIAO tool
specextract, we extracted radial profiles and spectra using 75 annuli
generated between 0–30” for the source and between 35–40” for the background,
after correcting for the ACIS readout streak. These resulting spectra were fit
with a simple model, tbabs$\times$powerlaw model. This model was input into
the ChaRT ray tracing simulation. The resulting rays were input into MARX,
where we simulated an observation with the same exposure time, with the ray
tracing file output by ChaRT. Initial fits to the spectrum of Swift
J1728.9$-$3613 in this observation are suggestive of moderate photon pile-up,
manifested in a hard power-law index and potentially flattened radial profile.
The source flux that was used to generate MARX simulations was based on fits
to the spectrum of Swift J1728.9$-$3613 including the pileup model in XSPEC
(Davis, 2001), which returned a slightly higher flux than models without this
component. The emission profiles that result are broadly consistent, and
indicate that any extended emission in the regions immediately surrounding
Swift J1728.9$-$3613 is due to photon pile-up. There is no evidence of
extension, especially not asymmetric, jet-like extension, that would indicate
an interaction with the remnant.

Figure 7: Anticipated absorbed Sedov emission from G351.9$-$0.9 versus
instrumental and astrophysical backgrounds. The data (in solid black)
represents the observed combination of instrumental and astrophysical
backgrounds observed with XMM-Newton. The models (dashed lines) show the
predicted absorbed Sedov emission from G351.9$-$0.9 as a function of distance
in 0.5 kpc steps, from 4 kpc (top) to 12 kpc (bottom). The model in red
depicts the absorbed emission from $d=7.5$ kpc; at this distance, the expected
emission is marginally detectable, though it peaks close to the instrumental
Al K$\alpha$ emission line at 1.49 keV. We therefore set this as a lower limit
on the distance to G351.9$-$0.9. Please see the text for additional details.

###  4.3 The Distance to G351.9-0.9

Observations have established a relationship between the radio surface
brightness of a SNR, $\Sigma$, and its diameter, $D$. The diameter and angular
size of a given remnant can then be used to determine its distance. Our
MeerKAT observtions yield a total flux density of 1.6–1.8 Jy; assuming major
and minor axes of 12’ and 9’, this translates into a surface brightness of
$\Sigma_{1.28\ GHz}=2.23-2.51\times\ 10^{-21}\rm W\ m^{-2}\ Hz^{-1}\ sr^{-1}$.
Recent statistical treatments of catalogs by Pavlovic et al. (2014) and
Vukotić et al. (2019) have derived improved $\Sigma-D$ relationships. Pavlovic
et al. (2014) employed a probability density function (PDF) to derive the
distance to G351.9$-$0.9, which has the advantage of retaining more
information about the calibration sample than a simple linear fit. Mean,
median, and mode distances of 8.5 kpc, 8.5 kpc, and 7.9 kpc were derived using
the PDF method, broadly consistent with a value of 9.0 kpc derived using an
orthogonal fitting method. For a derived diameter of $D=27.2$ pc, the results
of Xu et al. (2005) give a remnant age of $\tau=15,190$ years. It must be
noted that Pavlovic et al. (2014) do not quote errors for their diameter and
distance estimates, and the differences between the methods may be less than
the uncertainty in each. Vukotić et al. (2019) examine a number of selection
effects and potential sub-samples, and derive a $\Sigma-D$ relationship for
each. Their values are consistent with those derived by Pavlovic et al. (2014)
, but the errors allow for distances as large as 30 kpc.

The fact that G351.9$-$0.9 is likely in the Sedov phase, and therefore likely
dominated by relatively soft X-ray emission, means that the line of sight
column density and observed backgrounds could prevent its detection if it lies
at a sufficiently large distance. This can be translated into an effective
lower limit on the distance to G351.9$-$0.9 that is more constraining than the
$\Sigma-D$ relationship. With the assumption that the total explosion energy
is $10^{51}~{}{\rm erg}$, the Sedov solutions at a given distance are
completely determined by the fixed angular size of the remnant and the density
of the local ISM. We previously developed a characterization of the ISM
density as a function of distance along this line of sight; this also
determines the obscuring column density at a given distance. A set of
distinct, independent Sedov spectra, each determined by the unique conditions
at a given point along the line of sight, can therefore be calculated.

The XSPEC implementation of the Sedov model (Borkowski et al., 2001) captures
the emission in terms of the mean shock temperature, the electron temperature
behind the shock, the ionization age of the remnant ($\tau_{ion}=n\times\tau$,
where $n$ is the local ISM number density and $\tau$ is the remnant age), the
metal abundances (relative to solar), and a flux normalization that is simply
the emission measure diluted by distance. The age of the remnant can be
determined via the equation $r=12.9t_{4}^{2/5}(\epsilon/n)^{1/5}$ pc (where
$r$ is the effective radius at a given distance, $\epsilon$ is the explosion
energy in units of $10^{51}~{}{\rm erg}$, and $t_{4}$ is the remnant age in
$10^{4}$ years; see Cioffi et al. 1988); the ionization age is then also
determined. The mean shock temperature is given by $T_{\rm mean}=3\times
10^{6}t_{4}^{-6/5}(\epsilon/n)^{2/5}$ K (Cioffi et al., 1988), and it is
related to the electron and ion temperatures via
$T_{mean}=\frac{T_{e}+T_{i}}{2}=1.35\times 10^{5}v_{100}^{2}$ and
$\frac{T_{e}}{T_{i}}=\frac{4^{2}}{v_{100}^{2}}$ (where $v_{100}$ is the shock
speed in units of $10^{2}~{}{\rm km}~{}{\rm s}^{-1}$; see Ghavamian et al.
2007).

Within XSPEC, then, we constructed a set of absorbed Sedov spectra
corresponding to 0.5 kpc intervals between 4.0–12.0 kpc. Each spectrum was
simulated using the “fakeit” command, utilizing the observed MOS-2 background
described above and the corresponding MOS-2 response files. In each case, the
simulated source and background exposure times were set to be equivalent to
the total XMM-Newton exposure on G351.9$-$0.9 (184 ks).

Figure 7 shows the results of these simulations. The absorbed spectra are
shown as a sequence of dashed lines. Their emission is very soft, owing to
relatively low electron temperatures. Models corresponding to smaller
distances and less-absorbed Sedov emission spectra lie above the observed
MOS-2 background (the data plotted as connected black points). In contrast,
models corresponding to large distances and correspondingly higher levels of
obscuration lie below the MOS-2 background, and would not have been detected
in the combined 184 ks XMM-Newton exposure. The models that just exceed the
background only do so in the 1.1-1.4 keV band, with the column density
dominating at low energy and the instrumental line dominating at 1.49 keV. The
model in red corresponds to the expected absorption and emission at a distance
of $d=7.5$ kpc. In 184 ks, this model would nominally produce 15,000 counts
above background. On this basis, we suggest that the remnant would nominally
be detected at 7.5 kpc and smaller distances; so $d\geq 7.5$ kpc. At 7.5 kpc,
the remnant would be $\sim$20-25 pc in diameter, which is consistent with our
interpretation of the remnant’s morphology.

Although a slightly larger distance and higher column density might yield a
weaker detection, it is reasonable to be conservative in this circumstance.
The derived distance limit is based on an assumed explosion energy, and solar
abundances. It also hinges on our characterization of the MOS-2 background,
and therefore the efficacy of our efforts to characterize the soft proton
flaring and the accuracy of quiescent particle background calibrations.
Moreover, the accuracy with which emission can be separated from an
instrumental background – particularly a background with strong lines – is
necessarily over-estimated in simulated data. In practice, excluding a broader
range of energies around the instrumental lines would be pragmatic.

##  5 Probability of Coincidental Overlap

Figure 8: We show a histogram of the number of overlaps recorded from our
simulation of 10,000 MW-like galaxies. Orange bars denote runs with no
distance constraints, while purple bars refer to runs where we imposed the
distance constraints we derive in this work. Note that the y-scale is in
logarithmic units. It is clear that a distance constraint significantly
decreases the number of times an overlap occurs between a BH and SNR. On
average, there are 0.007 distance-constrained overlaps in each simulated
galaxy.

To accurately determine the probability of a chance alignment, we use Milky
Way supernova remnant (SNR) and black hole (BH) (Corral-Santana et al., 2016)
catalogs to get positions of all known SNRs and BHs in the Milky Way. Please
see Appendix B for details. We created distributions for galactic latitudes,
longitudes, distances, and angular sizes and sampled from MW catalogs
independently to “simulate” a galaxy. We sampled all parameters separately
(e.g., galactic latitudes for SNRs and then galactic latitudes for BHs) and
then used rejection sampling to draw a new set of values which were combined
with other new parameter values to create a simulated galaxy.

For example, we would take the MW distribution of SNRs Galactic latitudes, and
use rejection sampling to draw a new distribution of SNR latitudes. Then, we
do the same with SNR distances, and all other parameters (other than SNR
angular distance and longitude, which were sampled together), and combine all
the new values to “create” new SNRs and BHs that would have new distances,
longitudes, and latitudes. When looking at these parameters separately, e.g.,
new SNR distances, new SNR latitudes, etc., they would have probability
distributions that are similar to the distribution we see in our galaxy. The
end result was galaxies with populations of BHs and SNRs with l, b, distance,
and size values that had probability distributions that mimic the probability
distributions of these parameters observed in the Milky Way. In our
simulations, an overlap occurred if a BH is found within the central half of
the SNR, i.e. the BH had coordinates within half of the angular height and
width of the SNR on the sky, and it was within 90∘ of the GC (details in
Appendix B). We simulated 10,000 galaxies, and recorded the number of
overlaps. We can constrain the radial distance of the simulated BH and SNR
with the results of this work, that $d_{SNR}>7.5$ kpc (see Section 4.3) and
$d_{BH}=8.4\text{ kpc}\pm 0.8$ kpc (see Section 3.2). When we impose this
distance constraint, the number of chance overlaps decreases significantly, as
seen in Figure 8. Here, we show a histogram of the number of overlaps, with
(purple) and without (orange) our distance constraint (note that the y-axis is
in log-scale). The average number of overlaps per galaxy with constrained
radial distance is 0.007, while without a distance constraint this increases
to an average of 0.52 overlaps per galaxy. In other words, a given MW-type
galaxy has a 0.7% probability of having a BH overlap with the central half of
a SNR within 90∘ of the GC and with our distance constraints.

There are other ways of representing the histogram in Figure 8. If we divide
the number of overlaps by the number of BHs in every simulated galaxy, we can
calculate the probability of a given BH being located at a position covered by
the central half of a SNR (details in Appendix B). The probability of finding
a BH coincident with a SNR’s central half of $P<6.667\%$ without a distance
constraint and $P<1.695\%$ with a distance constraint, where we report the
3$\sigma$ upper limit from our simulations.

The number of known supernova remnants is limited by current observing
capabilities, but new telescopes such as MeerKAT are leading the way in the
search for supernova remnants. As we observe lower flux density supernova
remnants, the number of known supernova remnants will increase. Linearly
extending the current distribution (Green, 2019) of SNRs down to a radio flux
density of 0.001 Jy results in a two-fold increase in the number of SNRs.
Increasing the number of SNRs by a factor of 2 results in an increase in the
average number of overlaps to $P_{avg}\approx 2.5\%$. With these results, we
find it unlikely but not impossible that the overlap of Swift J1728.9$-$313
and G351.9$-$0.9 is coincidental.

##  6 Discussion

We have analyzed X-ray, radio, and IR observations of the transient X-ray
binary Swift J1728.9$-$3613, which lies within the interior of the known
supernova remnant G351.9$-$0.9 in the plane of the sky. Our analysis confirms
and extends independent examinations of Swift J1728.9$-$3613 that conclude
that the X-ray binary system harbors a black hole primary (Enoto et al., 2019;
Saha et al., 2022). We identified an IR counterpart in archival observations
that is spatially coincident with the best X-ray position of Swift
J1728.9$-$3613, and underwent a contemporaneous outburst. The quiescent colors
of this counterpart are broadly consistent with a relatively massive A or
B-type star, likely indicating that Swift J1728.9$-$3613 is an intermediate-
mass X-ray binary (IMXB) or high-mass X-ray binary (HMXB). Utilizing the
measured distribution of gas along the line of sight to Swift J1728.9$-$3613
and G351.9$-$0.9, we derived a distance to the X-ray binary system of
$d=8.4\pm 0.8$ kpc, and a lower limit on the distance to the SNR of $d\geq
7.5$ kpc. We ran extensive simulations to determine the probability of a
coincidental overlap between Swift J1728.9$-$3613 and G351.9$-$0.9, based on
the known numbers and positions of black holes and supernova remnants in the
Milky Way, the distribution of remnant sizes, and these distance constraints.
We find that the probability of coincidental overlap between a black hole and
a remnant is strictly $<1.7\%$. This suggests that the two are likely to be
physically connected.

If the same core-collapse supernova created the black hole in Swift
J1728.9$-$3613 and G351.9$-$0.9, the most fundamental implication is that not
all black holes are “born in the dark.” Rather, distinct formation channels
may exist for black holes. This would represent observational support for
numerical simulations that predict a subset of core-collapse supernovae can
leave black holes and remnants (Ugliano et al., 2012). It is worth repeating
here that SS 433 in W50, MAXI J1535$-$571, and W49B may also represent such
cases (Blundell et al., 2008; Kubota et al., 2010; Cherepashchuk et al., 2018;
Maxted et al., 2020; Lopez et al., 2013). Unlike SS 433 and W49B, however,
X-rays clearly detect a compact object and strongly favor a black hole primary
in Swift J1728.9$-$3613, and for these reasons it may represent the strongest
case for a distinct black hole formation channel. Another fundamental
implication is that the black hole in Swift J1728.9$-$3613 must be very young,
and its high spin ($a=0.86\pm 0.02$) must be set at birth rather than through
accretion. This is the subject of a companion paper; please see Draghis et al.
(2023).

Binary evolution models generally predict that much longer time scales are
required for accretion to commence, than is permitted by the potential youth
of Swift J1728.9$-$3613 and G351.9$-$0.9. For example, Podsiadlowski et al.
(2003) suggests that $\tau\simeq 1$ Gyr is required for mass transfer rates of
$\dot{M}\simeq 10^{-9}M_{\odot}$ per year. While Podsiadlowski et al. (2003)
accurately predicts the mass of the companion and the spin of the black hole
in GRS 1915$+$105, the same considerations predict a companion that is less
massive than the black hole in Cygnus X-1. In contrast, observations find that
the companion is about twice as massive as the black hole (Miller-Jones et
al., 2021). Moreover, Cir X-1 is transient neutron star X-ray binary in a
remnant with an age of $\tau\leq 4600$ years (Heinz et al., 2013), and it also
has an A- or B-type companion star. It may be the case that stronger binary
evolution is possible before supernova events in some instances, leading to
close, circularized binaries.

The binary evolution that is required if Swift J1728.9$-$3613 is physically
associated with G351.9$-$0.9 is broadly consistent with that described by
Maxted et al. (2020). Motivated by the possible association of MAXI
J1535$-$571 with G323.7$-1$.0, Maxted et al. (2020) employed the population
synthesis code “StarTrack” (Belczynski et al., 2008). This code has been used
to model X-ray binary populations in other galaxies; in NGC 1569, it predicts
a luminosity function that is within a factor of two of the observed
luminosity function (Belczynski et al., 2004). StarTrack also made key
predictions connecting short gamma-ray bursts to mergers of BH-BH/BH-NS/NS-NS
binaries that were confirmed by LIGO in 2017 (Abbott et al., 2017). Using
StarTrack, Maxted et al. (2020) found three different binary evolution paths
that could lead to a low-mass X-ray binary or IMXB in a SNR. The evolutionary
paths start with a $\sim 25~{}M_{\odot}$ primary and $\sim 3~{}M_{\odot}$
companion, and result in a system with a 5-11 $M_{\odot}$ black hole and a
companion with a mass of 2.6-4.5 $M_{\odot}$ (an A or B type star) in a very
close orbit, such that mass transfer begins approximately $\tau\simeq 24,000$
years after the supernova. Within the bounds of the data that are currently
available, these predictions broadly match Swift J1728.9$-$3613 and
G351.9$-$0.9.

New numerical simulations are needed to understand the relative rates of black
hole production through different channels, and understand binary evolution
mechanisms that allow for the formation of a highly spinning black hole (see
e.g., Qin et al. 2019). While the absolute number of supernova remnants with
evidence of black holes within the Milky Way is likely to be small, prompt
observations of supernovae in nearby galaxies are beginning to reveal aspects
of compact object formation (Pasham et al., 2021), and may soon provide
another window into black hole formation. In the longer run, the launch of the
Athena mission will deliver sensitive, high-resolution spectra of supernova
remnants in nearby galaxies; comparisons to remnants hosting black holes in
the Milky Way will help to understand black hole formation.

Especially given the broad implications of the potential association between
Swift J1728.9$-$3613 and G351.9$-$0.9, it is vital to critically examine the
evidence that Swift J1728.9$-$3613 harbors a black hole, and the likelihood
that it is physically associated with G351.9$-$0.9. It is equally vital to
examine alternative interpretations of the data. The rest of this section
examines caveats of our analysis, and discusses possible alternatives.

The 5.5 Hz QPOs observed in the NICER power spectra of Swift J1728.9$-$3613
are consistent with Type B QPOs that are typically found in the soft-
intermediate state of black hole X-ray binaries (Belloni et al., 2012).
However, observations and theory do not dictate that these oscillations cannot
be found in neutron stars. The evolution of the X-ray colors observed in NICER
observations may be a stronger indicator of a black hole primary in Swift
J1728.9$-$3613\. The source clearly populates a region that is exclusive to
black hole systems in a prior broad analysis, but two caveats must be attached
to our color analysis: (1) Swift J1728.9$-$3613 has a higher column density
than the ones upon which this diagnostic is based (Done & Gierliński, 2003),
and (2) the NICER data that we employed does not extend above 10 keV. In
principle, both concerns can be mitigated in the limit of very sensitive data;
moreover, we examined three different NICER background estimation methods in
order to check the validity of extrapolating more models to higher energies. A
more exhaustive analysis of X-ray binary colors than undertaken by Done &
Gierliński (2003) and/or the discovery of new transient neutron star X-ray
binaries may reveal a source or sources that populate the “black holes only”
region of the color-color diagram, invalidating this diagnostic. The value of
the spin parameter measured with NuSTAR, $a=0.86\pm 0.02$ (Draghis et al.
2023), may represent the most rigorous and quantitative evidence that Swift
J1728.9$-$3613 harbors a black hole. The error on the spin parameter excludes
the theoretical maximum spin for a neutron star at the $5\sigma$ level of
confidence. While the uncertainties given by Draghis et al. 2023 are simply
statistical, they perform a number of tests in order to assess any systematic
errors of the measurement, and no consideration points to a spin parameter
that is consistent with the range feasible for neutron stars. All measures of
black hole spin that utilize the disk depend upon the assumption that it
extends to the ISCO; NuSTAR observed Swift J1728.9$-$3613 at an Eddington
fraction where this assumption is likely valid (Draghis et al. 2022; also see
Shafee et al. 2008; Reynolds & Fabian 2008). Although the statistical mode of
the black hole spin distributions derived using relativistic disk reflection
and disk continuum spectroscopy are formally consistent
($a_{cont}=0.908^{+0.086}_{-0.558}$, $a_{refl}=0.968^{+0.028}_{-0.240}$,
Draghis et al. 2022), the continuum-derived distribution includes some low
spin values, and the two methods give statistically inconsistent results in a
few cases. It is currently not possible to test the relativistic reflection
results using the disk continuum of Swift J1728.9$-$3613, since the black hole
mass cannot be determined using optical techniques. We have also conducted
systematic searches for coherent pulsations and Type I X-ray bursts in Swift
J1728.9$-$3613, since these phenomena clearly indicate that a given X-ray
binary harbors a neutron star primary. No pulsations are detected in Swift
J1728.9$-$3613, and the upper limits we derive are comparable to the pulsation
amplitudes detected in a number of X-ray binaries. However, it is worth
recalling that a coherent pulsation at $\nu=550.27$ Hz was found in just 150 s
of data in the X-ray binary Aql X-1, amid a total exposure time of 1.3 Ms
(Casella et al., 2008). Similarly, we do not find any Type I X-ray bursts in
an extensive search of the light curves obtained from Swift J1728.9$-$3613,
but our observing campaign could have missed a burst from a source like Cir
X-1 that has a long and/or erratic recurrence time. (Tennant et al., 1986;
Linares et al., 2010). The null results of our searches do not conclusively
rule out a neutron star primary, and more exhaustive observing campaigns
during a future outburst of Swift J1728.9$-$3613 may be able to detect
pulsations or a Type I burst.

Our simulations indicate that a false association between Swift J1728.9$-$3613
and G351.9$-$0.9 is improbable, but not impossible. We have attempted to
anticipate the effect of a higher number of Galactic supernova remnants that
may be discovered in new radio surveys; the probability of a chance
association remains low, but it does grow accordingly. The radio image of
G351.9$-$0.9 is fully consistent with a supernova remnant in the Sedov phase,
and this is bolstered by its surface brightness. However, it is not possible
to completely exclude the possibility that the structure is actually a
superbubble resulting from multiple supernovae and/or massive stellar winds.
Thermal X-ray emission is seen in many superbubbles and some cases appear to
include non-thermal X-ray emission (e.g. Sasaki et al., 2022); if G351.9$-$0.9
is a superbubble, its non-detection with XMM-Newton indicates that it likely
does not lie in the foreground, but instead in the vicinity or background of
Swift J1728.9$-$3613.

We would like to thank the directors and teams at MeerKAT, XMM-Newton, NuSTAR,
Swift, Chandra, and NICER. NICER work at NRL is supported by NASA. The MeerKAT
telescope is operated by the South African Radio Astronomy Observatory, which
is a facility of the National Research Foundation, an agency of the Department
of Science and Innovation. AZ is supported by NASA under award number
80GSFC21M0002. We acknowledge helpful conversations with Teruaki Enoto, Duncan
Galloway, Marko Pavlovic, Steve Snowden, and Branislav Vukotic. The authors
would like to thank the anonymous referees for their comments and suggestions,
which greatly improved the quality of our manuscript.

##  7 Software & Data

The data and code used in our probability of overlap analysis can be found at
this link https://zenodo.org/badge/latestdoi/568961617.

Telescope | OBSID | Start Time (MJD) | Exposure Time (ks) | Target
---|---|---|---|---
NICER | 1200550101 | 58512 | 1.3 | Swift J1728.9$-$3613
NICER | 1200550102 | 58513 | 11.7 | Swift J1728.9$-$3613
NICER | 1200550103 | 58514 | 14.1 | Swift J1728.9$-$3613
NICER | 1200550104 | 58515 | 10.6 | Swift J1728.9$-$3613
NICER | 1200550105 | 58516 | 8.2 | Swift J1728.9$-$3613
NICER | 1200550106 | 58517 | 19.1 | Swift J1728.9$-$3613
NICER | 1200550107 | 58518 | 10.7 | Swift J1728.9$-$3613
NICER | 1200550108 | 58519 | 7.4 | Swift J1728.9$-$3613
NICER | 1200550109 | 58520 | 7.4 | Swift J1728.9$-$3613
NICER | 1200550110 | 58521 | 3.4 | Swift J1728.9$-$3613
NICER | 1200550111 | 58522 | 1.5 | Swift J1728.9$-$3613
NICER | 1200550112 | 58523 | 1.4 | Swift J1728.9$-$3613
NICER | 1200550113 | 58524 | 5.0 | Swift J1728.9$-$3613
NICER | 1200550114 | 58525 | 4.2 | Swift J1728.9$-$3613
NICER | 1200550115 | 58526 | 8.0 | Swift J1728.9$-$3613
NICER | 1200550116 | 58527 | 4.5 | Swift J1728.9$-$3613
NICER | 1200550117 | 58528 | 1.1 | Swift J1728.9$-$3613
NICER | 1200550118 | 58529 | 1.9 | Swift J1728.9$-$3613
NICER | 1200550119 | 58530 | 3.2 | Swift J1728.9$-$3613
NICER | 1200550120 | 58531 | 0.8 | Swift J1728.9$-$3613
NICER | 1200550121 | 58532 | 0.9 | Swift J1728.9$-$3613
NICER | 1200550122 | 58533 | 0.7 | Swift J1728.9$-$3613
NICER | 1200550123 | 58534 | 1.2 | Swift J1728.9$-$3613
NICER | 1200550124 | 58535 | 1.4 | Swift J1728.9$-$3613
NICER | 1200550125 | 58539 | 3.5 | Swift J1728.9$-$3613
NICER | 1200550126 | 58540 | 2.3 | Swift J1728.9$-$3613
NICER | 2200550104 | 58560 | 1.2 | Swift J1728.9$-$3613
NICER | 2200550105 | 58561 | 1.9 | Swift J1728.9$-$3613
NICER | 2200550108 | 58571 | 1.7 | Swift J1728.9$-$3613
NICER | 2200550110 | 58572 | 2.5 | Swift J1728.9$-$3613
NICER | 2200550111 | 58573 | 1.8 | Swift J1728.9$-$3613
NICER | 2200550112 | 58574 | 3.1 | Swift J1728.9$-$3613
NICER | 2200550113 | 58575 | 1.8 | Swift J1728.9$-$3613
NICER | 2200550115 | 58585 | 2.3 | Swift J1728.9$-$3613
NICER | 2200550117 | 58586 | 2.4 | Swift J1728.9$-$3613
NICER | 2200550120 | 58590 | 1.4 | Swift J1728.9$-$3613
NICER | 2200550122 | 58593 | 1.4 | Swift J1728.9$-$3613
NICER | 2200550123 | 58594 | 1.1 | Swift J1728.9$-$3613
NICER | 2200550124 | 58595 | 1.6 | Swift J1728.9$-$3613
NICER | 2200550130 | 58612 | 2.1 | Swift J1728.9$-$3613
NICER | 2200550133 | 58636 | 1.5 | Swift J1728.9$-$3613
NICER | 2200550137 | 58657 | 1.8 | Swift J1728.9$-$3613
NICER | 2200550138 | 58658 | 1.5 | Swift J1728.9$-$3613
NICER | 2200550139 | 58660 | 1.5 | Swift J1728.9$-$3613
NICER | 2200550140 | 58660 | 1.8 | Swift J1728.9$-$3613
NICER | 2200550142 | 58662 | 2.0 | Swift J1728.9$-$3613
NICER | 2200550145 | 58674 | 1.4 | Swift J1728.9$-$3613
NICER | 2200550147 | 58674 | 2.3 | Swift J1728.9$-$3613
NICER | 2200550151 | 58697 | 1.1 | Swift J1728.9$-$3613
NICER | 2200550152 | 58697 | 2.7 | Swift J1728.9$-$3613
Chandra | 21213 | 58526 | 5.1 | Swift J1728.9$-$3613
Chandra | 22289 | 58692 | 30.1 | Swift J1728.9$-$3613
NuSTAR | 90501303002 | 58518 | 16.8 | Swift J1728.9$-$3613
Swift | 00011201001 | 58554 | 1.4 | Swift J1728.9$-$3613
Swift | 00011201002 | 58694 | 2.1 | Swift J1728.9$-$3613
Swift | 00011201003 | 58697 | 3.0 | Swift J1728.9$-$3613
Swift | 00011201004 | 58746 | 2.9 | Swift J1728.9$-$3613
Swift | 00011201005 | 58748 | 1.8 | Swift J1728.9$-$3613
Swift | 00011201006 | 58750 | 1.9 | Swift J1728.9$-$3613
Swift | 00011201007 | 59844 | 1.5 | G351.9$-$0.9
XMM-Newton | 860140101 | 59116 | 18 | G351.9$-$0.9
XMM-Newton | 0901540101 | 59638 | 83 | G351.9$-$0.9
XMM-Newton | 0901540201 | 59642 | 83 | G351.9$-$0.9

Note. — Telescope, OBSIDs, starting days, exposure times, and targets for the
key observations in our analysis. We have omitted observations the remaining
Chandra observations for brevity and clarity.

Source Name | l (deg) | b (deg) |  $N_{H}$ $(10^{22}cm^{-2})$ | Distance from literature (kpc) | Distance from $N_{H}$ (kpc)
---|---|---|---|---|---
V404 Cyg | 73.1 | 2.1 | 6.4 (Kajava et al., 2020) | 2.4 $\pm$ 0.1 (Miller-Jones et al., 2009) | 2.1
GRS 1915+105 | 45 | -0.2 | 3.5 | 9 $\pm$ 3 (Chapuis & Corbel, 2004) | 9.0
SS 433 | 39.7 | -2.2 | 1.76 | 4.9 (Lockman et al., 2007) | 4.7
GRO J1744-28 | 0 | 0.3 | 5 | 8 (Mönkkönen et al., 2019) | 8.6
XTE J1739-285 | 359.7 | 1.3 | 1.73 $\pm$ 0.01 | 4 $\begin{subarray}{c}+4\\\ -2\end{subarray}$ (Bult et al., 2021) | 4.0
IGR J19140+0951 | 44.3 | -0.5 | 6 | 5.2 (Hannikainen et al., 2007) | 5.0
IGR J17544-2619 | 3.2 | -0.3 | 1.2 | 2 – 4 (Pellizza et al., 2006) | 4.2
SAX J1818.6-1703 | 16.9 | -1.3 | 1.17 (Kalberla et al., 2005) | 2.1 (Negueruela et al., 2006) | 3
4U 1728-34 | 354.3 | -0.2 | 2.6 (Worpel et al., 2013) | 5.2 (Galloway et al., 2008) | 5.0
GX 3+1 | 2.3 | 0.8 | 1.91 $\pm$ 0.05 (Rogantini et al., 2019) | 3 – 7 (Kuulkers & van der Klis, 2000) | 4.2
RX J1718.4-4029 | 347.3 | -1.7 | 3.1 $\pm$ 1.3 | 6.5 $\pm$ 0.5 (Cornelisse et al., 2004) | $5.9\begin{subarray}{c}+4.4\\\ -1.9\end{subarray}$
2S 1711-339 | 352.1 | 2.7 | 1.5 |  $<$ 7.5 (Cornelisse et al., 2002) | 3.6
XTE J1720-318 | 354.6 | 1 | 1.33 | 3 – 10 (Chaty & Bessolaz, 2006) | 3.2
XTE J1723-376 | 350.2 |  $-$0.9 | 1.2 (Kalberla et al., 2005) | 1.3 – 10 (Galloway et al., 2008) | 3
SAX J1713.4-4219 | 348.9 | 0.9 | 2 | 4.8 (Smith et al., 2012) | 4.5
4U 1630-47 | 336.9 | 0.3 | 5.8 $\pm$ 0.3 (Tomsick et al., 2014) | 4.7 – 11.5 (Kalemci et al., 2018) | 4.53 – 9
G350.1-0.3 | 350.1 | -0.3 | 3.2 $\pm$ 0.1 | 9 $\pm$ 3 (Yasumi et al., 2014) | $6.5\begin{subarray}{c}+0.5\\\ -0.4\end{subarray}$
Table 1: This table contains references and values for Figure 10, which shows
comparisons between our distance estimation technique and published values. If
no source is cited for the $N_{H}$ value, that value came from the same source
used for the distance. Figure 9: This figure shows how well the simulations in
our probability analysis reproduce known Milky Way distributions for different
populations. From top to bottom, left to right, the panels show plots of b vs.
l (BHs), b vs. l (SNRs), angular height vs. l (SNRs), angular width vs. l
(SNRs), and probability density histograms for BH distances and SNR distances.
In all panels, the original data from MW catalogs is shown in orange, and the
simulated populations are shown in purple. Figure 10: Comparisons of distances
calculated using our technique (purple) of integrating the galactic $n_{H}$
distribution from column densities in the literature with published distances
(orange). Note the breaks in the galactic longitude in order to better
visualize all the points. References are provided in Table 1. The red point
corresponds to distances calculated using this method for the Chandra
measurement of $N_{H}=4.22\pm 0.02\times\rm 10^{22}cm^{-2}$ for Swift J1728.

##  Appendix A Distance Estimation Technique

We tested our column density-based distance calculation on a set of well-known
sources within 5 degrees of the Galactic plane in latitude and within $\sim$
60 degrees of the Galactic Center in longitude. These restrictions are
appropriate for the hydrogen number density distributions that we used to
create a continuous distribution (Kalberla & Kerp, 2009; Marasco et al.,
2017). Our source sample draws heavily from the Galactic Center because
structure within that region, molecular clouds, gas dynamics, superbubbles,
and other factors are potentially important complications that could serve to
indicate the method is not robust.

Table 1 lists the position of each source, the reported column density along
the line of sight, an independent distance estimate, and the distance estimate
derived from the column density using our technique. Figure 10 compares our
technique with independent distance measurements from the literature. The
independent estimates are shown in orange, and the distances estimated using
our technique are shown in purple. The degree of correspondence is striking.
In many cases, there is formal agreement between the techniques; where the
techniques differ, the fractional disagreement is less than 10% in all cases
but one, and in that instance (SAX J1818.6$-$1701) the absolute difference is
less than 1 kpc ($d=2.1$ kpc, Negueruela et al. (2006), versus $d=3$ kpc using
our technique).

##  Appendix B Overlap Simulation

We used Milky Way supernova remnant (SNR) (Green, 2019) and black hole (BH)
(Corral-Santana et al., 2016) catalogs to get information for all known SNRs
and BHs in the galaxy. We omitted SNR G351.9$-$0.9 (our remnant) from the SNR
catalog.

For BHs, we used galactic latitudes, galactic longitudes, and distances from
the catalogs. For the SNRs, we used galactic latitudes, galactic longitudes,
distances, and angular sizes. For each parameter listed, a continuous
probability distribution was approximated by inputting the bin centers and
y-values on histograms into scipy.interpolate.interp1d. We sampled each
probability distribution in order to obtain a new distribution with
approximately the number of those objects in the galaxy using Monte Carlo
rejection sampling. Rejection sampling is a technique where a grid of points
is ’thrown’ on a histogram of a given distribution, and all points below the
curve of the histogram get accepted. A histogram of the new samples will
return a probability distribution very similar to that of the initial data.

It was important to simulate galaxies that had approximately the same number
of objects as the Milky Way: 294 SNRs and 62 BHs. In order to account for
this, we ran our rejection sampling code on each parameter distribution until
we reached an initial sample that would return, on average, the number of
those objects in the Milky Way. We note that many objects did not have
distance measurements, so we had to use more samples to simulate distance
distributions, but we ensured that approximately 294 SNRs and 62 BHs were
generated each time. The code used for this analysis is available at
https://zenodo.org/badge/latestdoi/568961617 with all relevant files.

The angular size of supernova remnants is closely tied to its galactic
longitude (and not as closely to its galactic latitude in our work). Remnants
with larger angular sizes are mostly observed towards the galactic anticenter.
This is physically reasonable; SNRs toward the anti-center are closer to us
than towards the Galactic Center, and will have larger angular sizes, and SNRs
in the galactic centre could also be physically smaller because of the higher
ambient gas density. This means that any coincidental overlaps are more likely
to occur toward the Galactic anti-center, where there are more remnants that
have large angular sizes. In order to account for the dependence between
angular size and Galactic longitude, we sampled from a 2-D probability
distribution, where 1-D rejection sampling was not sufficient. Here, we used
the Walker alias method (Walker, 1974) to randomly sample the distribution.
Using the alias method samples directly from the initial values, rather than
sampling from an initial probability distribution.

In Figure 9, we show an example of a simulated galaxy’s parameters when
compared to our galaxy. The simulation technique used in our code very
accurately samples from known values for SNRs, LMXBs, HMXBs, and BHs. Without
sampling SNR l and size together, we wouldn’t be able to recreate those
distributions when creating our simulated galaxies.

To calculate the probability of overlap in a simulated galaxy, BH galactic l,
b, and distance and SNR l, b, distance, and angular size, we iterated through
each simulated supernova remnant and treated it like a box. If a simulated
BH’s galactic latitude and longitude placed it within half of the supernova
remnant’s angular height and width, and that supernova remnant was located
within $90^{\circ}$ of the Galactic Center, it would count as an overlap. Our
distance constraint required $d_{SNR}>7.5$ kpc and $7.6\text{kpc }<d_{BH}<9.2$
kpc.

Our code simulates 10,000 Milky-Way like galaxies and returns the number of
overlaps, the number of simulated SNRs, and the number of simulated BHs. There
are several ways to convey these results. We can look at the average number of
overlaps (the sum of the overlaps divided by the number of galaxies), the
probability of an SNR being on top of a BH (number of overlaps divided by
number of BHs), or the probability of a BH being inside a SNR (number of
overlaps divided by number of SNRs).

With a distance constraint, the average number of overlaps tells us the
likelihood that a MW-type galaxy would have a BH in the central half of the
SNR where the SNR is at least 7.5 kpc from Earth and the BH is between 7.6 and
9.2 kpc from Earth. There are also other ways to statistically interpret our
results. For a given BH, what is the probability of it overlapping with a SNR?
Or for a given SNR, what is the probability of it appearing to contain a BH?

To demonstrate the difference between these two questions, we can introduce an
analogy. First, imagine a wall with all known SNRs. We take a dart, and throw
it at the wall. The probability that our dart would land within the central
half of a SNR on the wall is analogous to the probability of the central half
of a SNR being on top of a BH. This is not very dependent on the number of BHs
in our sample. However, if we increase the number of SNRs (and the area on the
wall taken up by SNRs), then this probability increases.

Now, imagine the wall has dots on it (locations of galactic BHs). Now, we take
a sticker and throw it at the wall. The probability that the central half of
the sticker would land on a dot on the wall is analogous to the probability of
a BH being seen in the central half of a SNR. If we increase the number of
SNRs in the galaxy, this probability would not change. If we increase the
number of BHs (dots on the wall), then this probability of overlap increases.

Based on our simulations, we can report the 3$\sigma$ upper limits on these
probabilities. The probability of finding a given BH inside the central half
of the SNR of $P<0.37\%$ and a probability on finding the central half of a
given SNR coincident with a BH is $P<1.695\%$. While we cannot definitively
rule out a coincidental overlap in the case of G351.9$-$0.9 and Swift
J1728.9$-$3613, we have determined that the probability of overlap is low. As
mentioned in Section 5, it is certain we have only observed a fraction of the
population of BHs and SNRs in the MW. However, any observational biases in our
simulations are reflected in the observational biases that affect the
detection of these objects.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Acero et al. (2016) Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 224, 8, doi: 10.3847/0067-0049/224/1/8
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Barthelmy et al. (2019) Barthelmy, S. D., D’Ai, A., D’Elia, V., et al. 2019, The Astronomer’s Telegram, 12436, 1
  * Bartlett (1948) Bartlett, M. S. 1948, Nature, 161, 686, doi: 10.1038/161686a0
  * Belczynski et al. (2008) Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223, doi: 10.1086/521026
  * Belczynski et al. (2004) Belczynski, K., Kalogera, V., Zezas, A., & Fabbiano, G. 2004, ApJ, 601, L147, doi: 10.1086/382131
  * Belloni et al. (2012) Belloni, T. M., Sanna, A., & Méndez, M. 2012, MNRAS, 426, 1701, doi: 10.1111/j.1365-2966.2012.21634.x
  * Blundell et al. (2008) Blundell, K. M., Bowler, M. G., & Schmidtobreick, L. 2008, ApJ, 678, L47, doi: 10.1086/588027
  * Borkowski et al. (2001) Borkowski, K. J., Lyerly, W. J., & Reynolds, S. P. 2001, ApJ, 548, 820, doi: 10.1086/319011
  * Bright et al. (2019) Bright, J., Fender, R., Woudt, P., & Miller-Jones, J. 2019, The Astronomer’s Telegram, 12522, 1
  * Bult et al. (2021) Bult, P., Altamirano, D., Arzoumanian, Z., et al. 2021, ApJ, 907, 79, doi: 10.3847/1538-4357/abd54b
  * Bykov et al. (2018) Bykov, A. M., Ellison, D. C., Marcowith, A., & Osipov, S. M. 2018, Space Sci. Rev., 214, 41, doi: 10.1007/s11214-018-0479-4
  * Casella et al. (2008) Casella, P., Altamirano, D., Patruno, A., Wijnands, R., & van der Klis, M. 2008, ApJ, 674, L41, doi: 10.1086/528982
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939, doi: 10.1086/156922
  * Chapuis & Corbel (2004) Chapuis, C., & Corbel, S. 2004, A&A, 414, 659, doi: 10.1051/0004-6361:20034028
  * Chaty & Bessolaz (2006) Chaty, S., & Bessolaz, N. 2006, A&A, 455, 639, doi: 10.1051/0004-6361:20054735
  * Cherepashchuk et al. (2018) Cherepashchuk, A. M., Postnov, K. A., & Belinski, A. A. 2018, MNRAS, 479, 4844, doi: 10.1093/mnras/sty1853
  * Cioffi et al. (1988) Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252, doi: 10.1086/166834
  * Comrie et al. (2021) Comrie, A., Wang, K.-S., Hsu, S.-C., et al. 2021, CARTA: The Cube Analysis and Rendering Tool for Astronomy, 2.0.0, Zenodo, doi: 10.5281/zenodo.4905459
  * Cornelisse et al. (2002) Cornelisse, R., Verbunt, F., in’t Zand, J. J. M., et al. 2002, A&A, 392, 885, doi: 10.1051/0004-6361:20020707
  * Cornelisse et al. (2004) Cornelisse, R., in’t Zand, J. J. M., Kuulkers, E., et al. 2004, Nuclear Physics B Proceedings Supplements, 132, 518, doi: 10.1016/j.nuclphysbps.2004.04.087
  * Corral-Santana et al. (2016) Corral-Santana, J. M., Casares, J., Muñoz-Darias, T., et al. 2016, A&A, 587, A61, doi: 10.1051/0004-6361/201527130
  * Cox (2000) Cox, A. N. 2000, Allen’s astrophysical quantities
  * Dauser et al. (2013) Dauser, T., Garcia, J., Wilms, J., et al. 2013, MNRAS, 430, 1694, doi: 10.1093/mnras/sts710
  * Davis (2001) Davis, J. E. 2001, ApJ, 562, 575, doi: 10.1086/323488
  * Done & Gierliński (2003) Done, C., & Gierliński, M. 2003, MNRAS, 342, 1041, doi: 10.1046/j.1365-8711.2003.06614.x
  * Draghis et al. (2022) Draghis, P. A., Miller, J. M., Zoghbi, A., et al. 2022, arXiv e-prints, arXiv:2210.02479.  https://arxiv.org/abs/2210.02479
  * Enoto et al. (2019) Enoto, T., Altamirano, D., Bult, P. M., et al. 2019, The Astronomer’s Telegram, 12455, 1
  * Fender et al. (2016) Fender, R., Woudt, P. A., Corbel, S., et al. 2016, in MeerKAT Science: On the Pathway to the SKA, 13.  https://arxiv.org/abs/1711.04132
  * Fesen et al. (2018) Fesen, R. A., Weil, K. E., Cisneros, I. A., Blair, W. P., & Raymond, J. C. 2018, MNRAS, 481, 1786, doi: 10.1093/mnras/sty2370
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. D. R. Silva & R. E. Doxsey, 62701V, doi: 10.1117/12.671760
  * Gabriel et al. (2004) Gabriel, C., Denby, M., Fyfe, D. J., et al. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 314, Astronomical Data Analysis Software and Systems (ADASS) XIII, ed. F. Ochsenbein, M. G. Allen, & D. Egret, 759
  * Galloway et al. (2008) Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360, doi: 10.1086/592044
  * García et al. (2013) García, J., Dauser, T., Reynolds, C. S., et al. 2013, ApJ, 768, 146, doi: 10.1088/0004-637X/768/2/146
  * Ghavamian et al. (2007) Ghavamian, P., Laming, J. M., & Rakowski, C. E. 2007, ApJ, 654, L69, doi: 10.1086/510740
  * Green (2019) Green, D. A. 2019, VizieR Online Data Catalog, VII/284
  * Hannikainen et al. (2007) Hannikainen, D. C., Rawlings, M. G., Muhli, P., et al. 2007, MNRAS, 380, 665, doi: 10.1111/j.1365-2966.2007.12092.x
  * Heinz et al. (2013) Heinz, S., Sell, P., Fender, R. P., et al. 2013, ApJ, 779, 171, doi: 10.1088/0004-637X/779/2/171
  * Heywood (2020) Heywood, I. 2020, oxkat: Semi-automated imaging of MeerKAT observations.  http://ascl.net/2009.003
  * Huppenkothen et al. (2019) Huppenkothen, D., Bachetti, M., Stevens, A. L., et al. 2019, apj, 881, 39, doi: 10.3847/1538-4357/ab258d
  * Huppenkothen et al. (2019) Huppenkothen, D., Bachetti, M., Stevens, A., et al. 2019, Journal of Open Source Software, 4, 1393, doi: 10.21105/joss.01393
  * Kaastra & Bleeker (2016) Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151, doi: 10.1051/0004-6361/201527395
  * Kajava et al. (2020) Kajava, J. J. E., Sánchez-Fernández, C., Alfonso-Garzón, J., Motta, S. E., & Veledina, A. 2020, A&A, 634, A94, doi: 10.1051/0004-6361/201937191
  * Kalberla et al. (2005) Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775, doi: 10.1051/0004-6361:20041864
  * Kalberla & Kerp (2009) Kalberla, P. M. W., & Kerp, J. 2009, ARA&A, 47, 27, doi: 10.1146/annurev-astro-082708-101823
  * Kalemci et al. (2018) Kalemci, E., Maccarone, T. J., & Tomsick, J. A. 2018, ApJ, 859, 88, doi: 10.3847/1538-4357/aabcd3
  * Kennea (2019) Kennea, J. A. 2019, The Astronomer’s Telegram, 12445, 1
  * Kenyon et al. (2018) Kenyon, J. S., Smirnov, O. M., Grobler, T. L., & Perkins, S. J. 2018, MNRAS, 478, 2399, doi: 10.1093/mnras/sty1221
  * Knowles et al. (2022) Knowles, K., Cotton, W. D., Rudnick, L., et al. 2022, A&A, 657, A56, doi: 10.1051/0004-6361/202141488
  * Kubota et al. (2010) Kubota, K., Ueda, Y., Kawai, N., et al. 2010, PASJ, 62, 323, doi: 10.1093/pasj/62.2.323
  * Kurtzer et al. (2017) Kurtzer, G. M., Sochat, V., & Bauer, M. W. 2017, PLoS ONE, 12, e0177459, doi: 10.1371/journal.pone.0177459
  * Kuulkers & van der Klis (2000) Kuulkers, E., & van der Klis, M. 2000, A&A, 356, L45.  https://arxiv.org/abs/astro-ph/0003180
  * Linares et al. (2010) Linares, M., Watts, A., Altamirano, D., et al. 2010, ApJ, 719, L84, doi: 10.1088/2041-8205/719/1/L84
  * Lo & Lin (2011) Lo, K.-W., & Lin, L.-M. 2011, ApJ, 728, 12, doi: 10.1088/0004-637X/728/1/12
  * Lockman et al. (2007) Lockman, F. J., Blundell, K. M., & Goss, W. M. 2007, MNRAS, 381, 881, doi: 10.1111/j.1365-2966.2007.12170.x
  * Lopez et al. (2013) Lopez, L. A., Ramirez-Ruiz, E., Castro, D., & Pearson, S. 2013, ApJ, 764, 50, doi: 10.1088/0004-637X/764/1/50
  * Makhathini (2018) Makhathini, S. 2018, PhD thesis, Rhodes University, South Africa
  * Marasco et al. (2017) Marasco, A., Fraternali, F., van der Hulst, J. M., & Oosterloo, T. 2017, A&A, 607, A106, doi: 10.1051/0004-6361/201731054
  * Mauch et al. (2020) Mauch, T., Cotton, W. D., Condon, J. J., et al. 2020, ApJ, 888, 61, doi: 10.3847/1538-4357/ab5d2d
  * Maxted et al. (2020) Maxted, N. I., Ruiter, A. J., Belczynski, K., Seitenzahl, I. R., & Crocker, R. M. 2020, arXiv e-prints, arXiv:2010.15341.  https://arxiv.org/abs/2010.15341
  * McMullin et al. (2007a) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007a, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * McMullin et al. (2007b) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007b, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Miller et al. (2017) Miller, J. M., Bautz, M. W., & McNamara, B. R. 2017, ApJ, 850, L3, doi: 10.3847/2041-8213/aa9566
  * Miller et al. (2018) Miller, J. M., Gendreau, K., Ludlam, R. M., et al. 2018, ApJ, 860, L28, doi: 10.3847/2041-8213/aacc61
  * Miller & Miller (2015) Miller, M. C., & Miller, J. M. 2015, Phys. Rep., 548, 1, doi: 10.1016/j.physrep.2014.09.003
  * Miller-Jones et al. (2009) Miller-Jones, J. C. A., Jonker, P. G., Dhawan, V., et al. 2009, ApJ, 706, L230, doi: 10.1088/0004-637X/706/2/L230
  * Miller-Jones et al. (2021) Miller-Jones, J. C. A., Bahramian, A., Orosz, J. A., et al. 2021, Science, 371, 1046, doi: 10.1126/science.abb3363
  * Mirabel & Rodrigues (2003) Mirabel, I. F., & Rodrigues, I. 2003, Science, 300, 1119, doi: 10.1126/science.1083451
  * Mitsuda et al. (1984) Mitsuda, K., Inoue, H., Koyama, K., et al. 1984, PASJ, 36, 741
  * Mönkkönen et al. (2019) Mönkkönen, J., Tsygankov, S. S., Mushtukov, A. A., et al. 2019, A&A, 626, A106, doi: 10.1051/0004-6361/201935507
  * Nasa High Energy Astrophysics Science Archive Research Center (2014) (Heasarc) Nasa High Energy Astrophysics Science Archive Research Center (Heasarc). 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, record ascl:1408.004.  http://ascl.net/1408.004
  * Negoro et al. (2019) Negoro, H., Nakahira, S., Nakajima, M., et al. 2019, The Astronomer’s Telegram, 12437, 1
  * Negueruela et al. (2006) Negueruela, I., Smith, D. M., Reig, P., Chaty, S., & Torrejón, J. M. 2006, in ESA Special Publication, Vol. 604, The X-ray Universe 2005, ed. A. Wilson, 165.  https://arxiv.org/abs/astro-ph/0511088
  * Offringa et al. (2014) Offringa, A. R., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606, doi: 10.1093/mnras/stu1368
  * Pasham et al. (2021) Pasham, D. R., Ho, W. C. G., Alston, W., et al. 2021, Nature Astronomy, doi: 10.1038/s41550-021-01524-8
  * Patruno & Watts (2021) Patruno, A., & Watts, A. L. 2021, in Astrophysics and Space Science Library, Vol. 461, Astrophysics and Space Science Library, ed. T. M. Belloni, M. Méndez, & C. Zhang, 143–208, doi: 10.1007/978-3-662-62110-3_4
  * Pavlovic et al. (2014) Pavlovic, M. Z., Dobardzic, A., Vukotic, B., & Urosevic, D. 2014, Serbian Astronomical Journal, 189, 25, doi: 10.2298/SAJ1489025P
  * Pellizza et al. (2006) Pellizza, L. J., Chaty, S., & Negueruela, I. 2006, A&A, 455, 653, doi: 10.1051/0004-6361:20054436
  * Podsiadlowski et al. (2003) Podsiadlowski, P., Rappaport, S., & Han, Z. 2003, MNRAS, 341, 385, doi: 10.1046/j.1365-8711.2003.06464.x
  * Qin et al. (2019) Qin, Y., Marchant, P., Fragos, T., Meynet, G., & Kalogera, V. 2019, ApJ, 870, L18, doi: 10.3847/2041-8213/aaf97b
  * Ransom et al. (2002) Ransom, S. M., Eikenberry, S. S., & Middleditch, J. 2002, AJ, 124, 1788, doi: 10.1086/342285
  * Raymond et al. (2020) Raymond, J. C., Chilingarian, I. V., Blair, W. P., et al. 2020, ApJ, 894, 108, doi: 10.3847/1538-4357/ab886d
  * Remillard et al. (2021) Remillard, R. A., Loewenstein, M., Steiner, J. F., et al. 2021, arXiv e-prints, arXiv:2105.09901.  https://arxiv.org/abs/2105.09901
  * Remillard et al. (2022) —. 2022, AJ, 163, 130, doi: 10.3847/1538-3881/ac4ae6
  * Reynolds (2021) Reynolds, C. S. 2021, ARA&A, 59, doi: 10.1146/annurev-astro-112420-035022
  * Reynolds & Fabian (2008) Reynolds, C. S., & Fabian, A. C. 2008, ApJ, 675, 1048, doi: 10.1086/527344
  * Rogantini et al. (2019) Rogantini, D., Costantini, E., Zeegers, S. T., et al. 2019, A&A, 630, A143, doi: 10.1051/0004-6361/201935883
  * Saha et al. (2022) Saha, D., Mandal, M., & Pal, S. 2022, arXiv e-prints, arXiv:2210.13748.  https://arxiv.org/abs/2210.13748
  * Saji et al. (2018) Saji, S., Matsumoto, H., Nobukawa, M., et al. 2018, PASJ, 70, 23, doi: 10.1093/pasj/psx158
  * Sasaki et al. (2022) Sasaki, M., Knies, J., Haberl, F., et al. 2022, A&A, 661, A37, doi: 10.1051/0004-6361/202141054
  * Shafee et al. (2008) Shafee, R., Narayan, R., & McClintock, J. E. 2008, ApJ, 676, 549, doi: 10.1086/527346
  * Shternin et al. (2011) Shternin, P. S., Yakovlev, D. G., Heinke, C. O., Ho, W. C. G., & Patnaude, D. J. 2011, MNRAS, 412, L108, doi: 10.1111/j.1745-3933.2011.01015.x
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Smith et al. (2012) Smith, D. M., Markwardt, C. B., Swank, J. H., & Negueruela, I. 2012, MNRAS, 422, 2661, doi: 10.1111/j.1365-2966.2012.20836.x
  * Snowden et al. (2004) Snowden, S. L., Collier, M. R., & Kuntz, K. D. 2004, ApJ, 610, 1182, doi: 10.1086/421841
  * Snowden et al. (2008) Snowden, S. L., Mushotzky, R. F., Kuntz, K. D., & Davis, D. S. 2008, A&A, 478, 615, doi: 10.1051/0004-6361:20077930
  * Tennant et al. (1986) Tennant, A. F., Fabian, A. C., & Shafer, R. A. 1986, MNRAS, 221, 27P, doi: 10.1093/mnras/221.1.27P
  * Tomsick et al. (2014) Tomsick, J. A., Yamaoka, K., Corbel, S., et al. 2014, ApJ, 791, 70, doi: 10.1088/0004-637X/791/1/70
  * Ugliano et al. (2012) Ugliano, M., Janka, H.-T., Marek, A., & Arcones, A. 2012, ApJ, 757, 69, doi: 10.1088/0004-637X/757/1/69
  * Vaughan et al. (1994) Vaughan, B. A., van der Klis, M., Wood, K. S., et al. 1994, ApJ, 435, 362, doi: 10.1086/174818
  * Veena et al. (2019) Veena, V. S., Vig, S., Sebastian, B., et al. 2019, MNRAS, 482, 4630, doi: 10.1093/mnras/sty3032
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vukotić et al. (2019) Vukotić, B., Ćiprijanović, A., Vučetić, M. M., Onić, D., & Urošević, D. 2019, Serbian Astronomical Journal, 199, 23, doi: 10.2298/SAJ1999023V
  * Walker (1974) Walker, A. J. 1974, Electronics Letters, 10, 127
  * Whiteoak & Green (1996) Whiteoak, J. B. Z., & Green, A. J. 1996, A&AS, 118, 329
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016
  * Worpel et al. (2013) Worpel, H., Galloway, D. K., & Price, D. J. 2013, ApJ, 772, 94, doi: 10.1088/0004-637X/772/2/94
  * Xu et al. (2005) Xu, J.-W., Zhang, X.-Z., & Han, J.-L. 2005, Chinese J. Astron. Astrophys., 5, 165, doi: 10.1088/1009-9271/5/2/007
  * Yasumi et al. (2014) Yasumi, M., Nobukawa, M., Nakashima, S., et al. 2014, PASJ, 66, 68, doi: 10.1093/pasj/psu043