Searching for Magnetar Binaries Disrupted by Core-Collapse Supernovae

Core-collapse Supernovae (CCSNe) are considered the primary magnetar formation channel, with 15 magnetars associated with supernova remnants (SNRs). A large fraction of these should occur in massive stellar binaries that are disrupted by the explosion, meaning that $\sim45\%$ of magnetars should be nearby high-velocity stars. Here we conduct a multi-wavelength search for unbound stars, magnetar binaries, and SNR shells using public optical ($uvgrizy-$bands), infrared ($J-$, $H-$, $K-$, and $K_s-$bands), and radio ($888$ MHz, $1.4$ GHz, and $3$ GHz) catalogs. We use Monte Carlo analyses of candidates to estimate the probability of association with a given magnetar based on their proximity, distance, proper motion, and magnitude. In addition to recovering a proposed magnetar binary, a proposed unbound binary, and 13 of 15 magnetar SNRs, we identify two new candidate unbound systems: an OB star from the Gaia catalog we associate with SGR J1822.3-1606, and an X-ray pulsar we associate with 3XMM J185246.6+003317. Using a Markov-Chain Monte Carlo simulation that assumes all magnetars descend from CCSNe, we constrain the fraction of magnetars with unbound companions to $5\lesssim f_u \lesssim 24\%$, which disagrees with neutron star population synthesis results. Alternate formation channels are unlikely to wholly account for the lack of unbound binaries as this would require $31\lesssim f_{nc} \lesssim 66\%$ of magnetars to descend from such channels. Our results support a high fraction ($48\lesssim f_m \lesssim 86\%$) of pre-CCSN mergers, which can amplify fossil magnetic fields to preferentially form magnetars.

It is unclear whether magnetars constitute a unique class of neutron star or instead that neutron stars exhibit a continuum of magnetar-like features depending on their magnetic fields and environments.The latter case would be supported by the magnetar-like activity (glitches, X-ray outbursts, and high magnetic field strengths) observed in some pulsars (e.g.Gavriil et al. 2008;Rea et al. 2016;Archibald et al. 2016; ★ E-mail: msherman@caltech.edu 1 While a number of Gamma Ray Bursts (GRBs) in nearby galaxies have been proposed as magnetar Giant Flares (e.g.Trigg et al. 2023;Burns et al. 2021;Svinkin et al. 2021;Mazets et al. 2008;Frederiks et al. 2007;Minaev et al. 2024;Mereghetti et al. 2024), we focus this work on magnetars and candidates within the Galaxy and Magellanic Clouds.Beniamini et al. 2023).In this work, we will consider magnetars to be a class of highly-magnetized neutron star distinguished by their magnetic field decay-powered emission.
Magnetars, like pulsars, are believed to be formed when massive stars undergo Core-collapse Supernovae (CCSNe, formed from Type-II, Type-Ib, and Type-Ic supernovae), with 15 magnetars or candidates having confirmed or potential associations with Supernova Remnants (SNRs) (e.g., Baade & Zwicky 1934;Wheeler 1966;Olausen & Kaspi 2014;Nomoto et al. 1995;Thielemann et al. 1996;Barker et al. 2023;Truelove & McKee 1999).However, what factors determine whether a CCSN leaves behind a magnetar or a normal neutron star remain unclear (e.g., Giacomazzo & Perna 2013;Giacomazzo et al. 2015;Fuller & Lu 2022;Revnivtsev & Mereghetti 2016;Popov 2020).Furthermore, the CCSN model may be inconsistent with trends in the growing neutron star population.Keane & Kramer (2008) discuss the birth rates of various species of galactic neutron star, including radio and X-ray pulsars, Rotating Radio Transients (RRATs), X-ray Dim Isolated Neutron Stars (XDINs), and magnetars.They find that while the total neutron star birth rate is between 5 − 10 century −1 , the CCSN rate is only 1.9 ± 1.1 century −1 , suggesting either our understanding of evolutionary links among neutron star species, our CCSN model, or our birth rate estimate is incomplete (e.g.Popov 2010;Keane & Kramer 2008;Kaspi 2010).CCSNe are also unlikely to account for the population of young pulsars in Galactic globular clusters (e.g.Kremer et al. 2023).
Magnetar birth rates  ≈  mag T −1 ∼ 0.3 century −1 , are dif-ficult to constrain due to less confident age T and population size  mag estimates.Most notably, the characteristic spindown age T rot ∝ (2 ) −1 ∼ 10 kyr−10 Myr used for radio pulsars (where  is the pulse period and  is the spindown rate) does not incorporate the magnetic field decay on ∼ 10 kyr timescales thought to drive magnetar emission2 .Therefore, T rot can significantly overestimate the true age (e.g.Nakano et al. 2015;Ferrario & Wickramasinghe 2008;Pons & Geppert 2007).While this is the major effect, 'glitches', or abrupt decreases in rotation period  (e.g.Tong & Huang 2020;Younes et al. 2023;Livingstone et al. 2011b), can also complicate magnetar age estimation when data is limited (e.g. for transient magnetars as discussed in Rea et al. 2016;Ibrahim et al. 2024;Lower et al. 2023;Scholz et al. 2017), as can 'anti-glitches' and 'glitch recovery' for which  increases (e.g.Scholz et al. 2014;Archibald et al. 2013;Olausen & Kaspi 2014).For these reasons, magnetic field decay ages T decay ∼ 1 − 100 kyr derived from X-ray luminosity and magnetic field strength may be preferred, although these are model-dependent, and may not account for rapid initial spindown models (see Appendix D1 in this work, and e.g.Ferrario & Wickramasinghe 2008;Mondal 2021;Arras et al. 2004, for examples of decay models).The ages of associated SNRs T SNR ∼ 1 kyr−100 kyr are more accurate, but can cover a wide range and are only known for a small fraction of magnetars (Suzuki et al. 2021;Borkowski & Reynolds 2017;Lyman et al. 2022;Gaensler et al. 1999a).
While Muno et al. (2008) estimate that  mag ≲ 540 magnetars reside in the Galaxy, a number of factors make this uncertain.First, it is unclear how magnetar beaming affects the observed population, particularly due to the small sample size (e.g.Keane & Kramer 2008).Furthermore, the detection of a globular cluster (M81) Fast Radio Burst (FRB), thought to originate in a magnetar's magnetosphere, offers a counterexample to our assumption that  mag scales with star formation, and supports the notion of multiple magnetar formation channels (Kirsten et al. 2022;Lu et al. 2022;Petroff et al. 2022).The evolutionary connections among magnetars, pulsars, XDINs, RRATs, and isolated neutron stars are also unclear, leading to potential double-counting errors in the Keane & Kramer (2008) birthrates.Thus, the continued search for and recovery of CCSNe and SNRs will help inform our models and improve age estimates.
In recent years, the presence (or absence) of binary magnetar systems has become a powerful probe of CCSN formation channels for the broader neutron star population.Chrimes et al. (2022b), Renzo et al. (2019), Kochanek et al. (2019), and Moe & Di Stefano (2017) conducted population synthesis simulations with massive stellar binaries to estimate what fraction of primary stars would remain bound to their secondaries as neutron stars following their CCSN.Moe & Di Stefano (2017) estimate that  0 ∼ 84% of massive stars should initially be in binaries, while Kochanek et al. (2019) predicts that   ∼ 48% of these merge prior to core-collapse3 ; the other 1 −   ∼ 52% undergo CCSNe as non-merged binaries.Renzo et al. (2019) suggest a lower merger rate   ∼ 22 +26 −9 %.They also find that ∼ 14 +22 −10 % of the binaries that do not merge remain bound following the CCSN.Thus only   ∼ 6 − 9% of all neutron stars should have massive binary companions.Chrimes et al. (2022b) affirm these results, obtaining a bound neutron star fraction   ∼ 5%.Notably, the binary statistics of magnetars may deviate from those above if they form from distinct progenitor populations, such as those with enhanced rotation or magnetic fields, though this has not been explored extensively through population synthesis (e.g.Popov & Prokhorov 2006)  4 .
Taking a census of young, non-degenerate neutron star-massive star binaries, 7 pulsars with bound OB star companions have been reported in the Australia Telescope National Facility 5 (ATNF) catalog, while SGR-0755-2933 is a tentatively confirmed magnetar binary (Chrimes et al. 2022b;Richardson et al. 2023;Doroshenko et al. 2021).A candidate stellar companion was also proposed for magnetar CXOU J171405.7-381031 using Infrared (IR) Hubble Space Telescope (HST) images, though the IR counterpart's stellar nature is unconfirmed (Chrimes et al. 2022b).The two proposed binaries out of 31 magnetars are marginally consistent with the   ∼ 5 − 9% expected to be in bound systems, while the 9 total neutron stars with massive companions out of the total catalog of ∼ 2950 neutron stars marginally underestimate the binary fraction due to selection effects (Lorimer 2008)  6 . 3Note this fraction includes only mergers where both stars are massive OB stars.For a discussion of CCSNe from mergers of intermediate mass stars, see Zapartas et al. (2017). 4We caution that the results of Popov & Prokhorov (2006) may be biased by their definition of 'magnetars', which they define as evolutionary paths that include "Coalescence prior to NS formation" or "Roche lobe overflow", rather than using an independent condition on their magnetic field strength.Therefore, the merger fractions and binary fractions from Tables 1 and 2 of Popov & Prokhorov (2006) should not be directly compared to those from Chrimes et al. (2022b), Renzo et al. (2019), Kochanek et al. (2019), Moe & Di Stefano (2017), or this work. 5https://www.atnf.csiro.au/research/pulsar/psrcat/ 6Our young neutron star census includes 2903 slow (period  > 30 ms) pulsars from the ATNF catalog, ∼ 10 thermally discovered isolated neutron stars including the Magnificent Seven ROSAT sample (Pires et al. 2009;Schwope et al. 1998;Haberl 2007;Turolla 2009;Tetzlaff et al. 2011;Avakyan et al. 2023), and 31 magnetars and magnetar candidates from the McGill magnetar catalog.We acknowledge this may be incomplete; for example, not all pulsars have sufficient timing solutions to rule out binary companions, and some may be old enough that their companions have undergone supernovae, as well.However, this provides a reasonable estimate of the number of known young neutron stars to date.We exclude older populations such as Highand Low-Mass X-ray Binaries (HMXBs, LMXBs) and millisecond pulsars However, one expects a much larger fraction (   ∼ 38 − 56%) of 'unbound' neutron star binaries disrupted by the CCSN, motivating a search for 'walkaway' or 'runaway' stellar companions (Renzo et al. 2019).Considering a physical scenario, a massive companion ejected at ∼ 30 km s −1 (Renzo et al. 2019) could travel ∼ 0.03 pc over ∼ 1 kyr.At a distance of 5−10 kpc, this should appear only ∼ 1.2−2.5 ′′ away from its previously bound magnetar.Kochanek et al. (2019) and Kochanek (2018Kochanek ( , 2021Kochanek ( , 2023) have used the Gaia optical catalog to search for both runaway and bound stars nearby 43 SNRs with compact objects, including central pulsars, HMXBs, LMXBs, radio sources, and gamma ray sources.Only the pulsar PSR J0538+2817 was confidently associated with a runaway companion star (Dinçel et al. 2015;Kochanek 2021).Clark et al. (2014) note that the Wolf-Rayet star WR 77F (Wd1-5) may be a runaway unbound companion of magnetar CXOU 164710.2-455216,but the absence of a magnetar proper motion measurement make this difficult to confirm.More general searches for runaway stars with high velocities relative to field stars have also been carried out, but haven't identified magnetar or pulsar companions (e.g.Carretero-Castrillo et al. 2023).With ∼ 2950 neutron stars and 31 magnetars, one expects ∼ 1100 − 1700 and ∼ 11 − 18 unbound companions, respectively 7 .
A comprehensive search of archived data can leverage catalogs like Gaia, the Two Micron All Sky Survey (2MASS), and Panoramic Survey Telescope and Rapid Response System (Pan-STARRS or PS1), as well as surveys such as the SkyMapper Southern Survey (SMSS), UKIRT Infrared Deep Sky Survey (UKIDSS) and the Visible and Infrared Survey Telescope for Astronomy (VISTA) Variables in the Vía Láctea (VVV) to find stellar companions (Flewelling et al. 2020;Skrutskie et al. 2006;Chambers et al. 2019;Brown et al. 2021;Lawrence et al. 2007a;Kaiser et al. 2002;Nikzat et al. 2022;Wolf et al. 2018;Onken et al. 2024).For example, Chrimes et al. (2023) predict the ability of Gaia Data Release 3 (DR3), as well as future releases from HST, the James Webb Space Telescope (JWST), Euclid, and the Nancy Grace Roman Space Telescope (NGRST) to detect unbound stellar objects.They find that while Gaia is limited by extinction (  ≲ 10), its absolute reference capability and proper motion measurements could increase the number of unbound companion discoveries by ∼ 8× when combined with next generation IR telescope data.Gaia data have also been used effectively to search for unbound companions of White Dwarfs and neutron stars (Kochanek et al. 2019;El-Badry et al. 2023;Shen et al. 2018).
In this paper, we conduct a search for unbound binary companions of magnetars using the Gaia catalog, supplemented by a search for bound companions and CCSN remnants.In Section 2 we describe the magnetar sample used for the search.Section 3 describes the search and trajectory analysis of sources in the Gaia archive for unbound magnetar companions and results.In Section 4 we describe the methods and results of the optical and IR search for bound magnetar companions.In Section 5 we summarize the radio image search using the RACS, VLASS, and NVSS catalogs.Section 6 discusses the sensitivity, completeness, and implications of our search.Section 7 concludes with a summary and outlook on next steps for analyzing magnetar formation.Section 8 provides information regarding data availability.Details on the statistical analysis are included in the Appendices.
since they are only observable in binary systems, and are therefore subject to selection effects (e.g.Lorimer 2008). 7See Table 4 for definitions of   ,   ,   , and other relevant fractions.These will be discussed in detail in Section 6.2.

MAGNETAR SAMPLE
In Table 1, we present the sample of 24 magnetars listed in the McGill magnetar catalog8 , as well as 7 unconfirmed magnetar candidates (Olausen & Kaspi 2014).For each we list the Right Ascension (), declination (), distance (), and proper motion () from the available literature.We also provide the -band dust extinction estimates along each magnetar's line-of-sight, which we derive using the Bayestar199 dust map for magnetars above declination  ≳ −30 • (Green et al. 2015;Schlafly et al. 2015;Green et al. 2018;Green et al. 2019).For magnetars outside this range we use the Predehl & Schmitt (1995) relation between the X-ray column density,   , and   if an   measurement is available.If not, the average photometric extinction from Gaia Data Release 3 (DR3) sources within 1 • is given (Gaia Collaboration et al. 2016;Collaboration et al. 2022).This will likely underestimate the true extinction, since most Gaia sources are in the foreground.12 magnetars and 3 magnetar candidates have secure or proposed SNR associations.

Gaia Data
The Gaia DR3 catalog was queried to search for unbound magnetar companions ejected by the CCSN (Gaia Collaboration et al. 2016;Collaboration et al. 2022).As is highlighted in Chrimes et al. (2023), the Gaia catalog is a unique tool given its absolute reference frame which allows for comparison of proper motion data.We build upon previous magnitude-based searches (e.g.Chrimes et al. 2022b) by analyzing the trajectories for each GDR3 source in reference to the magnetars'.Chrimes et al. (2023) note the Gaia survey depth  ,min = 20.7 and proper motion uncertainty ranging from 9.6 ≤  min ≤ 801.6 as yr −1 for 13 ≤   ≤ 20.7 limit the survey's detection of walkaway or runaway stars to those with low extinction (  ≲ 10; Lindegren et al. 2021) 10 .Despite this, our combination of magnitude analysis with proper motion data will still make a marked improvement over previous searches, seen for example in El-Badry et al. (2023) and inferred through simulations by Chrimes et al. (2023, Figure 8 row 1).
The Gaia query is limited to potential O and B type stars, as we expect many of the companions of magnetars to be massive stars, and removing dim stars significantly reduces contamination by false positives.To do this we impose an upper limit on the absolute -band magnitude   < −2.5log 10 (  / ⊙ )+4.68 = 1.185,where we use   ≈ 25 ⊙ , and color   −   < 2. Appendix A contains the Astronomical Data Query Language (ADQL) query used to compile the initial Gaia sample.
Following the initial query, extinction corrections and cutoffs were re-applied using the Bayestar1911 dust map (Green et al. 2018;Green et al. 2019) for any Gaia sources above declination  ≳ −30 • .From Bayestar19 we obtain the  −  color excess12 ,  − , which is    5.9 for details.
−band magnitudes,   and   .Details on conversions between photometric systems are in Appendix A.
For sources with  ≲ −30 • ,   is computed with Equations A5-A8 using   and    −  reported in Gaia, which are derived from   −   spectra in the GSP-Aeneas library (Vallenari et al. 2022).However, we proceed cautiously with the GSP extinction estimates, which are known to be unreliable within the Galactic Plane where the effective temperature is overestimated (see Andrae et al. 2022;Fouesneau et al. 2022, Section 3.4, and the Gaia DR3 Documentation Section 11.3.3 for details).
We use a Monte Carlo simulation to estimate p-values for the Gaia search.The criteria for association for a given sampled trajectory are below: (i) When the magnetar and source are extrapolated back in time over a sampled travel time (between 1 − 10 kyr), their initial angular separation Θ ≈ 0 • .
(ii) The absolute −band magnitude of the Gaia source is consistent with that of an OB (possibly A) star (−10 ≲   ≲ 3)14 .
(iii) The difference in 2D travel time for the magnetar and source to travel from their past intersection point to their present positions is Δ 2D ≈ 0.
(iv) The magnetar's physical travel time T 3D,mag is consistent with its age T .
(v) The source's parallax distance is consistent with magnetar's distance, such that Δ ≈ 0. which are used to define p-values  1 ,  2 ,  3 , and  4 respectively (the first two criteria are combined to define  1 ).Note that the Gaia radial velocity measurements are not used in our search since they have not been confirmed through spectroscopic follow-up and would overconstrain the sampled velocity vectors.Appendix B4 describes the definition of  1 ,  2 ,  3 , and  4 in more detail, which are combined as shown below: We choose this definition (as opposed to e.g. the Fisher statistic, Lancaster statistic, or a direct average) for its intuitive application to this search, as all four criteria must be met with   << 5% to obtain  < 5%.This also makes it straightforward to interpret p-values by tracing back each   to its criterion.Each magnetar's age is estimated as the age of its associated SNR T SNR (Suzuki et al. 2021).If no SNR is associated, its X-ray luminosity decay age T decay is used, which adopts the Ferrario & Wickramasinghe (2008) field model and is derived in Appendix D1.For magnetars with no persistent X-ray emission or magnetic field estimates, the characteristic spin-down age T rot is used (Olausen & Kaspi 2014).In Table 2 we provide the derived or reported timescales for each magnetar.

Gaia Search Test Using High Velocity Pulsars with Supernova Remnants
In order to test the search pipeline, we applied the Monte Carlo simulation to a sample of 15 pulsars with known SNR associations and proper motion measurements compiled from the ATNF catalog (van der Wateren et al. 2023;Andersen et al. 2023;Lyne et al. 2015;Stairs et al. 2001;Dinçel et al. 2015;Kochanek 2021).Within this sample, only PSR J0538+2817 is known to have a runaway companion B0.5V-type star, as discovered by Dinçel et al. (2015) and confirmed by Kochanek (2021).The upper limit on the travel time  is set to 100 kyr, the maximum age of SNR shells before they dissipate (e.g.Frail et al. 1994), given that pulsars can have ages > 100 kyr (e.g.Lorimer 2008;Kaspi & Beloborodov 2017).The initial Monte Carlo simulation (criteria 1 and 2) identifies between 2 and 800 potential candidates with  1 < 5% for each of the 15 pulsars.Following the trajectory analysis, only 1 pulsar, PSR B1951+32 remains with  < 5% and a sufficient fraction of valid trajectories15 .The companion of PSR J0538+2817 is not initially recovered.In this section we discuss the likelihood of PSR B1951+32's candidate companion, and the reasons for PSR J0538+2817's omission.Figure 2 shows the trajectory plots for both pulsars.

PSR J0538+2817
For PSR J0538+2817, we initially do not recover the runaway companion star, HD 37424, which was identified by Dinçel et al. (2015) and Kochanek (2021).The corresponding Gaia source, Gaia DR3 3441732292729818752, was assigned a p-value  > 99% 16 and therefore was not identified as a candidate.Upon fur-regions in {Θ,   } space than magnetars due to the higher age limit and lower distances, the high proper motions of PSR J0538+2817, PSR J1809-2332, and PSR B2334+61 further extend the region, resulting in unusually high false positive rates ≳ 99%, ≳ 91%, and ≳ 76%, respectively.We therefore increase the p-value thresholds for these three sources to match their minimum false positive rates.We use a threshold  1 < 100% for PSR J0538+2817 since the false positive rate is ≳ 99%.For PSR B2334+61 and PSR J1809-2332, no additional sources are identified with the new threshold, whereas HD 37424 is identified as a candidate for PSR J0538+2817 after increasing the threshold.This scenario also justifies our use of the search statistic f < f to identify candidates initially, as it implicitly accounts for this falsepositive rate bias.Our conditional recovery of PSR J0538+2817's runaway star demonstrates the accuracy and sensitivity of our search method, but indicates we must be wary of false positives and sources along unique sightlines.We address completeness of the magnetar search in detail in Section 6.

PSR B1951+32
One Gaia companion candidate is identified for PSR B1951+32 with  < 0.1%, Gaia DR3 2035109887486067456 (catalogued as TYC 2673-414-1; Høg et al. 2000a).The parallax distance 1.13 ± 0.02 kpc is only marginally consistent with the pulsar's ∼ 3 kpc distance, which has an uncertainty that is not well defined (Kulkarni et al. 1988).The pulsar's dispersion measure (DM) and more recent estimates of the associated SNR CTB 80 suggest a closer ∼ 2 kpc distance, which would be more consistent with TYC 2673-414-1 (Kramer et al. 2000;Migliazzo et al. 2002;Yao et al. 2017).From the 2MASS counterpart, 2MASS J19544429+3327472, we estimate absolute −band magnitude   ≈ −4.5 and  −  color  −  ≈ −0.1 when placed at the pulsar distance 3 kpc.While these are reasonable, they do fall outside of the expected ranges of IR magnitude and color from Chrimes et al. (2022b)'s simulations.From the valid trajectory samples, TYC 2673-414-1 is required to have velocity | − →  |≈ 75 +14 −6 km s −1 , which would be consistent with a more rare runaway star if ejected from the pulsar's CCSN (Renzo et al. 2019).While pulsar companions are not the focus of this search, we encourage follow-up of TYC 2673-414-1 as a possible unbound companion of PSR B1951+32.

Gaia Search Results
Through the Gaia Monte Carlo search, we identify potential candidates with  < 5% for 2 magnetars.In Figure 2, we show the trajectory plots and contour plots of the Monte Carlo samples for these, along with CXOU J164710.2-455216, whose proposed companion WR 77F we recover with  = 0.6% (Clark et al. 2014).We rule out one of the two candidates, 3XMM J185246.6+003317, as a false positive based on detailed analysis of its candidates' velocities and distances, which we discuss in Appendix C1.Having not recovered the proposed binaries SGR 0755-2933 and CXOU J171405.7-38131,we find this search method is relatively insensitive to bound stellar companions due to their position within magnetars' error ellipses (Chrimes et al. 2022b;Richardson et al. 2023;Doroshenko et al. 2021).However, it is a useful tool to identify which candidates are likely to be bound companions (see Section 4.3).In the following sections, we focus on the unbound companion candidates for SGR J1822.3-1606 and CXOU J164710.2-455216.

SGR J1822.3-1606
The primary candidate for an unbound system is SGR J1822.3-1606 and Gaia DR3 4097832458955594880, which is assigned a p-value  < 0.1%.Its distance,  = 1.00 +0.03  −0.02 kpc, is consistent with that of the magnetar, 1.6 ± 0.3 kpc, at the 2 level, and from the IR counterpart, 2MASS 18230515-1600296, we estimate absolute magnitudes   = −3.9± 0.5 and   = −4.6 +0.5 −0.4 , and color  −  = 0.7 ± 0.117 , which are both reasonable for a low-mass B-type star.SGR J1822.3-1606does not have a measured proper motion, and our simulation suggests   cos() = −29 +15 −24 mas yr −1 ,   = −25 +13 −19 mas yr −1 are required for the magnetar to have a valid trajectory 18 .A wide range of magnetar velocities are permitted for this association, with | − →  mag |= 413 +184 −175 km s −1 .Notably, the permitted stellar velocity | − →  |= 60 +7 −5 km s −1 would classify this as a rare 'runaway' companion star (Renzo et al. 2019).This is roughly consistent with the source's high transverse velocity, ∼ 55 km s −1 , making the association achievable without a significant radial component.While it is unlikely that such velocities can be reached from CCSN disruption alone, other mechanisms may contribute to the source's and magnetar's space motion, such as dynamical ejection from an OB association (e.g.Poveda et al. 1967;Tetzlaff et al. 2011).Therefore, we identify this as a promising candidate, and further characterization from follow-up observations would be of great interest.

CXOU J164710.2-455216
Next, we consider the Wolf-Rayet star WR 77F (or Wd1-5), which Clark et al. (2014) suggest as an unbound companion of CXOU 164710.2-455216based on the position, radial velocity and chemical content of the star.WR 77F is associated with Gaia DR3 5940106797374726784 which was not included in the initial Gaia query, possibly due to an underestimated Gaia extinction from GSP photometry.We add this candidate to the sample and run the simulation with an increased 1000 trials after finding a low fraction of valid trajectories with 100 samples.Given that the candidate is an evolved Wolf-Rayet star rather than a main-sequence OB star, and that magnetar CXOU J164710.2-455216lies along a high extinction sightline with   = 13.35 ± 0.28, we omit  1 from this analysis (which initially was  1 > 99.9%).With this adjustment, we recover a p-value  2,3,4 = 0.6%.The magnetar's derived physical travel time T 3D,mag ≈ 16.0 +25.6  −10.0 kyr, is consistent with its decay age T decay ≈ 0.2 Myr within ∼ 1.However, the localization of CXOU J164710.2-455216 and WR 77F to the Westerlund-1 massive star cluster implies that dynamical ejection may have contributed to the large observed velocities, as well.We conclude it is likely that CXOU J164710.2-455216 are WR 77F are an unbound magnetar binary disrupted by the former's CCSN, po-tentially supplemented by dynamical ejection, supporting the results of Clark et al. (2014).

OPTICAL AND IR SEARCH FOR BOUND COMPANIONS
In this section, we search for bound OB star-magnetar binaries using optical and near-IR images.As discussed in Section 1, bound companions should have negligible (≲ 1 ′ ) separation from the magnetar in high resolution (∼ 0.25 ′′ − 0.4 ′′ ) images.Chrimes et al. (2022a) conducted a similar search using HST Wide Field Camera 3 19 images with ∼ 27 magnitude depth and 0.1265 ′′ resolution.They place tight constraints with milliarcsecond precision on ∼ 15 magnetar NIR counterpart candidates using −band HST data, laying the foundation for the Chrimes et al. (2022b) search.
While Chrimes et al. (2022b) discuss that searching in NIR is more effective than optical due to the ∼ 7× lower extinction, searching optical data is still beneficial.The massive, un-evolved OB star companions we expect are intrinsically fainter in the NIR than optical wavelengths.The magnitude depth of the PS1 (∼ 24) and SkyMapper (∼ 18) catalogs therefore makes them powerful resources for this search (Maund et al. 2016;Kaiser 2004;Wolf et al. 2018;Onken et al. 2024).We also note that since the resolution (0.4 ′′ , 0.34 ′′ , and 3 ′′ , respectively) and depth (∼ 18, ∼ 15, and ∼ 19.3, respectively) of the UKIDSS, VVV, and 2MASS images are far exceeded by HST, any additional data to supplement the breadth of the search is valuable.
In addition to image cutouts, the PS1 and SkyMapper point source catalogs were queried to obtain known sources within a 30 ′′ search cone of each magnetar.We also include any bright sources identified 'by-eye' within the image that do not appear to coincide with catalogued point sources (see Appendix B3.1).We then use the catalogued magnitudes and position of each source to conduct a Monte Carlo simulation.Similarly to the Gaia search we define the association p-value using three criteria: (i) The PS1 source and magnetar are located at the same RA and declination.
The method above was applied only for magnetars with errors defined for the RA, declination and distance, and only for point sources with at least − and − band PSF and Kron magnitudes22 , and errors defined for RA and declination.Note that our evaluation of  2 contrasts with the Chrimes et al. (2022b) search, which uses the magnetar distance to estimate extinction values for nearby candidate stars.However, magnetar distances are not well-constrained, with errors of order ∼ kpc.In this work, we invert this process to obtain a distance estimate whose uncertainty is driven by the expected absolute magnitude range.Through this, we avoid making assumptions about the source distance beyond the requirement that it be an OB star.

UKIDSS, VVV, and 2MASS Data and IR Search Method
Infrared images were obtained from the UKIDSS survey catalog, which is a superset of the Large Area Survey (LAS), Galactic Plane Survey (GPS), Galactic Clusters Survey (GCS), Deep Extragalactic Survey (DES), and Ultra Deep Survey (UDS) catalogs (Dye et al. 2006).Collectively, this covers the declination range from −60 • <  < +40 • , including sightlines for 12 magnetars and 6 magnetar candidates from this sample.We obtained 1 ′ × 1 ′ image cutouts with 0.2 ′′ pixels in the −, −, and −bands (390 − 1220 nm) for each magnetar to search for binary companions.For 8 magnetars in the southern portion of the Galactic disk, the VISTA Variables in the Via Lactea (VVV) survey was queried for −, −, and   −band images with 0.2 ′′ pixels (Nikzat et al. 2022).Images for the five remaining magnetars were obtained from the Two-Micron All-Sky Survey (2MASS)23 , which surveys the full sky in the −, −, and   −bands at much lower ∼ 3 ′′ resolution (∼ 2 ′′ pixels; Skrutskie et al. 2006).The UKIDSS Galactic Plane Survey (UGPS)24 catalog, 2MASS point source catalog, and VVV Infrared Astrometric Catalogue (VIRAC) 25 were queried in a cone search for discrete objects within 0.5 ′ of each magnetar (Skrutskie et al. 2006;Nikzat et al. 2022;Lucas et al. 2008;Lawrence et al. 2012Lawrence et al. , 2013;;Smith et al. 2018).Additional sources were identified 'by-eye' in each image to include in the search.
The criteria used for the IR search are below:  (i) The source and magnetar are located at the same RA and declination.
(ii) If we assume the source's absolute −band magnitude and  −  color are stellar in nature and resemble that of a bound OB star (−2.76 ≲   ≲ 3.34, −0.14 ≲  −  ≲ 0.31 are 95% confidence ranges from Chrimes et al. (2022b) population synthesis), the distance is equal to the magnetar's distance.
for which we define p-values  1 and  2 , respectively.Kron magnitudes are not available from the included IR catalogs.The p-values are combined using: An explicit derivation of these p-values is provided in Appendix B2.Note the p-value is only evaluated for magnetars and sources with position errors, and only for 2MASS and VVV sources with at least − and − band magnitudes.Chrimes et al. (2022b) revealed that magnetar surface IR emission can be easily mistaken for a companion star, with IR magnitudes   ≳ 2.5 (see Chrimes et al. 2022b, Figure 5).Although interacting binaries such as SGR 0755-2933 can be far brighter (its companion CPD-29 2076 has   ∼ −4), passive binaries like PSR J0210+5845 are much fainter in the absence of mass transfer (its companion B0V star has   ∼ 0.5; Tan et al. 2020).Both magnitude and color must be considered to distinguish massive stellar companions from surface emission, which we account for in  2 .

Optical and IR Search Results
Using a 95% confidence level, we recover the proposed Be star bound companion of magnetar candidate SGR 0755-2933, CPD-29 2176, through the optical (PS1) and IR (UKIDSS) searches with  < 0.1%26 (Chrimes et al. 2022b;Richardson et al. 2023).In the following sections, we discuss the association of a PS1 source with the companion of SGR 0755-2933's error circle and the non-recovery of the proposed counterpart of CXOU J171405.7-38131.Figure 3 shows optical and IR images for these three magnetars.Candidate companions were identified with  < 5% for PSR J1622-4950 and SGR J1822.3-1606, but are most likely false positives as discussed in Appendix C3.

PS1 Counterpart of SGR 0755-2933's Be Star Companion
Five IR candidates and three optical candidates were recovered with  < 5% as potential companions of SGR 0755-2933.The abundance of sources is due to its large 3 ′ position error (Barthelmy et al. 2016).This sample includes the known 2MASS counterpart to the proposed companion star CPD-29 2176, 2MASS J07554248-2933535.We also recover a PS1 source within the 2MASS source's 60 mas error circle, PSO J118.9270-29.5649,and therefore conclude this is the optical counterpart of CPD-29 2176.However, we note that the association of SGR 0755-2933 with CPD-29 2176 is tenuous due to the former's lack of activity; only one burst localized to a 3 ′ region by the Swift telescope has been detected (Barthelmy et al. 2016).While a revised position uncertainty 0.4 ′′ was derived by Richardson et al. (2023) by assuming its association with CPD-29 2176, Doroshenko et al. (2021) identify 8 persistent X-ray sources within the 3 ′ error circle.This prevents a confident association of the burst with any persistent source, including CPD-29 2176.Although Richardson et al. (2023)'s spectral analysis confirms the presence of an HMXB containing CPD-29 2176 (formally 2SXPS J075542.5-293353), a second burst localized to the binary is required to confirm that it contains SGR 0755-2933.

CXOU J171405.7-38131
We do not recover the infrared source, VIRAC 384090103, proposed by Chrimes et al. (2022b) as a potential stellar companion to CXOU J171405.7-38131(Smith et al. 2018;Nikzat et al. 2022).If placed at the distance of the magnetar, 3.8 ± 0.5 kpc, the source has an −band magnitude   ≈ −0.9 ± 0.1 consistent with an OB star and its color  −  ≈ −0.5 ± 0.3 differs only marginally.We note that the magnetar lies along a highly extincted sightline outside of the Bayestar19 range, making its   ≈ 22.1 ± 0.8 from the   −   relation an upper limit (Predehl & Schmitt 1995).This biases the Monte Carlo p-value test since the extinction cannot be estimated as a function of distance.Therefore, the location of the magnetar and source make this analysis inconclusive, and we cannot rule out VIRACS 384090103 as a potential companion of CXOU J171405.7-38131.

RACS, VLASS, and NVSS Data and Search Method
While numerous continuum radio searches for SN shells have been conducted (see e.g. Green 2019;Lorimer et al. 1998;Dubner & Giacani 2015, for comprehensive reviews) as have searches for compact objects near SNRs (e.g.Kochanek et al. 2019;Kochanek 2018Kochanek , 2021Kochanek , 2023)), few targeted searches for SNRs near compact objects have been carried out beyond individual sources.Leveraging the young ages (≲ 1 Myr, e.g.Suzuki et al. 2021) we search for evidence of SNRs among the magnetar sample.
The RACS-low Data Release 1 (DR1) contains images with 15 ′′ resolution at 888 MHz, with a limiting sensitivity of 0.25 mJy beam −1 (McConnell et al. 2020;Hale et al. 2021).Public images are smoothed to a common pixel scale 25 ′′ .20 magnetars and 6 magnetar candidates without confirmed SNRs fall within the survey region −80 • <  < +30 • .We obtained 1 • × 1 • radio image cutouts via the RACS Data Access Portal27 .For two magnetars (SGR 0501+4516 and 1E 2259+586) with  > +30 • , images from the Karl G. Jansky Very Large Array Sky Survey (VLASS), which covers  > −40 • with 0.12 mJy sensitivity, were obtained using the Canadian Initiative for Radio Astronomy Data Analysis (CIRADA28 ) cutout service (Gordon et al. 2020(Gordon et al. , 2021(Gordon et al. , 2023a)).VLASS Quick-Look images29 are taken at 3 GHz with 1 ′′ pixel scale, which we bin to 25 ′′ to search for extended emission.For three magnetars covered by neither VLASS nor RACS (4U 0142+61, SGR 0418+5729, and SGR 2013+34), we obtain cutouts from the 1.4 GHz National Radio Astronomy and VLA Sky Survey (NVSS)30 , a precursor to VLASS (Condon et al. 1998).This catalog contains images with 15 ′′ pixels and RMS noise 0.5 mJy beam −1 , or ∼ 0.8 mJy beam −1 for 25 ′′ The proposed companion of CXOU J1714 is not recovered, but cannot be ruled out as a bound companion due to uncertainty in the extinction estimate.Note that the axes for CXOU J1714 differ between optical and infrared due to the use of different World Coordinate System (WCS) projections (ZPN) for VVV data (e.g.Greisen & Calabretta 2002).
MNRAS 000, 1-34 (2024) pixels 31 .One magnetar, SGR 1627-41, is not covered by any of the radio surveys discussed 32 .The Set of Identifications, Measurements and Bibliography for Astronomical Data (SIMBAD 33 ) database was queried for known SNRs, SNR candidates, and unidentified radio sources with 'ISM' object-types 34 within the same 1 • × 1 • region for each magnetar (Egret et al. 1991).If a previously proposed SNR candidate was not found in the initial SIMBAD query, it was added manually to the sample for the association test (for example, SNR HBH 9 is associated with SGR 0501+4516, but falls more than 1.3 • away and was missed in the query).Additional sources which may be SNRs are identified 'by eye' to include in the search (see Appendix B3.1).Sources identified as known radio galaxies or AGN were removed from the sample.While a radio image is not available for SGR 1627-41, we proceed with the association test using SIMBAD radio sources within 1 • .Through a Monte Carlo search, p-values are estimated with the criteria below: (i) The SNR and magnetar are located at the same RA and declination.
(ii) The magnetar's proper motion direction is consistent with the direction required to intersect the SNR.
(iii) The magnetar's proper motion and age are consistent with the observed angular separation of the magnetar and SNR.
(iv) The SNR's dynamical age is consistent with an SNR shell rather than a point source.
for which we define p-values  1 ,  2 ,  3 , and  4 , respectively.Appendix B3 contains a derivation of these p-values explicitly.Unlike in the binary companion searches, criteria (i)-(iii) relate to two independent cases: first that the magnetar is located within the SNR shell (i), and then that the magnetar has travelled some distance away from the SNR during its lifetime (ii and iii).Therefore, unless  1 < 0.1% or 1 − (1 −  2 )(1 −  3 ) < 0.1% 35 , we combine  1 ,  2 , and  3 using the Fisher statistic 36 : This follows a  2 distribution with 4 degrees of freedom, from which we compute the p-value   using the inverse survival function.We then combine this with  4 using: Note that  1 requires definition of the RA and declination errors for the magnetar and either the RA and declination errors or angular size and position angle for the SNR. 2 and  3 additionally require definition of the magnetar proper motion errors and the X-ray decay or characteristic age. 4 requires the angular size or position error to be defined.2005) distribution.Given that 9 of the 10 trials were successful in recovering the associated SNRs, we proceed with confidence that the image-plane and Monte Carlo search methods are effective.
Radio images for five magnetars with proposed SNR associations are shown in Figure 5.

Radio SNR Search Results
Among the 13 magnetars and magnetar candidates with no previous SNR associations, we identify two candidates which have  < 5%.
One of these, associated with SGR 0755-2933 is most likely a false positive which we discuss in Appendix C2.For SGR 1808-20, we find two viable SNR association candidates which require proper motion and distance measurements to properly assess.Of the eight magnetars and magnetar candidates with proposed but unconfirmed or rejected SNRs we recover five (and marginally recover a sixth with  = 5.6%), one of which is nearby the Galactic Center magnetar SGR J1745-2900 and another of which would suggest 3XMM J185246.6+003317 has an ejected pulsar companion.In this section, we discuss the SNR candidates for SGR 1808-20 and the recovered unconfirmed SNR candidates in detail.

SGR 1808-20
For SGR 1808-20, we find multiple candidate radio sources (29) with  < 0.1%, 10 of which were selected for analysis 'by-eye' rather than through the SIMBAD query.The large number of sources comes from the loose position ( ∼ 7 ′ error) and distance (no measurement) constraints, as only one burst has been detected from the magnetar since its discovery (Lamb et al. 2003). 2 ′ × 2 ′ VLASS cutouts at 1 ′′ scales were obtained for each candidate, revealing that 6 of the 29 sources were resolved and unlikely to be SNRs.One source is a known Young Stellar Object (YSO) while another is an unidentified sub-millimetric source, both of which are unlikely to be SNRs (Eden et al. 2017).14 are labelled as 'dense cores', or InfraRed Dark Clouds (IRDCs) which are pre-star forming regions that may be the birthplaces of massive stars (Jiao et al. 2023).IRDCs are distinguished from SNRs based on their high masses (∼ 100 ⊙ ) and hydrogen column densities (  ∼ 10 23 cm −2 ), essentially ruling out these candidates (e.g.Güver & Özel 2009).Six of the seven remaining sources are 'by-eye' candidates; we query the SIMBAD database within 1 ′ of each source to identify potentially missed sources of radio emission.One is nearby PSR J1808-2057 which has a radio flux 2.3 mJy and is likely responsible for the observed emission (Hobbs et al. 2004).Another is coincident with a ∼ 532 ⊙ molecular cloud discovered through sub-millimeter (∼ 0.1−1 THz) emission, making it a plausible radio source (Duarte-Cabral et al. 2021).One is coincident with an infrared-bright 'bubble'; these are star-forming sites in HII regions and are therefore also plausible radio sources (Simpson et al. 2012).Two others are coincident with stars, one an AGB star and the other an eclipsing binary system, both of which are potentially radio bright, though no radio flux has yet been catalogued (Heinze et al. 2018;Spangler et al. 1977;Cutri et al. 2003).While these five sources cannot be confidently ruled out, we will assume they are not associated given the possible contamination by nearby objects.This leaves two sources as candidates for association: SNR G009.7-00.0 and an unidentified radio source, which we formally label S24-RACS J180842-202849, with no SIMBAD objects within 1 ′ .For the remaining discussion, we will refer to the latter source as 'Source 1'.
Source 1 was identified by eye at 18 h 08 m 42 s − 20 • 28 ′ 45 ′′ and fit with an ellipse of approximate major and minor axes  ∼  ∼ 18 ′′ .This lies very near the position of another magnetar, SGR 1806-20, and the W31 HII region (e.g.Frail & Clifton 1989;Corbel & Eikenberry 2004;Moisés et al. 2011).The distance to W31 has been disputed in the literature; if we adopt the distance of SGR 1806-20, 8.7 +1.8  −1.5 kpc, Source 1 would have a radius ∼ 0.8 pc.At this size the SNR would be in the free-expansion phase, and assuming  0 ∼ 1 cm −3 , we estimate the age of the remnant to be between 0.3−1.2kyr (Sedov 2018).This is roughly consistent with T rot and T decay for magnetars, but the lack of data for SGR 1808-20 prevents a direct comparison.A detailed composition analysis is required to confirm Source 1 as a supernova remnant, though its morphological similarity to other false positives implies it is unlikely to be an SNR.SNR G009.7-00.0, on the contrary, was detected by Frail et al. (1994) and confirmed through its OH maser by Hewitt & Yusef-Zadeh (2009), which constrained its distance to ∼ 4.7 kpc.Yeung et al. (2016) find its distance is inconsistent with W31 and SGR 1806-20.At ∼ 4.7 kpc, the ∼ 11 − 15 ′ angular size corresponds to a radius ∼ 8 − 10 pc, suggesting a free-expansion age of ∼ 10 kyr.Both Source 1 and SNR G009.7-00.0appear to be viable candidates for association with SGR 1808-20, but a magnetar distance measure and more tightly constrained position are necessary to investigate further.The radio image is shown in Figure 5.
Considering this case, let the initially bound massive star progenitors of 3XMM J185246.6+003317 and CXOU J185238.6+004020be called stars 1 and 2, respectively.When the CCSN of star 1 disrupts the binary39 , the remnant magnetar (3XMM J185246.6+003317) and star 2 are ejected at ∼ 300 km s −1 and ∼ 6 km s −1 , respectively (Hobbs et al. 2005;Renzo et al. 2019).After travelling for ∼ 0.7 Myr based on the magnetar's age, the CCSN of star 2 would then kick the remnant pulsar (CXOU J185238.6+004020) at ∼ 300 km s −1 , leaving behind the remnant shell Kes 79.This sequence of events is shown in Figure 4. We investigate this scenario by applying a modified Gaia unbound companion search pipeline to the magnetar and CXOU J185238.6+004020'sGaia counterpart, Gaia DR3 4266508881354196736.We split the companion's trajectory into two parts with different velocities to account for the two CCSN events, assume the absolute magnitude of the companion is   > 0, and neglect  2 since the magnetar's and pulsar's travel times no longer coincide directly.10000 trials (rather than 100) are used since the magnetar has no proper motion measurment.
The resulting trajectory plot is shown in Figure 4; we find a pvalue  1,3 = 0.02%, suggesting that this scenario is likely.From the 497 valid samples out of 10000, the scenario requires the magnetar to have velocity −5.8 mas yr −1 ,   ≈ −9.8 +4.9 −20.8 mas yr −1 ) and a travel time T 3D,mag ≈ 40.1 +56.1 −21.8 kyr, which is consistent with expected magnetar ages, but disagrees with Rea et al. (2014), Zhou et al. (2014), and T decay .Despite the robustness of the X-ray timing solution, magnetars are known to glitch, anti-glitch, and slowly recover from glitches during outbursts such as the one used for timing analysis (e.g.Tong & Huang 2020;Younes et al. 2023;Scholz et al. 2014;Archibald et al. 2013).Coupled with Rea et al. (2014) observing the outburst at a late time, this implies that 3XMM J185245.6+003317 may have experienced a glitch event leading to its undetectable spindown.If we impose T 3D,mag as an age constraint, then using the period  = 11.55871346(6)s we derive a revised spindown rate of  ≈ 4.6 +5.4  −2.7 × 10 −12 s s −1 , and dipolar magnetic field  dip ≈ 2.3 +1.1 −0.2 × 10 14 G, both of which are more consistent with the magnetar sample.
The pulsar requires a velocity | − →  psr,1 |≈ 8 +11 −4 km s −1 following the first CCSN, and | − →  psr,2 |≈ 422 +192 −172 km s −1 following the second.The most supportive evidence for this scenario is the position of the pulsar's CCSN, estimated from valid samples as 18 h 52 m 38.s 60 + 00 • 40 ′ 20.′′ 09 with 1 errors of 0.2 ′ and 0.4 ′ in RA and declination respectively.This is roughly 8 away from the centroid of Kes 79, 18 h 52 m 48 s +00 • 41 ′ 00 ′′ , but is well within its ∼ 5 ′ radius.We obtain an estimate of 18 h 52 m 56.s 5 + 00 • 43 ′ 57.′′ 5 with errors   ∼ 3 ′ and   ∼ 6 ′ for the position of the first CCSN.This also lies within the extent of Kes 79, and suggests that its unique multi-shell structure may result from interactions between the magnetar and pulsar SNR shock fronts.Other single-SNR scenarios have been presented to explain this structure as, for example, the interaction of the SNR with a molecular cloud or the progenitor's stellar wind (Velusamy et al. 1991;Giacani et al. 2009;Sun et al. 2004;Sato et al. 2016;Zhou et al. 2016).We conclude that it is likely that the progenitors of 3XMM J185246.6+003317 and CXOU J185238.6+004020 began in a bound system prior to the former's CCSN, and a proper motion measurement for 3XMM J185246.6+003317 would further constrain their formation history.

SGR J1745-2900, the Galactic Center Magnetar
SGR J1745-2900 is a unique magnetar at the Galactic Center (GC) of the Milky Way (Mori et al. 2013b).The radio emission from the supermassive black hole Sag A * and high density of radio sources in this region make it difficult to identify any associations with confidence.The black hole would also distort an SNR shell from SGR J1745-2900 (e.g.Yalinewich et al. 2017).However, Ponti et al. (2015) used high resolution X-ray observations to catalog extended regions near the GC, identifying the GC X-ray/radio lobes as a warped and mirrored SNR.SGR J1745-2900 lies near the intersection of the two lobes, and Yalinewich et al. (2017) use numerical simulations to demonstrate that an 'engulfing' SNR explosion at this location could produce the observed SNR structure.
Given this analysis, the radio lobes were added to the simulation using the centroid position (17 h 45 m 32.s 3 − 29 • 2 ′ 32.′′ 5) and angular diameter ∼ 5.8 ′ from Ponti et al. (2015).Through the Monte Carlo simulation, we recover the radio lobes, as well as a gamma-ray source 0FGL J1746.0-2900found not to be an SNR (Acero et al. 2016) and an HII region (Nguyen et al. 2015), with  < 0.1%.In our simplified model of the lobes, we do not account for the irregular structure of the SNR, but instead assume the SNR is centered at the lobes' centroid with angular size equal to a single lobe.This simplification should not greatly impact the results since only a lower limit is imposed on the SNR's dynamical age based on the angular size.The position of the radio lobes' centroid is marginally consistent with the trajectory of SGR J1745+1900 (see Figure 5).With a proper motion magnitude ∼ 6.4 mas yr −1 , the ∼ 2.7 ′ separation suggests an SNR age ∼ 25 kyr, which is slightly lower than the decay age estimate T decay ≈ 64.9 kyr.Further interpretation requires a more detailed analysis of the GC lobe geometry, but our results support a scenario where the bipolar radio lobes are the SNR of SGR J1745+1900.

Other Magnetars with Unconfirmed SNRs
Five other magnetars in the sample were identified in previous searches as potential associations of nearby supernovae, but have not yet been confirmed: SGR 1900+14, SGR 1806-20, SGR 0526-66, SGR 0501+4516, and SGR 1801-23.SGR 1900+14's association with SNR G042.8+00.6 (for which we find  = 5.6%) was rejected by Vrba et al. (2000) given the former's more likely association with a star cluster.We identify one additional by-eye source, S24-RACS J190707+091737, as an SNR candidate with  < 0.1%, but rule it a false positive as it disagrees with the proper motion direction of SGR 1900+14.SGR 1806-20's association with AX J1808.6-2024(SNR G10.0-0.3) was also rejected by Hurley et al. (1999), suggesting it is not an SNR but a radio nebula powered by its central LBV star.While there is insufficient position data to estimate a p-value for the nebula, we do recover it when using the angular size of its potentially associated gamma ray emission, Fermi J1808.2-2029(H.E. S. S. Collaboration et al. 2018).Finally, we do not recover XTE J1810-197's association with SNR G11.0-0.0,supporting the age considerations which ruled out the candidate in Ding et al. (2020b).
SGR 0526-66 is located within the Large Magellanic Cloud (LMC) near SNR N49 (SNR J052559-660453; Klose et al. 2004).However, its association with SNR N49 could not be confirmed as the magnetar's proper motion is unknown, which is necessary to localize its birth place to SNR N49's host star cluster, SL 463.With our pipeline, we obtain a p-value  < 0.1% for this association.Notably, five other 'by-eye' radio sources were recovered with  < 5%; however, each of these fall within 1 ′ of possible known radio sources that are not within the LMC.While SGR 0526-66 does not have a proper motion measurement, its position error of 0.6 ′′ compared to the 1.4 ′ angular size result in a highly confident association.The large distance 53.6 ± 1.2 kpc suggests N49 has a tranverse size of roughly 22 pc, while the 0.6 ′′ position error for the magnetar corresponds to ∼ 0.1 pc.Therefore it is likely that N49 contains SGR 0526-66, supporting Klose et al. (2004)'s claim of association.
For SGR 0501+4516, Gaensler & Chatterjee (2008) localized the magnetar to the outer rim of SNR HBH 9 (HB9 or SNR G160.9+2.6) at a position ∼ 15 ′ from the next nearest SNR.The large angular separation ∼ 1.3 • , however, initially excluded SNR HBH 9 from the test sample.After adding this as a bye-eye radio source, we obtain a p-value  = 7.0%, providing insufficient evidence to confirm the association.This is most likely due to the large angular separation ∼ 1.3 • ; at the magnetar's distance 2 ± 0.3 kpc, this corresponds to a transverse separation of ∼ 45 pc.For the magnetar's decay age T decay ≈ 59.4 kyr, this would require a velocity ∼ 730 km s −1 for the magnetar, which agrees well with the Hobbs et al. ( 2005) distribution expected for neutron stars.Therefore, given the marginal p-value recovered and the achievable kick velocity, we find this association remains possible, and a proper motion measurement is required to confirm or reject it.

DISCUSSION
Through this multiwavelength search of 31 magnetars and magnetar candidates, we identify one magnetar with a potential unbound OB star companion, SGR J1822.3-1606, but require proper motion data to explore further.We note that magnetar 3XMM J185246.6+003317 may be an unbound system whose progenitor was bound to that of an X-ray pulsar, CXOU J185238.6+004020.While we recover the Be star companion of SGR 0755-2933 proposed by Richardson et al. ( 2023) and the unbound Wolf-Rayet star companion of CXOU J164710.2-455216 proposed by Clark et al. (2014), we are unable to confidently associate the IR counterpart of CXOU J171405.7-381031 proposed as a bound companion by Chrimes et al. (2022b).Through the radio search, we recover nine of ten confirmed SNR associations, propose scenarios supporting four of five unconfirmed SNR associations, and identify a new candidate SNR associated with SGR 1808-20 which requires a proper motion and distance measurement for confirmation.
Table 3 summarizes our results, which prompt a discussion of the CCSN model for magnetar formation.The recovery of only two unbound binary candidates is inconsistent with Chrimes et al. (2022b) and Renzo et al. (2019) simulations predicting that   ∼ 38 − 56% of magnetars should be nearby unbound companion stars.We begin this section by analyzing our search's sensitivity before exploring the implications of the search on magnetars' massive star progenitors.We then consider potential limitations in our work and investigate the possibility of more exotic magnetar formation channels.MNRAS 000, 1-34 (2024)

Completeness and Sensitivity Considerations
For each search, we estimate the expected ranges of intrinsic parameters and derive the ranges of their corresponding observables in order to compare to our sensitivity thresholds.For the bound companion search, we search for OB stars with intrinsic magnitudes in the following ranges: (i) Radio luminosity  ≈ 10 42 erg s −1 in the Rayleigh-Jeans limit (Barker et al. 2022(Barker et al. , 2023;;Hamuy 2003) (ii) Angular size 25 ′′ ≲ ΔΘ ≲ 1 • (from observed SNR sizes in the Green (1987) catalog and considerations of free-and Sedov-Taylor-expansion in Truelove & McKee (1999)) By combining these constraints with the distance, extinction, and proper motion of each magnetar (see Table 1), we estimate the expected range of apparent magnitudes (−, −, and −band), SNR flux at 1 GHz, and relative proper motion Δ between the magnetar and an unbound companion.Figure 6 and Table 3 summarize the results.Four magnetars have no distance measurements and are not included in this analysis.
For the Gaia search we observe that the relative proper motion of each magnetar and a walk/runaway companion should be detectable down to   ∼ 21 given the proper motion uncertainty   ≈ 1.4 mas yr −1 (Collaboration et al. 2022).Only the SMC magnetar CXOU J010043.1-721134 is partially incomplete for faint stars due to its large distance.However, four magnetars lie along high-extinction sightlines with   ≳ 17 (CXOU J171405.7-381031,1E 1547.0-5408,PSR J1622-4950, SGR 1627-41) making unbound companions undetectable at Gaia's   < 20.7 sensitivity.19 magnetars are partially complete, meaning that our search may miss fainter B-star companions.These results are more promising than indicated by the simulations of Chrimes et al. (2023), which determined that without the aid of next generation deep-field surveys, Gaia may recover only a small fraction of unbound neutron star companions.This demonstrates that incorporating a trajectory search into our Monte Carlo analysis can make significant improvements to our search sensitivity.
The bound search is complete or partially complete for most magnetars covered by the IR survey; only SGR 1627-41, which has a large extinction   ≈ 55.87 ± 11.17, and the two magnetars in Magellanic Clouds, CXOU J010043.1-721134 and SGR 0526-66 are incomplete.The optical search, on the other hand, is incomplete for six magnetars due to their high extinction.This is in agreement with Chrimes et al. (2022b), which concluded that IR searches are more effective for OB stars due to the ∼ 4× lower extinction in −band compared to −band (e.g. Green 2019).
Finally, the SNR search is complete in all cases for typical radio luminosities and angular sizes in the free-or Sedov-Taylor expansion 40 See Section 5.5 of Gaia DR3 Documentation (https://gea.esac.esa. int/archive/documentation/GDR3/).
phases (e.g.Truelove & McKee 1999;Salpeter 1955).We note, however, that a bright synchrotron radio shell may not be produced if the magnetar inhabits a low-density local environment or lacks a significant internal energy source (e.g.Frail et al. 1995;Chevalier & Clegg 1985;Nomoto 1987;Leung et al. 2020).Future work may expand on this sensitivity analysis with detailed forward modelling of SNRs, but this is beyond the scope of this work.Table 3 indicates whether each magnetar's search is considered 'partially or fully complete' (✓) or 'incomplete' (✗) in the above analysis.We conclude that our search for companion stars and SNRs is reasonably sensitive along sightlines with typical extinction and distances < 10 kpc.

Markov-Chain Monte Carlo Analysis of Magnetar Binary Fractions
With more well defined estimates of our search completeness, we next determine the implications of our results on the CCSN model.We adapt the method of Kochanek (2023) to constrain the binary parameters of the magnetar progenitor population using a Markov-Chain Monte Carlo (MCMC) simulation.Uniform priors are chosen to estimate the fraction of systems which were not binaries at death,   ; the fraction that were unbound at death,   ; and the fraction of interacting binaries remaining bound after the CCSN,   .We relate these to the total fraction of binaries remaining bound after the CCSN,   , and the fraction of non-interacting binaries remaining bound after the CCSN,   , by The likelihood is the joint multinomial probability of the observed samples; we define the sample for this analysis incorporating the completeness ratios in Table 3.Of the 31 magnetars and candidates, three systems, CXOU J164710.2-455216,3XMM J185246.6+003317, and SGR J1822.3-1606, were proposed as unbound binary systems.The candidate companion of CXOU J171405.7-381031 from Chrimes et al. (2022b) is unconfirmed based on the optical and IR search results.For the purpose of our analysis, we will accept SGR 0755-2933 and CPD-29 2176 as a bound interacting binary system based on the results of Richardson et al. (2023).Of the remaining sample, three magnetars are constrained to be not currently bound, two are constrained not to have unbound companions, and the remaining 21 are considered not to be binaries at the time of the CCSN.With this, the likelihood can be written as: where (1 −   ) is the fraction that are not unbound systems and (   +   ) is the fraction that are not bound systems.The sampled posteriors are used to estimate the median values and 90% confidence intervals for   ,   , and   .
In keeping with the notation of Kochanek (2023), we define the case above as 'SNR incomplete', where we assume that all magnetars in the sample originate in CCSNe.While this is likely the primary formation channel, it is possible that other processes may form magnetars at lower rates.Considering this, we define the 'SNR complete' case, in which magnetars without confirmed or proposed SNR associations are not formed via CCSNe, and therefore do not contribute to the binary fractions above.In this case, the likelihood becomes: where   is the fraction of magnetars which do not originate in CCSNe.This modifies the definition for   to   = 1 −   −   −   .✓(S,V)  For each magnetar (written with its abbreviated name), the SNR candidates, bound companion candidates, and unbound companion candidates with  < 5% are given, with p-values in brackets.Candidates suspected to be false positives are given in italics.Sources identified 'by-eye' which are not confidently associated with any known object are labelled using the format, 'S24-survey-Jhhmmss±ddmmss', where survey is the survey used and the italicized portion is the J2000 RA and declination of the centroid.For each search, the completeness along each magnetar's sightline is indicated by a ✓, meaning the search is either partially or fully sensitive to stars or remnants, or by an ✗, meaning that even the brightest expected stars or remnants would not be detectable.(2000).Detailed discussion for the latter is provided in Section 5.3.* CXOU J185238.6+004020 is an X-ray pulsar associated with Kes 79.We explore the scenario that the progenitors of 3XMM J1852 and this pulsar were bound in Section 5.3.
Note this case is more restrictive, as SGR 0755+2933, SGR J1822.3-1606, and 13 other magnetars for which the radio search is > 75% complete have no confirmed SNR in this or any other work.A final case, 'SNR/unbound complete' is defined if we also take unbound magnetar binaries as evidence of CCSNe; in this case, we include CXOU J164710.2-455216,SGR J1822.3-1606, and 3XMM J185246.6+003317 as unbound systems, and the likelihood becomes: Results for the three cases are summarized in Table 4.In the 'SNR incomplete' case, we find a median bound binary fraction   ≈ 9 +6 −4 %, and unbound binary fraction   ≈ 12 +6 −5 %, where the errors are 1 uncertainties.While the bound fraction is marginally consistent with the predictions from Chrimes et al. (2022b), Kochanek et al. (2019), andRenzo et al. (2019) population synthesis results, the unbound fraction significantly disagrees with their prediction of   ∼ 38 − 56%.Furthermore, we find that a median of only 1 −   ≈ 22 +8 −7 % of magnetars were binaries at death, which disagrees with the estimate of ∼ 50 − 60% from population synthesis with 90% confidence.
In the 'SNR complete' case, we estimate a median unbound fraction   ≈ 3 +3 −2 %, and a 90% upper limit on the bound fraction of   ≲ 12%, the former of which also disagrees with the Chrimes et al. (2022b), Kochanek et al. (2019), andRenzo et al. (2019) population synthesis estimates.We find a median 1 −   ≈ 63 +9 −8 % of magnetars undergo CCSNe while in binaries, which more closely agrees with population synthesis than the 'SNR incomplete' case.Most significantly, we constrain the fraction of magnetars formed by channels other than CCSNe to   ≈ 53 +8 −10 %, with a 90% confidence interval from 37 − 66%.This would contradict the view that CCSNe are the dominant formation pathway (e.g.Schneider et al. 2019;White et al. 2022;Keane & Kramer 2008).Furthermore, no magnetars have been confirmed to form from non-CCSN channels.
A similar result is obtained in the 'SNR/unbound complete' case; notably, the larger unbound fraction   ≈ 8 +5 −3 % allows a lower non-CCSN fraction   ≈ 46 ± 9%, while the fraction of magnetars from binaries 1 −   ≈ 61 ± 9% still agrees with population synthesis results.While the 'SNR/unbound complete' and 'SNR complete' cases do not differ in their underlying assumptions about alternate formation channels, they demonstrate that the detection of unbound magnetar binaries allows for tighter constraints to be placed on   .Furthermore, these cases motivate discussion of alternate formation channels, such as Accretion Induced Collapse (AIC) of magnetized White Dwarfs, electron-capture supernovae of less massive stars, or binary Neutron Star mergers, to which we return in Section 6.3 (Ruiter et al. 2019;Fryer et al. 1999;Giacomazzo & Perna 2013;Giacomazzo et al. 2015;Miyaji et al. 1980).

Implications for Massive OB Star Binary and Merger Fractions
Many models of magnetar formation attribute their strong magnetic fields to a dynamo process driven by accretion from ablated stellar companions (e.g.Popov 2020;Revnivtsev & Mereghetti 2016;Fuller & Lu 2022, and references therein).The low   and   fractions in Table 4 may indicate that other magnetic field growth mechanisms are more common, such as the merger of massive stars.In this scenario, the magnetic field of the primary star is first amplified up to ≲ 10 8 G through accretion before further amplification through merger of the cores, producing a highly magnetized massive star such as  Sco and the Wolf-Rayet star in the HD 45166 binary (Schneider et al. 2019;Shenar et al. 2023).If magnetic flux is conserved during the CCSN event, a neutron star remnant with a magnetar-like surface magnetic field (∼ 10 14 G) could be achieved (e.g.Popov 2020;Schneider et al. 2019;Obergaulinger et al. 2014, and references therein).A high frequency of such pre-CCSN mergers could potentially explain the lack of magnetars with bound or unbound companions.The fraction of massive stellar binaries that merge prior to CCSNe,   , can be inferred by assuming an initial massive star companion fraction  0 ≈ 84 ± 9%, derived in Moe & Di Stefano (2017), for early B-type stars. 41While some simulations suggest that CCSNe  can form from mergers of intermediate mass stars (5 − 8 ⊙ ), we focus this discussion on binaries where both stars are > 7.5 ⊙ (e.g.Zapartas et al. 2017).Relating the initial companion fraction to the fraction of systems in binaries at death: where   = 0.426 is the fraction of binaries with a secondary star less massive than the primary 42 .  effectively estimates the fraction of CCSNe from explosions of the secondary star, which are excluded from this analysis since we consider CCSN of stars with stellar companions 43 (see Kochanek 2009;Kochanek et al. 2019;Kobulnicky & Fryer 2007, for detailed discussions).If we adopt   and   from Table 4 and use a Gaussian distribution 44 for  0 , we obtain median values   ≈ 70 +10 −12 % for the 'SNR complete' case,   ≈ 88 +5 −8 % for the 'SNR incomplete' case, and   ≈ 82 +7 −9 % for the 'SNR/unbound complete' case.In all three cases,   is significantly higher than the merger fractions estimated for neutron star progenitors in population synthesis from Kochanek et al. (2019) (∼ 48%) and Renzo et al. (2019)

(∼ 22 +26
−9 %).It is also possible that magnetar progenitors have an above average merger rate compared to other neutron star species.Corner plots summarizing this analysis are shown in the left-hand plot of Figure 7.
The 'SNR incomplete' case does not allow for alternate formation channels, suggesting instead that a large fraction   ≈ 80% were formed from single stars.For a fixed binary fraction  0 , this requires an inflated merger rate   ≈ 75%.The 'SNR/unbound complete' 42 This is derived by integrating a Salpeter-like ( = −2.35)initial mass function (IMF;Salpeter 1955;Moe & Di Stefano 2017); see Section 1 and Equation 1 of Kochanek et al. (2019), and references therein.Note that in Kochanek et al. (2019) Equation 1,   is defined as the fraction of CCSNe that occur in stellar binaries, whereas we here define   as the fraction of systems for which the binary remains bound following the CCSN. 43We assume that magnetar 3XMM J185246.6+003317 was formed from the first explosion, as is implied from our analysis in Section 4. 44 More rigorously,  0 follows a binomial distribution, but we find this makes little difference to the results.The Gaussian method is adopted to more easily incorporate the uncertainty on  0 .
case is one step removed, taking unbound companions as evidence of a disruption event (CCSN) but also allowing for  nc ≈ 45% to form from other formation channels (e.g.AIC, neutron star merger).For this reason,   ≈ 40% is lowered, but the unbound fraction   ≈ 10% is relatively unchanged and the merger rate   ≈ 80% increases only marginally (< 1 change).The 'SNR complete' case, however, excludes both SGR 1822.3-1606 and CXOU J164710.2-455216, which do not have SNR detections, from the unbound sample.SGR 0755-2933, the only proposed bound system, is excluded for the same reason.Therefore, both the merger fraction   ≈ 90% and the fraction from alternate channels  nc ≈ 55% increase to accommodate the lowered unbound fraction   ≈ 5%.In all three cases, the merger fraction   is larger than expected, requiring us to reconsider the conditions that lead to massive OB star mergers.However, a more strict interpretation ('complete' cases) is consistent with  nc ≈ 31 − 66% forming from alternate channels.
Alternatively, we can allow both  0 and   to vary, leaving the progenitor mass and stellar type unconstrained.Under this assumption, the right-hand plot of Figure 7 shows the 90% confidence regions for  0 and   in each case.This demonstrates that a wider range of merger rates can be tolerated if magnetars can form from the mergers of less massive stars (Moe & Di Stefano 2017).This would marginally support the analysis of Zapartas et al. (2017), which identified a population of CCSNe from the mergers of intermediate mass stars (5 − 7.5 ⊙ ).However, a reliable comparison requires a broader search that targets both massive and intermediate mass stars, which we leave to future work.Figure 7 summarizes the results of this analysis in a flowchart showing the evolution of stellar progenitors to magnetar remnants.

Limitations in Comparison to Population Synthesis Results
While we compare our resulting binary, merger, and CCSN fractions to population synthesis simulations, a discussion of the limitations of such analysis is warranted.We focus this discussion on the Chrimes et al. (2022b) simulations, though similar comparisons can be made with the Kochanek et al. (2019) and Renzo et al. (2019) simulations.Chrimes et al. (2022b) utilized Binary Population and Spectral Synthesis (BPASS) models and imposed constraints on the total mass (> 2.0 ⊙ ), CO core mass (> 1.38 ⊙ ),    (2022b, teal).Observational constraints from this work are shown in red for the 'SNR incomplete' case and blue for the 'SNR/unbound complete' case.A lower binary fraction would allow for a lower merger fraction for magnetar progenitors, but would require lower mass progenitors.MNRAS 000, 1-34 ( 2024) and remnant mass (1.38 <  rem < 3.0) to identify stellar models that evolve into magnetars.One potential source of error is the Initial Mass Function (IMF), for which Chrimes et al. (2022b) use the two-component Kroupa et al. (1993) broken power law model with index  = −1.30for masses 0.1 − 0.5  ⊙ and  = −2.35for masses 0.5 − 300  ⊙ .Alternatively, some analyses have proposed a smoothly varying IMF with a mass-dependent  which is steeper for high mass stars (e.g.Kroupa 2002;Kroupa & Weidner 2003;Miller & Scalo 1979;Chabrier 2001;Weidner & Kroupa 2005).This could lead to a larger ratio of low mass to high mass stars.The relationship of the IMF to the binary fraction also has a complex dependency on the mass-ratio ( 2 / 1 ) distribution.While a uniform distribution is assumed to compute   ≈ 0.426 (e.g.Kochanek 2018;Kochanek et al. 2019;Kouwenhoven et al. 2005), Renzo et al. (2019) find that a negative power law distribution results in a larger merger fracton   .This demonstrates that both the IMF and mass ratio distribution can have a significant impact on the massive star binary fraction, and motivates detailed consideration in future studies.An additional caveat is that relatively loose constraints were used to gather the magnetar sample for this work in comparison to previous analyses.For example, Chrimes et al. (2022b) excludes magnetars in Magellanic Clouds, as well as those without distance estimates, imaging data, or photometric data.With the small 31-magnetar sample, we have chosen not to place arbitrary limits beyond the completeness constraints discussed in Section 6.1.The reported p-values provide quantitative measures of confidence that incorporate large distances, errors, and accuracy of imaging data.Furthermore, the irregularity of transient magnetar outbursts could have marginal effects on our analysis; for example, both SGR 1935+2154 and PSR J1846-0258 have recently entered X-ray outburst phases which may tighten constraints on their ages (Ibrahim et al. 2024;Laha et al. 2020;Krimm et al. 2020;Sathyaprakash et al. 2024).SGR 0755-2933 and SGR 1808-20, contrarily have only one burst detection, each (Lamb et al. 2003;Barthelmy et al. 2016).Future tests may benefit from partitioning the sample based on their distance, proper motion, extinction properties, or outburst phase to characterize its affect on the search.

Comparing Magnetar Binary Statistics to the Neutron Star Population
The apparent differences between the magnetar and neutron star formation processes limit the generalization of our results to the full neutron star population.Magneto-hydrodynamic (MHD) plasma simulations of CCSNe demonstrate that the magneto-rotational instability (MRI) or a dynamo action can amplify the progenitor stars' fossil fields during collapse (e.g.Akiyama et al. 2003;Obergaulinger et al. 2014Obergaulinger et al. , 2018;;Raynaud et al. 2020;Obergaulinger et al. 2009).However both MRI and dynamos require either high progenitor magnetic fields (≳ 10 8−11 G core magnetic field; Obergaulinger et al. 2018;Schneider et al. 2019) or high specific angular momentum (≳ 9 × 10 14 cm 2 s −1 ; Obergaulinger et al. 2018).Such conditions were considered at a high level by Popov & Prokhorov (2006), whose population synthesis experiment estimates that ∼ 8 − 9% of neutron stars would be magnetars, though the details of MRI and dynamo evolution were not incorporated.Thus magnetars and neutron stars probe different progenitor star populations that are not well-defined.Such a partition of the progenitor star population would also disfavor models presenting magnetar emission as a continuum among neutron stars.Furthermore, some evidence suggests that highly magnetized massive stars can be formed through the merger of stars with standard surface magnetic field strengths (∼ 1 − 100 G; e.g.Schneider et al. 2019).This is presented as a possible origin of blue straggler stars with ∼ 1 kG surface fields (e.g. Sco; Schneider et al. 2016;Donati et al. 2006), which are thought to be magnetar progenitors.The high merger rate found in this work,   ∼ 70−80%, supports this scenario since it implies neutron stars with typical field strengths may have a smaller progenitor merger rate.This may also explain the discrepancy between our work and the ∼ 22 +26 −9 % merger rate suggested by Renzo et al. (2019)'s population synthesis.Without a more detailed understanding of the partitioning of magnetar and pulsar progenitors, our results are confined to the magnetar population and cannot be confidently applied to neutron stars as a whole.

Alternate Magnetar Formation Channels
In this section, we consider the implications of the 'SNR complete' and 'SNR/unbound complete' cases which require large contributions   ∼ 45 − 55% from non-CCSN formation channels.First, consider how this alters the birth rate assumptions for the neutron star population.Keane & Kramer (2008) predicted that the Galactic CCSN rate (1.9 ± 1.1 century −1 ) is too low to source the observed population of magnetars, pulsars, XDINS, and RRATs (total 10.8 +7.0 −5.0 century −1 ).If we incorporate the results of the 'SNR complete' and 'SNR/unbound complete' cases, which predict   ≈ 55% and   ≈ 45% of magnetars are formed by alternate channels, respectively, then the total neutron star birth rate from CCSNe is decreased to ∼ 5 century −1 and ∼ 6 century −1 .This is still inconsistent with the CCSN rate, and suggests that alternate channels are unlikely to wholly explain the observed binary fractions.
It is possible that alternate formation channels are partially responsible for our observations.One such channel is Accretion Induced Collapse (AIC), in which ONe White Dwarfs accrete from Helium burning giant companions before collapsing (e.g., Ruiter et al. 2019;Fryer et al. 1999).Following AIC, either a bound neutron star-white dwarf binary is formed, or the donor star is ejected as a runaway.While the ejected mass from AIC is expected to be lower than that of a CCSN (∼ 10 −3 − 0.05 ⊙ ), low mass white dwarf companions would require less energy to become unbound than massive stars (Ruiter et al. 2019;El-Badry et al. 2023;Shen et al. 2018).
Direct evidence for AIC has not yet been observed, making them unlikely to make up the required   ∼ 45 − 55%; however, Sharma et al. (2023) observed that long delay times and high-mass host galaxies for two Fast Radio Bursts (FRBs) support origins from neutron stars or magnetars formed through AIC (see also Woodland et al. 2023;Gordon et al. 2023b, for comparisons of offset and host galaxy flux distributions between FRBs and CCSNe).More unique formation processes such as neutron star mergers require gravitational wave detection and followup, but are unlikely to supercede CCSNe in magnetar production (e.g.Giacomazzo & Perna 2013;Giacomazzo et al. 2015).Electron-capture supernovae, which occur in less massive (8 − 11 ⊙ ) stars, also are unlikely to produce magnetars as they impart weaker explosion kicks to their companions, resulting in a higher surviving binary fraction,   (Poelarends et al. 2017;Miyaji et al. 1980;Nomoto 1987;Leung et al. 2020).Therefore, the results of this search continue to support CCSNe as the most common magnetar formation channel, and we expect future improvements on radio survey sensitivity to recover any remaining SNRs with sufficient local density and explosion energy (e.g.Frail et al. 1995).
In this work, we have searched optical, IR, and radio image and point source catalogs for evidence of CCSNe and massive OB stars associated with magnetars.Out of the 31 known magnetars and candidates, we identify two candidates for unbound magnetar binaries at > 95% confidence: a massive OB star in association with SGR 1822.3-1606, and an X-ray pulsar in association with 3XMM J185246.6+003317.We also recover the proposed Be star companion of SGR 0755-2933 and marginally detect the proposed unbound companion of CXOU J164710.2-455216.Through the radio search we recover nine out of ten confirmed SNR associations and four of five unconfirmed associations; additional SNR candidates for SGR 1808-20 are identified but require a distance and proper motion measurement for the magnetar to confirm.Through an MCMC analysis, we find only 5 − 24% of magnetars have unbound massive star companions with 90% confidence.This is much lower than is predicted by neutron star population synthesis simulations from Chrimes et al. (2022b), Kochanek et al. (2019), andRenzo et al. (2019).
Our results support a high merger rate for massive OB star progenitors between 48 − 86%.It is unclear if only magnetars or neutron stars as a whole preferentially form from mergers.While alternate formation channels such as AIC and neutron star mergers may account for a small fraction of magnetars, it is unlikely that they are solely responsible for the lack of observed unbound binaries, as this would require ∼ 45 − 55% of magnetars to follow such channels and would not entirely solve the neutron star birth rate problem.
Significant improvements to this search can be made with proper motion measurements for the 23 magnetars without constraints.This may be achieved with next-generation, high sensitivity infrared telescopes like the JWST and NGRST.Similarly, Very Long Baseline Interferometry (VLBI) parallax measurements are available for only two magnetars; targeting the remaining sample will provide independent distance measurements to better inform the association tests.Next generation radio telescopes like the Deep Synoptic Array (DSA-2000) and the SKA will also detect a large sample of FRBs to investigate their origins, potentially constraining the magnetar birth rate from an extragalactic perspective.

DATA AVAILABILITY
The data underlying this article, including the magnetar parameters in Tables 1 and Tables 2 as well as derived values and p-values for each optical, IR, and radio source included in the search, are available in machine-readable format in the CaltechDATA repository, at https://doi.org/10.22002/cfj2d-zbc71.
ity (https://ror.org/05qajvd42).Operation of ASKAP is funded by the Australian Government with support from the National Collaborative Research Infrastructure Strategy.ASKAP uses the resources of the Pawsey Supercomputing Research Centre.Establishment of ASKAP, Inyarrimanha Ilgari Bundara, the CSIRO Murchison Radio-astronomy Observatory and the Pawsey Supercomputing Research Centre are initiatives of the Australian Government, with support from the Government of Western Australia and the Science and Industry Endowment Fund.This paper includes archived data obtained through the CSIRO ASKAP Science Data Archive, CASDA (https://data.csiro.au).This work has utilized data from the LoTSS survey and image archive, which is described in detail in Shimwell et al. (2022).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.This research has made use of the CIRADA cutout service at cutouts.cirada.ca,operated by the Canadian Initiative for Radio Astronomy Data Analysis (CIRADA).CIRADA is funded by a grant from the Canada Foundation for Innovation 2017 Innovation Fund (Project 35999), as well as by the Provinces of Ontario, British Columbia, Alberta, Manitoba and Quebec, in collaboration with the National Research Council of Canada, the US National Radio Astronomy Observatory and Australia's Commonwealth Scientific and Industrial Research Organisation.
We proceed as in the PS1 search to compute the search statistic ĝ1 and its p-value  1 , utilizing the 95 th percentile contour as a threshold and simulating the H0 distribution of angular separation, Δ, using mean separations between 1−5 pos {1, 1} for RA and 1−5 pos {0, 0} for declination.
For the second condition we follow the same method as for PS1 sources using instead the − band and − band magnitudes,   and   .Following Chrimes et al. (2022b), we expect OB stars to have −band magnitudes in the range −2.76 ≲   ≲ 3.34 and  −  colors between −0.14 ≲  − ≲ 0.31.After drawing samples of   and  −  from uniform distributions in these ranges, we compute   = ( − ) +   .The Bayestar19 extinction library is again queried, using   ≈ 0.4690 − and   ≈ 0.7927 − , and the inferred distance is estimated from both the − and − bands,  IR,J and  IR,H .The search statistic is taken as the fraction of sources for which either  IR,J or  IR,H falls outside the 95 th percentile of the H1 distribution.

B3 Radio Search for SNR Shells
The radio search uses the position, proper motion, and angular size to define association p-values, which cover two independent cases: first, that the magnetar resides within the SNR shell, and second, that the magnetar was kicked from the center of the SNR.H0 in both cases is that the magnetar did not originate at the center of the SNR. 1 is derived similarly to the IR case; when the angular size of the SNR is known, we assume the semi-major ( SNR /2) and semi-minor axes ( SNR /2) are the 3 half-widths of the error ellipse.If not, the error ellipse is used to derive the binormal distribution.We again derive the search statistic ĝ1 using the 95 th percentile contour, and the p-value  1 by comparing to the H0 distribution of ĝ1 , which is derived using mean angular separations between 1 − 5 /  of the origin.
The second criteria is similar to that above, but instead compares the angle of the line connecting the magnetar location and the SNR,  SNR = tan −1 (Δ/Δ), and the direction of the magnetar's proper motion,  mag = tan −1 (  /(  cos())).The H0 distribution of ĝ2 is computed using mean angular differences Δ = | SNR −  mag | between the 95 th percentile of the observed distribution and 180 • . 3 uses the magnetar proper motion and the median magnetic field decay age T decay computed from the X-ray luminosity (see Appendix D1 and Ferrario & Wickramasinghe 2008) to derive samples of the expected angular separation between magnetar and SNR, Δ mag = | − →  |/T decay .If the magnetar has no X-ray luminosity or magnetic field estimate, the spin-down age T rot from the McGill magnetar catalog is used instead.Samples of the true angular separation, Δ SNR are derived using samples of the RA and declination of the magnetar and SNR.The search statistic ĝ2 is the fraction of samples for which the difference in angular separation, Δ(Δ) = |Δ SNR − Δ mag | is outside the 95 th percentile from the H1 case. 3 is then derived from the H0 distribution of ĝ3 , which assumes mean angular separations between the 95 th percentile and 1 ′ .
Finally, we estimate  4 based on the angular size Δ of the SNR; we impose a lower limit on the angular size required to be associated with the magnetar given the magnetar's distance and age T decay/rot .By sampling the magnetar distance  we compute samples of the physical radius of the SNR as  SNR =  mag Δ and estimate the corresponding dynamical age T dyn from either free-expansion or Sedov-Taylor phase expansion as described in Appendix D2 (Sedov 2018).For this we also sample the progenitor envelope mass from a uniform distribution between 2.75 − 120 ⊙47 and use an explosion energy 10 51 ergs48 (Morozova et al. 2017).We then draw samples of the magnetar age T decay/rot and compute the statistic ĝ3 as the fraction of samples for which T decay/rot > T dyn .We impose only a loose lower limit since magnetar spindown and decay ages are uncertain.The p-value  3 is taken as the percentile of ĝ3 within the H0 distribution which assumes mean ages between T dyn and the 95  ℎ percentile of the H1 distribution. 1 ,  2 , and  3 are combined using Equations 5 and 6.

B3.1 'By-Eye' Identification and Fitting of Candidates
Optical, infrared, and radio sources identified 'by eye' which do not have known counterparts in public catalogs considered here were fit with ellipses parameterized by the RA and declination of the center (  ,   ), major and minor axes (, ), and position angle measured North of East ().Shifting the center of the ellipse to 0, the fit equation for the radius  as a function of polar angle  is given by: (B5) () = (( 2 cos 2 () +  2 sin 2 ())cos 2 () + ( 2 sin 2 () +  2 cos 2 ())sin 2 () + 2( 2 −  2 )cos()sin()cos()sin()) −1/2 For each source, the contour at one-half the maximum value within one arcminute (for radio sources; one arcsecond for optical and IR sources) of the identified center.The mean  and  were recomputed from this contour prior to fitting with an ellipse using the equation above.The major and minor axes and position angle were then used identically to those of SNRs for the association test (see Section B3.Note that apparent magnitudes were not estimated for these sources, and p-values could not be estimated for the optical and IR searches.Any point sources of interest were queried in SIMBAD and public catalogs to determine the likelihood of association.

B4 Gaia Search for Unbound Companions
Criteria for association of unbound Gaia sources with magnetars is based on their current positions, proper motion magnitudes and directions, and the consistency of the magnetar travel time with age estimates.H0 is that the magnetar and source positions are not consistent with having an origin at the same location.For  1 , we assume the Gaia DR3 proper motion, parallaxes, and positions are Gaussian (or 2-sided Gaussian, where applicable) distributed, with standard deviation equal to the uncertainty.To account for bias in the conversion from parallax to distance ( ≈  −1 ) we adopt the photo-geometric distance from the Bailer-Jones (2023) catalog.
For each magnetar, samples of the distance , RA () and declination () are drawn from their respective distributions to get samples of the position vector, − →  mag,i (the subscript  indicates the quantity is Monte Carlo sampled).Magnetar proper motion measurements are often unavailable due to their reliance on Very Long Baseline Interferometry (VLBI) (e.g Tendulkar et al. 2012).In such cases, we draw samples of the magnetar velocity magnitude from a Maxwellian distribution with scale parameter  = 265 km s −1 as described in Hobbs et al. (2005).The velocity direction is drawn from a 3D uniform distribution to obtain − →  mag,i .If proper motions are known, the total velocity is still drawn from the truncated Hobbs et al. (2005) distribution, and the radial velocities are then computed from the velocity and proper motion samples with the sign chosen from a uniform distribution.Samples of the age   , defined here as the time since the magnetar supernova would have occurred, are drawn from a uniform distribution from 1 − 10 kyr.A sub-sample of Gaia sources within 5 • of the magnetar-under-test are selected from the full sample to test for association.In a similar method to that above, samples of position − →   and velocity − →   are drawn for each Gaia source; here we sample the empirical walk/runaway velocity distribution from the fiducial simulation of Renzo et al. (2019) for which ∼ 95% of simulated ejected massive stars have velocities ≲ 30 km s −1 (see also Blaauw 1961;Tetzlaff et al. 2010;Sayers et al. 1996;Feast & Shuttleworth 1965).We choose not to utilize the Gaia radial velocity measurements since they have not been confirmed through spectroscopic follow-up, and many (∼ 10%) of sources within our sample do not have radial velocity measurements.This also prevents us from over-constraining the velocity samples while we impose the Renzo et al. ( 2019) distribution.
Next, the magnetar and source trajectories are traced for each sample from their current positions − →  mag,i and − →   to their locations at time   in the past: We have utilized the astropy library's apply_space_coord method to perform this propagation, and use the astropy.SkyCoord object for conversion between Cartesian and Spherical coordinate systems.
The initial RA and declination of source and magnetar,  , , , and  ,mag, , ,mag, , are derived from − →  , and − →  o,mag,i , and are used to compute the angular separation at the time of the supernova, Θ  .We compute this using SkyCoord.angular_separationmethod, which implements the Vincenty formula for the Great-circle distance between two points on a spherical surface (Vincenty 1975).Additionally, using the G-band magnitude for each source, the  − color derived through Johnson-Cousins relations, and the distance samples which are converted to distance moduli, we also derive samples of the V-band absolute magnitude   , : where ΔΦ is the angular separation between the intersection point and the current coordinate.Since we expect  2D,mag, =  2D, , we estimate a  2 statistic comparing the  2D,mag, vs.  2D, curve to a line with slope 1, taking the p-value  2 as another measure of association.Similarly, we compute the physical travel time for the magnetar assuming it is at the Gaia source distance: We expect each magnetar's travel time to be less than its age; if an SNR candidate with an age measurement is confidently associated with the magnetar, we use the SNR age, T SNR .Five magnetars in our sample have T SNR estimates, four of which are measured assuming Sedov-Taylor plasma dynamics and the other which is computed kinematically using the magnetar proper motion (Sedov 2018; Dohm- rate (TPR) for trial thresholds from 0.1 to 1, 10000 injection sources were simulated with randomly chosen ages between 1-10 kyr and velocity drawn from the fiducial simulation described in Renzo et al. (2019).A test magnetar was simulated with velocity drawn from the Hobbs et al. (2005), and both magnetar and source were propagated to present-day positions, starting from the same start position ( 0 =  0,mag ).Each simulated source was given an OB-star-like absolute magnitude   = −1.5.The Monte Carlo search was conducted with 1000 trials for each simulated source and the simulated magnetar.Search statistics f were computed for each source.The TPR for each trial f was computed as the fraction of sources with f < f .The false negative rate (FNR) is computed as FNR = 1-TPR.
To compute the false positive rate (FPR) for a given magnetar location, we tile the region nearby with 10 ′ pointings, starting at 10 ′ from the magnetar.For each tile, we pass all sources through the Monte Carlo association test, using the tile's center as the coordinates of the 'fake' magnetar.Since no magnetar is in the query region, these test samples are used as true negative injections.The FPR for each trial f was computed as the fraction of total sources (from all tiles) with f < f .
We expect differences in each magnetar's local environment, both in source density and, potentially, stellar types.Therefore, the ROC analysis is carried out for each magnetar individually.Note that in most cases, the FPR is very low; the FPR is higher for SGR 1801-23 because of its lack of distance measure and large position error.We choose the threshold for each magnetar where the TPR is maximized and the FPR is minimized; in cases where the FPR is consistently 0%, the threshold is f = 100%.In these cases, a Gaia source is a candidate for association if any Monte Carlo sample falls within the 90 th percentile contour in {Θ,   } space.Similarly, in cases where the FNR is consistently 0% (which occurs when the proper motion, position, and distance are well-constrained), the threshold is f = 0%.In these cases, a Gaia source is a candidate for association only if all Monte Carlo samples fall within the 90 th percentile contour in {Θ,   } space.Overall, we find this has little effect on the results, as most Gaia sources have either all or none of their samples within the contour region.We list f for each magnetar in Table B1.

APPENDIX C: FALSE POSITIVE ANALYSES C1 Gaia Search False Positives
In this section we discuss in detail the reasons that the candidates for unbound stellar companions of 3XMM J185246.6+003317 were rejected as false positives.
C1.1 3XMM J185246.6+003317−10 km s −1 , respectively.Without uncertainties on the distance of 3XMM J185246.6+00317, it is difficult to evaluate if each source is consistent with the magnetar's distance.However, the magnetar's estimated distance ∼ 7.1 kpc far exceeds the four Sources' parallax distances, which are between 1 − 3 kpc.3XMM J185246.6+00317'sdistance was estimated from the nearby SNR Kes 79 based on their similar hydrogen column densities   ≈ 1.5 × 10 22 cm −2 (Zhou et al. 2014).This is questionable given that Kes 79 is confirmed in association with a separate X-ray pulsar, CXOU 185238.6+004020(Seward et al. 2003), although we explore a scenario in which both the pulsar and magnetar may be associated with the SNR in Section 3. It is also known that some X-ray sources have significant local extinction that can bias the foreground   estimate (e.g.Chaty & Rahoui 2012;Predehl & Schmitt 1995).However, the lack of notable IR emission from 3XMM J185246.6+003317distinguish it from sources with significant local   such as HMXBs.Therefore, it is likely that   is indeed from the foreground and its similarity to that of Kes 79 suggest they are at the same distance ∼ 7.1 kpc.If this is the case, then all five sources are in the foreground and cannot be associated with 3XMM J185246.6+003317.A precise proper motion measurement and independent distance measurement would allow more confident conclusions to be drawn.

C2 Radio SNR Search False Positives
In this section we justify why candidates for SNR association with SGR 0755-2933 were rejected as false positives.

C2.1 SGR 0755-2933
Eight radio sources were identified with  < 5% as potential SNRs associated with SGR 0755-2933 Richardson et al. (2023) and Doroshenko et al. (2021) (see also Surnis et al. 2016;Archibald et al. 2013;Barthelmy et al. 2016).VLASS image cutouts with 1 ′′ pixels were obtained for each candidate, and four are ruled out as point sources.Two other sources are found to be coincident with known stellar objects within 1 ′ (an eclipsing binary and a Carbon star).The remaining two sources are each ∼ 14 ′′ across, or ∼ 0.22 pc at SGR 0755-2933's distance of 3.5 ± 2 kpc.Assuming the free-expansion phase, this corresponds to an SNR age range of 80 − 300 years, which is marginally consistent with magnetar ages in general.If we accept the Doroshenko et al. (2021) association of SGR 0755-2933 with the X-ray source 2SXPS J075542.5-293353,we can use the X-ray luminosity   ≈ 10 34 erg s −1 to estimate T decay for a typical magnetic field strength  ∼ 10 14 G (see Appendix D1).From this method we estimate T decay ≈ 60 kyr which is much larger than the free-expansion ages for the two SNR candidates.The detection of a subsequent burst resulting in more precise localization could confidently rule out one source which lies outside the 3 error region of SGR 0755-2933.With the information at hand, we conclude that both sources are likely to be false positives based on the large position uncertainty and decay age estimate.

C3 Bound Companion Search False Positives
In this section we justify why the bound companion candidates for PSR J1622-4950 and SGR J1822.3-1606 were rejected as false positives.

C3.1 PSR J1622-4950
Through the optical search, the SkyMapper source SMSS J162244.99-495055.7 was recovered with p-value  = 1% as a potential bound companion of PSR J1622-4950.This magnetar lies along a highly extincted sightline, with   ≈ 30.17 +8.94  −7.82 from the Predehl & Schmitt (1995)   −   relation.From the − and −band apparent magnitudes, we therefore estimate an absolute magnitude   ∼ −23.4,which would be extremely bright for an OB star.From the mass-luminosity relation, this would correspond to a ∼ 10 13  ⊙ star, which is unphysical.We also note that the candidate has a large angular separation ∼ 3.31 ′′ , corresponding to 0.14 pc at the magnetar's 9 ± 1.35 kpc distance.This is impractical for the orbit of a neutron star (∼ 1.4 ⊙ ) and massive star (≲ 120 ⊙ ), suggesting an orbital period of ∼ 10 18 years.Therefore, we conclude that this source is most likely a false positive.Chrimes et al. (2022a) identify a potential infrared counterpart for PSR J1622-4950 using HST data, but this is not recovered in either the optical or IR search in this work.This is consistent with their claim in Chrimes et al. (2022b) that this is more likely to be surface emission from the neutron star than a companion star.
The second source, 2MASS 18221794-1604259, is much closer to the magnetar at ∼ 0.01 ′′ ± 0.7 ′′ , or ∼ 20 au at the magnetar's distance ∼ 1.6 ± 0.3 kpc, though this is still too distant for a bound system.This source was similarly identified as evidence that SGR J1822.3-1606 may be a Be star X-ray binary, like SGR 0755-2933, when initially discovered by e.g.Halpern (2011), Gogus & Gorosabel (2011), Gorosabel et al. (2011).The source's Gaia counterpart, Gaia DR3 4097826476059022848, has a negative parallax  = −1.2± 0.6 mas which excluded it from the Gaia search (see Section 3 and Appendix A).Its parallax distance is  ≈ 7.0 +3.1 −2.8 kpc, which agrees with the magnetar's distance within 1.At the magnetar's distance, the intrinsic magnitude and color are   ≈ 2.5 ± 0.1 and  −  ≈ 1.5 ± 0.1, which are stellar in nature.However, the color is redder than expected for OB stars, and is more characteristic of an evolved supergiant star (Chrimes et al. 2022b).Notably, the star's proper motion from Gaia is   cos() = −5.6 ± 0.7 mas yr −1 ,   = −8.2± 0.6 mas yr.Applying the Gaia Monte Carlo test results in a p-value  = 79% due to the significant travel time difference Δ 2D ≈ 0.3 +0.7  −0.2 kyr required to trace trajectories to a common origin.While a magnetar proper motion measurement would motivate re-analysis, we conclude that the close proximity and companion's significant proper motion make it unlikely that the candidate is associated with SGR 1822.3-1606.

APPENDIX D: MAGNETAR AGE ESTIMATION
Magnetars' magnetic field decay and irregular outbursting and 'glitching' behavior make their characteristic spindown ages T rot = (2 ) −1 overestimates of their true age (Nakano et al. 2015;Kaspi & Beloborodov 2017;Olausen & Kaspi 2014;Ferrario & Wickramasinghe 2008).In this section we provide a brief overview of two alternative estimates based on the X-ray luminosity and the dynamical plasma age of an associated SNR.

D1 Derivation of Magnetar Age X-ray Luminosity Decay
We define T decay as the age of the magnetar estimated from the decaying crustal X-ray luminosity,   .This model is described in detail in Ferrario & Wickramasinghe (2008) and Ferrario & Wickramasinghe (2006) based on Pons & Geppert (2007).In this model, the dipole magnetic field component is responsible for spindown energy loss, but does not decay on timescales below ∼ 100 Myr.The toroidal component decay on timescales ∼ 1 − 10 kyr drives energy loss through quiescient crustal X-ray emission   .We assume the luminosity is modelled by an exponential decay following the initial cooling phase: where   is the dipole field strength,  0 = 10 33 erg s −1 , and   is the decay timescale.The decay timescale is inversely proportional to the magnetic field, and is modelled as a power law dependence with  = 1.3: where  0 = 500 kyr.Plugging this into Equation D2 and solving for  = T decay , which we take to be the present-day: where   (T decay ) is the X-ray luminosity at present-day.The constants  0 ,  0 , and  were determined through simulations by Ferrario & Wickramasinghe (2008).Figure D1 shows each magnetar's T decay as a function of   and   .

D2 Derivation of Free-Expansion and Sedov-Taylor Expansion SNR Ages
In this section, we provide a brief description of the free-expansion and Sedov-Taylor expansion phases for SNRs as are relevant to this paper.We do not offer a detailed derivation, but refer the reader to Sedov (2018, Chapters IV. and V.), Truelove & McKee (1999), and Suzuki et al. (2021) from which this section is summarized.See also Taylor (1950) and Sedov (1946) for additional details.In this work, we are concerned with the age-radius relations for SNRs, and focus the discussion on this.
The earliest stage of a supernova explosion, which we refer to as 'free expansion', is dominated by the momentum of the ejecta mass,  ej , rather than mass swept up from the surrounding medium (Truelove & McKee 1999).This, coupled with the negligible radiative energy losses during this phase, allow the ejecta from the supernova to expand freely into the surrounding medium.During this phase, the radius of the expanding ejecta front is approximately linear with time: FE 1 km ≈ (7090 km/s) √︄  10 51 ergs where  is the energy of the supernova and T is the age of the remnant.The expansion velocity is lower than the sound speed in the surrounding medium, causing a 'blast wave' shock to precede the ejecta shell and inducing a reverse shock to oppose the ejecta.This marks the start of the transition to the Sedov-Taylor expansion phase, and occurs roughly on the characteristic timescale: T ch =  5/6  1/2 ( ,0 ) 1/3 (D5) where  ≈ 1.4  is the reduced mass of the ejecta and  ,0 is the initial electron density, which increases with time through interaction with the ambient medium.The corresponding characteristic radius is  ch =   ,0 1/3 (D6) (Truelove & McKee 1999;Sedov 2018).
Beyond this timescale and radius, the system enters the Sedov-Taylor expansion phase, in which the expansion is damped by the ambient medium.On late timescales, the ejecta expands as a power law with  ∝ T 2/5 .To first order, we can form a piecewise model of (T ) by joining the free expansion and Sedov-Taylor solutions at  ch (T ch ).Imposing this condition, we can derive the coefficient for the Sedov-Taylor solution: which we note is independent of mass (Truelove & McKee 1999;Sedov 2018).We use this age estimate as an upper limit on the required age of a candidate SNR when deriving  3 in the radio SNR search (see Section B3) and when considering the associations of SNRs in Section 5.3.We do not use this approximation for the SNR age estimates in Table 2, and opt only to use  SNR for those with age estimates in the literature.

Figure 1 .
Figure 1.(Top) Locations of Gaia search sample in Celestial coordinates.The tangential velocity magnitude in km s −1 is indicated by the colorbar.The Bayestar19 Dust Map is used for sources at  > −30 • (red outline), and the Gaia photometric extinction is used for  < −30 • (black outline; Vallenari et al. 2022; Green 2019).The 5 • radius cones queried around each magnetar are shown in red, and those around each pulsar used for the test sample (see Section 3.2) are shown in violet.Lower transverse velocities are found in the Galactic Plane (GP), meaning that the unbound companions of GP magnetars can be found at smaller angular separations.(Bottom) Color-Magnitude Diagrams for the Gaia search sample.The left plot shows the apparent −band magnitude   plotted against the apparent    −    color.The right plot shows absolute  −band magnitude   for each source against the absolute color   −   .  is derived using the Johnson-Cousins relationships in Equations A1-A8 (see Appendix A).Sources with  > −30 • are again outlined in black while those with  ≥ −30 • are outlined in red.The parallax of each source is indicated by the colorbar in milliarcseconds (mas).Sources from the initial query (see Appendix A) which were cut by the magnitude (black) and color (red) limits are shown in grey.GSP-Aeneas extinction estimates may be unreliable in the Galactic plane according to the Gaia DR3 Documentation (Section 11.3.3),resulting in a rough color cutoff near  BP −  RP ≈ −0.5(Vallenari et al. 2022;Andrae et al. 2022;Fouesneau et al. 2022).

Figure 2 .
Figure 2. (Left) Monte Carlo sampled trajectories for candidate Gaia DR3 sources with p-values  < 5% possibly associated with 3 magnetars (right, bottom), PSR J0538+2817 (top left), and PSR B1951+32 (middle left).Each figure includes only Gaia OB star candidates with association p-value  < 5%, including the unbound companion of CXOU J164710.2-455216 proposed by Clark et al. (2014), WR 77F (top right).Present-day magnetar/pulsar and GDR3 source locations are shown as red stars and circular markers, respectively; the latter has colorbar corresponding to the total p-value , except for PSR 0538+2817 and CXOU J164710.2-455216, for which we exclude  1 (see Section 3).Red lines show the initial magnetar/pulsar trajectory samples while thin lines extrapolate the magnetar/pulsar and source trajectories to their past intersection.Note for these initial trajectories extend farther for pulsars (top left, middle left) than for magnetars (right, bottom) since travel time samples were drawn from a uniform distribution up to 1 Myr for pulsars in contrast to 10 kyr for magnetars.If proper motion measurements are available, the proper motion direction is indicated by a purple arrow.Each GDR3 source is shown with its source ID, 2D travel time difference Δ 2D , the magnetar's/pulsar's 3D travel time T 3D,mag/psr , an upper limit age estimate T SNR/decay/rot , and the relevant p-value.(Right) Contour plot of Monte Carlo samples in {Θ −   } space.The 95% contour region is in beige, and candidates with search statistic f less than the critical threshold, f , have colorbar corresponding to f (see Appendix B4).f is indicated by a red horizontal line on each colorbar.We recover the proposed companions of CXOU J164710.2-455216(WR 77F; upper right) and PSR J0538+2817 (HD 37424; top left) with  < 5%.We also identify potential runaway companions of SGR 1822.3-1606 and PSR B1951+32 (middle row) with  < 5%.The companions shown for 3XMM J185246.6+003317(bottom) are most likely false positives as discussed in Appendix C1.

Figure 3 .
Figure 3. Optical (Left) and Infrared (Right) 1 ′ × 1 ′ image cutouts for 2 magnetars and 1 magnetar candidate.Each image is in grayscale where RGB channels are combined using 0.2989 + 0.5870 + 0.1140.Optical images use the −, −, and −band filters for PS1 images (CXOU J1714) and the −,  −, and −band filters for SkyMapper images (SGR 0755 and SGR J1822).Infrared images use the  −,  −, and  −band filters for UKIDSS images (SGR 1822) and the  −,  −, and   −band filters for VVV (CXOU J1714) and 2MASS (SGR 0755) images.The magnetar location is indicated by red cursors and the 3 error ellipse is shown in red.Note the error ellipse is much larger than the image for SGR 0755-2933.Each marker indicates a point source from the optical or infrared catalogs with colorbar corresponding to the association p-value, .Sources with insufficient data to determine a p-value are grey.Sources with  < 5% are outlined in green, while those with  > 5% are outlined in red.(Left) PS1 sources are shown as circles, SkyMapper sources are shown as triangles, and sources identified 'by-eye' are shown as squares.(Right) 2MASS sources are shown as circles, VVV sources are shown as triangles, UKIDSS sources are shown as diamonds, and sources idenified 'by-eye' are shown as squares.Sources with no corresponding marker were not identified in the catalogs considered, and could not be well-fit with an ellipse.The Be-star companion of SGR 0755 is recovered in the optical and infrared searches, and candidates are identified for SGR 1822.The proposed companion of CXOU J1714 is not recovered, but cannot be ruled out as a bound companion due to uncertainty in the extinction estimate.Note that the axes for CXOU J1714 differ between optical and infrared due to the use of different World Coordinate System (WCS) projections (ZPN) for VVV data (e.g.Greisen & Calabretta 2002).MNRAS 000, 1-34 (2024)

Figure 4 .Figure 5 .
Figure 4. Radio Image and Trajectory Analysis of 3XMM J185246.6+00317.(Top left) A RACS DR1 888 MHz continuum 1• × 1 • image cutout is shown for 3XMM J185246.6+00317 with 25 ′′ × 25 ′′ pixels.The magnetar and pulsar locations are indicated by red and blue pointers, while circular and square markers show known and 'by-eye' identified radio sources, respectively.Note that no markers are shown for catalogued Active Galactic Nuclei or radio galaxies, as they were not included in this search.Other extended sources without markers were not identified in SIMBAD and could not be well-fit with an ellipse.The colorbar corresponds to the association p-value.The purple shaded region indicates the 1, 2, and 3 past travel distances based on theHobbs et al. (2005) velocity distribution.The angular extent of Kes 79 (labelled) is shown as a blue circle.(Bottom) The cartoon demonstrates how a bound massive stellar binary (orange) may evolve into an unbound magnetar (red) and pulsar (blue) following two CCSNe.(Top right) The trajectory plot for 3XMM J185246.6+00317 and the X-ray pulsar CXOU J185238.6+004020 is shown with 468 valid trajectories (out of 10000 trials).Blue traces show the pulsar's trajectory prior to the second CCSN, and yellow traces show its trajectory after the CCSN.The blue circle indicates the location and extent of Kes 79.

Figure 6 .
Figure 6.Completeness plots versus distance modulus for (rows 1 and 2) the Gaia search, (row 3) the optical and (row 4) IR searches, and (row 5) the radio search.Circular, square, and diamond markers with errorbars indicate known companions or remnants of magnetars, pulsars, and HMXBs, respectively.Proper motion, apparent magnitude, and 1 GHz flux ranges, are shown as vertical lines, with color relating to the survey used for the search.The sensitivity limits of each survey are indicated by horizontal lines.We expect most magnetars with SNR shells or stellar companions would have been detected if present.Sensitivity limits for past searches are shown as black dotted lines for comparison (see Appendix E for references).
Figure 7. (Upper Left:) Corner Plots from the Magnetar Binary Fraction Simulation in Section 6.2.The 'SNR complete', 'SNR incomplete', and 'SNR/unbound complete' cases as defined in Section 6 and Table 4 are shown in black, red, and blue, respectively.Vertical lines indicate the median values from Table4in each case.The 'SNR complete' and 'SNR/unbound complete' cases would suggest that   ≈ 45 − 55% of magnetars to descend from non-CCSN channels, whereas the 'SNR incomplete' case suggests a higher stellar merger rate   ≈ 70% among OB stars is sufficient to produce the observed sample.The merger rates displayed assume a massive star binary fraction  0 = 84 ± 9%(Moe & Di Stefano 2017).(Upper Right:) 90% Confidence Regions for the merger fraction   and OB star binary fraction  0 .The colors for each region correspond to the cases in the contour plot, and vertical purple lines show  0 predicted for different stellar types through population synthesis(Moe & Di Stefano 2017).The green and yellow dashed lines indicate the median merger fractions from Renzo et al. (2019) and Kochanek et al. (2019); the teal dashed line indicates the   vs.  0 curve inferred from   and   in Chrimes et al. (2022b).The green and teal shaded regions indicate the 1 errors for the Renzo et al. (2019) and Chrimes et al. (2022b) cases.Each of these appear to underestimate the true fraction for early B-type stars.(Bottom:) Flowchart showing the progression of massive stars into magnetars with constraints from population synthesis and observations in this work.CCSNe are indicated by orange arrows while alternate formation channels are indicated by blue arrows.Constraints on   ,   ,   ,   , and  0 (see Section 3.3) are shown from simulations by Moe & Di Stefano (2017, purple), Renzo et al. (2019, green), Kochanek et al. (2019, yellow), and Chrimes et al.(2022b, teal).Observational constraints from this work are shown in red for the 'SNR incomplete' case and blue for the 'SNR/unbound complete' case.A lower binary fraction would allow for a lower merger fraction for magnetar progenitors, but would require lower mass progenitors.MNRAS 000, 1-34 (2024)

Table 1 .
Magnetars and Magnetar Candidates from the McGill Catalog (revised position errors are boldfaced for clarity) Predehl & Schmitt (1995)es are given with full names in parentheses, which are used throughout this work.RA (), declination ( ), and distance are taken from the McGill Magnetar Catalog(Olausen & Kaspi 2014), while proper motions (  cos( ) and   ) were found in the available literature.Italicized entries are unconfirmed magnetar candidates and boldfaced entries have confirmed or proposed supernova remnant associations.Proposed associations are indicated by a * .extinctionestimates are from the Bayestar19 map for magnetars with  ≳ −30 • .Magnetars outside this range have   from the X-ray column density   (indicated by a ‡ ) using thePredehl & Schmitt (1995)relation:   ≈ (5.59 × 10 −22 cm 2 )  .If no   measurement is available, we list the median   from Gaia GSP photometry for massive sources within 1

2 Radio SNR Search Test with Known Magnetar-SNR Associations 8
Corbel et al. 1999)en in the B configuration which is sensitive to angular scales between 2.1−58 ′′ .NVSS data is taken in the D and DnC configurations which are sensitive to angular scales between 14 − 970 ′′ and 46 − 970 ′′ , respectively.See https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/resolutionfor deails.magnetars and 2 magnetar candidates with confirmed SNR associations are used to test the radio image search method.With the analysis above, we accurately recover the associated SNRs with  < 5% for 7 magnetars and both magnetar candidates.SGR 1627-41 is the only confirmed SNR association that we fail to recover in the radio search (SNR G337.0-00.1;Corbeletal. 1999), which is likely due to the large angular offset ∼ 15.5 ′ .No proper motion is known for SGR 1627-41, but a velocity ∼ 1275 km s −1 would be required given its T decay ≈ 37.0 kyr, which is in marginal agreement with the Hobbs et al. ( 32While this magnetar does fall in the survey range of RACS-DR1, no continuum images are available that contain SGR 1627-41 in its field-of-view. 33tps://simbad.u-strasbg.fr/simbad/sim-fbasic 34e https://simbad.cds.unistra.fr/guide/otypes.htxfordetails on SIMBAD object-types.35Theminimum p-value is set by the 1000 trials used for the SNR search. 36https://en.wikipedia.org/wiki/Fisher%27s_method 5. Zhou et al. (2014)erive a Sedov-Taylor expansion age of ∼ 5 kyr for the SNR which contrasts with both the decay age T decay ≈ 0.7 Myr 37 and spindown age T rot ≳ 1.3 Myr of 3XMM J185246.6+003317.38Magnetarsaretypicallyexpectedtobe unobservable beyond ∼ 10 4 years; however,the spindown age T rot is a well-constrained lower limit derived byRea et al. (2014)from an X-ray timing solution.They find a pulse period  = 11.55871346(6)s with no discernible spindown, leading to upper limits of  < 1.4 × 10 −13 s s −1 and  dip ≲ 10 13 G. on the spindown rate and magnetic field.While these paramters are reminiscent of typical radio pulsars,Zhou et al. (2014)noted in their (Gotthelf & Halpern 2008;Halpern et al. 2007;Gotthelf et al. 2013simulations byde Lima et al. (2024)suggest T decay ≳ 1.2 Myr, though they caution not to over-interpret their initial estimates.38While thpindown age of CXOU J185238.6+004020 is ∼ 192 Myr, 'antimagnetars' or 'compact central objects' (CCO) like this are known to have overestimated spindown ages due to the complex relation of their low magnetic fields to a dynamo process(Gotthelf & Halpern 2008;Halpern et al. 2007;Gotthelf et al. 2013).

Table 3 .
Search Results and Completeness Ratios

Table 4 .
Constraints on magnetar binaries from this work (Moe & Di Stefano 2017)es based on the completeness of the RACS SNR search: 'SNR incomplete' assumes that all magnetars originate in CCSNe, whether one is detected or not.'SNRcomplete'assumes that magnetars with RACS search completeness ratios > 75% which had no SNR detected did not form through CCSNe.We also consider the 'SNR/unbound complete' case in which unbound magnetar binaries are evidence of the CCSN progenitor channel as well.Median values are reported with 1 errors.†,,   , and   are defined such that   +   +   = 1 in the 'SNR incomplete' case, and   +   +   +   = 1 in the 'SNR complete' and 'SNR/unbound complete' cases.andsum to the total binary fraction   +   =   .The merger rate   is the fraction of binary stars that merge, while the other fractions are taken out of the total number of CCSNe.‡   is estimated assuming the progenitor stars are early B-type stars (9 ⊙ <  * < 16 ⊙ ) with binary fraction  0 = 84 ± 9%(Moe & Di Stefano 2017).

CCSN Non-merged Binaries Non-Merger Merged Stars CCSN Alternate Channel Alternate Channel Alternate Channel Alternate Channel
Corner Plots from the Magnetar Binary Fraction Simulation in Section 6.2.The 'SNR complete', 'SNR incomplete', and 'SNR/unbound complete' cases as defined in Section 6 and Table4are shown in black, red, and blue, respectively.Vertical lines indicate the median values from Table4in each case.The 'SNR complete' and 'SNR/unbound complete' cases would suggest that   ≈ 45 − 55% of magnetars to descend from non-CCSN channels, whereas the 'SNR incomplete' case suggests a higher stellar merger rate   ≈ 70% among OB stars is sufficient to produce the observed sample.
+ .1 ≈ $* ± &% B8)Θ  and   , are used as the parameters of interest and compared to a threshold contour in {Θ,   } space.The contour is computed by injecting 100 positively-associated simulated sources with an OBstar-like absolute magnitude, −10 ≤   ≤ +3.The 90 th percentile contour of the 2D histogram is used as a threshold, and we define the search statistic f as the fraction of Monte Carlo samples falling outside the contour.A separate threshold contour is determined for each magnetar, using the position, distance, and proper motion of the magnetar as available.Critical values f were derived through a Receiver Operating Characteristic (ROC) analysis for each magnetar, which is summarized in Appendix B4.1.By comparing the final f to samples in the H1 hypothesis (obtained during the ROC test), we estimate a p-value,  1 , as the percentile of f within the H1 sample.To further narrow down the candidate Gaia sources with f < f , we iteratively compute the true intersection points − →  , of the sampled trajectories with those of the magnetars, eliminating any samples with non-physical solutions.We then use the Monte Carlo sampled proper motions and positions to estimate the 2-dimensional travel times for each magnetar-candidate pair:

Table B1 .
Critical Thresholds for Magnetar Association Tests