11institutetext: Universidade do Estado de Santa Catarina, Joinville,
89219-710, SC, Brazil
11email<EMAIL_ADDRESS>22institutetext: Nicolaus Copernicus
Astronomical Center, Polish Academy of Sciences, Warsaw, 00-716, Poland
22email<EMAIL_ADDRESS>33institutetext: Núcleo de Astrofísica e
Cosmologia (Cosmo-Ufes) & Departamento de Física, Universidade Federal do
Espírito Santo, Vitória, 29075-910, ES, Brazil
33email<EMAIL_ADDRESS>44institutetext: Divisão de Astrofísica,
Instituto Nacional de Pesquisas Espaciais, São José dos Campos, 12227-010, SP,
Brazil 55institutetext: Instituto de Física, Universidade Federal do Rio
Grande do Sul, Porto Alegre, 91501-970, RS, Brazil 66institutetext: Instituto
de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo,
São Paulo, 05508–900, SP, Brazil 77institutetext: RECOD.ai, Institute of
Computing, University of Campinas, Campinas, 13083-852, SP, Brazil
88institutetext: INFN Sezione di Ferrara, Ferrara, 44122, Italy

# Evidence for 3XMM J185246.6+003317 as a massive magnetar with a low magnetic
field

Rafael C. R. de Lima 11 Jonas P. Pereira 2233 Jaziel G. Coelho 3344 Rafael C.
Nunes 5544 Paulo E. Stecchini 4466 Manuel Castro 77 Pierre Gomes 11 Rodrigo R.
da Silva 44 Claudia V. Rodrigues 44 José C. N. de Araujo 44 Michał Bejger 2288
Paweł Haensel 22 J. Leszek Zdunik 22

3XMM J185246.6+003317 is a slowly rotating soft-gamma repeater (neutron star)
in the vicinity of the supernova remnant Kes 79. So far, observations have
only set upper limits to its surface magnetic field and spindown, and there is
no estimate for its mass and radius. Using ray-tracing modelling and Bayesian
inference for the analysis of several light curves spanning a period of around
three weeks, we have found that it may be one of the most massive neutron
stars to date. In addition, our analysis suggests a multipolar magnetic field
structure with a subcritical field strength and a carbon atmosphere
composition. Due to the time-resolution limitation of the available light
curves, we estimate the surface magnetic field and the mass to be
$\log_{10}(B/{\rm G})=11.89^{+0.19}_{-0.93}$ and $M=2.09^{+0.16}_{-0.09}$
$M_{\odot}$ at $1\sigma$ confidence level, while the radius is estimated to be
$R=12.02^{+1.44}_{-1.42}$ km at $2\sigma$ confidence level. They were verified
by simulations, i.e., data injections with known model parameters, and their
subsequent recovery. The best-fitting model has three small hot spots, two of
them in the southern hemisphere. These are, however, just back-of-the-envelope
estimates and conclusions, based on a simple ray-tracing model with
anisotropic emission; we also estimate the impact of modelling on the
parameter uncertainties and the relevant phenomena on which to focus in more
precise analyses. We interpret the above best-fitting results as due to
accretion of supernova layers/interstellar medium onto 3XMM J185246.6+003317
leading to burying and a subsequent re-emergence of the magnetic field, and a
carbon atmosphere being formed possibly due to hydrogen/helium diffusive
nuclear burning. Finally, we briefly discuss some consequences of our findings
for superdense matter constraints.

##  1 Introduction

The X-ray pulsar 3XMM J185246.6+003317 (hereafter 3XMM J1852+0033) was
discovered in the field-of-view of an XMM-Newton observation of the supernova
remnant Kes 79, which hosts a central compact object (CCO) (Seward et al.,
2003). These observations were independently analysed by Zhou et al. (2014)
and Rea et al. (2014), who reported a bright point-like source, located 7.4′
away from the CCO just outside the southern boundary of Kes 79, having a very
prominent periodic modulation of $P{\approx}11.6$ s. The source is increasing
its period at a rate $\dot{P}<1.4\times 10^{-13}$ s/s, and its X-ray
luminosity is higher than its spin-down luminosity, ruling out a rotation-
powered nature (see, e.g., Coelho et al., 2017). This implies a surface
dipolar magnetic field strength $B<4.1\times 10^{13}$ G, and a characteristic
age $\tau_{\mathrm{age}}>1.3$ Myr. The foreground absorption ($N_{H}$) toward
3XMM J1852+0033 is similar to that of Kes 79, suggesting a similar distance of
${\sim}7.1$ kpc (see Sun et al., 2004; Zhou et al., 2014, 2016, for details).
3XMM J1852+0033 has been classified within the Soft Gamma Repeaters (SGRs) and
the Anomalous X-ray Pulsars (AXPs) class, usually called magnetars, which are
neutron stars (NSs) characterized by a quiescent soft X-ray ($2-10$ keV)
luminosity of the order of $10^{30}-10^{35}$ erg/s, spin period in the range
$2-12$ s, and a spindown rate from $10^{-15}$ to $10^{-10}$ s/s (see, e.g.,
Olausen & Kaspi, 2014; Turolla et al., 2015; Kaspi & Beloborodov, 2017). In
particular, the low dipolar magnetic field inferred from the spindown suggests
that this source is a SGR with low-B (Rea et al., 2014).

Since SGRs/AXPs are usually isolated NSs, inferring their macroscopic
properties is a complicated task; it is usually assumed that they have a
canonical mass of $1.4\,M_{\odot}$. However, the existence of high-mass NSs is
also well-known (in binaries): PSR J1614–2230 has a mass $M=1.97\pm 0.04$
$M_{\odot}$ (Demorest et al., 2010); the mass of PSR J0348+0432 is $2.01\pm
0.04$ $M_{\odot}$ (Antoniadis et al., 2013), and PSR J0740+6620 has an
estimated mass of $2.08\pm 0.07$ $M_{\odot}$ (Fonseca et al., 2021).
Observationally, the probability distribution for known NS masses presents a
two component Gaussian mixture model, with mean values around
$1.34\,M_{\odot}$ and $1.8\,M_{\odot}$, and standard deviations of
${\sim}0.1\,M_{\odot}$ (see Alsing et al., 2018, for details).

This paper reports on our analysis of the XMM-Newton X-ray data of 3XMM
J1852+0033 using ray-tracing modelling. We assume that its X-ray pulse profile
is due to the emission of hot spots on its surface. Our main motivations are:
(i) magnetars are usually slowly rotating NSs in which ray-tracing models are
simple and relativistic phenomena are important; (ii) as far as we are aware
of, there are no inferences of masses and radii of magnetars by means of ray-
tracing techniques. Therefore, ray-tracing allows us to learn more about the
physics taking place around magnetars (and any other class of NSs) and
independently complement what other phenomena/models already tell us about
them.

In Sect. 2 we describe aspects of the data selection and its reduction.
Section 3 is devoted to the description of the model and the parameter
estimation. The best-fitting results for the pulse profile of 3XMM J1852+0033
are presented and discussed in Sect. 4. Finally, Sec. 5 summarizes our main
findings. Appendix A contains details about the pulse profile model. Further
details about the atmosphere models for NSs are given in Appendix B.

##  2 Data selection, reduction, and preparation

###  2.1 The Observations

The field around 3XMM J1852+0033 was monitored by XMM-Newton on several
occasions from 2004 to 2009. The source entered a bright state at some time
before 2009; it is unclear whether it went into a quiescent stage in 2009 (see
Zhou et al., 2014, and references therein). We chose to retrieve five
observations during the bright state, namely ObsIDs 0550670201 (2008 Sep. 19),
0550670301 (2008 Sep. 21), 0550670401 (2008 Sep. 23), 0550670501 (2008 Sep.
29) and 0550670601 (2008 Oct. 10), which we refer to as epochs A, B, C, D and
E, respectively. The above choice is mainly due to data quality and due to
reported characteristic timescales for significant hot spot motions. With
relation to data quality, we have that for epochs away from the outburst the
net source count rate is approximately one order of magnitude lower (than
during it). Regarding hot spot motions, we made use of knowledge stemming from
observations of other magnetars. Careful monitoring of SGR 1830-0645 by NICER,
for instance, has shown that characteristic timescales are on the order of a
month (Younes et al., 2022), suggesting that fixed hot spots for smaller
timescales would be a reasonable approximation. Hence, observations spanning
no more than a couple of weeks may be combined, which increases the quantity
of data to be fit, and improves the quality of the statistics.

###  2.2 Science products extraction

During the observations, the instrument EPIC-pn (Strüder et al., 2001) was
operating in small window mode and did not have 3XMM J1852+0033 within its
field-of-view. EPIC-MOS (Turner et al., 2001) cameras were operating in full
window mode and did contain 3XMM J1852+0033; however, for MOS 1, the source
fell into a CCD that was switched off in two of the observations. Thus, we
only make use of data from MOS 2. Standard data reduction and filtering
procedures were conducted with the XMM-Newton Science Analysis System (SAS,
v.19.1.0).

Source photons were extracted, for all observations, from circular regions of
40′′ centred around the object’s position. The background regions were chosen
with the aid of the SAS task ebkgreg. Because this task indicates the optimal
background region based solely on the detector geometry, the regions suggested
were slightly shrunk to avoid the inclusion of source photons (see Figure 1).

The barycen and epiclccorr tasks were applied to extract the light curves. The
former converts the photon arrival registered time into the solar system
barycenter time reference and the latter performs a series of corrections to
minimize effects that may impact the detection efficiency (e.g. dead time,
chip gaps, point-spread-function variation111
https://heasarc.gsfc.nasa.gov/docs/xmm/sas/help/epiclccorr/) before producing
a background-subtracted light curve. The timing resolution was limited by the
camera’s operation mode, that is 2.6 seconds. Light curves were extracted in
two energy bands: 0.3–10 and 3–8 keV. More details are given in the next
sections.

Although spectral analysis is not in the scope of this study, in order to
obtain information on the source’s flux during each observation, we also
extracted source and background spectra for the same aforementioned energy-
band regions. Standard tasks rmfgen and arfgen were used to create the
redistribution matrix file (RMF) and the ancillary response file (ARF).

Figure 1: MOS 2 raw image for observation 0550670201. The image LUT is in log-
scale; units are counts per pixel. Source and background extraction regions
for 3XMM J1852+0033 are indicated.

###  2.3 Phase folded light curves production

We folded the light curves (pulse profiles) at the periods provided by Lomb-
Scargle periodograms computed for each observations’ light curves. The period
values (P $\sim$ 11.56 s) found for each set differ only beyond the fourth
decimal place, as already shown by Zhou et al. (2014) and Rea et al. (2014).
They also pointed out that there is no significant time derivative amongst
observations. This is true for either energy bands (0.3–10 and 3–8 keV). The
folded light curves were binned to have 50 bins/phase, which will be the main
input to our model. For comparative testing, we also produced folded light
curves with 16 bins/phase.

To derive the flux emitted by 3XMM J1852+0033 during each observation, we
fitted the 0.3–10 keV spectra with a simple absorbed blackbody model
(phabs*bbody). We used the X-ray spectral fitting package XSPEC (Arnaud,
1996), version 12.11.1. The overall values of the equivalent hydrogen column
($N_{\text{H}}$ $\sim$ 1.2–1.5 $\times$ 1022 cm-2) and blackbody temperature
($kT$ $\sim$ 70–80 eV) provided by the best fits are in agreement with those
reported by Zhou et al. (2014). The unabsorbed fluxes (0.3–10 keV) for each
observation were computed by the flux command after setting the absorption
component to zero. To explore how the choice of model may affect the
unabsorbed flux obtained, we added a phenomenological power-law to the
previous model, i.e., the combination phabs*(bbody+powerlaw). The blackbody
model parameters do not change, considering the confidence range. The hydrogen
column tends to assume slightly higher values (from 5% to 15%), but it is also
less constrained, agreeing within 90% error with the blackbody model alone;
the powerlaw index is not well constrained. These variations in the absorption
component value impact the 0.3–10 keV unabsorbed flux obtained accordingly.
For further comparative testing, we have also computed the unabsorbed fluxes
in the 3–8 keV band. In this band, the role of the hydrogen column is much
smaller, and the unabsorbed fluxes derived for the two models agree within 2%.

Our modelling of the folded light curves relies on the shape of each curve,
meaning that it is able to fit the curves regardless of the choice of
normalization. Nevertheless, as it will be clarified in the following
sections, using a sort of representative value of the source’s emission, such
as the unabsorbed fluxes, may help improve the quality of the results.

##  3 Model description and parameter estimation

In our previous work (de Lima et al., 2020), as well as in the original
reference that introduced the method (Turolla & Nobili, 2013), the
calculations were performed using the normalized flux, defined as
$(F_{max}+F_{min})/2$. However, in the present work, we make use of a modified
normalization, $(F_{max}+F_{min})/2\bar{F}$, where $\bar{F}$ denotes the mean
value of the observed unabsorbed flux in the period, derived from the spectral
fitting. This modification is of great significance because normalizing the
theoretical pulse to the correct level improves the statistical sampling and
prevents the acceptance of parameters that describe unrealistic situations,
which do not correspond to the actual flux. Details are given in Appendix A.

The spacetime outside the star is well-described by the Schwarzschild metric,
i.e., we neglect rotational effects, as the source is a slow rotator. We also
make use of the spectral tables for models of magnetic NS atmospheres
calculated by Ho et al. (2007, 2008) and Mori & Ho (2007), which have been
made public and implemented in the XSPEC package (see Arnaud, 1996; Ho &
Heinke, 2009; Ho, 2013). These neutron star atmosphere models were obtained
using the most up-to-date equation of state (EOS) and opacity results for a
partially ionized, strongly magnetized hydrogen or mid-Z element plasma (e.g.,
carbon, oxygen, and neon). The associated spectra come from the solution to
the coupled radiative transfer equations for the two photon polarization modes
in a magnetized medium. The atmosphere is assumed to be in radiative and
hydrostatic equilibrium. The atmosphere models depend on the NS surface
temperature, the NS surface magnetic field strength $B$ and its inclination
relative to the radial direction. Moreover, there is a dependence on the NS
surface gravity $g=(1+z_{g})GM/R^{2}$, where the gravitational redshift is
$1+z_{g}=1/\sqrt{1-2GM/Rc^{2}}$. More precisely, the tables used for the
spectra of the neutron star atmospheres were calculated using discrete
magnetic fields: $[10^{10},\leavevmode\nobreak\ 10^{11},\leavevmode\nobreak\
10^{12},\leavevmode\nobreak\ 3.0\times 10^{13}]$ G (nsmaxg, for an H
atmosphere), and $[10^{12}$, $10^{13}]$ G (nsmaxg, for a C atmosphere). For
low magnetic fields ($B\leq 10^{10}$ G), one can use nsx for a non-magnetic
atmosphere (H, He, and Ca). The tables also depend on discrete values of the
surface gravity $g$. For a C atmosphere, we adopted
$\log_{10}(g\leavevmode\nobreak\ [\rm{cm/s^{2}}])=[13.6-15.3]$, and for a H
atmosphere, $\log_{10}(g\leavevmode\nobreak\ [\rm{cm/s^{2}}])=[13.6-15]$. To
be able to consider continuous values of the model parameters in the Bayesian
analysis, the model emission has been interpolated using the tabulated
emission values. Further details about the atmosphere models for NSs are given
in Appendix B.

We used the Markov Chain Monte Carlo (MCMC) method to determine the set of
parameters $\bar{\theta}_{i}$ that best describes the 3XMM J1852+0033 light
curves. For the case of three spots, which is the maximum number of spots we
considered in the modelling, they are

$\displaystyle\bar{\theta}_{i}=\Big{\\{}$ $\displaystyle
M,\,R,\,\log(B/G),\,i,\,\theta_{1},\,\theta_{2},\,\phi_{2},\,\theta_{3},\,\phi_{3},$
(1)
$\displaystyle\theta_{c1},\,T_{1},\,\theta_{c2},\,T_{2},\,\theta_{c3},\,T_{3}\Big{\\}},$

with the parameters separated into two categories: the fixed and the secular
parameters. The fixed parameters are presented in the first row and do not
change for the epochs investigated (NS mass $M$ and radius $R$, its magnetic
field $B$, position of the hot spots $\theta_{i}$ and $\phi_{i}$ $(i=1,2,3)$
and, for reference, $\phi_{1}\equiv 0$). The secular parameters may have
different values for each epoch (sizes and temperatures of hot spots,
$\theta_{ci}$ and $T_{i}$). Table 1 summarizes their definitions and prior
range distributions. The posterior probability distribution function is

$p({\cal D}\mid\bar{\theta})\propto\exp\Big{(}-\frac{1}{2}\chi^{2}\Big{)}\,,$
(2)

with

$\chi^{2}=\sum_{i=1}^{N}(F^{obs}-F^{th}(\bar{\theta}_{i}))^{T}\Sigma^{-1}(F^{obs}-F^{th}(\bar{\theta}_{i})),$
(3)

where $F^{obs}$, $F^{th}$ and $\Sigma^{-1}$ are the observational flux as
explained in Appendix A, the total theoretical flux produced by the spots, and
the square N-by-N covariance matrix associated with the observed flux,
respectively. The goal of any MCMC approach is to draw $N$ samples
$\bar{\theta}_{i}$ from the general posterior probability density

$p(\bar{\theta}_{i},\alpha\mid{\cal
D})=\frac{1}{Z}p(\bar{\theta},\alpha)p({\cal D}\mid\bar{\theta},\alpha)\,,$
(4)

where $p(\bar{\theta},\alpha)$ and $p({\cal D}\mid\bar{\theta},\alpha)$ are
the prior distribution and the likelihood function, respectively. Here, ${\cal
D}$ and $\alpha$ are the set of observations and possible nuisance parameters
and $Z$ is the normalization factor.

We perform the statistical analysis based on the emcee algorithm (Foreman-
Mackey et al., 2013), assuming the theoretical model described in Appendix A
and the uniform priors summarized in Table 1. During our analysis, we
discarded the first 20% steps of the chain as burn-in. We follow the Gelman-
Rubin convergence criterion (Gelman & Rubin, 1992), checking that all
parameters in our chains have ${\cal R}-1<0.01$, with ${\cal R}$ quantifying
the Gelman-Rubin statistics, also known as the potential scale reduction
factor.

In order to better quantify the agreement (and/or disagreement) among models,
we have adopted the Akaike information criterion (AIC) (Akaike, 1974), defined
as

$\text{AIC}\equiv-2\ln\mathcal{L}_{\rm max}+2p,$ (5)

where $\mathcal{L}_{\rm max}$ is the maximum likelihood value and $p$ is the
total number of free parameters in the model. The AIC is derived from an
approximate minimization of the Kullback-Leibler divergence (Kullback &
Leibler, 1951) between the data-fit distribution and the true distribution.
For statistical comparison, the AIC difference between the model under study
and the reference model is calculated. It has been argued in Tan & Biswas
(2011) that a model is preferred over another one if their AIC difference is
larger than a threshold value $\Delta_{\rm thr}$. As a rule of thumb,
$\Delta_{\rm thr}=5$ can be considered the minimum value to assert a strong
support in favour of the model with a smaller AIC value, regardless of the
properties of the models under comparison (Liddle, 2007).

Table 1: Ranges of the model parameters. The table shows the parameters used
in the model, which are divided into fixed parameters, whose values remain
constant between epochs, and secular parameters, which can evolve over time.
The magnetic field range depends on the chosen atmosphere, such as H
($10^{10}-10^{13}$ G), O ($10^{12}-10^{13}$ G), and C ($10^{12}-10^{13}$ G).
Although extrapolations beyond this range may be allowed for numerical
purposes, any walker chains outside the acceptable interval are discarded to
ensure full coverage. Fixed parameters
---
$M(\text{M}_{\odot}$) | stellar mass | $0.1-2.3$
$R$(km) | stellar radius | $7.5-20.0$
$T$ (keV) | stellar surface temperature | $0.0-0.2$
$\log_{10}(B/G)$ | magnetic field strength | $10.0-13.0$
$i$ | angle between the LOS and the rotation axis | $0-90^{\circ}$
$\theta_{1}$ | spot one’s colatitude | $0-90^{\circ}$
$\theta_{2}$ | spot two’s colatitude | $0-90^{\circ}$
$\phi_{2}$ | spot two’s longitude | $0-360^{\circ}$
$\theta_{3}$ | spot three’s colatitude | $0-180^{\circ}$
$\phi_{3}$ | spot three’s longitude | $0-360^{\circ}$
Secular parameters
$\theta_{c1}$ | spot one’s semi-aperture | $0-15^{\circ}$
$T_{1}$ (keV) | spot one’s temperature | $0.0-1.0$
$\theta_{c2}$ | spot two’s semi-aperture | $0-15^{\circ}$
$T_{2}$ (keV) | spot two’s temperature | $0.0-1.0$
$\theta_{c3}$ | spot three’s semi-aperture | $0-15^{\circ}$
$T_{3}$ (keV) | spot three’s temperature | $0.0-1.0$
Figure 2: Parametric space at 1$\sigma$ (dark blue) and $2\sigma$ (light blue)
confidence levels (CLs) and one-dimensional marginalized distribution of $M$,
$R$ and $\log(B)$ for three hot spots on the surface of 3XMM J1852+0033. For
completeness, $M$, $R$ and $B$ are measured in units of $M_{\odot}$, km, and
$G$, respectively. The quoted values of the mass, radius and surface magnetic
field are at 1$\sigma$ CL. Our control simulations indicate that the radius of
3XMM J1852+0033 is reliable at $2\sigma$ CL
$\left(R=12.02^{+1.44}_{-1.42}\rm{km}\right)$, while its mass and surface
magnetic field are reliable at $1\sigma$ CL.  Figure 3: $M-R$ relations for
(purely hadronic) soft, intermediate and stiff cEFT EOS models and also for
the SLy4 EOS. For completeness, selected polytropic EOS models for hybrid
stars with different energy density jumps ($\eta$ values) are presented (for
details on them, see the main text). The cusps in the associated curves mark
the phase transition masses beyond which an exotic (e.g., quark) phase
appears. The dark blue and the light blue areas are the posterior constraints
for 1$\sigma$ (68$\%$) and 2$\sigma$ ($95\%$) CL, respectively. The black
cross marks the median of our best-fitting stellar parameters.

In order to minimise the effects from the light curve’s time resolution, as
well as from the binning choice when building the folded light curves (pulse
profiles), two precautions were taken prior to the fitting. First, our pulse
profile model was calculated taking into account the data time resolution,
i.e., the model was convoluted with a top-hat function with 2.6 s width.
Secondly, synthetic pulse profiles with different binning were produced to
test if after convergence – according to our Bayesian analysis – the injection
data could be recovered. It turns out that the synthetic pulse profiles
produced with the same time resolution of our folded light curve data allow
recovery of the injected mass and surface magnetic field of the star within
$1\sigma$ ($68\%$) CL. The injected radius of the star, though, could only be
retrieved at $2\sigma$ ($95\%$) CL. The number of hot spots has been correctly
recovered in all cases, while the positions and temperatures of the spots only
within $95\%$ CL. Thus, to be conservative given our data quality, we only
interpret radius related observables, geometrical properties and temperatures
of hot spots at $2\sigma$ CL.

##  4 Results and Discussion

We have investigated models with different numbers of hot spots and atmosphere
compositions. Details about them can be found on Table 2. The best model has
been found to have three spots and a carbon atmosphere. Its AIC was 362.36
($\chi^{2}/{\rm d.o.f.}=280.36/230=1.22$), while for two hot spots, it was
647.41 ($\chi^{2}/{\rm d.o.f.}=581.41/242=2.40$). A model with a single spot
not even presented a convergence of the chains, and it has been discarded. For
a model with four hot spots, we have not obtained a decrease in $\chi^{2}$
larger than $1\sigma$ with respect to the three spots’ case. Given this
penalty, it is disfavoured over a model with three spots. Note that in the
case of three hot spots, we are working on a 16D reference parameter space.
However, the total number of free parameters depend on the number of epochs
analyzed and also on the number of fixed and secular parameters. In our case
of 5 epochs, 6 secular and 10 fixed parameters, one has in total $10+6\times
5=40$ free parameters. Only for reference, the total number of data-points we
have for the 5 epochs is 270 (thus, ${\rm d.o.f.=230}$). The best fit to the
data (done for all the epochs simultaneously) and the associated hot spot
geometry on the star are shown in Figs. 4 and 5, respectively.

Figure 4: Phase-folded pulse profiles in counts per second for each epoch we
have investigated. Observational datapoints (black histogram) and the best
fits (red solid curves) are plotted. The energy band is $0.3-10$ keV. The
dashed blue curves correspond to the contributions to the phase-folded pulse
profiles due to spot 1 (identified by $\theta_{1}$), the dot-dashed green
curves due to spot 2 (identified by $\theta_{2}$), and the dotted magenta
curves due to spot 3 (identified by $\theta_{3}$).

Figure 5: Geometry of the best fit: 3D representation (left) and stereographic projection (right). On the left (right) are shown the spots as red circles (red dots), and the rotation axis and the line of sight with the blue arrow (blue dot) and the black arrow (black dot), respectively.  Table 2:  Comparison of AICs and $\chi^{2}$s for models with one, two, and three hot spots and different atmosphere compositions for the energy band 0.3–10keV. Atmosphere | number of spots | AIC | $\mathbf{\chi^{2}}$ | d.o.f.
---|---|---|---|---
Hydrogen | one | 6251.00 | 7852.00 | 254
two | 1312.11 | 1902.11 | 242
three | 2735.00 | 2223.00 | 230
Oxygen | one | 3630.27 | 3540.27 | 254
two | 1630.15 | 1540.15 | 242
three | 1029.94 | 939.94 | 230
Carbon | one | 2581.00 | 2766.00 | 254
two | 647.41 | 581.41 | 242
three | 362.36 | 280.36 | 230

###  4.1 Best-fitting stellar parameters

Table 3 summarizes the main results of our statistical analysis at 1$\sigma$
(68% CL). Throughout the text, we also give results for the 95% CL. Due to our
data quality and the results from our data injections, physical outcomes
involving the radius of 3XMM J1852+0033, are physically interpreted only at
the $95\%$ CL. Figure 2 summarizes the parameter space for some key NS
parameters, namely $M$, $R$ and $B$. The all-epoch fit for 3XMM J1852+0033
data suggests that $M=2.09^{+0.16}_{-0.09}$ $M_{\odot}$,
$R=12.02^{+0.72}_{-0.35}$ km and $\log_{10}(B/{\rm G})=11.89^{+0.19}_{-0.93}$
at $1\sigma$; for $95\%$ CL, we have that $M=2.09^{+0.21}_{-0.49}$
$M_{\odot}$, $R=12.02^{+1.44}_{-1.42}$ km and $\log_{10}(B/{\rm
G})=11.89^{+1.51}_{-1.89}$. Our results are in full agreement with the recent
NICER estimates for the equatorial radius of PSR J0740+6620 (see Miller et
al., 2021; Riley et al., 2021). In terms of compactness, ${\cal C}\equiv
MG/(Rc^{2})$, one has that $0.246\leq{\cal C}\leq 0.268$. For the median
parameters, it follows that ${\cal C}=0.256$.

One can also assess the detectability of the NS parameters found by means of
the universal relations connecting the compactness to the dimensionless tidal
deformability $\Lambda$, the later obtained by the gravitational-wave (GW)
detectors during the inspiral of binary NS systems (Abbott et al., 2018):

${\cal C}=\sum_{k=0}^{2}a_{k}(\ln\Lambda)^{k},$ (6)

with $a_{0}=0.360$, $a_{1}=-0.0355$, $a_{2}=0.000705$ (Yagi & Yunes, 2017). By
numerically solving Eq. 6, one finds the associated tidal deformations for the
possible combinations of masses and radii. For the median $M$ and $R$ values,
it follows that $\Lambda_{{\cal C}=0.256}=23$, which is at the threshold of
detectability by third generation GW detectors (Chatziioannou, 2022).

We compare the results against selected purely hadronic reference equations of
state: the chiral effective field theory (cEFT) models (see, eg., Hebeler et
al., 2013; Greif et al., 2020), which are in agreement with nuclear physics
experiments and astrophysics, and also the benchmark nucleonic SLy4 EOS
(Douchin & Haensel, 2001). For completeness, we also use some polytropic sharp
phase transition EOSs for hybrid stars as presented in Pereira et al. (2022)
with different energy density jump $\eta$ values222The polytropic parameters
for the EOS with $\eta=0.18$ are: $(n_{0}=0.475\,{\rm
fm}^{-3},n_{02}=0.413\,{\rm
fm}^{-3},\gamma_{1}=3.5,\gamma_{2}=6.0,n_{cc}=0.20\,{\rm fm}^{-3})$. The
parameters for $\eta=0.69$ are: $(n_{0}=0.50\,{\rm fm}^{-3},n_{02}=0.322\,{\rm
fm}^{-3},\gamma_{1}=3.75,\gamma_{2}=6.5,n_{cc}=0.20\,{\rm fm}^{-3})$. For the
precise meaning of each parameter, see Pereira et al. (2022). . As shown in
Fig. 3, the intermediate stiffness chiral EOS is privileged to the detriment
of (very) soft and stiff EOSs. Additionally, within $95\%$ CL, the SLy4 EOS is
a viable possibility. Rotation may change the picture in general, increasing
the equatorial radii and masses for any EOS (Breu & Rezzolla, 2016). But this
is not relevant for 3XMM J1852+0033 due to its large period (${\sim}11.6$ s);
the effect of rotation on spacetime is of the order of $10^{-8}$ (Haensel et
al., 2007). More precise analysis should be performed to check for specific
impact of our findings on current multi-messenger constraints for the EOS.

Table 3: Estimated fixed and secular parameters of 3XMM J1852+0033 obtained
from the modelling of XMM-Newton light curves.

Fixed parameters
---
$M(\text{M}_{\odot}$) | $R(\text{km})$ | log$(B/G)$ | $T(\text{keV})$ | $i$ | $\theta_{1}$ | $\theta_{2}$ | $\phi_{2}$ | $\theta_{3}$ | $\phi_{3}$
$2.10^{+0.16}_{-0.09}$ | $12.02^{+0.72}_{-0.35}$ | $11.89^{+0.19}_{-0.93}$ | $0.16^{+0.14}_{-0.06}$ | $41.77^{+0.81}_{-1.25}$ | $86.10^{+1.08}_{-0.91}$ | $176.0^{+1.7}_{-2.3}$ | $76.6^{+3.4}_{-3.0}$ | $140.5^{+2.1}_{-2.1}$ | $80.8^{+3.4}_{-4.2}$

Secular parameters
---
Epoch | ObsID | $\theta_{c1}$ | $T_{1}(\text{keV})$ | $\theta_{c2}$ | $T_{2}(\text{keV})$ | $\theta_{c3}$ | $T_{3}(\text{keV})$
A | 0550670201 | $7.75^{+0.56}_{-0.52}$ | $0.462^{+0.024}_{-0.031}$ | $10.17^{+0.48}_{-0.41}$ | $0.801^{+0.028}_{-0.031}$ | $9.47^{+0.22}_{-0.18}$ | $0.779^{+0.034}_{-0.024}$
B | 0550670301 | $5.37^{+0.45}_{-0.49}$ | $0.684^{+0.030}_{-0.035}$ | $12.26^{+0.46}_{-0.38}$ | $0.683^{+0.035}_{-0.011}$ | $7.71^{+0.23}_{-0.19}$ | $0.036^{+0.018}_{-0.030}$
C | 0550670401 | $7.99^{+0.45}_{-0.54}$ | $0.550^{+0.030}_{-0.024}$ | $10.68^{+0.44}_{-0.63}$ | $0.644^{+0.024}_{-0.027}$ | $9.42^{+0.18}_{-0.28}$ | $0.65^{+0.023}_{-0.035}$
D | 0550670501 | $9.31^{+0.56}_{-0.43}$ | $0.459^{+0.024}_{-0.028}$ | $10.24^{+0.46}_{-0.41}$ | $0.955^{+0.024}_{-0.029}$ | $9.41^{+0.19}_{-0.27}$ | $0.958^{+0.028}_{-0.022}$
E | 0550670601 | $7.78^{+0.68}_{-0.33}$ | $0.331^{+0.029}_{-0.027}$ | $7.83^{+0.36}_{-0.45}$ | $0.905^{+0.026}_{-0.030}$ | $7.51^{+0.16}_{-0.27}$ | $0.897^{+0.025}_{-0.031}$

22footnotetext: All angles are in degrees.

###  4.2 Magnetar with standard pulsar magnetic field

We have found that the best-fitting value for the magnetic field strength
($B\sim 10^{11}-10^{12}$ G) is consistent with the magnetic field inferred
from the star’s spindown ($B<4.1\times 10^{13}$ G) (see Rea et al., 2014).
However, in this latter case, only the dipolar field is taken into account.
The presence of three hot spots suggests that the magnetic field should
actually be multipolar (see e.g., de Lima et al., 2020, and references
therein). The subcritical value inferred, comparable to ordinary pulsars,
suggests that 3XMM J1852+0033 has experienced burying due to accretion.
Possible origins to the accreted material are supernova layers or the
interstellar medium. Indeed, a small accreted mass amount (when compared to
the stellar mass) of around $5\times 10^{-4}M_{\odot}$ would already be enough
to bury magnetic fields of the order of $10^{14}$ G and lead to surface fields
around $10^{11}-10^{12}$ G for around $10^{3}-10^{4}$ yrs (Viganò & Pons,
2012). Considering that supernova remnant Kes 79 can be a physically near
object and that its estimated age is between 4.4 and 6.7 kyr (see, e.g. Zhou
et al., 2016; He et al., 2022), our results suggest that we may be just seeing
the re-emergence of $B$.

###  4.3 High mass and the possibility of dense-matter phase transitions

A relevant consequence of our best-fitting results is the possibility of a
high mass for 3XMM J1852+0033. The values obtained are in the same range of
the most massive pulsars (Demorest et al., 2010; Antoniadis et al., 2013;
Fonseca et al., 2016; Cromartie et al., 2020), and below the estimated maximum
mass (for nonrotating stars) inferred from GW and kilonova analysis
(${\sim}2.16^{+0.17}_{-0.15}\,M_{\odot}$) (Rezzolla et al., 2018). Therefore,
if confirmed, 3XMM J1852+0033 would constitute one of the most massive neutron
stars to date. Interestingly, the best-fitting value for $M$ using three spots
is very close to the inferred mass for two hot spots (at $1\sigma$ CL,
$M=2.00^{+0.19}_{-0.15}\,M_{\odot}$), meaning that a large mass persists for
different models, strengthening the case of 3XMM J1852+0033 being a massive
magnetar.

The possibility of a high mass for 3XMM J1852+0033 also makes it natural to
speculate if it may have an exotic (e.g., quark) core. Statistical studies
suggest that the phase transition mass should be high, around $2M_{\odot}$
(Annala et al., 2020), and this would render 3XMM J1852+0033 another strong
candidate for probing exotic cores. For this, as a first approach, we have
investigated a large sample of polytropic sharp phase transition EOSs (for
further details, see Pereira et al. (2022)). With our posteriors, one can
restrict the exotic-hadron density jumps in stars. We find, by using dozens of
thousands of parametric EOSs fulfilling known astrophysical constraints that
at $1\sigma$ CL, $\eta\equiv\Delta\rho/\rho_{\rm
trans}\equiv\rho_{-}/\rho_{+}-1=0.129^{+0.211}_{-0.099}$ where
$\rho_{-}(\rho_{+})$ is the density at the top of the exotic phase (bottom of
the hadronic phase). For $95\%$ CL, it follows that
$\eta=0.129^{+0.861}_{-0.129}$. These results are in broad agreement with more
rigorous Bayesian inferences combining multi-messenger astronomy, e.g.,
$0.40^{+0.20}_{-0.15}$ (Xie & Li, 2021) and $0.28^{+0.39}_{-0.25}$ (Li et al.,
2021), while slightly favoring weak phase transitions (small energy density
jumps; see, e.g., Sieniawska et al. (2019) and references therein). At $95\%$
CL, strong phase transitions are however also possible. For an example, see
Fig. 3 ($\eta=0.69$). Our polytropic EOSs and posteriors also allow us to
estimate the phase transition mass $M_{\rm trans}$, beyond which an exotic
core appears. At $1\sigma$ CL, we get $M_{\rm
trans}=1.39^{+0.36}_{-0.31}M_{\odot}$. Within $95\%$ CL, $M_{\rm
trans}=1.39^{+0.75}_{-0.56}M_{\odot}$. This means that exotic matter cores may
be present in NS with masses around the Chandrasekhar mass (Alsing et al.,
2018) all the way to massive NSs, in agreement with other statistical studies
(Annala et al., 2020). For comparison, we also estimate $\eta$ and
$M_{\rm{trans}}$ with our set of polytropic EOSs using the inferred parameters
of PSR J0030+0451 ($M=1.34^{+0.15}_{-0.14}M_{\odot}$,
$R=12.71^{+1.14}_{-1.19}$ km) (Miller et al., 2019; Riley et al., 2019). At
$1\sigma$ CL, we have $\eta=0.243^{+0.377}_{-0.183}$ and $M_{\rm
trans}=1.50^{+0.29}_{-0.26}M_{\odot}$. For $95\%$ CL,
$\eta=0.243^{+0.987}_{-0.240}$ and $M_{\rm
trans}=1.50^{+0.63}_{-0.48}M_{\odot}$. Clearly, 3XMM J1852+0033 and PSR
J0030+0451 (and PSR J0740+6620) lead to consistent phase-transition parameter
estimates.

###  4.4 Carbon atmosphere

Another aspect from our best-fitting result concerns the atmosphere of 3XMM
J1852+0033. We find that the best fit points to the carbon composition, and
disfavours (ionized) hydrogen or helium atmospheres. Blackbody models (no
atmosphere and beaming) are not favored by our analysis either. While 3XMM
J1852+0033 is definitely not the CCO in the Kes 79 supernova remnant (Ho,
2011), its carbon atmosphere renders it a CCO-like object, as it shares this
feature with nearly half of the young CCOs, e.g., the youngest of them, the
Cas A remnant (Ho et al., 2021; Klochkov et al., 2013; Ho & Heinke, 2009). As
for CCOs, a carbon atmosphere associated with a buried magnetic field could be
explained by the depletion of hydrogen due to its diffusion into the denser
and hotter layers of a young NS, where protons are captured by heavier
elements. This is the diffusive nuclear burning of hydrogen acting in young
NSs with surface temperatures $T_{\rm eff}>1.5$ MK (see e.g. Chang & Bildsten,
2003, 2004). This process generates a hydrogen current flowing downwards into
the hydrogen-burning layer, and it leads eventually to the emergence of the
carbon atmosphere. A similar mechanism can work for helium atmospheres, but
for somewhat higher temperatures ($T_{\rm eff}>2$ MK, see Chang et al.
(2010)). The effective temperatures of the hot spots of 3XMM J1852+0033, shown
in Tab. 3, are higher than 5 MK ($0.43{\rm keV}$), and therefore satisfy the
condition for an efficient diffusive nuclear burning. Specifically, at such
$T_{\rm eff}$, a hydrogen/helium atmosphere is completely burned into carbon
in less than a few days (Chang & Bildsten, 2004; Chang et al., 2010). Any
hydrogen/helium atmosphere after an outburst of duration of at least 200 days
(the estimated span of the bright phase of 3XMM J1852+0033) could only result
from the accretion of these elements from its environment. For completeness,
carbon debris from the supernova remnant Kes 79 falling onto 3XMM J1852+0033
might be another possibility for its carbon atmosphere. However, in this case,
carbon should reemerge after accretion of hydrogen and helium from the
interstellar medium, which could possibly happen via the DNB of the freshly
accreted plasma.

###  4.5 Impact of the modelling on the parameter uncertainties

It may seem surprising that the uncertainties in mass and radius are only
around $10\%-15\%$, given that the XMM-Newton MOS2 data has a time resolution
that is only about 1/4 of the rotational period of the magnetar. However, we
have taken this into account in our analysis by convolving the flux
integration with a step function that has a temporal width equivalent to the
satellite resolution. Importantly, our precise results are due to the use of
data from five epochs and their simultaneous fit. Additionally, a considerable
number of our parameters are fixed, which significantly increases the
precision of our results. Finally, although the number of bins per phase is
relatively low (50), we obtained precise results for the temperature and hot
spot maps on the star because the pulse profiles have an approximate
sinusoidal reconstruction (meaning that the precise number of points is not
crucial) and because we are fitting the parameters of a specific ray-tracing
model. The price to pay for this pulse-profile reconstruction is the
difficulty of finding fine-structure details on the star, which we have indeed
ignored. Although what we obtain are basically “smoothed out” results, they
still tell us interesting aspects about the star. For example, models with one
or two hot spots are not as good as three ones, and this can already give us
interesting insights about the properties of the magnetic fields in magnetars.
In addition, three hot spots are needed due to the slight deviation of the
phase-folded pulse profiles from sinusoids. It is known that smaller
uncertainties for the stellar parameters arise in the case of sinusoidal
waveforms when hot spots are around the stellar equator (Lo et al., 2013),
which is the case for two hot spots in our best-fitting results. This could
also help explain the relatively small mass and radius uncertainties we found.

We also stress that the choice of model to compute the unabsorbed fluxes for
the normalization leads to systematic uncertainties in our parameter
estimations. As mentioned, these values generally agree within $10\%-15\%$ for
different models in the whole band (0.3–10 keV). As we discussed in our
previous work (de Lima et al., 2020), these systematics roughly translate into
uncertainties of up to $5\%$ for the NS mass and $7\%$ for the NS radius.
However, these uncertainties are still smaller than the reliable ones we
estimated for the mass and radius of 3XMM J1852+0033.

In highly magnetized neutron star atmospheres, such as those of magnetars, the
orientation of the magnetic field lines with respect to the star surface’s
normal can affect the observed flux and spectra (Pavlov et al., 1994; Ho &
Lai, 2001). In order to assess the relevance of that for small hot spots, one
should look at the local fluxes (Pavlov et al., 1994). Figures 2 and 3 of
Pavlov et al. (1994) help us estimate the expected relative changes for
different angles $\@vec{B}$ makes with the normal to the surface ($0$ degrees
and approximately $40$ degrees) when taking the fixed value $B=4.7\times
10^{12}$ G. One roughly has that relative differences for the spectral
intensity are up to $40\%-50\%$ for the peak transparency angles and magnetic
field poles, and up to around $20\%-30\%$ elsewhere. (The numbers change
slightly depending on the photon energies.) One would expect smaller
differences for smaller magnetic field strengths, such as the ones coming from
our best-fitting results. For larger hot spots, such as ours, the resultant
effect on the flux of the orientation of the magnetic field with respect to
the normal to the NS surface may also decrease because it should be integrated
over the whole spot, which smooths out the effect (Ho & Lai, 2001). However,
the above flux differences are not negligible in general and would also impact
the mass and radius inferences by around $5\%-10\%$. Although smaller than our
estimated errors for the macroscopic parameters of 3XMM J1852+0033, the
general conclusion here is that more precise analysis should not ignore the
angle the magnetic field makes with the NS’s surface normal.

For the modelling of magnetars’ hot spot emissions, beaming should not be
overlooked either (Hu et al., 2019). While it should arise self-consistently
from a given microscopic atmospheric model, its description is not simple in
general. A way to overcome that when one wants to test several atmosphere
models is to assume an empirical beaming function whose free parameters are
inferred simultaneously with other stellar parameters by Bayesian analyses for
each atmosphere (Salmi et al., 2023). Here, since we were interested in back-
of-the-envelope estimates, we have considered the $I\propto\cos^{n}\alpha$
phenomenological beaming model (DeDeo et al., 2001) for all possible
atmospheric compositions, where $\alpha$ is the angle the ray makes with the
normal to the surface (for further details, see Appendix A). Due to the
absence of accretion columns, which favor $n=2-3$ choices (DeDeo et al.,
2001), we took $n=1$. Since this effective beaming model does not take into
account specific atmosphere models, the carbon atmosphere that we have
obtained as the best-fitting result should be seen as indicative. However,
fits of free beaming parameters of NICER’s PSR J0030+0451 show significant
overlaps for different atmosphere compositions and little impact on the
originally inferred PSR J0030+0451 parameters, based on a fully ionized
hydrogen atmosphere model (Salmi et al., 2023). This suggests that in some
cases the atmospheric composition does not largely affect the NS parameters
and the beaming function. Indeed, we found that different choices of the
beaming function, such as the Hopf one (DeDeo et al., 2001) and the one
associated with $n=3$, did not result in acceptable fits for the folded pulse
profiles in any atmosphere case. Thus, based on all the above, uncertainties
for the mass and radius of 3XMM J1852+0033 might be small when associated with
variations of the (correct) beaming functions for different atmospheric
compositions. We leave further investigation on this important point for
future work.

###  4.6 Additional tests

In order to increase confidence in our main results, we conducted a few tests
using subsets of our data. The main goal is to assess the impact of data
treatment choices and model-dependent factors on the analysis. Specifically,
how the binning of the pulse profile and/or the energy band used to extract
the light curve affect the results. First, we re-performed the modelling for
the same bins/phase as before (50), but for data within a more restricted
energy band (from 3 to 8 keV). The lower limit chosen is so that the
computation of unabsorbed flux is less sensitive between models (as mentioned,
in this band the fluxes agree within $2\%$), and the upper limit is just to
avoid a high background portion of the data. We briefly talk about the case
with 16 bins/phase at the end of the section.

We proceeded with the Bayesian analysis exactly as in the case of the 0.3–10
keV energy band and obtained consistent results. Indeed, the best-fitting
model in the 3–8 keV energy band was also found to be the one with three hot
spots and a carbon atmosphere. Moreover, the inferred values for the mass,
radius, and magnetic field were all in agreement (within the uncertainties)
with those obtained from the broader energy band. Specifically, we obtained
$2.20^{+0.10}_{-0.05}M_{\odot}$, $11.80^{+0.56}_{-0.52}$ km, and ${\rm
log_{10}}(B/G)=12.05^{+0.27}_{-0.03}$ for the mass, radius, and surface
magnetic field, respectively (within 1$\sigma$ CL). The difference between the
results from the more restricted energy band and the broader one can be
entirely explained in terms of flux uncertainties due to the modelling of
unabsorbed flux (which contributes to approximately $5\%-10\%$ uncertainty in
the mass and radius of 3XMM J1852+0033). Within 1$\sigma$ CL, the
uncertainties for the macroscopic stellar parameters in the more restricted
energy band are smaller due to reduced flux fluctuations. However, we decided
to report the results from the broader energy band since they encompass a more
representative dataset, in terms of the source’s emission, and try to account
for systematic modelling uncertainties that we did not explore in detail.

For completeness, we did another Bayesian analysis also for the 3–8 keV energy
range with 16 bins/phase. We obtained that the main effect of decreasing the
binning was the increase in the uncertainties of the observables with respect
to the case with 50 bins/phase. Median values only changed slightly. For
instance, radius uncertainties at $1\sigma$ CL were around 1.3 km, while mass
uncertainties were around 0.4 $M_{\odot}$. Larger uncertainties are very much
expected due to the reduction of data.

###  4.7 Future improvements

Our work represents an initial step towards systematic studies of the light-
curve predictions of isolated NSs. To obtain more accurate NS parameter
extractions, further improvements to the modelling and data reduction are
needed. In future studies, we plan to (i) investigate atmospheric models that
consider a wider range of elements (beyond hydrogen, oxygen and carbon) and
magnetic field orientations (e.g., cases where the magnetic field is not
parallel to the NS surface’s normal), (ii) examine a more diverse range of hot
spot geometries, beyond circular ones, (iii) do not ignore the orientation of
magnetic field lines with the surface’s normal, (iv) implement the beaming
function self-consistently in our models, and (v) allow more parameters to
vary when fitting the observations.

##  5 Summary

X-ray pulse profile analysis of 3XMM J1852+0033 suggests that it is a low-B
($\sim 10^{11}-10^{12}$ G) magnetar with a large mass ($\sim
2.0-2.2\,M_{\odot}$) and a carbon atmosphere. Diffuse nuclear burning of
hydrogen/helium would be a possible mechanism for the latter. Posterior $M(R)$
distributions favor intermediate stiffness EOSs for purely hadronic stars. If
exotic (e.g., quark) phases are present, phase transitions in which energy
density jumps are small enough (weak phase transitions) are favored. Our best-
fitting models feature three small hot spots, which suggest a multipolar
magnetic field structure. The estimated magnetic field strength is subcritical
(i.e., $B<4.4\times 10^{13}$ G). The above is in agreement with magnetars if
their surfaces have been buried and their fields are just re-emerging.
Finally, our results can be improved if better time resolution data is
available, so future observations with 3XMM J1852+0033 as the main target are
greatly encouraged.

###### Acknowledgements.

R.C.R.L. acknowledges the support of Fundação de Amparo à Pesquisa e Inovação
do Estado de Santa Catarina (FAPESC) under grant No. 2021TR912. J.G.C. is
grateful for the support of FAPES (1020/2022, 1081/2022, 976/2022, 332/2023),
CNPq (311758/2021-5), FAPESP (2021/01089-1), and NAPI “Fenômenos Extremos do
Universo” of Fundação Araucária. R.C.N. thanks the CNPq for partial financial
support under the project No. 304306/2022-3 and the FAPERGS for partial
financial support under the project No. 23/2551-0000848-3. P.E.S. acknowledges
PCI/INPE/CNPq for financial support under grant No. 300320/2022-1. J.P.P.,
M.B., P.H. and J.L.Z. gratefully acknowledge the financial support of the
National Science Center Poland grants no. 2016/22/E/ST9/00037 and
2018/29/B/ST9/02013. C.V.R. thanks the support from CNPq (310930/2021-9).
J.C.N.A. thanks CNPq (308367/2019-7) for partial financial support.

## References

  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101
  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
  * Alsing et al. (2018) Alsing, J., Silva, H. O., & Berti, E. 2018, MNRAS, 478, 1377
  * Annala et al. (2020) Annala, E., Gorda, T., Kurkela, A., Nättilä, J., & Vuorinen, A. 2020, Nature Physics, 16, 907
  * Antoniadis et al. (2013) Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448
  * 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
  * Beloborodov (2002) Beloborodov, A. M. 2002, ApJ, 566, L85
  * Breu & Rezzolla (2016) Breu, C. & Rezzolla, L. 2016, MNRAS, 459, 646
  * Chang & Bildsten (2003) Chang, P. & Bildsten, L. 2003, ApJ, 585, 464
  * Chang & Bildsten (2004) Chang, P. & Bildsten, L. 2004, ApJ, 605, 830
  * Chang et al. (2010) Chang, P., Bildsten, L., & Arras, P. 2010, ApJ, 723, 719
  * Chatziioannou (2022) Chatziioannou, K. 2022, Phys. Rev. D, 105, 084021
  * Coelho et al. (2017) Coelho, J. G., Cáceres, D. L., de Lima, R. C. R., et al. 2017, A&A, 599, A87
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nature Astronomy, 4, 72
  * de Lima et al. (2020) de Lima, R. C. R., Coelho, J. G., Pereira, J. P., Rodrigues, C. V., & Rueda, J. A. 2020, ApJ, 889, 165
  * DeDeo et al. (2001) DeDeo, S., Psaltis, D., & Narayan, R. 2001, ApJ, 559, 346
  * Demorest et al. (2010) Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081
  * Douchin & Haensel (2001) Douchin, F. & Haensel, P. 2001, A&A, 380, 151
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJ, 915, L12
  * Fonseca et al. (2016) Fonseca, E., Pennucci, T. T., Ellis, J. A., et al. 2016, ApJ, 832, 167
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gelman & Rubin (1992) Gelman, A. & Rubin, D. B. 1992, Statistical Science, 7, 457
  * Greif et al. (2020) Greif, S. K., Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2020, ApJ, 901, 155
  * Haensel et al. (2007) Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007, Astrophysics and Space Science Library, 326
  * He et al. (2022) He, X., Cui, Y., Yeung, P. K. H., et al. 2022, ApJ, 928, 89
  * Hebeler et al. (2013) Hebeler, K., Lattimer, J., Pethick, C., & Schwenk, A. 2013, Astrophys. J., 773, 11
  * Ho et al. (2007) Ho, W. C., Chang, P., Kaplan, D. L., et al. 2007, Advances in Space Research, 40, 1432
  * Ho (2011) Ho, W. C. G. 2011, MNRAS, 414, 2567
  * Ho (2013) Ho, W. C. G. 2013, Proceedings of the International Astronomical Union, 9, 435–438
  * Ho (2014) Ho, W. C. G. 2014, in Magnetic Fields throughout Stellar Evolution, ed. P. Petit, M. Jardine, & H. C. Spruit, Vol. 302, 435–438
  * Ho & Heinke (2009) Ho, W. C. G. & Heinke, C. O. 2009, Nature, 462, 71
  * Ho & Lai (2001) Ho, W. C. G. & Lai, D. 2001, MNRAS, 327, 1081
  * Ho et al. (2008) Ho, W. C. G., Potekhin, A. Y., & Chabrier, G. 2008, The Astrophysical Journal Supplement Series, 178, 102
  * Ho et al. (2021) Ho, W. C. G., Zhao, Y., Heinke, C. O., et al. 2021, MNRAS, 506, 5015
  * Hu et al. (2019) Hu, C.-P., Ng, C. Y., & Ho, W. C. G. 2019, MNRAS, 485, 4274
  * Kaspi & Beloborodov (2017) Kaspi, V. M. & Beloborodov, A. M. 2017, ARA&A, 55, 261
  * Klochkov et al. (2013) Klochkov, D., Pühlhofer, G., Suleimanov, V., et al. 2013, A&A, 556, A41
  * Kullback & Leibler (1951) Kullback, S. & Leibler, R. A. 1951, The Annals of Mathematical Statistics, 22, 79
  * Li et al. (2021) Li, A., Miao, Z., Han, S., & Zhang, B. 2021, ApJ, 913, 27
  * Liddle (2007) Liddle, A. R. 2007, Monthly Notices of the Royal Astronomical Society: Letters, 377, L74–L78
  * Lo et al. (2013) Lo, K. H., Miller, M. C., Bhattacharyya, S., & Lamb, F. K. 2013, ApJ, 776, 19
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJ, 918, L28
  * Mori & Ho (2007) Mori, K. & Ho, W. C. G. 2007, MNRAS, 377, 905
  * Olausen & Kaspi (2014) Olausen, S. A. & Kaspi, V. M. 2014, ApJS, 212, 6
  * Pavlov et al. (1994) Pavlov, G. G., Shibanov, Y. A., Ventura, J., & Zavlin, V. E. 1994, A&A, 289, 837
  * Pereira et al. (2022) Pereira, J. P., Bejger, M., Leszek Zdunik, J., & Haensel, P. 2022, arXiv e-prints, arXiv:2201.01217
  * Rea et al. (2014) Rea, N., Viganò, D., Israel, G. L., Pons, J. A., & Torres, D. F. 2014, ApJ, 781, L17
  * Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, ApJ, 852, L25
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, 918, L27
  * Salmi et al. (2023) Salmi, T., Vinciguerra, S., Choudhury, D., et al. 2023, arXiv e-prints, arXiv:2308.09319
  * Seward et al. (2003) Seward, F. D., Slane, P. O., Smith, R. K., & Sun, M. 2003, ApJ, 584, 414
  * Sieniawska et al. (2019) Sieniawska, M., Turczański, W., Bejger, M., & Zdunik, J. L. 2019, A&A, 622, A174
  * Strüder et al. (2001) Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18
  * Sun et al. (2004) Sun, M., Seward, F. D., Smith, R. K., & Slane, P. O. 2004, ApJ, 605, 742
  * Tan & Biswas (2011) Tan, M. Y. J. & Biswas, R. 2011, MNRAS, 419, 3292–3303
  * Turner et al. (2001) Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27
  * Turolla & Nobili (2013) Turolla, R. & Nobili, L. 2013, ApJ, 768, 147
  * Turolla et al. (2015) Turolla, R., Zane, S., & Watts, A. L. 2015, Reports on Progress in Physics, 78, 116901
  * Viganò & Pons (2012) Viganò, D. & Pons, J. A. 2012, MNRAS, 425, 2487
  * Xie & Li (2021) Xie, W.-J. & Li, B.-A. 2021, Phys. Rev. C, 103, 035802
  * Yagi & Yunes (2017) Yagi, K. & Yunes, N. 2017, Phys. Rep., 681, 1
  * Younes et al. (2022) Younes, G., Lander, S. K., Baring, M. G., et al. 2022, ApJ, 924, L27
  * Zhou et al. (2014) Zhou, P., Chen, Y., Li, X.-D., et al. 2014, ApJ, 781, L16
  * Zhou et al. (2016) Zhou, P., Chen, Y., Safi-Harb, S., et al. 2016, ApJ, 831, 192

##  Appendix A Pulse Profile model

Following our previous approach (see de Lima et al. (2020) for details), we
consider an observer at $r\rightarrow\infty$ and a photon that arises from the
stellar surface at $dS=R^{2}\sin\theta d\theta d\phi$, where $R$ is the
stellar radius, making an angle $\alpha$ with the local normal to the surface
($0\leq\alpha\leq\pi/2$). The photon path is bent by an additional angle
$\beta$ owing to the spacetime curvature, and the effective emission angle as
seen by the observer is $\psi=\alpha+\beta$. The geometry is symmetric
relative to $\phi$. Beloborodov (2002) has shown that a simple approximate
formula can be used to relate the emission angle $\alpha$ to $\theta$:

$1-\cos\alpha=(1-\cos\theta)\left(1-\frac{R_{s}}{R}\right),$ (7)

where $R_{s}=2GM/c^{2}$ is the Schwarzschild radius, $G$ is the gravitational
constant and $M$ is the stellar mass. We note that Eq. (7) is a very good
approximation for $R>3R_{s}$ since it typically leads to very small errors
($\lesssim 1\%$). For the range of masses and corresponding radii of interest
here, errors would be up to a few percent.

We assume that the stellar emission follows a local atmospheric spectrum, and
that the observed flux comes mainly from hot spots. The flux from the stellar
photosphere is also considered, which is an improvement relative to de Lima et
al. (2020) that only consider the emission from hot spots.

The intensity $I_{\nu}(T)$ at a given point on the stellar surface ($\nu$ is
the photon frequency) depends on the temperature $T$ there, but it could also
be related to an “out of spot” atmospheric region. The flux is proportional to
the visible area of the emitting region ($S_{V}$) plus a relativistic
correction, and it is given by (Beloborodov 2002; Turolla & Nobili 2013)

$\displaystyle F_{\nu}$
$\displaystyle=\left(1-\frac{R_{s}}{R}\right)\int_{S_{V}}I_{\nu}(T)\cos\alpha\frac{d\cos\alpha}{d(\cos\theta)}ds\
,$ (8)

The model calculates the absolute count rates observed at an infinite
distance.

In polar coordinates, the circular hot spot has its center at $\theta_{0}$ and
a semi-aperture $\theta_{c}$. The spot is bound by the function
$\phi_{b}(\theta)$, where $0\leq\phi_{b}\leq\pi$, given by (Turolla & Nobili
2013)

$\phi_{b}=\left[\frac{\cos\theta_{c}-\cos\theta_{0}\cos\theta}{\sin\theta_{0}\sin\theta}\right]\
.$ (9)

Only the visible part of the star must be considered, therefore the spot must
be also limited by a constant $\theta_{F}$ defined by

$\theta_{F}=\arccos{\left(1-\frac{c^{2}R}{2GM}\right)^{-1}}\ .$ (10)

For a given bending angle $\beta$, that diverges from the emission angle
$\alpha$, $\theta_{F}$ occurs for the maximum emission $\alpha$, i.e.
$\alpha=\pi/2$ [see Fig. 1 of de Lima et al. (2020)]. In Newtonian gravity,
where $\beta=0$, the maximum visible angle is $\theta_{F}=\pi/2$, meaning that
half of the stellar surface is visible. However, for a relativistic star,
$\theta_{F}>\pi/2$.

By using the approximation given by Eq. 7 into Eq. 8, and integrating in
$\phi$, one obtains

$\displaystyle F_{\nu}$
$\displaystyle=2\left(1-\frac{R_{s}}{R}\right)^{2}\int_{\theta_{min}}^{\theta_{max}}I_{\nu}(T,\theta)A(\theta)\sin\theta\
\phi_{b}(\theta)d\theta\ ,$ (11)

where $I_{\nu}=I_{\nu}(T)\cos\theta$, which takes into account non-isotropic
emission (DeDeo et al. 2001) (e.g., due to the magnetic field lines), and

$\displaystyle A(\theta)$
$\displaystyle=\left[\frac{R_{S}}{R}+(1-\frac{R_{S}}{R})\cos\theta\right],$
(12) $\displaystyle\theta_{\text{min}}$
$\displaystyle=\text{min}[\text{max}(0,\theta_{0}-\theta_{c}),\theta_{F}]\ ,$
(13) $\displaystyle\theta_{\text{max}}$
$\displaystyle=\text{min}[\theta_{0}+\theta_{c},\theta_{F}]\ .$ (14)

The total flux produced by $N_{\sigma}$ spots, where the $\sigma$-th spot has
a semi-aperture $\theta_{c\sigma}$ and a temperature $T_{\sigma}$, can be
calculated by adding up each contribution, and so we have

$\displaystyle F_{\nu}^{TOT}$ $\displaystyle=$
$\displaystyle\sum_{\sigma}F_{\nu}(T_{\sigma})\ .$ (15)

The pulse profile in a given energy band $[\nu_{1},\nu_{2}]$ for a given spot
$\sigma$ is

$F_{\sigma}(\nu_{1},\nu_{2})=\int_{\nu_{1}}^{\nu_{2}}F_{\nu}(T_{\sigma})d\nu\
.$ (16)

Therefore, one can rewrite Eq. (15) for a given energy band, and it becomes

$\displaystyle F^{TOT}=\sum_{\sigma}F_{\sigma}(\nu_{1},\nu_{2})\ .$ (17)

We define the Line of Sight (LOS) as the unity vector $\hat{\mathbf{l}}$
pointing from the neutron stellar center to the observer in the infinity.
Also, we define by $\hat{\mathbf{r}}$ the unit vector parallel to the rotation
axis of the star, whose angular velocity is $\Omega=2\pi/P$.

We take $i$ as the angle between the LOS and the rotation axis ($\cos
i=\hat{\mathbf{r}}\cdot\hat{\mathbf{l}}$). As the star rotates, the polar
coordinate of the spot’s center, $\theta_{0}$, changes. Let $\gamma(t)=\Omega
t$ be the star’s rotational phase. Thus, from a geometrical reasoning we have
that

$\cos\theta_{0}(t)=\cos i\cos\theta-\sin i\sin\theta\cos\gamma(t)\ ,$ (18)

where we have taken that $i$ and $\theta$ do not change with time.

##  Appendix B Atmosphere model table

Table 4 gives the parameters for which the Nsmaxg model is tabulated. The
energy range for the model is 0.05–10 keV (with 117 tabulated points), and it
assumes hydrogen, carbon, or oxygen as the main atmospheric constituent. The
magnetic field strength, local gravity, and temperature are the key inputs for
computing the flux. Note that the atmosphere tables assume that the flux
emerges at a zero angle with respect to the normal to the star surface. The
number of tabulated points for $\log T_{eff}$ varies between 12 and 14 for a
given $B$, surface gravitational field and energy value. For further details,
see (Ho et al. 2008; Ho 2014).

Table 4: Parameters of the Nsmaxg model that are tabulated in XSPEC. Element |  $B$ ($10^{12}$ G) |  $\log g$ (cm/s2) | $\log T_{eff}$
---|---|---|---
H | 0.01 | 2.4 | 5.5 - 6.7
H | 0.04 | 2.4 | 5.5 - 6.7
H | 0.07 | 2.4 | 5.5 - 6.7
H | 0.1 | 2.4 | 5.5 - 6.7
H | 1.0 | 0.4 - 2.5 | 5.5 - 6.8
H | 2.0 | 2.4 | 5.5 - 6.8
H | 4.0 | 2.4 | 5.5 - 6.8
H | 7.0 | 2.4 | 5.5 - 6.8
H | 10.0 | 2.4 | 5.6 - 6.8
H | 30.0 | 2.4 | 5.7 - 6.8
H | 1.26 | 1.6 | 5.5 - 6.8
H | 7.0 | 1.6 | 5.5 - 6.8
O | 1.0 | 2.4 | 5.8 - 6.9
O | 10.0 | 2.4 | 5.8 - 6.9
C | 1.0 | 2.4 | 5.8 - 6.9
C | 10.0 | 2.4 | 5.8 - 6.9


454 | 013417.1671+304843.974 | 0.1 | 0.4 | 0.3 | 0.7 | – | 1.2${}^{+0.14}_{-0.13}$ | 0.21
460 | 013418.0482+304208.638 | 0.5 | -0.2 | 0.1 | 1.0 | – | 0.4${}^{+0.15}_{-0.13}$ | 0.26
463 | 013418.5457+304709.500 | 0.5 | 0.9 | -0.7 | 0.7 | – | 0.2${}^{+0.16}_{-0.19}$ | 0.31
471 | 013421.1706+303931.110 | -0.0 | -0.1 | 0.2 | 0.5 | 2.0${}^{+5.2}_{-2.0}$ | 1.6${}^{+0.14}_{-0.14}$ | 0.11
472 | 013421.2721+304335.749 | 0.4 | 0.7 | 0.0 | 0.7 | – | 0.8${}^{+0.0}_{-0.33}$ | 0.21
477 | 013422.6599+304623.249 | -0.0 | 0.7 | -0.0 | 1.1 | – | 1.8${}^{+0.52}_{-0.11}$ | 0.21
477 | 013422.6275+304624.150 | -0.0 | 0.7 | 0.6 | 1.1 | – | 0.2${}^{+0.12}_{-0.2}$ | 0.21
497 | 013426.4962+304446.801 | -0.1 | -0.1 | 0.1 | 1.3 | 0.2${}^{+2.0}_{-0.2}$ | 0.2${}^{+0.15}_{-0.13}$ | 0.31
497 | 013426.5357+304446.729 | -0.1 | -0.1 | 1.0 | 1.3 | 0.2${}^{+2.0}_{-0.2}$ | 2.8${}^{+0.14}_{-0.14}$ | 0.31
502 | 013427.0123+304313.236 | -0.1 | -0.1 | 0.3 | 0.4 | 0.4${}^{+1.3}_{-0.4}$ | 0.4${}^{+0.14}_{-0.19}$ | 0.06
508 | 013428.2930+304319.763 | 0.0 | 0.5 | -0.2 | 0.5 | – | 0.6${}^{+0.05}_{-0.58}$ | 0.11
523 | 013430.8829+304510.032 | 0.1 | -0.2 | -0.1 | 1.3 | – | 0.4${}^{+0.13}_{-0.16}$ | 0.21
541 | 013433.9145+304652.378 | -0.2 | -0.3 | 0.1 | 5.0 | – | 1.4${}^{+0.18}_{-0.12}$ | 0.01
542 | 013434.3598+304701.946 | -0.9 | -0.2 | 0.0 | 8.7 | – | 0.4${}^{+0.23}_{-0.1}$ | 0.01
545 | 013435.1355+304711.653 | -0.4 | 0.3 | 0.1 | 0.6 | – | 0.6${}^{+0.14}_{-0.19}$ | 0.26
545 | 013435.0362+304712.012 | -0.4 | 0.3 | 0.2 | 0.6 | – | 0.4${}^{+0.18}_{-0.16}$ | 0.26

Note. — We present the multi-wavelength properties used to evaluate the
quality of each HMXB candidate. For a full discussion of how we determined the
highest quality HMXB candidate sample, see Section 3.4. The Chandra hardness
ratios listed here were calculated using count rates from Tüllmann et al.
(2011). The F110W-F160W color is calculated using the HST photometry for each
HMXB companion star candidate, listed in Table 3. We include the mean SFR over
the last 80 Myr, calculated from the recent SFH maps of M33 generated with the
PHATTER optical photometry, presented in Lazzarini et al. (2022). Spectral
fits were performed on X-ray sources in the ChASeM33 catalog with sufficient
counts, and we present the resulting best-fit NH for these sources from
Tüllmann et al. (2011). The SED-derived AV value for each HMXB companion star
candidate comes from our SED fits, and is also listed in Table 1. The SFH-
derived AV value comes from Lazzarini et al. (2022).

Table 6: Flags used to determine best HMXB candidate sample

T11 ID | PHATTER Pos. | FLAG: Chandra | FLAG: FG | FLAG: Flat | FLAG: No | FLAG: NH | FLAG: Spectral | FLAG
---|---|---|---|---|---|---|---|---
|  | HRs | IR Color | SED shape | SF $<$ 80 Myr | AV mismatch | Type | SUM
192 | 013330.2309+304255.014 | 0 | 0 | 0 | 0 | 0 | 0 | 0
192 | 013330.1734+304254.931 | 0 | 0 | 0 | 0 | 0 | 0 | 0
200 | 013331.2921+303402.538 | 0 | 1 | 0 | 0 | 0 | 0 | 1
204 | 013331.7587+303555.916 | 0 | 0 | 0 | 0 | 0 | 0 | 0
205 | 013331.8256+304011.190 | 0 | 0 | 0 | 1 | 0 | 0 | 1
205 | 013331.8551+304012.079 | 0 | 0 | 0 | 1 | 0 | 0 | 1
207 | 013332.1735+303656.933 | 0 | 1 | 1 | 0 | 1 | 0 | 3
212 | 013332.5421+303618.159 | 0 | 0 | 0 | 0 | 0 | 0 | 0
213 | 013332.6992+303338.903 | 0 | 0 | 0 | 0 | 0 | 0 | 0
221 | 013333.7188+303109.771 | 0 | 0 | 0 | 0 | 0 | 0 | 0
225 | 013334.1333+303211.355 | 0 | 0 | 0 | 0 | 0 | 0 | 0
225 | 013334.1489+303211.229 | 0 | 0 | 0 | 0 | 0 | 1 | 1
225 | 013334.1388+303211.457 | 0 | 0 | 0 | 0 | 0 | 0 | 0
225 | 013334.1134+303211.388 | 0 | 1 | 0 | 0 | 0 | 0 | 1
225 | 013334.0783+303211.376 | 0 | 0 | 0 | 0 | 0 | 0 | 0
228 | 013334.5519+303555.630 | 0 | 0 | 0 | 0 | 0 | 0 | 0
233 | 013335.5034+303729.430 | 0 | 0 | 0 | 0 | 0 | 0 | 0
237 | 013336.0402+303333.097 | 0 | 0 | 0 | 0 | 0 | 1 | 1
250 | 013337.8783+303837.497 | 0 | 1 | 0 | 0 | 0 | 0 | 1
252 | 013337.9678+304035.531 | 0 | 0 | 0 | 0 | 0 | 0 | 0
261 | 013339.1501+303216.771 | 0 | 0 | 0 | 0 | 0 | 0 | 0
266 | 013339.8512+304349.810 | 0 | 0 | 0 | 0 | 0 | 0 | 0
268 | 013340.1282+304322.598 | 0 | 0 | 0 | 0 | 0 | 0 | 0
272 | 013340.8370+303524.269 | 0 | 0 | 0 | 0 | 0 | 0 | 0
272 | 013340.8230+303523.898 | 0 | 0 | 0 | 0 | 0 | 0 | 0
272 | 013340.7570+303523.918 | 0 | 0 | 0 | 0 | 0 | 0 | 0
272 | 013340.8203+303524.111 | 0 | 1 | 0 | 0 | 0 | 0 | 1
274 | 013341.2625+303213.735 | 0 | 0 | 0 | 0 | 0 | 0 | 0
274 | 013341.2839+303213.275 | 0 | 0 | 0 | 0 | 0 | 0 | 0
276 | 013341.4248+303815.823 | 0 | 0 | 0 | 0 | 0 | 0 | 0
276 | 013341.4588+303815.839 | 0 | 0 | 0 | 0 | 0 | 0 | 0
277 | 013341.5638+304136.445 | 0 | 1 | 0 | 0 | 0 | 0 | 1
280 | 013342.1069+304852.999 | 0 | 0 | 0 | 0 | 0 | 0 | 0
281 | 013342.4948+304253.566 | 0 | 0 | 0 | 0 | 0 | 0 | 0
294 | 013345.3022+304134.902 | 0 | 0 | 0 | 0 | 0 | 0 | 0
294 | 013345.1904+304134.824 | 0 | 0 | 0 | 0 | 0 | 0 | 0
299 | 013346.5634+303748.811 | 0 | 0 | 0 | 0 | 0 | 0 | 0
304 | 013347.5081+304042.825 | 0 | 0 | 0 | 0 | 0 | 0 | 0
315 | 013350.4962+304858.789 | 0 | 0 | 1 | 0 | 0 | 0 | 1
316 | 013350.5300+303822.076 | 0 | 0 | 0 | 0 | 0 | 0 | 0
316 | 013350.4516+303821.461 | 0 | 0 | 1 | 0 | 0 | 0 | 1
316 | 013350.5026+303821.603 | 0 | 0 | 1 | 0 | 0 | 0 | 1
319 | 013351.1466+303823.370 | 0 | 0 | 0 | 0 | 0 | 0 | 0
319 | 013351.1192+303823.231 | 0 | 0 | 0 | 0 | 0 | 0 | 0
321 | 013352.1445+303844.851 | 0 | 0 | 0 | 0 | 0 | 0 | 0
329 | 013354.4304+303414.753 | 0 | 0 | 0 | 0 | 0 | 0 | 0
340 | 013355.7346+303924.919 | 0 | 1 | 0 | 0 | 0 | 0 | 1
347 | 013356.7712+303729.789 | 0 | 0 | 0 | 0 | 0 | 1 | 1
347 | 013356.7851+303729.343 | 0 | 1 | 0 | 0 | 0 | 0 | 1
348 | 013356.8378+303706.541 | 0 | 0 | 0 | 0 | 0 | 1 | 1
348 | 013356.7901+303706.260 | 0 | 0 | 0 | 0 | 0 | 0 | 0
358 | 013358.2439+303437.894 | 0 | 0 | 0 | 0 | 0 | 0 | 0
358 | 013358.2340+303438.742 | 0 | 1 | 0 | 0 | 0 | 0 | 1
363 | 013358.7032+304228.773 | 1 | 0 | 0 | 0 | 0 | 0 | 1
370 | 013400.2685+303057.540 | 0 | 0 | 0 | 0 | 0 | 0 | 0
377 | 013400.7323+303945.581 | 0 | 0 | 0 | 0 | 0 | 0 | 0
377 | 013400.7047+303944.987 | 0 | 1 | 0 | 0 | 0 | 0 | 1
382 | 013401.0581+303709.917 | 0 | 0 | 0 | 0 | 0 | 0 | 0
382 | 013401.1408+303711.803 | 0 | 0 | 0 | 0 | 0 | 0 | 0
384 | 013401.1653+303242.555 | 0 | 1 | 1 | 0 | 0 | 0 | 2
384 | 013401.1405+303242.573 | 0 | 0 | 0 | 0 | 0 | 0 | 0
390 | 013402.0897+304702.750 | 0 | 0 | 0 | 0 | 0 | 0 | 0
390 | 013402.1260+304702.936 | 0 | 0 | 0 | 0 | 0 | 0 | 0
391 | 013402.3852+303136.836 | 0 | 0 | 0 | 0 | 0 | 0 | 0
396 | 013402.8161+304941.003 | 0 | 0 | 1 | 0 | 0 | 0 | 1
397 | 013402.8397+303439.553 | 0 | 0 | 0 | 0 | 0 | 0 | 0
398 | 013402.8833+304151.341 | 0 | 0 | 0 | 0 | 0 | 0 | 0
398 | 013402.8632+304151.772 | 0 | 0 | 0 | 0 | 0 | 0 | 0
398 | 013402.8386+304151.547 | 0 | 1 | 0 | 0 | 0 | 0 | 1
406 | 013404.4372+304521.215 | 0 | 0 | 0 | 0 | 0 | 0 | 0
414 | 013407.6334+303901.935 | 0 | 0 | 0 | 0 | 0 | 0 | 0
414 | 013407.5945+303902.037 | 0 | 0 | 0 | 0 | 0 | 0 | 0
414 | 013407.6468+303902.541 | 0 | 0 | 0 | 0 | 0 | 0 | 0
415 | 013407.8274+303554.131 | 0 | 1 | 1 | 0 | 0 | 0 | 2
416 | 013408.3441+303852.435 | 0 | 0 | 0 | 0 | 0 | 0 | 0
424 | 013410.5196+303946.307 | 0 | 1 | 1 | 0 | 0 | 0 | 2
424 | 013410.5483+303946.919 | 0 | 0 | 0 | 0 | 0 | 0 | 0
439 | 013414.4690+303303.881 | 0 | 0 | 0 | 0 | 0 | 0 | 0
442 | 013415.0522+304748.936 | 0 | 0 | 0 | 0 | 0 | 0 | 0
452 | 013417.1086+303426.560 | 0 | 0 | 1 | 0 | 0 | 0 | 1
452 | 013417.0590+303426.483 | 0 | 0 | 0 | 0 | 0 | 0 | 0
452 | 013417.1327+303426.957 | 0 | 0 | 0 | 0 | 0 | 0 | 0
454 | 013417.1671+304843.974 | 0 | 0 | 0 | 0 | 0 | 0 | 0
460 | 013418.0482+304208.638 | 0 | 0 | 0 | 0 | 0 | 0 | 0
463 | 013418.5457+304709.500 | 0 | 0 | 0 | 0 | 0 | 0 | 0
471 | 013421.1706+303931.110 | 0 | 0 | 0 | 0 | 0 | 0 | 0
472 | 013421.2721+304335.749 | 0 | 0 | 0 | 0 | 0 | 0 | 0
477 | 013422.6599+304623.249 | 0 | 0 | 0 | 0 | 0 | 0 | 0
477 | 013422.6275+304624.150 | 0 | 1 | 0 | 0 | 0 | 1 | 2
497 | 013426.4962+304446.801 | 0 | 0 | 0 | 0 | 0 | 0 | 0
497 | 013426.5357+304446.729 | 0 | 0 | 0 | 0 | 0 | 0 | 0
502 | 013427.0123+304313.236 | 0 | 0 | 0 | 0 | 0 | 0 | 0
508 | 013428.2930+304319.763 | 0 | 0 | 0 | 0 | 0 | 0 | 0
523 | 013430.8829+304510.032 | 0 | 0 | 0 | 0 | 0 | 0 | 0
541 | 013433.9145+304652.378 | 0 | 0 | 0 | 0 | 0 | 0 | 0
542 | 013434.3598+304701.946 | 1 | 0 | 0 | 0 | 0 | 0 | 1
545 | 013435.1355+304711.653 | 0 | 0 | 0 | 0 | 0 | 0 | 0
545 | 013435.0362+304712.012 | 0 | 0 | 0 | 0 | 0 | 0 | 0

Note. — For each identified HMXB companion star candidate, we list whether it
raised each of our flags used to evaluate the quality of HMXB candidates. The
values used to determine the flags are listed in Table 5 and a discussion of
criterion is in Section 3.4.

Table 7: M33 HMXB Production Rates per Unit SFR

| HMXB Prod. Rate | HMXB Prod. Rate
---|---|---
Sample | 0$-$50 Myr ago | 0$-$80 Myr ago
| [HMXBs/(M⊙ yr-1)] | [HMXBs/(M⊙ yr-1)]
Full sample | 120$-$153 | 181$-$240
$\leq$1 Flag sample | 117$-$150 | 175$-$233
0 Flag sample | 107$-$136 | 150$-$199

Note. — HMXB production rates calculated per unit SFR for our HMXB candidate
sample. The rows correspond to subsamples based on our evaluation system
described in more detail in Section 3.4. The samples used for the calculations
in each row match the samples in Figure 4.

Table 8: Mean SFR in M31 and M33 over different timescales

Galaxy | 50 Myr mean SFR [M⊙ yr-1] | 80 Myr mean SFR [M⊙ yr-1] | 100 Myr mean SFR [M⊙ yr-1]
---|---|---|---
M31 | 0.23$\pm$0.01 | 0.24$\pm$0.01 | 0.24$\pm$0.01
M33 | 0.29$\pm$0.04 | 0.31$\pm$0.04 | 0.32$\pm$0.02

Note. — Mean SFR over various timescales in M31 and M33 in the areas of each
galaxy covered by overlapping HST and Chandra observations. The SFR values
come from spatially resolved recent SFHs of each galaxy (Lewis et al., 2015;
Lazzarini et al., 2022).

## References

  * Antoniou et al. (2009) Antoniou, V., Hatzidimitriou, D., Zezas, A., & Reig, P. 2009, ApJ, 707, 1080, doi: 10.1088/0004-637X/707/2/1080
  * Antoniou & Zezas (2016) Antoniou, V., & Zezas, A. 2016, MNRAS, 459, 528, doi: 10.1093/mnras/stw167
  * Antoniou et al. (2010) Antoniou, V., Zezas, A., Hatzidimitriou, D., & Kalogera, V. 2010, ApJ, 716, L140, doi: 10.1088/2041-8205/716/2/L140
  * Antoniou et al. (2019) Antoniou, V., Zezas, A., Drake, J. J., et al. 2019, ApJ, 887, 20, doi: 10.3847/1538-4357/ab4a7a
  * Basu-Zych et al. (2016) Basu-Zych, A. R., Lehmer, B., Fragos, T., et al. 2016, ApJ, 818, 140, doi: 10.3847/0004-637X/818/2/140
  * Basu-Zych et al. (2013) Basu-Zych, A. R., Lehmer, B. D., Hornschemeier, A. E., et al. 2013, ApJ, 774, 152, doi: 10.1088/0004-637X/774/2/152
  * Beasley et al. (2015) Beasley, M. A., San Roman, I., Gallart, C., Sarajedini, A., & Aparicio, A. 2015, MNRAS, 451, 3400, doi: 10.1093/mnras/stv943
  * Bodaghee et al. (2012) Bodaghee, A., Tomsick, J. A., Rodriguez, J., & James, J. B. 2012, ApJ, 744, 108, doi: 10.1088/0004-637X/744/2/108
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Brorby et al. (2016) Brorby, M., Kaaret, P., Prestwich, A., & Mirabel, I. F. 2016, MNRAS, 457, 4081, doi: 10.1093/mnras/stw284
  * Brown et al. (2019) Brown, M. J. I., Duncan, K. J., Landt, H., et al. 2019, Monthly Notices of the Royal Astronomical Society, 489, 3351, doi: 10.1093/mnras/stz2324
  * Carrera et al. (2008) Carrera, R., Gallart, C., Aparicio, A., et al. 2008, AJ, 136, 1039, doi: 10.1088/0004-6256/136/3/1039
  * Cioni (2009) Cioni, M. R. L. 2009, A&A, 506, 1137, doi: 10.1051/0004-6361/200912138
  * Coe & Kirk (2015) Coe, M. J., & Kirk, J. 2015, MNRAS, 452, 969, doi: 10.1093/mnras/stv1283
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Dalgarno & McCray (1972) Dalgarno, A., & McCray, R. A. 1972, ARA&A, 10, 375, doi: 10.1146/annurev.aa.10.090172.002111
  * de Grijs & Bono (2014) de Grijs, R., & Bono, G. 2014, AJ, 148, 17, doi: 10.1088/0004-6256/148/1/17
  * Dolphin (2002) Dolphin, A. E. 2002, MNRAS, 332, 91, doi: 10.1046/j.1365-8711.2002.05271.x
  * Garofali et al. (2019) Garofali, K., Levesque, E. M., Massey, P., & Williams, B. F. 2019, ApJ, 880, 8, doi: 10.3847/1538-4357/ab286e
  * Garofali et al. (2018) Garofali, K., Williams, B. F., Hillis, T., et al. 2018, MNRAS, doi: 10.1093/mnras/sty1612
  * Garofali et al. (2017) Garofali, K., Williams, B. F., Plucinsky, P. P., et al. 2017, MNRAS, 472, 308, doi: 10.1093/mnras/stx1905
  * Gilfanov et al. (2004) Gilfanov, M., Grimm, H.-J., & Sunyaev, R. 2004, MNRAS, 351, 1365, doi: 10.1111/j.1365-2966.2004.07874.x
  * Gordon et al. (2016) Gordon, K. D., Fouesneau, M., Arab, H., et al. 2016, ApJ, 826, 104, doi: 10.3847/0004-637X/826/2/104
  * Greig & Mesinger (2018) Greig, B., & Mesinger, A. 2018, MNRAS, 477, 3217, doi: 10.1093/mnras/sty796
  * Grimm et al. (2005) Grimm, H. J., McDowell, J., Zezas, A., Kim, D. W., & Fabbiano, G. 2005, ApJS, 161, 271, doi: 10.1086/468185
  * Güver & Özel (2009) Güver, T., & Özel, F. 2009, MNRAS, 400, 2050, doi: 10.1111/j.1365-2966.2009.15598.x
  * Haberl et al. (2022) Haberl, F., Maitra, C., Vasilopoulos, G., et al. 2022, A&A, 662, A22, doi: 10.1051/0004-6361/202243301
  * Haberl & Pietsch (2001) Haberl, F., & Pietsch, W. 2001, A&A, 373, 438, doi: 10.1051/0004-6361:20010595
  * Haberl & Sturm (2016) Haberl, F., & Sturm, R. 2016, A&A, 586, A81, doi: 10.1051/0004-6361/201527326
  * Harris & Zaritsky (2004) Harris, J., & Zaritsky, D. 2004, AJ, 127, 1531, doi: 10.1086/381953
  * Harris & Zaritsky (2009a) —. 2009a, AJ, 138, 1243, doi: 10.1088/0004-6256/138/5/1243
  * Harris & Zaritsky (2009b) —. 2009b, AJ, 138, 1243, doi: 10.1088/0004-6256/138/5/1243
  * Johnson et al. (2022) Johnson, L. C., Wainer, T. M., Torresvillanueva, E. E., et al. 2022, ApJ, 938, 81, doi: 10.3847/1538-4357/ac8def
  * Justham & Schawinski (2012) Justham, S., & Schawinski, K. 2012, MNRAS, 423, 1641, doi: 10.1111/j.1365-2966.2012.20985.x
  * Koplitz et al. (2023) Koplitz, B., Johnson, J., Williams, B. F., et al. 2023, arXiv e-prints, arXiv:2303.07318, doi: 10.48550/arXiv.2303.07318
  * Krivonos et al. (2018) Krivonos, R., Sazonov, S., Tsygankov, S. S., & Poutanen, J. 2018, MNRAS, 480, 2357, doi: 10.1093/mnras/sty1995
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Lamers & Levesque (2017) Lamers, H. J. G. L. M., & Levesque, E. M. 2017, Understanding Stellar Evolution, doi: 10.1088/978-0-7503-1278-3
  * Lazzarini et al. (2018) Lazzarini, M., Hornschemeier, A. E., Williams, B. F., et al. 2018, ApJ, 862, 28, doi: 10.3847/1538-4357/aacb2a
  * Lazzarini et al. (2019) Lazzarini, M., Williams, B. F., Hornschemeier, A. E., et al. 2019, ApJ, 884, 2, doi: 10.3847/1538-4357/ab3f32
  * Lazzarini et al. (2021) Lazzarini, M., Williams, B. F., Durbin, M., et al. 2021, ApJ, 906, 120, doi: 10.3847/1538-4357/abccca
  * Lazzarini et al. (2022) Lazzarini, M., Williams, B. F., Durbin, M. J., et al. 2022, ApJ, 934, 76, doi: 10.3847/1538-4357/ac7568
  * Lebouteiller et al. (2006) Lebouteiller, V., Kunth, D., Lequeux, J., et al. 2006, A&A, 459, 161, doi: 10.1051/0004-6361:20053161
  * Lehmer et al. (2010) Lehmer, B. D., Alexander, D. M., Bauer, F. E., et al. 2010, ApJ, 724, 559, doi: 10.1088/0004-637X/724/1/559
  * Lehmer et al. (2019) Lehmer, B. D., Eufrasio, R. T., Tzanavaris, P., et al. 2019, arXiv e-prints, arXiv:1905.05197.  https://arxiv.org/abs/1905.05197
  * Lehmer et al. (2021) Lehmer, B. D., Eufrasio, R. T., Basu-Zych, A., et al. 2021, ApJ, 907, 17, doi: 10.3847/1538-4357/abcec1
  * Lewis et al. (2015) Lewis, A. R., Dolphin, A. E., Dalcanton, J. J., et al. 2015, ApJ, 805, 183, doi: 10.1088/0004-637X/805/2/183
  * Lin et al. (2017) Lin, Z., Hu, N., Kong, X., et al. 2017, ApJ, 842, 97, doi: 10.3847/1538-4357/aa6f14
  * Linden et al. (2010) Linden, T., Kalogera, V., Sepinsky, J. F., et al. 2010, ApJ, 725, 1984, doi: 10.1088/0004-637X/725/2/1984
  * Linden et al. (2009) Linden, T., Sepinsky, J. F., Kalogera, V., & Belczynski, K. 2009, ApJ, 699, 1573, doi: 10.1088/0004-637X/699/2/1573
  * Liu et al. (2006) Liu, Q. Z., van Paradijs, J., & van den Heuvel, E. P. J. 2006, A&A, 455, 1165, doi: 10.1051/0004-6361:20064987
  * Long et al. (1996) Long, K. S., Charles, P. A., Blair, W. P., & Gordon, S. M. 1996, ApJ, 466, 750, doi: 10.1086/177548
  * Long et al. (2002) Long, K. S., Charles, P. A., & Dubus, G. 2002, ApJ, 569, 204, doi: 10.1086/339172
  * Long et al. (1981) Long, K. S., Dodorico, S., Charles, P. A., & Dopita, M. A. 1981, ApJ, 246, L61, doi: 10.1086/183553
  * Long et al. (2010) Long, K. S., Blair, W. P., Winkler, P. F., et al. 2010, ApJS, 187, 495, doi: 10.1088/0067-0049/187/2/495
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, ApJ, 840, 39, doi: 10.3847/1538-4357/aa6af9
  * Magrini et al. (2010) Magrini, L., Stanghellini, L., Corbelli, E., Galli, D., & Villaver, E. 2010, A&A, 512, A63, doi: 10.1051/0004-6361/200913564
  * Magrini et al. (2009) Magrini, L., Stanghellini, L., & Villaver, E. 2009, ApJ, 696, 729, doi: 10.1088/0004-637X/696/1/729
  * Maravelias et al. (2014) Maravelias, G., Zezas, A., Antoniou, V., & Hatzidimitriou, D. 2014, MNRAS, 438, 2005, doi: 10.1093/mnras/stt2302
  * Marigo et al. (2008) Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883, doi: 10.1051/0004-6361:20078467
  * Marigo et al. (2017) —. 2017, ApJ, 835, 77, doi: 10.3847/1538-4357/835/1/77
  * Markert & Rallis (1983) Markert, T. H., & Rallis, A. D. 1983, ApJ, 275, 571, doi: 10.1086/161555
  * McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4, doi: 10.1088/0004-6256/144/1/4
  * Mesinger et al. (2014) Mesinger, A., Ewall-Wice, A., & Hewitt, J. 2014, MNRAS, 439, 3262, doi: 10.1093/mnras/stu125
  * Mineo et al. (2012) Mineo, S., Gilfanov, M., & Sunyaev, R. 2012, MNRAS, 419, 2095, doi: 10.1111/j.1365-2966.2011.19862.x
  * Misanovic et al. (2006) Misanovic, Z., Pietsch, W., Haberl, F., et al. 2006, A&A, 448, 1247, doi: 10.1051/0004-6361:20053745
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Oey (1999) Oey, M. S. 1999, in IAU Symposium, Vol. 193, Wolf-Rayet Phenomena in Massive Stars and Starburst Galaxies, ed. K. A. van der Hucht, G. Koenigsberger, & P. R. J. Eenens, 627
  * Orosz et al. (2011) Orosz, J. A., McClintock, J. E., Aufdenberg, J. P., et al. 2011, ApJ, 742, 84, doi: 10.1088/0004-637X/742/2/84
  * Orosz et al. (2007) Orosz, J. A., McClintock, J. E., Narayan, R., et al. 2007, Nature, 449, 872, doi: 10.1038/nature06218
  * Orosz et al. (2009) Orosz, J. A., Steeghs, D., McClintock, J. E., et al. 2009, ApJ, 697, 573, doi: 10.1088/0004-637X/697/1/573
  * Pietsch et al. (2006) Pietsch, W., Haberl, F., Sasaki, M., et al. 2006, ApJ, 646, 420, doi: 10.1086/504704
  * Pietsch et al. (2004) Pietsch, W., Misanovic, Z., Haberl, F., et al. 2004, A&A, 426, 11, doi: 10.1051/0004-6361:20047068
  * Pietsch et al. (2009) Pietsch, W., Haberl, F., Gaetz, T. J., et al. 2009, ApJ, 694, 449, doi: 10.1088/0004-637X/694/1/449
  * Plucinsky et al. (2008) Plucinsky, P. P., Williams, B., Long, K. S., et al. 2008, ApJS, 174, 366, doi: 10.1086/522942
  * Politakis et al. (2020) Politakis, B., Zezas, A., Andrews, J. J., & Williams, S. J. 2020, MNRAS, 493, 5369, doi: 10.1093/mnras/staa561
  * Predehl & Schmitt (1995) Predehl, P., & Schmitt, J. H. M. M. 1995, A&A, 500, 459
  * Quirk et al. (2022) Quirk, A. C. N., Guhathakurta, P., Gilbert, K. M., et al. 2022, AJ, 163, 166, doi: 10.3847/1538-3881/ac5324
  * Ramachandran et al. (2022) Ramachandran, V., Oskinova, L. M., Hamann, W. R., et al. 2022, arXiv e-prints, arXiv:2208.07773.  https://arxiv.org/abs/2208.07773
  * Ranalli et al. (2003) Ranalli, P., Comastri, A., & Setti, G. 2003, A&A, 399, 39, doi: 10.1051/0004-6361:20021600
  * Reig (2011) Reig, P. 2011, Ap&SS, 332, 1, doi: 10.1007/s10509-010-0575-8
  * Reig et al. (2016) Reig, P., Nersesian, A., Zezas, A., Gkouvelis, L., & Coe, M. J. 2016, A&A, 590, A122, doi: 10.1051/0004-6361/201628271
  * Rivinius et al. (2013) Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&A Rev., 21, 69, doi: 10.1007/s00159-013-0069-0
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schulman & Bregman (1995) Schulman, E., & Bregman, J. N. 1995, ApJ, 441, 568, doi: 10.1086/175383
  * Shao & Li (2014) Shao, Y., & Li, X.-D. 2014, ApJ, 796, 37, doi: 10.1088/0004-637X/796/1/37
  * Sidoli & Paizis (2018) Sidoli, L., & Paizis, A. 2018, MNRAS, 481, 2779, doi: 10.1093/mnras/sty2428
  * Tauris et al. (2017) Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170, doi: 10.3847/1538-4357/aa7e89
  * Telford et al. (2023) Telford, O. G., McQuinn, K. B. W., Chisholm, J., & Berg, D. A. 2023, ApJ, 943, 65, doi: 10.3847/1538-4357/aca896
  * Toribio San Cipriano et al. (2016) Toribio San Cipriano, L., García-Rojas, J., Esteban, C., Bresolin, F., & Peimbert, M. 2016, MNRAS, 458, 1866, doi: 10.1093/mnras/stw397
  * Trinchieri et al. (1988) Trinchieri, G., Fabbiano, G., & Peres, G. 1988, ApJ, 325, 531, doi: 10.1086/166025
  * Tüllmann et al. (2011) Tüllmann, R., Gaetz, T. J., Plucinsky, P. P., et al. 2011, ApJS, 193, 31, doi: 10.1088/0067-0049/193/2/31
  * Verley et al. (2009) Verley, S., Corbelli, E., Giovanardi, C., & Hunt, L. K. 2009, A&A, 493, 453, doi: 10.1051/0004-6361:200810566
  * Vulic et al. (2018) Vulic, N., Hornschemeier, A. E., Wik, D. R., et al. 2018, ApJ, 864, 150, doi: 10.3847/1538-4357/aad500
  * Wainer et al. (2022) Wainer, T. M., Johnson, L. C., Seth, A. C., et al. 2022, ApJ, 928, 15, doi: 10.3847/1538-4357/ac51cf
  * Walter et al. (2015) Walter, R., Lutovinov, A. A., Bozzo, E., & Tsygankov, S. S. 2015, A&A Rev., 23, 2, doi: 10.1007/s00159-015-0082-6
  * West et al. (2018) West, L. A., Lehmer, B. D., Wik, D., et al. 2018, ApJ, 869, 111, doi: 10.3847/1538-4357/aaec6b
  * Wik et al. (2014) Wik, D. R., Lehmer, B. D., Hornschemeier, A. E., et al. 2014, ApJ, 797, 79, doi: 10.1088/0004-637X/797/2/79
  * Williams et al. (2013) Williams, B. F., Binder, B. A., Dalcanton, J. J., Eracleous, M., & Dolphin, A. 2013, ApJ, 772, 12, doi: 10.1088/0004-637X/772/1/12
  * Williams et al. (2014) Williams, B. F., Lang, D., Dalcanton, J. J., et al. 2014, ApJS, 215, 9, doi: 10.1088/0067-0049/215/1/9
  * Williams et al. (2015) Williams, B. F., Wold, B., Haberl, F., et al. 2015, ApJS, 218, 9, doi: 10.1088/0067-0049/218/1/9
  * Williams et al. (2018) Williams, B. F., Lazzarini, M., Plucinsky, P. P., et al. 2018, ApJS, 239, 13, doi: 10.3847/1538-4365/aae37d
  * Williams et al. (2021) Williams, B. F., Durbin, M. J., Dalcanton, J. J., et al. 2021, arXiv e-prints, arXiv:2101.01293.  https://arxiv.org/abs/2101.01293
  * Yang et al. (2017a) Yang, J., Laycock, S., Christodoulou, D., et al. 2017a, The Astrophysical Journal, 839, 119
  * Yang et al. (2017b) Yang, J., Laycock, S., Drake, J., et al. 2017b, Astron. Nachr./AN, 338, 220, doi: 10.1002/asna. 201713334
  * Yang et al. (2022) Yang, J., Wik, D. R., Lehmer, B. D., et al. 2022, ApJ, 930, 64, doi: 10.3847/1538-4357/ac6351
  * Yukita et al. (2016) Yukita, M., Hornschemeier, A. E., Lehmer, B. D., et al. 2016, ApJ, 824, 107, doi: 10.3847/0004-637X/824/2/107
  * Zapartas et al. (2017) Zapartas, E., de Mink, S. E., Izzard, R. G., et al. 2017, A&A, 601, A29, doi: 10.1051/0004-6361/201629685