# Where are the magnetar binary companions? Candidates from a comparison with
binary population synthesis predictions

A. A. Chrimes,1 A. J. Levan1,2, A. S. Fruchter3, P. J. Groot1,4,5, P. G.
Jonker1,6, C. Kouveliotou7,8,J. D. Lyman2, E. R. Stanway2, N. R. Tanvir9 and
K. Wiersema10
1Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL,
The Netherlands
2Department of Physics, University of Warwick, Gibbet Hill Road, Coventry, UK
3Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD
21218, USA
4Inter-University Institute for Data Intensive Astronomy, Department of
Astronomy, University of Cape Town, Private Bag X3, Rondebosch 7701, South
Africa
5South African Astronomical Observatory, P.O. Box 9, 7935 Observatory, South
Africa
6SRON, Netherlands Institute for Space Research, Niels Bohrweg 4, 2333 CA,
Leiden, The Netherlands
7Department of Physics, The George Washington University, Corcoran Hall, 725
21st St NW, Washington, DC 20052, USA
8GWU/Astronomy, Physics and Statistics Institute of Sciences (APSIS)
9School of Physics and Astronomy, University of Leicester, University Road,
Leicester, LE1 7RH, UK
10Physics Department, Lancaster University, Lancaster, LA1 4YB, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

It is well established that magnetars are neutron stars with extreme magnetic
fields and young ages, but the evolutionary pathways to their creation are
still uncertain. Since most massive stars are in binaries, if magnetars are a
frequent result of core-collapse supernovae, some fraction are expected to
have a bound companion at the time of observation. In this paper, we utilise
literature constraints, including deep Hubble Space Telescope imaging, to
search for bound stellar companions to magnetars. The magnitude and colour
measurements are interpreted in the context of binary population synthesis
predictions. We find two candidates for stellar companions associated with
CXOU J171405.7–381031 and SGR 0755–2933, based on their $J$-$H$ colours and
$H$-band absolute magnitudes. Overall, the proportion of the Galactic magnetar
population with a plausibly stellar near-infrared counterpart candidate, based
on their magnitudes and colours, is between 5 and 10 per cent. This is
consistent with a population synthesis prediction of 5 per cent, for the
fraction of core-collapse neutron stars arising from primaries which remain
bound to their companion after the supernova. These results are therefore
consistent with magnetars being drawn in an unbiased way from the natal core-
collapse neutron star population, but some contribution from alternative
progenitor channels cannot be ruled out.

###### keywords:

stars: neutron – stars: magnetars – binaries: general

††pubyear: 2022††pagerange: Where are the magnetar binary companions?
Candidates from a comparison with binary population synthesis predictions–A

##  1 Introduction

Magnetars comprise a class of highly magnetised young neutron stars, with
magnetic field strengths up to 1015 G (Kouveliotou et al., 1998). They are
typically discovered when they produce ${\gamma}$-ray or X-ray bursts during
active episodes; their nature is also confirmed through the detection of their
persistent X-ray emission (Kaspi & Beloborodov, 2017). Some were originally
classified as soft gamma repeaters (SGRs), others as anomalous X-ray pulsars
(AXPs), but both populations are now believed to be manifestations of the same
phenomenon (Thompson & Duncan, 1996). Approximately 30 are currently known;
two are in the Magellanic clouds (one each, Olausen & Kaspi, 2014), while the
rest reside in our Galaxy. Their very rare Giant Flares are bright enough to
be seen at extragalactic distances (Hurley et al., 2005; Palmer et al., 2005);
two originated from magnetars in the Milky Way and a handful have been
identified in other galaxies (e.g. Hurley et al., 2010; Burns et al., 2021).

Measurements of magnetar spin periods $P$ and period derivatives $\dot{P}$
imply extremely young ages of 103-105 yr, under the assumption of a surface
dipole field and magnetic braking. These estimates are unreliable for
magnetars, since the fundamental assumption of a constant magnetic field (and
hence braking index) is shown to vary significantly (Younes et al., 2017).
However, these are upper limits for their decaying magnetic fields, leading to
overestimates of their true ages. The young magnetar ages are further backed
up by their associations with star forming regions and supernova remnants
(e.g. Muno et al., 2006; Gelfand & Gaensler, 2007; Nakano et al., 2015).
Moreover, current sophisticated methods now exist that account for their
magnetic field decay, in addition to particle winds and gravitational wave
emission (Mondal, 2021).

Despite some associations with H II regions and supernova remnants, around
half of the magnetar population have not yet been associated with a high-mass
progenitor (Olausen & Kaspi, 2014). This may be attributed to detection biases
due to e.g., dust, source crowding, and high column densities. There are,
however, magnetars for which there is no obvious association with recent star
formation despite a relatively clear sightline (see e.g., SGR 0501$+$4516
which lies in the Galactic anti-centre direction). Putative non-core-collapse
production channels include binary neutron star mergers and the accretion or
merger induced collapse of ONe white dwarfs resulting in magnetars (Duncan &
Thompson, 1992; Levan et al., 2006; Ablimit, 2021).

The rarity of magnetars (Kouveliotou et al., 1994; Beniamini et al., 2019)
raises questions about their origin. Are specific conditions required to
generate neutron stars with such strong magnetic fields? If their magnetic
fields decay fast enough, it may be that a large fraction of core-collapse
supernovae result in magnetars (Nakano et al., 2015; Beniamini et al., 2019),
which then may decay into X-ray dim isolated neutron stars (XDINS, Jawor &
Tauris, 2021) before fading beyond detection limits. Correspondingly, if
magnetars are intrinsically rare, particular evolutionary pathways may be
required. Suggestions include core-collapse following main sequence mergers
(Schneider et al., 2020), rapid core rotation pre-collapse (White et al.,
2021), and non-core collapse channels, as mentioned above.

Owing to their large energy reservoirs, magnetars have also been invoked as
the central engines in a variety of transients, from compact object mergers,
which produce short gamma-ray bursts (GRBs, Bucciantini et al., 2012; Gompertz
et al., 2013; Rowlinson et al., 2013; Gompertz et al., 2014; Metzger & Piro,
2014; Zou et al., 2021), core-collapse long GRBs (Wheeler et al., 2000;
Komissarov & Barkov, 2007; Metzger et al., 2011), superluminous supernovae
(Kasen & Bildsten, 2010; Woosley, 2010; Metzger et al., 2015) and fast blue
optical transients (Prentice et al., 2018; Mohan et al., 2020). Many fast
radio burst models also implicate magnetars (e.g. Lyubarsky, 2014; Metzger et
al., 2019), a connection strengthened by the discovery of low-luminosity FRBs
from the Galactic magnetar SGR 1935$+$2154 (CHIME/FRB Collaboration et al.,
2020; Bochenek et al., 2020; Younes et al., 2021). The host environments of
extragalactic FRBs are broadly representative of where we expect to find young
neutron stars (Heintz et al., 2020; Safarzadeh et al., 2020; Bochenek et al.,
2021; Mannings et al., 2021; Chrimes et al., 2021; Bhandari et al., 2021;
Katz, 2021). However, there is also evidence that FRBs do not follow the
cosmic star formation rate density (James et al., 2020; Zhang & Zhang, 2021),
and there are other objects which seem inconsistent with core-collapse
magnetars. For example, assuming that magnetars do produce FRBs, Tendulkar et
al. (2021) find that FRB 20180916B is too offset from a nearby star forming
region to be a young ejected magnetar. Subsequently, FRB 20200120E was located
in a globular cluster outside M81 (Kirsten et al., 2021). Since globular
clusters have old stellar populations, a young core-collapse magnetar
explanation for this source is challenging (Lu et al., 2021; Kremer et al.,
2021).

It remains unclear which progenitor channels contribute to magnetar
production, or if any particular channel dominates their formation rate. If we
assume that magnetars are drawn from the natal core-collapse neutron star
population in an unbiased way, then some fraction should have a bound
companion star. The majority of massive stars are in binaries (e.g. Sana et
al., 2014), and not all supernovae in binaries unbind the system, as evidenced
by the existence of binary neutron star and black hole mergers, X-ray binaries
and pulsars in binaries.

In this paper, we use published photometry of counterpart candidates in
conjunction with distance and extinction estimates to determine whether any of
the Galactic magnetar near-infrared (NIR) counterparts could be a bound
companion star. We determine the fraction of the population that could
plausibly have a bound companion, and compare these results to population
synthesis predictions for the fraction of neutron stars born with a bound
secondary companion.

The paper begins with population synthesis predictions in Section 2, and a
exploration of the major sources of uncertainty in these predictions. We
describe the sample selection in Section 3, before applying distance modulus
and extinction corrections in Section 4, allowing us to compare population
synthesis predictions with observations in Section 5. We summarise the results
in Section 6, discuss our findings in Section 7 in the context of magnetar
progenitor and evolutionary models, and conclude in Section 8. All magnitudes
are quoted in the Vega system, where conversions from AB magnitudes (Oke &
Gunn, 1983) have been made using
stsynphot111https://github.com/spacetelescope/stsynphot_refactor for Hubble
Space Telescope (HST ) data. AB magnitudes are obtained by adding 0.8227,
0.9204, 1.0973, and 1.2741 to Vega magnitudes for F110W, F125W, F140W and
F160W respectively. For $J$, $H$ and $K$-band filters, the conversions of
Blanton & Roweis (2007) are used.

Figure 1: The predicted orbital period, mass, $H$-band absolute magnitude and
$J$–$H$ colour distributions for bound companions to natal neutron stars. The
median absolute $H$-band magnitude is -1.2, and the median colour is -0.08.
The axes of the top two panels are on a log scale to clearly show the
distribution shape across the full range of parameter space.

##  2 Binary Population Synthesis

###  2.1 Predictions

To establish predictions for the fraction of neutron stars born in primary
supernovae, and the fraction of those with a surviving bound companion, we
make use of the population synthesis code BPASS (Binary Population and
Spectral Synthesis, Eldridge et al., 2017; Stanway & Eldridge, 2018, v2.2.1).
BPASS consists of thousands of binary and single stellar evolution models
across a grid of masses, mass ratios, orbital periods, and metallicities. Each
model is assigned a weighting, corresponding to how often it is predicted to
occur in a 106 M⊙ stellar population, with the initial weightings based on
observations of Galactic populations (Moe & Di Stefano, 2017).

For a fiducial estimate, we use the default BPASS model weightings and broken
power law initial mass function (IMF), with a maximum zero age main sequence
mass of 300 M⊙, and gradients of –1.30 (0.1-0.5 M⊙) and –2.35 (0.5-300 M⊙,
based on Kroupa et al., 1993). Solar metallicity (a metal mass fraction of
0.02) is assumed. To identify models which produce neutron stars in core-
collapse, we require in the final time step of each model (i) a total mass
$>2.0$ M⊙, (ii) a CO core mass $>1.38$ M⊙, (iii) a non-zero ONe core mass, and
(iv) a remnant mass in the range $1.38<M_{\rm rem}<3.0$ (see also Eldridge et
al., 2017, 2019; Chrimes et al., 2020; Briel et al., 2021). We then sum the
weightings for all models which satisfy these criteria, and examine the
contributions from different model types (single, primary, bound secondary,
unbound secondary). To identify how many primary neutron stars are born with a
bound companion, the weightings for secondary star models whose companion is a
neutron star at $t=0$ are summed. We find that around half of all neutron
stars are born in the supernovae of primary stars in binaries. Another
${\sim}$25 per cent arise from single stars or mergers, and the remaining
${\sim}$25 per cent from secondary stars.

Around ${\sim}$10 per cent of neutron stars born from primaries remain bound
after the supernova. We therefore predict that the fraction of natal neutron
stars with a bound, pre-supernova companion is f${}_{\rm bound}=0.05$. Figure
1 shows the predicted Log10(P/days), Log10(M/M⊙), $H$-band absolute magnitude
and $J$–$H$ colour distributions for bound companions to neutron stars, in the
first time step of secondary BPASS models. The bulk of the colours correspond
to massive main sequence companions, with the smaller, redder bump arising
from a minority of models where the secondary has already evolved off the main
sequence when the primary goes supernova. This is in qualitative agreement
with previous predictions for surviving supernova companions (e.g. Kochanek,
2009; Zapartas et al., 2017).

We focus on NIR colours because of the available observational constraints on
Galactic magnetar counterparts. Since the population lies predominantly in the
Galactic plane, and out to large distances ($>10$ kpc), high extinctions
typically render the optical magnitudes of magnetar counterparts extremely
faint (with only a few detections reported). Conversely, there $\sim$15
(depending on the confidence of the association) counterpart candidates in the
$H$-band (see Section 3). Since the companions to natal neutron stars are
primarily expected to be massive, un-evolved main sequence stars, other works
in this area - for example searches for surviving companions to supernovae -
typically use bluer optical bands (Maund et al., 2016). While such companions
will be intrinsically fainter in the NIR, typically by $\sim$2 magnitudes
compared to the optical, just $A_{\rm V}\sim 2$ of extinction negates this
advantage. Approximately two magnitudes of optical extinction is the lowest
value inferred for any Galactic magnetar (see Section 4), and as extinction
increases, the NIR becomes relatively more favourable. Combined with many more
detections of magnetar counterparts in the NIR compared with the optical, this
justifies our focus on the NIR as the best wavelength range in which to search
for companions.

Figure 2: How the prediction for fbound varies with 60 different realisations
of the input distributions, drawn through bootstrap re-sampling (as in Stanway
et al., 2020). Default assumptions yield f${}_{\rm bound}=0.05$. The upper
left distribution arises from varying the orbital periods (Log10(P/days),
labelled LogP), mass ratios ($q$) and binary fraction (BF) simultaneously.
Orbital periods are the largest individual contributor to the uncertainty.

###  2.2 Uncertainties

Previous estimates for fbound have varied, in part due to uncertainties on the
population synthesis input distributions (for details of how BPASS handles
binary dynamics following a supernova, see Tauris & Takens, 1998; Eldridge et
al., 2011). For example, Eldridge et al. (2011) predict that the number of
neutron stars born from primaries which remain bound is a few per cent. The
choice of binary fraction, initial period and mass ratio distributions will
change fbound, as will the assumed natal kick distribution.

Stanway et al. (2020) investigated the impact of uncertainties in the BPASS
input parameters on the integrated light of stellar populations. Using the
same set of iterations across orbital period, mass ratio and binary fraction,
and drawing from the observed distribution of these inputs in the Galactic
stellar populations, we determine the impact on fbound. The default period
distribution (LogP), binary fraction (BF) and mass ratio (q) distributions are
based on table 13 of Moe & Di Stefano (2017) as follows,

  1. 1.

LogP: A period distribution is applied in each of the five mass bins studied
by Moe & Di Stefano (2017). For example, 0.095$\pm$0.018 of Solar mass stars
have a companion with LogP$=$4.5–5.5 (P in days), whereas 0.30$\pm$0.09 of
stars above 16 M⊙ have a companion in this range,

  2. 2.

BF: defined as one minus the single star fraction in each mass bin. The BF in
each mass bin has an associated observational uncertainty, e.g. the single
star fraction for Solar-type stars is 0.60$\pm$0.04, and among the most
massive stars it is 0.06$\pm$0.06,

  3. 3.

q: the mass ratio distribution is a broken power law following $p_{q}\propto
q^{\gamma}$, which varies in each LogP and (primary star) mass bin. For
example, the probability $p_{q}$ of a mass ratio in the range 0.3–1.0 for
Solar mass primaries in a binary with LogP$=$5 is given by $\gamma=-0.5\pm
0.3$. In the mass ratio range 0.1–0.3, $\gamma=0.3\pm 0.4$.

Finally, the elevated fraction of twin systems is taken into consideration,
and the distributions are interpolated across the BPASS model grid. Randomly
sampling these input distributions within the quoted uncertainties (assuming
Gaussian uncertainties, as fully described in Stanway et al., 2020) yields the
results shown in Figure 2. Uncertainty on the initial orbital periods of
massive binaries is the largest source of binary parameter uncertainty in
fbound. The result when all uncertainties are varied together is
0.049$\pm$0.005, the wider spread showing that there are correlated
uncertainties between the parameters.

We also investigate the impact of assuming different neutron star natal kick
velocities - the kick distributions of Hobbs et al. (2005) and Igoshev (2020)
yield very similar fbound values of 0.049 and 0.050. Igoshev et al. (2021)
also found little difference between these kicks in terms of the Be star mass
distribution in X-ray binaries. The kick of Bray & Eldridge (2016, 2018),
however, gives a substantially larger value of 0.23. For an observational
comparison, we refer to Kochanek (2021) who performed a Gaia search for
unbound companions to 10 supernova remnants. Combining these results with a
previous search for surviving bound companions in supernova remnants, they
find that 0.036–0.197 (90 per cent confidence, excluding triple systems) of
core-collapse neutron stars are born in primaries which remain bound post
supernova. This figure also agrees well with the results of the StarTrack
binary population synthesis code (yielding 0.05, assuming an 84 per cent
binary fraction; Kochanek et al., 2019) and binary$\\_$c
($0.11^{+0.07}_{-0.09}$, assuming 22 per cent of massive binaries merge before
the first supernova, Renzo et al., 2019). Therefore, the prediction of
f${}_{\rm bound}=0.05$ is consistent with other population synthesis codes and
observational constraints.

##  3 Sample Selection

Our sample consists of the 23 Galactic magnetars (excluding the Magellanic
clouds) which have imaging detections or limits, and a distance estimate. Of
the 31 objects listed in the McGill magnetar catalogue (Olausen & Kaspi,
2014)222http://www.physics.mcgill.ca/~pulsar/magnetar/main.html, CXOU
J010043.1–721134 and SGR 0526–66 are in the Large and Small Magellanic clouds
respectively, while SGRs 1801–23 and 1808–20 have neither imaging nor a
distance estimate. Swift J1818.0-1607 and PSR J1846–0258 have no photometric
data while SGR 1833–0832 and AX J1818.8–1559 do not have a distance estimate.
The full list of magnetars in the sample, and the photometry used, is provided
in Table 1. Most sources have HST F160W ($H$) and F125W ($J$) observations, as
reported in Chrimes et al. (2022). The NIR counterparts are localised, where
possible, through astrometric alignment of a Chandra X-ray localisation with
the HST images, by finding common sources in the field or doing so through an
intermediate image. The RMS residuals of this translation, combined with the
source cetroiding uncertainty, determines the size of the error circle in the
HST frame. In other cases, a radio localisation is used (for which the
absolute astrometric uncertainty of HST dominates), or simply the absolute
astrometric uncertainty of the X-ray observations, if no common sources could
be found and no alternative positional measurements were available. Shortened
magnetar names, as indicated in this table, are used hereafter. Given the
population synthesis predictions, we would expect 1$\pm$1 (Poisson uncertainty
only) of the 23 objects in the sample to have a bound companion. This
prediction holds under the assumptions that (i) the young age estimates are
broadly correct (since the probability of a companion going supernova within
the age of the magnetar is extremely low) and (ii) they are drawn from the
natal core-collapse neutron star population in an unbiased away.

Table 1: The full magnetar sample considered in this paper (23 objects).
Photometry (Vega magnitudes) for the most likely counterpart(s) is listed in
each case. Shortened magnetar names, as in the first column without brackets,
are used throughout this paper. The criteria for selection are a distance
estimate and imaging. 1RXS J1708 has three candidate counterparts in (or just
outside) the X-ray error circle (see Chrimes et al., 2022), the candidates A
and B are labelled following Israel et al. (2003) and Durant & van Kerkwijk
(2006a). The magnitudes listed here are the primary ones used (i.e. for the
colour-magnitude diagram), other values from the literature are also plotted
later in Figure 4, following references in the McGill catalogue (Olausen &
Kaspi, 2014)2.

Magnetar |  $H$-band (Vega) |  $J$-band (Vega) | Instrument/survey | Filter(s) | Reference
---|---|---|---|---|---
4U 0142($+$61) | 20.80$\pm$0.01 | 21.72$\pm$0.01 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
SGR 0418($+$5729) |  $>$25.12 |  $>$26.08 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
SGR 0501($+$4516) | 22.56${}^{+0.06}_{-0.07}$ | 23.33$\pm$0.07 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
1E 1048(.1–5937) | 22.15${\pm}$0.20 |  $>$24.08 |  HST /NICMOS | F160W$+$F110W | Tam et al. (2008)
1E 1547(.0–5408) |  $>$20.22 |  $>$22.45 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
PSR J1622(–4950) | 22.39$\pm 0.05$ | 23.95$\pm 0.08$ |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
SGR 1627(–41) | 20.48$\pm$0.01 | 21.97$\pm 0.01$ |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
CXOU J1647(10.2–455216)$\dagger$ | 21.0$\pm 0.1$ | 23.5$\pm 0.2$ | VLT/NACO |  $H$,$J$ | Testa et al. (2018)
1RXS J1708(49.0–400910) A | 18.77${\pm}$0.01 | 20.61${\pm}$0.01 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
1RXS J1708(49.0–400910) B | 20.58${\pm}$0.01 | 22.37$\pm$0.02 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
1RXS J1708(49.0–400910) C | 21.80${\pm}$0.08 | 23.71${}^{+0.14}_{-0.16}$ |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
CXOU J1714(05.7–381031) | 17.45$\pm$0.01 | 19.01$\pm$0.01 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
SGR 1745(–2900) |  $>$15.97 |  $>$18.98 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
SGR 1806(–20) | 21.75$\pm$0.75 |  $>$21.1 | Keck/NIRC2, VLT/NAOS |  $K_{p}$,$J$ |  Tendulkar et al. (2012), Israel et al. (2005)
XTE J1810(–197)$\dagger$ | 21.67$\pm$0.12 | 22.92$\pm$0.22 | Gemini/NIRI |  $H$,$J$ | Testa et al. (2008)
Swift J1822(.3–1606) | 19.76$\pm$0.01 | 21.04$\pm 0.02$ |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
Swift J1834(.9–0846) | 21.74$\pm 0.02$ | 23.21$\pm 0.03$ |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
1E 1841(–045) | 20.80${\pm}$0.4 |  $>$22.10 | Magellan /Panic |  $H$,$J$ | Durant (2005)
3XMM J1852(46.6$+$003317)$\star$ | 18.85$\pm$0.01 | 21.97$\pm$0.02 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
SGR 1900($+$14)$\dagger$ | 21.17$\pm$0.50 | - | Keck/NIRC2 | $H$ | Tendulkar et al. (2012)
SGR 1935($+$2154) |  $\sim$24 | - |  HST /WFC3 | F140W |  Levan et al. (2018), Lyman et al. (2022)
1E 2259($+$586) | 22.52$\pm 0.04$ | 23.61${\pm}$0.05 |  HST /WFC3 | F160W$+$F125W | Chrimes et al. (2022)
SGR 0755(–2933)$\star\star$ | 9.52$\pm$0.01 | 9.69$\pm$0.01 | 2MASS |  $H$,$J$ | Doroshenko et al. (2021)
AX J1845(.0–0258) |  $>$21 | - | NTT | $H$ | Israel et al. (2004)
SGR 2013($+$34) |  $>$18.5 |  $>$19.3 | PAIRITEL |  $H$,$J$ | Bloom (2005)

$\star$ \- Reddest object in large error circle, unlikely to be the
counterpart (Chrimes et al., 2022), $\star\star$ \- a possibly unrelated HMXB
(Doroshenko et al., 2021)
$\dagger$ \- has HST imaging in Chrimes et al. (2022), not used due to non-
detections or single-band imaging

##  4 Extinctions and distances

Since the photometry in Table 1 is reported in terms of dust attenuated
apparent magnitudes, we must now correct the observed photometry for the
distance and dust extinction along the magnetar sight-lines, in order to
compare with the un-attenuated, absolute magnitude predictions of BPASS.

We adopt the primary distance estimates listed in the McGill catalogue
(Olausen & Kaspi, 2014, and reference therein)2. These are derived using a
variety of methods, including neutral hydrogen column densities/red clump
stars (e.g., Durant & van Kerkwijk, 2006b), dispersion measures (e.g., Levin
et al., 2010), and associations with objects that have a well-measured
distance (such as supernova remnants, molecular gas clouds, or young clusters
Corbel et al., 1999; Bibby et al., 2008). In a few cases the distance is
derived by assuming that the magnetar is in a spiral arm along that line of
sight. In this instance, and any others where the distance uncertainty is not
quantified, we assume a 15 per cent error (the mean uncertainty of those that
are quantified, and large enough to account for spiral arm width). For SGR
1935, we adopt the estimate of Bailes et al. (2021), which is derived from the
dispersion measure and is constrained to $1.5<d<6.5$ kpc. All distance
estimates adopted and their sources are provided in Table 2.

To infer the reddening due to dust, $E(B-V)$, we try two methods in the
following order of preference,

  1. 1.

If the magnetar lies above declination –30$\degr$, we check the extinction
estimate from the Bayestar 3D dust map of Green et al. (2019, G19), using the
dustmaps python package (Green, 2018). These maps do not have 100 per cent
coverage of the Galaxy, and are only reliable out to a certain distance along
each sight-line because stars (in the input surveys Gaia, Pan-STARRS and
2MASS) need to be detected in order to infer an extinction. To determine
whether the extinction estimates are reliable, we compare the total cumulative
extinction along each sight-line to the 2D extinction estimate from the
Schlafly & Finkbeiner (2011, SF11) dust maps. Where these totals are
comparable, we adopt the 3D map estimate at the distance of the magnetar. If
there is a large discrepancy (greater than a 50 per cent difference in $A_{\rm
V}$) we move on to the next method. An extreme example is SGR 1745 towards the
Galactic centre, where we have $A_{\rm V}\sim 3$ from Bayestar and $A_{\rm
V}\sim 300$ from the 2D map.

  2. 2.

The neutral hydrogen column density $N_{\rm H}$ is known to be proportional to
visual extinction $A_{\rm V}$. Predehl & Schmitt (1995, PS95) find $A_{\rm
V}=N_{\rm H}$/$[(1.79\pm 0.03)\times 10^{21}\rm{cm}^{2}]$. Many magnetars have
an X-ray derived $N_{\rm H}$ estimate, but the uncertainties on this number
are typically large, and the $N_{\rm H}$-$A_{\rm V}$ relation itself has some
scatter. This method therefore more accurately probes the direct sight-line to
the magnetar, but is typically less precise.

For the first method, the uncertainty on the extinction is correlated with the
uncertainty on the distance. For example, if a 3D estimate is used, the upper
bound on the extinction is the dust map value at $d+{\delta}d$. The $N_{\rm
H}$ inferred $A_{\rm V}$ values depend on X-ray measurements and we treat them
independently from the distance estimate.

In Figure 3 we plot $N_{\rm H}$ versus $A_{\rm V}$ for each magnetar with a
$N_{\rm H}$ measurement. If a 3D dust map or cluster-based extinction is used,
that value is plotted. For example, we adopt a value based on observations of
the Wd1 cluster for CXOU J1647 (Clark et al., 2014). If an $A_{\rm V}$-$N_{\rm
H}$ value is adopted, the 2D value for that sight-line is indicated. The
Predehl & Schmitt (1995) relation is also shown. The quoted scatter on the
$A_{\rm V}$-$N_{\rm H}$ relation is a significant underestimate, if the
differences between the 3D dust map values and $N_{\rm H}$ derived values are
representative. The difference might partly arise due to gas and dust local to
the magnetars, not captured by the dust map spatial resolution (variable, but
between 3-14 arcmin for Bayestar), explaining the tendency for the 3D dust map
values to be lower than the $N_{\rm H}$ estimates. The standard deviation of
the relative differences between the 3D dust map and NH derived extinctions,
i.e. $\Delta A_{\rm V}$/$A_{\rm V}$, is 0.34. This is indicated in Figure 3 by
the grey shaded region on the $A_{\rm V}$-$N_{\rm H}$ relation. We plot this
region equally either side of the relation, although there is a clear bias
towards $A_{\rm V}$-$N_{\rm H}$ values exceeding dust map values.

The total 2D extinctions greatly exceed the NH estimates as expected, bar one
case at low $A_{\rm V}$, putting an upper limit on the degree to which local,
unresolved dust/gas enhancement is contributing to the higher NH derived
estimates. Since the NH and 3D map methods agree to within $\sim$30 per cent,
this also implies that the distance estimates are reasonable, since the
extinction is a function of distance for the 3D map estimates. In each case,
we choose the most reliable method as described above, where the 3D map values
are correlated with distance and both the quoted scatter and NH uncertainties
are included if Predehl & Schmitt (1995) is used. We note that in both cases,
the extinction uncertainties may be underestimated, but that they will be no
more than $\sim$30 per cent in any case.

Finally, we convert visual $A_{\rm V}$ extinctions into $H$,$J$, F160W, F140W,
F125W and F110W extinctions using the effective wavelengths of these filters
and the python extinction module (Barbary, 2016), assuming R${}_{\rm V}=3.1$
and a Fitzpatrick extinction curve (Fitzpatrick, 1999). The extinction values
using each method, and the adopted values, are listed in Table 2.

Figure 3: Extinction as a function of NH derived from (i) magnetar distance
estimates and the 3D maps of Green et al. (2019, G19), or a cluster
association-based extinction estimate if available, (ii) the NH-$A_{\rm V}$
relation of Predehl & Schmitt (1995, PS95) and (iii) the total integrated line
of sight extinction out of the Galaxy Schlafly & Finkbeiner (2011, SF11). We
adopt 3D dust map estimates if they are reliable, if not and no cluster-based
estimate is available, then the $N_{\rm H}-A_{\rm V}$ relation is used. The
shaded region around the AV-NH relation corresponds to the typical $\sim$0.3
relative difference between 3D dust map and NH derived estimates.

##  5 Comparison with population synthesis predictions

We now compare the observed colours and magnitudes of the candidate bound
companions to the predictions from Section 2.

###  5.1 H-band magnitude comparisons

Using the distance estimates in Table 2, the absolute magnitudes of NIR
counterpart candidates can be calculated. As a first step in comparing the
predictions of Section 2 to observations, we instead correct the $H$-band
absolute magnitude distribution for bound companions (Figure 1) for the
distance and extinction of each magnetar. These dust attenuated apparent
magnitudes are compared to observations in Figure 4. Dashed vertical lines
enclose 95 per cent of the $H$-band magnitude probability. Uncertainties on
the correction to the BPASS outputs are indicated by error-bars on the lower
(brighter) percentile. In calculating these uncertainties (provided in Table
2), we assume for dust-map derived extinctions that the distance and
extinction uncertainties are correlated, i.e. the upper distance error is used
with the upper extinction error to derive the maximum correction within the
uncertainties, and vice versa. In every other case we assume that the
uncertainties are uncorrelated, taking the average of the upper and lower
distance and extinction uncertainties if they are asymmetric. The observed
magnitudes are either from HST measurements (F140W is drawn on this figure as
F160W, F110W as F125W) or other observations from the literature (if not HST
they are assigned to $H$,$J$ or $K$, whichever is the closest match) as
indicated in Table 1. Although SGR 1806 does not have a $H$-band detection, if
we assume it has a similar intrinsic $H$–$K=1$ colour to other counterparts,
combined with 2 magnitudes more extinction in $H$ than $K$ ($A_{\rm V}=$29),
this places it a few magnitudes outside the lower (in terms of luminosity)
percentile.

Table 2: Details for each of the 23 magnetars with a distance estimate and
imaging, including distance estimates, and extinctions calculated using the
two methods outlined in the text. The majority of this information is as
provided in the McGill magnetar catalogue and references therein Olausen &
Kaspi (2014)2, with some additions and updates. The extinction used for each
magnetar is in bold, and listed in the final column with the uncertainty. For
dust-map derived extinctions, the upper extinction uncertainty corresponds to
the upper distance uncertainty, and vice versa.

Magnetar | Dist. | Distance reference | Total $A_{\rm V}$ |  $A_{\rm V}$ at $d$ |  $A_{\rm V}$-$N_{\rm H}$ |  $A_{\rm V}$ used
---|---|---|---|---|---|---
| [kpc] |  | SF11 (tot. G19) | G19 | PS95 | with uncertainty
4U 0142 | 3.6$\pm$0.4 | Durant & van Kerkwijk (2006b) | 4.90 (4.93) | 3.79 | 5.59 | 3.79${}_{-0.25}^{+0.40}$
SGR 0418 | 2$\pm$0.3${\dagger}$ | Van der Horst et al. (2010) | 2.25 (2.17) | 1.55 | 0.64 | 1.55$\pm$0.01
SGR 0501 | 2$\pm$0.3${\dagger}$ | Lin et al. (2011) | 3.45 (3.62) | 2.73 | 4.92 | 3.79${}_{-0.39}^{+0.15}$
IE 1048 | 9.0$\pm$1.7 | Durant & van Kerkwijk (2006b) | 3.75 | - | 5.42 | 5.42$\pm$0.11
1E 1547 | 4.5$\pm$0.5 | Tiengo et al. (2010) | 43.42 | - | 17.88 | 17.88${}_{-0.11}^{+1.15}$
PSR J1622 | 9$\pm$1.35${\dagger}$ | Levin et al. (2010) | 53.22 | - | 30.17 | 30.17${}_{-7.84}^{+8.94}$
SGR 1627 | 11$\pm$0.3 | Corbel et al. (1999) | 125.4 | - | 55.87 | 55.87${}_{-11.22}^{+11.21}$
CXOU J1647 | 3.9$\pm$0.7 | Kothes & Dougherty (2007) | 30.61 | 10.66 $[1]$ | 13.35 | 10.66$\pm$0.40
1RXS J1708 | 3.8$\pm$0.5 | Durant & van Kerkwijk (2006b) | 53.84 | - | 7.60 | 7.60${}_{-0.26}^{+0.25}$
CXOU J1714 | 13.2$\pm$1.98${\dagger}$ | Tian & Leahy (2012) | 31.36 | - | 22.07 | 22.07${}_{-0.87}^{+0.91}$
SGR J1745 | 8.3$\pm$0.3 | Bower et al. (2014) | 301.66 (3.34) | unreliable | 75.42 | 75.42${}_{-3.53}^{+3.58}$
SGR 1806 | 8.7${}^{+1.8}_{-1.5}$ | Bibby et al. (2008) | 45.5 (6.85) | 29 [2] | 38.55 | 29$\pm$2
XTE J1810 | 2.5${}^{+0.4}_{-0.3}$ | Ding et al. (2020) | 60.05 (5.43) | unreliable | 3.52 | 3.52$\pm$0.29
Swift J1822 | 1.6$\pm$0.3 | Scholz et al. (2012) | 12.75 (9.73) | 2.21 | 2.53 | 2.21${}_{-1.32}^{+1.18}$
Swift J1834 | 4.2$\pm$0.3 | Leahy & Tian (2008) | 96.78 (4.55) | unreliable | 98.27 | 98.27${}_{-20.13}^{+19.11}$
IE 1841 | 8.5${}^{+1.3}_{-1.0}$ | Tian & Leahy (2008) | 52.86 (5.61) | unreliable | 12.29 | 12.29${}_{-0.59}^{+0.60}$
3XMM J1852${\ddagger}$ | 7.1$\pm$1.07${\dagger}$ | Zhou et al. (2014) | 65.17 (9.89) | unreliable | 7.63 | 7.63${}_{-0.56}^{+0.05}$
SGR 1900 | 12.5$\pm$1.7 | Davies et al. (2009) | 11.05 (10.94) | 10.60 | 11.84 | 10.60${}_{-0.03}^{+0.07}$
SGR 1935 | 4.0$\pm 2.5$ | Bailes et al. (2021) | 13.13 (8.58) | 6.11 | 8.94 | 6.11${}_{-0.31}^{+0.28}$
IE 2259 | 3.2$\pm$0.2 | Kothes & Foster (2012) | 3.82 (3.22) | 2.54 | 5.64 | 2.54${}_{-0.03}^{+0.04}$
SGR 0755 | 3.5$\pm$0.2 | Doroshenko et al. (2021) | 2.6 (2.39) | 1.77 | 0.73 | 1.77${}_{-0.03}^{+0.09}$
AX J1845 | 8.5$\pm$1.28 | Torii et al. (1998) | 50.35 (5.58) | unreliable | 43.58 | 43.58${}_{-10.09}^{+12.87}$
SGR 2013 | 8.5$\pm$1.32 | Sakamoto et al. (2011) | 5.64 (6.20) | 5.36 | 6.13 | 5.36${}_{-1.24}^{+0.56}$

${\dagger}$ \- distance uncertainty unquantified, 15 per cent assumed,
${\ddagger}$ \- Reddest source in error circle
$[1,2]$ \- extinctions from cluster associations (Clark et al., 2014;
Eikenberry et al., 2004).

Figure 4: Photometry for the 23 magnetars with a distance estimate, compared
to BPASS predictions for the apparent magnitudes of bound companions
(histogram), assuming the distance and extinction estimates in Table 2. The
2.5 and 97.5 percentiles (enclosing 95 per cent of the distribution) are
indicated by vertical lines, the lower percentile has an errorbar indicating
the uncertainty on the distribution due to distance and $A_{\rm V}$
uncertainties. Previously reported apparent magnitudes (uncorrected) are
indicated by the black points ($K$-band, arbitrarily placed at $y=0.25$), red
($H$-band, 0.2) and blue ($J$-band, 0.05). If the source is variable and a
range of values have been reported, the brightest and faintest magnitudes
previously reported are joined by a line. If two separate sources have
previously been measured as candidates, they are plotted separately. Triangles
represent upper limits. The HST F160W (F125W) photometry listed in Table 2 is
plotted in orange (cyan) at $y=0.15$ ($y=0.01$).

Figure 5: A NIR colour magnitude diagram for a Solar metallicity stellar
population at 10 Myr. The right hand panel zooms out to show SGR 1806, PSR
J1622, SGR 1627 and Swift J1834 (whose extinctions are large and uncertain).
Extinctions derived from 3D dust map/cluster extinctions, or the $A_{\rm
V}$-NH relation, are labelled with orange circles or cyan squares
respectively. Immediately post primary supernova, 95 per cent of all bound
companion $H$-band magnitudes lie in the region bounded by the horizontal dot-
dashed lines, and 95 per cent of their $J$–$H$ colours lie between the
vertical dashed lines. The reddening vector is indicated in the upper right -
for a fixed apparent magnitude and observed colour, the inferred absolute
magnitude and intrinsic colour move in this direction with increasing
extinction. Objects that lie away from the shaded regions either have
intrinsically non-stellar colours, or the incorrect distance and extinction
applied (i.e., they are chance alignments). The underlying density contours
(an order of magnitude in number density per level) were plotted with the
python package for BPASS, HOKI (Stevance et al., 2020).

###  5.2 Colour-magnitude comparisons

While informative, single-band magnitudes are only part of the information
available. To add further constraints, we now examine the colours of the
counterpart and candidates. For those sources which have observations in $H$
and $J$ and a detection in at least one, we apply the distance and extinction
corrections to the magnetar photometry and place this on a colour-absolute
magnitude diagram in Figure 5. Here we show BPASS predictions for a Solar
metallicity stellar population at 107 yrs. Ages of a fews tens of Myr
correspond to the expected lifetimes of supernova progenitors, hence our
choice of 107 yr, but comparisons with populations at 106 yr and 108 yr are
also provided in Appendix A.

The corrected magnetar measurements are coloured by the method used to derive
the extinction. Horizontal dashed lines mark the upper and lower luminosity
percentiles enclosing 95 per cent of the distribution, previously adjusted in
each case in Figure 4. Vertical dashed lines represent the same for the
predicted $J$–$H$ colour distribution of bound companions. The box enclosed by
these 4 lines is therefore where companion stars are predicted. However, we
plot the whole stellar population at the given age underneath for two reasons.
Firstly, this is in order to be agnostic about the properties of any bound
companions, which represent a subset of the wider stellar population. For
example, as companions to massive stars, they will be biased to higher masses
since very high mass ratio systems are disfavoured (Moe & Di Stefano, 2017).
However, by plotting the whole population, we are open to the possibility of
companions outside of the predicted range, either because the predictions are
inaccurate, or because magnetar progenitor systems are atypical. Secondly, if
a source does not fall on the shaded regions (where each shading level
represents an order of magnitude difference in the number density of sources),
this is indicative of either,

  * •

Genuinely non-stellar colours, assuming the adopted $d$ and $A_{\rm V}$ are at
least approximately correct, or,

  * •

A chance alignment, where the source in the localisation error circle is
unrelated to the magnetar, and hence the wrong distance and extinction is
being applied,

where ‘non-stellar’ is defined as colour-magnitude properties which are
inconsistent with stars undergoing nuclear fusion (any phase from the main
sequence through to core-collapse or white dwarf formation), or white dwarfs.
Examples of non-stellar emission in this context include magnetospheric or
debris disc emission. We now discuss the placement of the counterparts in this
parameter space after distance and extinction correction, and the implications
in each case.

There are three cases where the wrong corrections have likely been applied for
the sources in the error circle. These are SGR 1627, PSR J1622 and Swift
J1834, which have very large uncertainties due to large and uncertain NH
measurements. None the less, they are unphysically blue (the right panel of
Fig. 5 is zoomed out to show this). Given the high extinction on these sight-
lines, it is likely that these are foreground chance alignments, and therefore
the distance and extinction corrections are much too high. This is also seen
by the direction of the extinction vector, which traces back to the main
sequence - the region of highest probability density - for these sources.

At first glance, the three sources in (or just outside, we have added source B
of Israel et al. (2003) and Durant & van Kerkwijk (2006a) to compare) the 1RXS
J1708 error circle also look like chance alignments, since they cannot all be
counterparts. However, Israel et al. (2003) claim that their source A (the
brightest of the three) is non-stellar in $J$–$H$, $J$–$K$ colour space.
Durant & van Kerkwijk (2006a) showed that source B can only be shifted into
the cloud of stellar sources if an $A_{\rm V}$ of ${\sim}$20 is assumed,
inconsistent with the column density extinction estimate (and consistent with
the main sequence offset shown in Figure 5). It is therefore unclear which is
associated with 1RXS J1708, but it is likely that two are anomalously red
because they are background chance alignments with an insufficient
distance/extinction applied, whereas the other is the magnetar counterpart,
with a non-stellar SED. AX J1818.0-1607 is a similar case, for which this was
shown to be a plausible scenario by Chrimes et al. (2022). In any case, all of
the 1RXS J1708 candidates are outside the bound companion magnitude range.

There are another nine sources below the population synthesis bound companion
limit. 4U 0142, 1E 2259 and SGR 0501 are unambiguously the counterpart based
on precise localisations and a lack of other sources in the immediate
vicinity. 1E 2259 stands out as being unusually red, whereas 4U 0142 does not,
despite previous claims of a dusty debris disc. We attribute this to the $J$
and $H$ bands used, which are not as sensitive in this case to the upturn
which becomes more significant from the $K$-band into the mid-infrared (e.g.
Muñoz-Darias et al., 2016). The source associated with Swift J1822 is also in
this faint, red area of parameter space.

1E 1841, 1E 1048 and SGR 1806 all have $J$-band non-detections and hence their
colours are lower limits. For 1E 1048 this is a constraining limit which
places it in the same faint, red region of parameter space as the other
counterparts. 1E 1841 has $H$ and $K$ imaging presented by Testa et al. (2008)
and the likely counterpart is identified as their source 9, due to 4$\sigma$
$K$-band variability. We therefore adopt the $H$ and $J$-band measurements of
the same object in Durant (2005), labelled there as source B. The shallow
$J$-band limit means that the $J$–$H$ colour is relatively poorly constrained.
The source placement is outside of the predicted companion magnitude range and
marginally consistent with the predicted colours. For SGR 1806, which has a
$K$-band detection (Tendulkar et al., 2012) and limits in $H$ and $J$ (Israel
et al., 2005), we first assume a $H$–$K$ colour of one. This is typical of the
intrinsic colours of known counterparts (see e.g. the least extincted examples
in Figure 4) and is consistent with the $H>19.5$ limit. Correcting for the
distance and extinction then places it at M${}_{H}\sim$5\. Given the $J>21.1$
limit, the high extinction ($A_{\rm V}=29$) yields an unconstraining
$J$–$H>-3.6$ colour limit.

Although there is HST imaging for CXOU 1647, this is only in one band (F140W,
Chrimes et al. (2022)), so we adopt the $J=23.5\pm 0.2$ and $H=21.0\pm 0.1$
measurements of (Testa et al., 2018) instead. The last of the eight sources
with colour information to lie below the population synthesis prediction for
bound companion $H$-band magnitudes is XTE J1810. The counterpart is not
detected in 2018 HST F160W and F125W imaging (Chrimes et al., 2022), so we use
measurements from Testa et al. (2008) instead. Because the source is variable,
we chose just the first epoch (MJD 52900) of observations presented
($H=21.67\pm 0.12$, $J=22.92\pm 0.22$).

Finally, there are two sources which are notably different from the cloud of
objects in the lower right of Figure 4. Not only do they lie in the magnitude
range predicted for companions, but they also have stellar colours. These are
SGR 0755 and CXOU J1714. The star associated with SGR 0755, as previously
discussed, is too bright to measure from the HST images, so 2MASS (Skrutskie
et al., 2006) magnitudes are adopted. This system, 2SXPS J075542.5–293353, is
an HMXB (Doroshenko et al., 2021). However, it remains an open question
whether the accretor could be a magnetar (Popov, 2022), or whether the
magnetar burst came from a different source in the BAT error circle. The
counterpart candidate for CXOU J1714 is unambiguous, well localised, and has
stellar colours, albeit slightly bluer than predicted for neutron star
companions. The source also stands apart from previously established
counterparts in terms of its X-ray to NIR power law index (Chrimes et al.,
2022). Unless the distance and extinction are drastically overestimated for
this magnetar, the source is clearly distinct from previously confirmed non-
stellar counterparts.

##  6 Summary of results

Out of the 23 sources considered in this sample, 3 are unconstrained in terms
of their $H$-band magnitude - AX J1845 only has a limit that does not preclude
a bound companion, SGR 1745 also has an unconstraining limit, and we deem 3XMM
J1852 as unconstrained given the poor localisation. The overall results,
including colour information where available, are summarised in Table 3. Here
we list whether the observations for each source can rule out or allow for the
possibility of a bound companion. Of the 20 that have a magnitude constraint,
2 are plausibly bound companions based on their absolute magnitude and colour.
These are SGR 0755 and CXOU J1714. CXOU J1714 is slightly outside the colour
predictions but is otherwise well within the magnitude range expected. SGR
0755 satisfies both criteria. IE 1841 is not far ($\sim 2\sigma$) outside the
lower magnitude bound, but the colour limit also precludes most of the $J$–$H$
distribution, so we deem it inconsistent with predictions.

Of the two candidates, which appear distinctly bright and blue in this
parameter space, SGR 0755 is an X-ray binary, and it is possible that the
magnetar is actually elsewhere in the BAT error circle (Doroshenko et al.,
2021). CXOU J1714 is perhaps the most promising bound companion candidate. We
can therefore constrain the fraction of magnetars with a bound companion to
0.05–0.10 (1–2 out of 20). This range is quoted at 95 per cent confidence
since both the localisation uncertainties and bound companion magnitude
predictions are at 95 per cent. This bound companion fraction is consistent
with population synthesis predictions for a regular core-collapse origin, as
discussed in Section 2. If there is a contribution amongst this population
from non-core-collapse channels, it cannot be very large before the bound
companion fraction among the remainder disagrees with predictions. For
instance, if 5/20 do not have a core-collapse origin, the bound fraction
amongst the rest would be 7–13 per cent, which would imply a bias towards
magnetar production in primary stars, assuming that $f_{\rm bound}=0.05$ is
accurate for the wider neutron star population.

Table 3: A summary of the results. Whether each source is plausibly a bound
companion based on the $H$-band or location on the colour-magnitude diagram
(CMD), is indicated. The tolerance for being consistent in Figure 4 is
extended beyond the plotted percentiles according to the indicated uncertainty
on the distance and extinction corrections. 3/23 are unconstrained. Of the
remaining 20, two have apparently stellar properties broadly consistent with
expectations for bound companions (CXOU J1714 and SGR 0755, for which the
bright star is a confirmed HMXB, but may but a chance alignment with the
magnetar in question).

Magnetar | Plausible | Plausible | Classification
---|---|---|---
| from H mag? | from CMD? |
4U 0142 | N | N | Non-stellar
SGR 0418 | N | - | Non-stellar
SGR 0501 | N | N | Non-stellar
1E 1547 | N | - | Non-stellar
IE 1048 | N | N | Non-stellar
PSR J1622 | Y | N | Non-stellar${\star\star}$
SGR 1627 | Y | N | Non-stellar${\star\star}$
CXOU J1647 | N | N | Non-stellar
1RXS J1708 | N | N | Non-stellar
CXOU J1714 | Y | Y | Possibly stellar
SGR J1745 | Unconstrained | - | Unconstrained
SGR 1806 | N⋆ | N | Non-stellar
XTE J1810 | N | N | Non-stellar
Swift J1822 | N | N | Non-stellar
Swift J1834 | N | N | Non-stellar⋆⋆
IE 1841 | N | N | Non-stellar
3XMM J1852 | Unconstrained | - | Unconstrained
SGR 1900 | N | - | Non-stellar
SGR 1935 | N | - | Non-stellar
IE 2259 | N | N | Non-stellar
SGR 0755 | Y | Y | Stellar
AX J1845 | Unconstrained | - | Unconstrained
SGR 2013 | N | - | Non-stellar

${\star}$ \- Assuming an intrinsic $H$–$K=1$
${\star\star}$ \- extremely blue unattenuated colours, high likelihood of
chance alignment

##  7 Discussion

###  7.1 Bound companion candidates

####  7.1.1 SGR0755

A system of particular interest here is 2SXPS J075542.5–293353, a HMXB which
may (or may not) be harbouring SGR 0755 (Doroshenko et al., 2021). If this is
a magnetar accreting from a high mass companion, it would be the first such
system identified (some works claim that such systems are highly unlikely to
form, since accretion would dampen the magnetic field, or at least that
accretion does not power magnetar emission e.g. King & Lasota, 2019;
Doroshenko et al., 2020). While it may be a chance coincidence that this HMXB
lies in the BAT error circle derived from the SGR 0755 discovery burst, it
could also be that the burst came from this system even if does not contain a
magnetar (there are several X-ray pulsars known in HMXBs e.g. Salganik et al.,
2021). There are also models which predict magnetar survival (i.e. avoiding
rapid magnetic field decay) in certain accreting systems (Igoshev & Popov,
2018; Bozzo et al., 2022). For example, the orbital and spin periods measured
for 2SXPS J075542.5–293353 agree with the magnetar X-ray binary model
predictions of Xu et al. (2021). Unfortunately, we cannot compare the age of
the magnetar in this system with an estimate age of the star, because the
$P/\dot{P}$ age becomes unreliable in an accreting system. Furthermore, the
lack of an independent magnetar distance estimate precludes a comparison with
the Gaia distance of 3.5$\pm$0.2 kpc (EDR3, Gaia Collaboration et al., 2016,
2021).

####  7.1.2 CXOUJ1714

The other magnetar with a coincident NIR source consistent with being a
massive stellar companion is CXOU J1714. It has an inferred absolute magnitude
consistent with BPASS predictions for a bound companion, and although bluer
than the predicted colours, it is still consistent with being stellar in
nature. This source is much fainter than the donor in 2SXPS J075542.5–293353,
and has not yet been unambiguously identified as stellar in nature. In
principle, there might be an orbital period signature in the X-ray lightcurve
of the magnetar, if the emission is somehow coupled to interactions with the
companion. Figure 1 shows the BPASS prediction for the orbital periods of
bound neutron star-massive companion systems. X-ray observations of CXOU J1714
are somewhat sparse, with 11 epochs since its identification as a magnetar in
2010 (Sato et al., 2010). There are also a handful of preceding epochs, where
the target was the surrounding supernova remnant CTB 37B (Gotthelf et al.,
2019). The 3.8s period of the magnetar is clear, and there is clear evidence
for variability over timescales of $\sim$years within these data (Sato et al.,
2010). However, whether there is day-week timescale periodicity in this long
term variation, as predicted in a bound system, remains to be seen.

We also search for catalogues containing the NIR source associated with CXOU
J1714. A Gaia star lies 2.7 arcsec away from the centre of the CXO
localisation of CXOU J1714, and appears in the Gaia \- neutron star
association catalogue of Antoniadis (2021). No proper motion or parallax is
available for this source, due to its faint 20.9 $G$-band magnitude. Given the
good CXO localisation of the magnetar (ruling the star out as a bound
companion), and given the extinction of $A_{\rm V}\sim 22$ along the magnetar
sight-line, it would have an unphysically luminous absolute $G$-band magnitude
of $\sim$–10 if placed at the same distance and extinction as the magnetar
(i.e., as a companion).

The only optical/NIR catalogue in which the NIR source at the position of CXOU
J1714 is detected, is the VVV Infrared Astrometric Catalogue (Smith et al.,
2018). It has VVV magnitudes of $H=17.3\pm 0.1$ and $J=19.1\pm 0.3$, agreeing
well with our measurements (see Table 2), given the bandpass differences. The
source has a VVV proper motion measurement of $6.3\pm 3.1$ mas yr-1. Following
Verbunt et al. (2017) and using the same model parameters as Lyman et al.
(2022), we correct for the relative motion due to Galactic rotation and the
peculiar velocity of the Sun. The IAU standard of 220 kms-1 is assumed for the
Galactic rotation at both the Solar location and the source (reasonable given
that all locations considered are outside the central few kpc of the Galaxy,
Clemens, 1985). The resultant transverse velocity in the local standard of
rest (LSR) of the source is $349\pm 151$ kms-1, where the 1 $\sigma$
uncertainty is derived from 105 Monte Carlo trials with Gaussian uncertainties
on U, V and W (the Solar peculiar motion components), $\mu_{\alpha}$ and
$\mu_{\delta}$ (the observed proper motion components) and distance $d$. This
is a high value which is more typical of the Galactic pulsar population, but
the lower limit is consistent with the upper end of compact remnant-massive
star bound systems, specifically high-mass X-ray binaries (Mirabel et al.,
2002; Atri et al., 2019; Igoshev et al., 2021). If we instead adopt a distance
of 5 kpc - the lowest distance estimated for CXOU J1714 before it was revised
upwards by association with the CTB 37 supernova remnant complex (Tian &
Leahy, 2012) \- the transverse velocity estimate drops to $140\pm 67$ kms-1,
entirely consistent with bound systems. If we change the distance, the
inferred absolute magnitude is only affected by a change in distance modulus
because the extinction is derived from the column density. Consequently, at 5
kpc, the source is 2.1 magnitudes fainter, which still places it firmly in the
stellar region of parameter space in Figure 5. Although a lower distance is
preferred given the VVV proper motion measurement of the NIR source and the
typical systemic velocities of X-ray binaries, the large uncertainties mean a
reasonable velocity cannot be ruled out even at $d=13.2$ kpc.

###  7.2 Pre-main sequence stars

Another possibility, if the companion is significantly less massive than the
magnetar progenitor, is that it may still be in a pre-main sequence (PMS)
period of evolution. The NIR variability of PMS stars is well-established
(e.g. Carpenter et al., 2001; Eiroa et al., 2002; Alves de Oliveira & Casali,
2008), and several magnetar counterparts are variable. The variations in PMS
stars are thought to be due to changes in the disc accretion rate, structure,
and clumpy circumstellar dust distributions. Although the time-scales of PMS
evolution of low-mass stars are indeed comparable to the entire life-times of
massive O stars (tens of Myr; Palla & Stahler, 1993) – with the magnetar’s
life-time adding an insignificant amount of time – this scenario would invoke
the need for a somewhat extreme mass ratio binary. This scenario is
disfavoured due to the initial binary parameter distributions (Sana et al.,
2012), and the preference for these systems to unbind upon primary supernova.
However, we know such systems can occur, due to the existence of low-mass
X-ray binaries, so the pre-main sequence scenario cannot be ruled out.

###  7.3 Comparison with pulsars

In addition to the population synthesis predictions of Section 2, we can also
compare our results to the pulsar population. Seven per cent of pulsars in the
ANTF catalogue are listed as having a main sequence companion (Manchester et
al., 2005)333https://www.atnf.csiro.au/research/pulsar/psrcat/. However, this
number is highly uncertain due to detection and follow-up biases - some may be
obscured by dust, while others have not had deep observations to establish the
presence of a companion, or lack thereof. Antoniadis (2021) performed a cross-
match between the Gaia DR3 catalogue and 1534 pulsars. They find that the
binary fraction of young pulsars (defined as being un-recycled and with
magnetic field strength $>10^{10}$ G) is $<5$ per cent, under the assumption
that all pulsar stellar companions with orbital periods of $<50$ yr and
brighter than 20.5 mag should have been identified. However, some systems with
bound companions will manifest as X-ray binaries. We note that the ages of
HMXBs, typically tens of Myrs, (inferred from kinematics and cluster
associations, e.g. Coleiro & Chaty, 2013) are more typical of (young) pulsars
(Manchester et al., 2005) than magnetars (Olausen & Kaspi, 2014). This might
suggest that more pulsar bound companions will be ‘missing’ (because these
systems are instead identified as XRBs) than magnetar bound companions. The
lifetimes of bound companions also suggest that they will leave the main
sequence on a similar time-scale to the ages of HMXBs, consistent with this
scenario. Overall, the population synthesis, pulsar and magnetar estimates for
$f_{\mathrm{bound}}$ are all broadly consistent with ${\sim}$5 per cent.

###  7.4 Other constraints and future progress

The origin of magnetars has been probed in a variety of ways, from
investigating their colours/spectral energy distributions (SEDs) as in this
work, their proper motions, and associations with recent star
formation/supernovae. The results presented here do not require non-core-
collapse channels to explain the Galactic population, and this is backed up by
previous studies which point towards core-collapse as the dominant channel
(e.g. Beniamini et al., 2019; Burns et al., 2021). However, some contribution
from other channels cannot be excluded. For example, accretion or merger
induced collapse magnetars are predicted on theoretical grounds (Duncan &
Thompson, 1992; Levan et al., 2006; Margalit et al., 2019; Gautam et al.,
2021; Ablimit, 2021), and are a possible explanation for FRB202001E, which
resides in a globular cluster (Kirsten et al., 2021; Lu et al., 2021). If the
majority of magnetars are born in core-collapse events, it is still uncertain
whether this is a typical outcome, or if it requires specific pre-collapse
conditions (this is largely dependent on how rapidly their magnetic fields
decay, Nakano et al., 2015; Beniamini et al., 2019).

Another avenue to compare magnetars with other neutron stars is their proper
motions. Thus far a handful have measured proper motions, (Tendulkar et al.,
2012; Deller et al., 2012; Tendulkar et al., 2013; Ding et al., 2020). Lyman
et al. (2022) also add the proper motion of SGR 1935 to this distribution,
making use of the stability and spatial resolution of HST. The sample remains
small, but they show that there is currently no statistically significant
difference between magnetar and pulsar proper motions.

Potential further tests include searching for unbound (runaway) companions. A
runaway companion to magnetar CXOU J1647 in the young star cluster
Westerlund-1 has been claimed (Wd1-5, a blue supergiant, Clark et al., 2014).
Due to the carbon rich composition and over-luminous nature of Wd1-5, they
infer an evolutionary history in which the (initially slightly less massive)
secondary evolves more rapidly following accretion, polluting the primary with
carbon rich Wolf-Rayet winds, before going supernova and unbinding the system.
Wd1-5 has H${\sim}$8.5, at the upper end of the population synthesis
predictions for bound companions (Fig. 4). We predict that 90 per cent of
primary neutron stars should have unbound from their companion following the
supernova (47 per cent of core-collapse neutron stars overall, including
single stars and mergers). The unbound companions will have a range of
velocities, and some will be runaways with velocites of hundreds of kms-1.
These numbers are less prone to low number statistics and may provide a more
robust constraint on magnetar progenitors, but the large number of possible
matches means that finding magnetar-secondary pairs is not trivial. Wd1-5 was
selected due to its high radial velocity from spectroscopy, but tracing back
proper motions may still be a promising route to identifying associations
(e.g. Fraser & Boubert, 2019). We leave searches for unbound runaway
companions, which requires significant additional analyses beyond the scope of
this paper, for further studies. However, we qualitatively note that the bound
companion search presented here also rules out a fraction of this parameter
space, namely unbound companions with small natal kicks.

Another constraint is the masses of the progenitors. For magnetars in a
stellar cluster, these can be inferred by ageing the cluster with the red
supergiant luminosity distribution. Assuming (i) that the magnetar age is
(much) less than the progenitor evolutionary time, and (ii) that the cluster
is coeval, the age of the cluster corresponds to a progenitor mass range (e.g.
Davies et al., 2009). Progenitor mass estimates range from $>40$ M⊙
(CXOUJ1647, SGR1806 Muno et al., 2006; Bibby et al., 2008) down to
${\sim}$15–20M⊙ (SGR1900,1E1841 Davies et al., 2009; Borkowski & Reynolds,
2017). Supernova remnant associations can also be used to estimate the
progenitor mass, based on modelling of the remnant itself. For example, Zhou
et al. (2019) use spatially resolved X-ray observations of supernova remnants
to infer progenitor masses of less than 20M⊙ for magnetars 1E 1841, SGR
0526–66 and 1E 161348–5055. If these progenitor mass distributions can be
expanded, comparisons can be made to the overall massive star IMF and pre-SN
masses, to see if magnetars represent a biased subset.

In this work, we have assumed Solar metallicity for the progenitors. If they
are biased to lower metallicities, this could be reflected in the progenitor
masses and delay times distributions. Quantifying the metallicity of the
parent stellar populations is plausible aim, particularly for sources which
have cluster associations, e.g. CXOU 1647 in Westerlund 1 (approximately
Solar, Crowther et al., 2006).

Finally, spectra or SEDs of magnetar counterparts can conclusively determine
their nature. They can distinguish stellar from non-stellar, particularly if
the SED extends into the $K$-band and beyond, where dusty debris discs start
to dominate the spectra (Wang et al., 2006). The nature of any stellar
companion can also be distinguished. The colours of O-stars and lower-mass
pre-main sequence stars differ considerably in the NIR, and also differ from
non-thermal/disc models. Characterisation of the NIR SED, as well as a better
understanding of the time-scales and amplitude of variability in these
sources, would therefore be able to better determine if the emission can be
explained by a companion. In principle, neutron star-mass companions could be
inferred from radial velocity measurements of bound stellar candidates, but
this becomes increasingly difficult at higher mass ratios. Given the locations
of magnetars in the Galactic plane, and the associated dusty sight-lines,
magnetar spectroscopy and even SED construction is also challenging, but the
possibilities will broaden with the advent of JWST (James Webb Space
Telescope, Gardner et al., 2006).

##  8 Conclusions

In this paper, we have employed photometry of candidate magnetar counterparts
to identify any bound companion stars among the Galactic population. Distance
and extinction estimates are used to infer intrinsic colours and absolute
magnitudes, which are compared to predictions from binary population
synthesis. We find that 5–10 per cent of the Galactic population could
plausibly have a bound companion, corresponding to one or two NIR counterparts
which appear stellar in nature. One of these is associated with SGR 0755–2933
and has been identified as a high mass X-ray binary, but it remains unclear
whether the accretor in this system is the magnetar, or if this is a chance
alignment. The other is the counterpart to CXOU J171405.7–381031. These
numbers are consistent with binary population synthesis predictions, assuming
that the dominant magnetar formation channel is through core-collapse. Further
observations are needed to determine if magnetars arise ubiquitously from
core-collapse events or a biased subset of this population, and to what extent
alternative progenitor channels contribute.

## Acknowledgements

AAC is supported by the Radboud Excellence Initiative. AJL has received
funding from the European Research Council (ERC) under the European Union’s
Seventh Framework Programme (FP7-2007- 2013) (Grant agreement No. 725246). PJG
acknowledges support from the NRF SARChI program under grant number 111692.
JDL acknowledges support from a UK Research and Innovation Fellowship
(MR/T020784/1). CK and ASF acknowledge support for this research, provided by
NASA through a grant from the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc. We
thank Rens Waters and Gijs Nelemans for useful discussions.

This work made use of v2.2.1 of the Binary Population and Spectral Synthesis
(BPASS) models as described in Eldridge et al. (2017) and Stanway & Eldridge
(2018). This work has made use of ipython (Perez & Granger, 2007), numpy
(Harris et al., 2020), scipy (Virtanen et al., 2020); matplotlib (Hunter,
2007), Seaborn packages (Waskom, 2021), astropy,444https://www.astropy.org a
community-developed core Python package for Astronomy (Astropy Collaboration
et al., 2013; Price-Whelan et al., 2018), and the python package for BPASS,
HOKI (Stevance et al., 2020). We have also made use of the python modules
statsmodels (Seabold & Perktold, 2010), dustmaps (Green, 2018), extinction
(Barbary, 2016) and photutils, an Astropy package for detection and photometry
of astronomical sources (Bradley et al., 2021). This research has made use of
the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/)
supported from the Spanish MINECO through grant AYA2017-84089 (Rodrigo et al.,
2012; Rodrigo & Solano, 2020).

This publication makes use of data products from the Two Micron All Sky
Survey, which is a joint project of the University of Massachusetts and the
Infrared Processing and Analysis Center/California Institute of Technology,
funded by the National Aeronautics and Space Administration and the National
Science Foundation. This work has made use of data from the European Space
Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the
Gaia Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

Finally, we thank the referee for their constructive comments on this
manuscript.

## Data Availability

Based on observations made with the NASA/ESA Hubble Space Telescope, obtained
from the data archive at the Space Telescope Science Institute. STScI is
operated by the Association of Universities for Research in Astronomy, Inc.
under NASA contract NAS 5–26555. These observations are associated with
programs 14805, 15348 and 16019 (Levan). The scientific results reported in
this article are based on data obtained from the Chandra Data Archive.

## References

  * Ablimit (2021) Ablimit I., 2021, arXiv e-prints, p. arXiv:2110.12140
  * Alves de Oliveira & Casali (2008) Alves de Oliveira C., Casali M., 2008, A&A, 485, 155
  * Antoniadis (2021) Antoniadis J., 2021, MNRAS, 501, 1116
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Atri et al. (2019) Atri P., et al., 2019, MNRAS, 489, 3116
  * Bailes et al. (2021) Bailes M., et al., 2021, MNRAS,
  * Barbary (2016) Barbary K., 2016, extinction v0.3.0, doi:10.5281/zenodo.804967, https://doi.org/10.5281/zenodo.804967
  * Beniamini et al. (2019) Beniamini P., Hotokezaka K., van der Horst A., Kouveliotou C., 2019, MNRAS, 487, 1426
  * Bhandari et al. (2021) Bhandari S., et al., 2021, arXiv e-prints, p. arXiv:2108.01282
  * Bibby et al. (2008) Bibby J. L., Crowther P. A., Furness J. P., Clark J. S., 2008, MNRAS, 386, L23
  * Blanton & Roweis (2007) Blanton M. R., Roweis S., 2007, AJ, 133, 734
  * Bloom (2005) Bloom J. S., 2005, GRB Coordinates Network, 4042, 1
  * Bochenek et al. (2020) Bochenek C. D., Ravi V., Belov K. V., Hallinan G., Kocz J., Kulkarni S. R., McKenna D. L., 2020, Nature, 587, 59
  * Bochenek et al. (2021) Bochenek C. D., Ravi V., Dong D., 2021, ApJ, 907, L31
  * Borkowski & Reynolds (2017) Borkowski K. J., Reynolds S. P., 2017, ApJ, 846, 13
  * Bower et al. (2014) Bower G. C., et al., 2014, ApJ, 780, L2
  * Bozzo et al. (2022) Bozzo E., Ferrigno C., Oskinova L., Ducci L., 2022, MNRAS, 510, 4645
  * Bradley et al. (2021) Bradley L., et al., 2021, astropy/photutils: 1.1.0, doi:10.5281/zenodo.4624996, https://doi.org/10.5281/zenodo.4624996
  * Bray & Eldridge (2016) Bray J. C., Eldridge J. J., 2016, MNRAS, 461, 3747
  * Bray & Eldridge (2018) Bray J. C., Eldridge J. J., 2018, MNRAS, 480, 5657
  * Briel et al. (2021) Briel M. M., Eldridge J. J., Stanway E. R., Stevance H. F., Chrimes A. A., 2021, arXiv e-prints, p. arXiv:2111.08124
  * Bucciantini et al. (2012) Bucciantini N., Metzger B. D., Thompson T. A., Quataert E., 2012, MNRAS, 419, 1537
  * Burns et al. (2021) Burns E., et al., 2021, ApJ, 907, L28
  * CHIME/FRB Collaboration et al. (2020) CHIME/FRB Collaboration et al., 2020, Nature, 587, 54
  * Carpenter et al. (2001) Carpenter J. M., Hillenbrand L. A., Skrutskie M. F., 2001, AJ, 121, 3160
  * Chrimes et al. (2020) Chrimes A. A., Stanway E. R., Eldridge J. J., 2020, MNRAS, 491, 3479
  * Chrimes et al. (2021) Chrimes A. A., Levan A. J., Groot P. J., Lyman J. D., Nelemans G., 2021, MNRAS, 508, 1929
  * Chrimes et al. (2022) Chrimes A. A., Levan A. J., Fruchter A. S., Groot P. J., Kouveliotou C., Lyman J. D., Tanvir N. R., Wiersema K., 2022, MNRAS,
  * Clark et al. (2014) Clark J. S., Ritchie B. W., Najarro F., Langer N., Negueruela I., 2014, A&A, 565, A90
  * Clemens (1985) Clemens D. P., 1985, ApJ, 295, 422
  * Coleiro & Chaty (2013) Coleiro A., Chaty S., 2013, ApJ, 764, 185
  * Corbel et al. (1999) Corbel S., Chapuis C., Dame T. M., Durouchoux P., 1999, ApJ, 526, L29
  * Crowther et al. (2006) Crowther P. A., Hadfield L. J., Clark J. S., Negueruela I., Vacca W. D., 2006, MNRAS, 372, 1407
  * Davies et al. (2009) Davies B., Figer D. F., Kudritzki R.-P., Trombley C., Kouveliotou C., Wachter S., 2009, ApJ, 707, 844
  * Deller et al. (2012) Deller A. T., Camilo F., Reynolds J. E., Halpern J. P., 2012, ApJ, 748, L1
  * Ding et al. (2020) Ding H., et al., 2020, MNRAS, 498, 3736
  * Doroshenko et al. (2020) Doroshenko V., Santangelo A., Suleimanov V. F., Tsygankov S. S., 2020, A&A, 643, A173
  * Doroshenko et al. (2021) Doroshenko V., Santangelo A., Tsygankov S. S., Ji L., 2021, A&A, 647, A165
  * Duncan & Thompson (1992) Duncan R. C., Thompson C., 1992, ApJ, 392, L9
  * Durant (2005) Durant M., 2005, ApJ, 632, 563
  * Durant & van Kerkwijk (2006a) Durant M., van Kerkwijk M. H., 2006a, ApJ, 648, 534
  * Durant & van Kerkwijk (2006b) Durant M., van Kerkwijk M. H., 2006b, ApJ, 650, 1070
  * Eikenberry et al. (2004) Eikenberry S. S., et al., 2004, ApJ, 616, 506
  * Eiroa et al. (2002) Eiroa C., et al., 2002, A&A, 384, 1038
  * Eldridge et al. (2011) Eldridge J. J., Langer N., Tout C. A., 2011, MNRAS, 414, 3501
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Australia, 34, e058
  * Eldridge et al. (2019) Eldridge J. J., Stanway E. R., Tang P. N., 2019, MNRAS, 482, 870
  * Fitzpatrick (1999) Fitzpatrick E. L., 1999, PASP, 111, 63
  * Fraser & Boubert (2019) Fraser M., Boubert D., 2019, ApJ, 871, 92
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Gardner et al. (2006) Gardner J. P., et al., 2006, Space Sci. Rev., 123, 485
  * Gautam et al. (2021) Gautam A., Crocker R. M., Ferrario L., Ruiter A. J., Ploeg H., Gordon C., Macias O., 2021, arXiv e-prints, p. arXiv:2106.00222
  * Gelfand & Gaensler (2007) Gelfand J. D., Gaensler B. M., 2007, ApJ, 667, 1111
  * Gompertz et al. (2013) Gompertz B. P., O’Brien P. T., Wynn G. A., Rowlinson A., 2013, MNRAS, 431, 1745
  * Gompertz et al. (2014) Gompertz B. P., O’Brien P. T., Wynn G. A., 2014, MNRAS, 438, 240
  * Gotthelf et al. (2019) Gotthelf E. V., Halpern J. P., Mori K., Beloborodov A. M., 2019, ApJ, 882, 173
  * Green (2018) Green G., 2018, The Journal of Open Source Software, 3, 695
  * Green et al. (2019) Green G. M., Schlafly E., Zucker C., Speagle J. S., Finkbeiner D., 2019, ApJ, 887, 93
  * Harris et al. (2020) Harris C. R., et al., 2020, arXiv e-prints, p. arXiv:2006.10256
  * Heintz et al. (2020) Heintz K. E., et al., 2020, ApJ, 903, 152
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Hurley et al. (2005) Hurley K., et al., 2005, Nature, 434, 1098
  * Hurley et al. (2010) Hurley K., et al., 2010, MNRAS, 403, 342
  * Igoshev (2020) Igoshev A. P., 2020, MNRAS, 494, 3663
  * Igoshev & Popov (2018) Igoshev A. P., Popov S. B., 2018, MNRAS, 473, 3204
  * Igoshev et al. (2021) Igoshev A. P., Chruslinska M., Dorozsmai A., Toonen S., 2021, MNRAS, 508, 3345
  * Israel et al. (2003) Israel G. L., et al., 2003, ApJ, 589, L93
  * Israel et al. (2004) Israel G., et al., 2004, in Camilo F., Gaensler B. M., eds, IAU Symposium Conference Proceedings Vol. 218, Young Neutron Stars and Their Environments. p. 247 (arXiv:astro-ph/0310482)
  * Israel et al. (2005) Israel G., et al., 2005, A&A, 438, L1
  * James et al. (2020) James C. W., et al., 2020, ApJ, 895, L22
  * Jawor & Tauris (2021) Jawor J. A., Tauris T. M., 2021, MNRAS,
  * Kasen & Bildsten (2010) Kasen D., Bildsten L., 2010, ApJ, 717, 245
  * Kaspi & Beloborodov (2017) Kaspi V. M., Beloborodov A. M., 2017, ARA&A, 55, 261
  * Katz (2021) Katz J. I., 2021, arXiv e-prints, p. arXiv:2110.10847
  * King & Lasota (2019) King A., Lasota J.-P., 2019, MNRAS, 485, 3588
  * Kirsten et al. (2021) Kirsten F., et al., 2021, arXiv e-prints, p. arXiv:2105.11445
  * Kochanek (2009) Kochanek C. S., 2009, ApJ, 707, 1578
  * Kochanek (2021) Kochanek C. S., 2021, MNRAS, 507, 5832
  * Kochanek et al. (2019) Kochanek C. S., Auchettl K., Belczynski K., 2019, MNRAS, 485, 5394
  * Komissarov & Barkov (2007) Komissarov S. S., Barkov M. V., 2007, MNRAS, 382, 1029
  * Kothes & Dougherty (2007) Kothes R., Dougherty S. M., 2007, A&A, 468, 993
  * Kothes & Foster (2012) Kothes R., Foster T., 2012, ApJ, 746, L4
  * Kouveliotou et al. (1994) Kouveliotou C., et al., 1994, Nature, 368, 125
  * Kouveliotou et al. (1998) Kouveliotou C., et al., 1998, Nature, 393, 235
  * Kremer et al. (2021) Kremer K., Piro A. L., Li D., 2021, ApJ, 917, L11
  * Kroupa et al. (1993) Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545
  * Leahy & Tian (2008) Leahy D. A., Tian W. W., 2008, AJ, 135, 167
  * Levan et al. (2006) Levan A. J., Wynn G. A., Chapman R., Davies M. B., King A. R., Priddey R. S., Tanvir N. R., 2006, MNRAS, 368, L1
  * Levan et al. (2018) Levan A., Kouveliotou C., Fruchter A., 2018, ApJ, 854, 161
  * Levin et al. (2010) Levin L., et al., 2010, ApJ, 721, L33
  * Lin et al. (2011) Lin L., et al., 2011, ApJ, 739, 87
  * Lu et al. (2021) Lu W., Beniamini P., Kumar P., 2021, arXiv e-prints, p. arXiv:2107.04059
  * Lyman et al. (2022) Lyman J. D., Levan A. J., Wiersema K., Kouveliotou C., Chrimes A. A., Fruchter A. S., 2022, ApJ, 926, 121
  * Lyubarsky (2014) Lyubarsky Y., 2014, MNRAS, 442, L9
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993
  * Mannings et al. (2021) Mannings A. G., et al., 2021, ApJ, 917, 75
  * Margalit et al. (2019) Margalit B., Berger E., Metzger B. D., 2019, ApJ, 886, 110
  * Maund et al. (2016) Maund J. R., Pastorello A., Mattila S., Itagaki K., Boles T., 2016, ApJ, 833, 128
  * Metzger & Piro (2014) Metzger B. D., Piro A. L., 2014, MNRAS, 439, 3916
  * Metzger et al. (2011) Metzger B. D., Giannios D., Thompson T. A., Bucciantini N., Quataert E., 2011, MNRAS, 413, 2031
  * Metzger et al. (2015) Metzger B. D., Margalit B., Kasen D., Quataert E., 2015, MNRAS, 454, 3311
  * Metzger et al. (2019) Metzger B. D., Margalit B., Sironi L., 2019, MNRAS, 485, 4091
  * Mirabel et al. (2002) Mirabel I. F., Mignani R., Rodrigues I., Combi J. A., Rodríguez L. F., Guglielmetti F., 2002, A&A, 395, 595
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Mohan et al. (2020) Mohan P., An T., Yang J., 2020, ApJ, 888, L24
  * Mondal (2021) Mondal T., 2021, ApJ, 913, L12
  * Muñoz-Darias et al. (2016) Muñoz-Darias T., de Ugarte Postigo A., Casares J., 2016, MNRAS, 458, L114
  * Muno et al. (2006) Muno M. P., et al., 2006, ApJ, 636, L41
  * Nakano et al. (2015) Nakano T., Murakami H., Makishima K., Hiraga J. S., Uchiyama H., Kaneda H., Enoto T., 2015, PASJ, 67, 9
  * Oke & Gunn (1983) Oke J. B., Gunn J. E., 1983, ApJ, 266, 713
  * Olausen & Kaspi (2014) Olausen S. A., Kaspi V. M., 2014, ApJS, 212, 6
  * Palla & Stahler (1993) Palla F., Stahler S. W., 1993, ApJ, 418, 414
  * Palmer et al. (2005) Palmer D. M., et al., 2005, Nature, 434, 1107
  * Perez & Granger (2007) Perez F., Granger B. E., 2007, Computing in Science and Engineering, 9, 21
  * Popov (2022) Popov S. B., 2022, arXiv e-prints, p. arXiv:2201.07507
  * Predehl & Schmitt (1995) Predehl P., Schmitt J. H. M. M., 1995, A&A, 500, 459
  * Prentice et al. (2018) Prentice S. J., et al., 2018, ApJ, 865, L3
  * Price-Whelan et al. (2018) Price-Whelan A. M., et al., 2018, AJ, 156, 123
  * Renzo et al. (2019) Renzo M., et al., 2019, A&A, 624, A66
  * Rodrigo & Solano (2020) Rodrigo C., Solano E., 2020, in XIV.0 Scientific Meeting (virtual) of the Spanish Astronomical Society. p. 182
  * Rodrigo et al. (2012) Rodrigo C., Solano E., Bayo A., 2012, SVO Filter Profile Service Version 1.0, IVOA Working Draft 15 October 2012, doi:10.5479/ADS/bib/2012ivoa.rept.1015R
  * Rowlinson et al. (2013) Rowlinson A., O’Brien P. T., Metzger B. D., Tanvir N. R., Levan A. J., 2013, MNRAS, 430, 1061
  * Safarzadeh et al. (2020) Safarzadeh M., Prochaska J. X., Heintz K. E., Fong W.-f., 2020, ApJ, 905, L30
  * Sakamoto et al. (2011) Sakamoto T., et al., 2011, Advances in Space Research, 47, 1346
  * Salganik et al. (2021) Salganik A., Tsygankov S. S., Djupvik A. A., Karasev D. I., Lutovinov A. A., Buckley D. A. H., Gromadzki M., Poutanen J., 2021, MNRAS,
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Sana et al. (2014) Sana H., et al., 2014, ApJS, 215, 15
  * Sato et al. (2010) Sato T., Bamba A., Nakamura R., Ishida M., 2010, PASJ, 62, L33
  * Schlafly & Finkbeiner (2011) Schlafly E. F., Finkbeiner D. P., 2011, ApJ, 737, 103
  * Schneider et al. (2020) Schneider F. R. N., Ohlmann S. T., Podsiadlowski P., Röpke F. K., Balbus S. A., Pakmor R., 2020, MNRAS, 495, 2796
  * Scholz et al. (2012) Scholz P., Ng C. Y., Livingstone M. A., Kaspi V. M., Cumming A., Archibald R. F., 2012, ApJ, 761, 66
  * Seabold & Perktold (2010) Seabold S., Perktold J., 2010, in 9th Python in Science Conference.
  * Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
  * Smith et al. (2018) Smith L. C., et al., 2018, MNRAS, 474, 1826
  * Stanway & Eldridge (2018) Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75
  * Stanway et al. (2020) Stanway E. R., Chrimes A. A., Eldridge J. J., Stevance H. F., 2020, MNRAS, 495, 4605
  * Stevance et al. (2020) Stevance H., Eldridge J., Stanway E., 2020, The Journal of Open Source Software, 5, 1987
  * Tam et al. (2008) Tam C. R., Gavriil F. P., Dib R., Kaspi V. M., Woods P. M., Bassa C., 2008, ApJ, 677, 503
  * Tauris & Takens (1998) Tauris T. M., Takens R. J., 1998, A&A, 330, 1047
  * Tendulkar et al. (2012) Tendulkar S. P., Cameron P. B., Kulkarni S. R., 2012, ApJ, 761, 76
  * Tendulkar et al. (2013) Tendulkar S. P., Cameron P. B., Kulkarni S. R., 2013, ApJ, 772, 31
  * Tendulkar et al. (2021) Tendulkar S. P., et al., 2021, ApJ, 908, L12
  * Testa et al. (2008) Testa V., et al., 2008, A&A, 482, 607
  * Testa et al. (2018) Testa V., Mignani R. P., Hummel W., Rea N., Israel G. L., 2018, MNRAS, 473, 3180
  * Thompson & Duncan (1996) Thompson C., Duncan R. C., 1996, ApJ, 473, 322
  * Tian & Leahy (2008) Tian W. W., Leahy D. A., 2008, ApJ, 677, 292
  * Tian & Leahy (2012) Tian W. W., Leahy D. A., 2012, MNRAS, 421, 2593
  * Tiengo et al. (2010) Tiengo A., et al., 2010, ApJ, 710, 227
  * Torii et al. (1998) Torii K., Kinugasa K., Katayama K., Tsunemi H., Yamauchi S., 1998, ApJ, 503, 843
  * Van der Horst et al. (2010) Van der Horst A. J., et al., 2010, ApJ, 711, L1
  * Verbunt et al. (2017) Verbunt F., Igoshev A., Cator E., 2017, A&A, 608, A57
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wang et al. (2006) Wang Z., Chakrabarty D., Kaplan D. L., 2006, Nature, 440, 772
  * Waskom (2021) Waskom M. L., 2021, Journal of Open Source Software, 6, 3021
  * Wheeler et al. (2000) Wheeler J. C., Yi I., Höflich P., Wang L., 2000, ApJ, 537, 810
  * White et al. (2021) White C. J., Burrows A., Coleman M. S. B., Vartanyan D., 2021, arXiv e-prints, p. arXiv:2111.01814
  * Woosley (2010) Woosley S. E., 2010, ApJ, 719, L204
  * Xu et al. (2021) Xu K., Li X.-D., Cui Z., Li Q.-C., Shao Y., Liang X., Liu J., 2021, arXiv e-prints, p. arXiv:2110.10438
  * Younes et al. (2017) Younes G., Baring M. G., Kouveliotou C., Harding A., Donovan S., Göğüş E., Kaspi V., Granot J., 2017, ApJ, 851, 17
  * Younes et al. (2021) Younes G., et al., 2021, Nature Astronomy, 5, 408
  * Zapartas et al. (2017) Zapartas E., et al., 2017, ApJ, 842, 125
  * Zhang & Zhang (2021) Zhang R. C., Zhang B., 2021, arXiv e-prints, p. arXiv:2109.07558
  * Zhou et al. (2014) Zhou P., Chen Y., Li X.-D., Safi-Harb S., Mendez M., Terada Y., Sun W., Ge M.-Y., 2014, ApJ, 781, L16
  * Zhou et al. (2019) Zhou P., Vink J., Safi-Harb S., Miceli M., 2019, A&A, 629, A51
  * Zou et al. (2021) Zou L., et al., 2021, MNRAS, 508, 2505

##  Appendix A Colour-magnitude diagrams at other ages

In this appendix, we show colour-magnitude diagrams at two other ages,
$10^{6}$ yr and $10^{8}$ yr, covering the majority of the age range where
core-collapse supernovae are expected (Figure 6). Along with Figure 5, we
therefore show the range of stellar population colours and magnitudes at the
expected age of bound companions to recently born core-collapse neutron stars,
agnostic to any prior assumptions about the expected properties of bound
companions.

Figure 6: Left: a BPASS NIR colour-magnitude diagram made with HOKI, as in
Figure 5, at $10^{6}$ yr rather than $10^{7}$ yr. Right: the same plot at
$10^{8}$ yr. The age range of 1 to 100 Myr represents the expected lifetimes
of supernova progenitors; the stellar populations shown therefore represent
the possible $J$ and $H$ photometric properties of all stellar objects at the
time of primary star supernovae, without restricting ourselves specifically to
the predicted properties of bound companions.