Corecollapse 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 sim 45 rm per cent of magnetars should be nearby highvelocity stars. Here we conduct a multiwavelength search for unbound stars magnetar binaries and SNR shells using public optical uvgrizy bands infrared J H K and KSUBsSUB bands and radio 888 MHz 1.4 GHz and 3 GHz catalogues. 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 catalogue we associate with SGR J1822.31606 and an Xray 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 5lesssim fu lesssim 24 rm per cent 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 31lesssim fnc lesssim 66 rm per cent of magnetars to descend from such channels. Our results support a high fraction 48lesssim fm lesssim 86 rm per cent of preCCSN mergers which can amplify fossil magnetic fields to preferentially form magnetars.
['Astrophysics - High Energy Astrophysical Phenomena']
Searching for magnetar binaries disrupted by corecollapse supernovae
{'Searching for Magnetar Binaries Disrupted by Core-Collapse Supernovae': 'Myles B. Sherman, 1 ★ Vikram Ravi, 1 Kareem El-Badry, 1 Kritti Sharma, 1 Stella Koch Ocker, 1 , 2 Nikita Kosogorov, 1 Liam Connor, 1 Jakob T. Faber 1 1 Cahill Center for Astronomy and Astrophysics, MC 249-17 California Institute of Technology, Pasadena CA 91125, USA.  \n2 The Observatories of the Carnegie Institution of Washington, Pasadena, CA 91101, USA.  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': '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 ∼ 45% 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 ( 𝑢𝑣𝑔𝑟𝑖𝑧𝑦 -bands), infrared ( 𝐽 -, 𝐻 -, 𝐾 -, and 𝐾 𝑠 -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 3XMMJ185246.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 ≲ 𝑓 𝑢 ≲ 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 ≲ 𝑓 𝑛𝑐 ≲ 66% of magnetars to descend from such channels. Our results support a high fraction (48 ≲ 𝑓 𝑚 ≲ 86%) of pre-CCSN mergers, which can amplify fossil magnetic fields to preferentially form magnetars.  \nKey words: stars: magnetars, stars: neutron, binaries: general, stars: massive, ISM: supernova remnants, proper motions', '1 INTRODUCTION': "Magnetars are highly magnetized ( 𝐵 ≳ 10 14 G) neutron stars which emit X-ray, gamma ray, or radio emission driven by their magnetic fields (Kaspi & Beloborodov 2017). At the time of submission, there are 24 confirmed magnetars and 7 magnetar candidates identified based on their high-energy X-ray and gamma ray bursts (Richardson et al. 2023; Cline et al. 2000; Lamb et al. 2003; Mereghetti et al. 2012; Torii et al. 1998a; Sakamoto et al. 2011a; Leahy & Tian 2007) 1 . With burst events initially categorized as either Anomalous X-ray Pulsars (AXPs) or Soft Gamma Repeaters (SGRs), magnetars were later identified as the common source of their emission (Gavriil et al. 2002).  \nIt 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;  \n- ★ E-mail: msherman@caltech.edu  \n1 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.  \nBeniamini 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.  \nMagnetars, 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).  \nMagnetar birth rates 𝛽 ≈ 𝑁 mag T -1 ∼ 0 . 3 century -1 , are dif-  \nconstrain 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 emission 2 . 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).  \nWhile 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.  \nSNR radio emission surrounding a magnetar is the most direct indicator of a CCSN progenitor (e.g., Kaspi 1996). Supernova shells are detectable in radio and X-rays expanding from a young (1 -10 kyr) magnetar's birth site at its center; most of the known magnetar-SNR associations were identified by their remnant shell radio emission (e.g. Gelfand & Gaensler 2007; Anderson et al. 2012a). The significant increase in data from radio all-sky surveys will boost the sensitivity considerably. The Rapid ASKAP Continuum Survey (RACS), for example, uses the 36-antenna Australian Square Kilometre Array Pathfinder (ASKAP) telescope to survey the southern sky (declination range -80 · < 𝛿 < +30 · ) from 700-1800 MHz (McConnell et al. 2020). The Karl G. Jansky Very Large Array (VLA) Sky Survey (VLASS) and National Radio Astronomy Observatory (NRAO) VLA Sky Survey (NVSS) similarly cover the northern sky above 𝛿 > -40 · at 1.4 and 3 GHz (Gordon et al. 2020, 2021, 2023a; Condon et al. 1998), while the Low Frequency Array (LOFAR) covers lower frequencies near 135 MHz (Shimwell et al. 2022; Shimwell 2016; Condon et al. 1998; Gordon et al. 2020, 2021, 2023a). Identifying supernova shells through radio emission relies on the magnetar  \nbeing young enough that the shell has not yet dissipated ; for expansion rates up to ∼ 1000 km s -1 , SNRs with ages ≲ 100 kyr may be detected (e.g., Anderson et al. 2012a; Truelove & McKee 1999; Sedov 2018).  \nIn 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-collapse 3 ; 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 .  \nTaking 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 CXOUJ171405.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 .  \n- 5 https://www.atnf.csiro.au/research/pulsar/psrcat/  \n6 Our 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  \nHowever, 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 (2018, 2021, 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 WolfRayet 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 .  \nAcomprehensive 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).  \nIn 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.  \nsince they are only observable in binary systems, and are therefore subject to selection effects (e.g. Lorimer 2008).", '2 MAGNETAR SAMPLE': "In Table 1, we present the sample of 24 magnetars listed in the McGill magnetar catalog 8 , 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 Bayestar19 9 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.", '3.1 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).  \nThe 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.  \nFollowing the initial query, extinction corrections and cutoffs were re-applied using the Bayestar19 11 dust map (Green et al. 2018; Green et al. 2019) for any Gaia sources above declination 𝛿 ≳ -30 · . From Bayestar19 we obtain the 𝐵 -𝑉 color excess 12 , 𝐸 𝐵 -𝑉 , which is  \n- 8 http://www.physics.mcgill.ca/~pulsar/magnetar/main.html\n- 9 https://dustmaps.readthedocs.io/en/latest/modules.html 10\n- https://www.cosmos.esa.int/web/gaia/science-performance Section 1\n- 11 https://dustmaps.readthedocs.io/en/latest/modules.html\n- 12 Note that Bayestar19 reports reddening in arbitrary units and is converted to 𝐸 𝐵 -𝑉 by multiplying by 0 . 995 as described at http://argonaut. skymaps.info/usage (Schlafly & Finkbeiner 2011). We proceed assuming this factor is approximately 1, which is accurate to the first decimal place.  \nTable 1. Magnetars and Magnetar Candidates from the McGill Catalog (revised position errors are boldfaced for clarity)  \nAbbreviated magnetar names 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 ∗ . 𝐴 𝑉 extinction estimates are from the Bayestar19 map for magnetars with 𝛿 ≳ -30 · . Magnetars outside this range have 𝐴 𝑉 from the X-ray column density 𝑁 𝐻 (indicated by a ‡ ) using the Predehl & 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 · (indicated by a † ; see Section 3 for details).  \nRelevant citations referenced in the rightmost column are below: (1) Haschke et al. (2012) (2) Lamb et al. (2002) (3) Durant & van Kerkwĳk (2006) (4) Hulleman et al. (2004) (5)Tendulkar et al. (2013) (6) van der Horst et al. (2010) (7) Lin et al. (2011) (8) Göˇgüş et al. (2010b) (9) Gaensler & Chatterjee (2008) (10) Kulkarni et al. (2003) (11) Klose et al. (2004) (12) Wang & Chakrabarty (2002) (13) Gaensler et al. (2005) (14) Tiengo et al. (2010) (15) Deller et al. (2012) (16) Gelfand & Gaensler (2007) (17) Levin et al. (2010) (18) Anderson et al. (2012a) (19) Corbel et al. (1999) (20) Wachter et al. (2004) (21) Kothes & Dougherty (2008) (22) Muno et al. (2006) (23) Israel et al. (2003) (24) Tian & Leahy (2012) (25) Halpern & Gotthelf (2010a) (26) Halpern & Gotthelf (2010b) (27) Bower et al. (2014) (28) Mori et al. (2013a) (29) Shannon & Johnston (2013) (30) Bower et al. (2015) (31) Bibby et al. (2008) (32) Israel et al. (2005) (33) Tendulkar et al. (2012) (34) Minter et al. (2008) (35) Helfand et al. (2007) (36) Karuppusamy et al. (2020) (37) Stamatikos et al. (2020) (38) Ding et al. (2020a) (39) Scholz et al. (2012) (40) Pagani et al. (2011) (41) Göˇgüş et al. (2010a) (42) Leahy & Tian (2008a) (43) Kargaltsev et al. (2012) (44) Tian & Leahy (2008) (45) Wachter et al. (2004) (46) Vasisht & Gotthelf (1997) (47) Zhou et al. (2014) (48) Davies et al. (2009) (49) Frail et al. (1999) (50) Gaensler (2014) (51) Israel et al. (2016) (52) Zhong et al. (2020) (53) Kothes & Foster (2012) (54) Fahlman & Gregory (1981) (55) Hulleman et al. (2001) (56) Barthelmy et al. (2016) (57) Doroshenko et al. (2021) (58) Richardson et al. (2023) (59) Cline et al. (2000) (60) Lamb et al. (2003) (61) Mereghetti et al. (2012) (62) Torii et al. (1998b) (63) Gaensler et al. (1999b) (64) Tam et al. (2006) (65) Sakamoto et al. (2011b) (66) Leahy & Tian (2007) (67) Gotthelf et al. (2000) (68) Helfand et al. (2003) (69) Blumer et al. (2019) (70) Lin et al. (2012) (71) Roming et al. (2008) (72) Barthelmy et al. (2016) (73) Jing et al. (2023) (74) Ponti et al. (2015) (75) Yalinewich et al. (2017) (76) Cline et al. (1982) (77) Park et al. (2012) (78) Park et al. (2020) (79) Sano et al. (2023) (80) Borkowski & Reynolds (2017) (81) Zhou et al. (2019) (82) Kothes et al. (2018) (83) Zhou et al. (2020) (84) Nakano et al. (2017) (85) Leahy & Tian (2008b) (86) Straal et al. (2023) (87) Ding et al. (2020b)  \nconverted to 𝑉 -band extinction with 𝐴 𝑉 ≈ 3 . 2 𝐸 𝐵 -𝑉 and the 𝐺 -𝑉 excess using the Johnson-Cousins Relations 13 . These are used along with the parallax distance 𝑑 from the Bailer-Jones (2023) catalog (see also Bailer-Jones et al. 2021) to obtain the absolute 𝐺 -and  \n𝑉 -band magnitudes, 𝑀 𝐺 and 𝑀 𝑉 . Details on conversions between photometric systems are in Appendix A.  \nFor sources with 𝛿 ≲ -30 · , 𝑀 𝑉 is computed with Equations A5A8 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  \net al. 2022; Fouesneau et al. 2022, Section 3.4, and the Gaia DR3 Documentation Section 11.3.3 for details).  \nSources were constrained to have parallax 0 . 008 < 𝜔 ± 𝜎 𝜔 < 1 . 25 mas based on the distances of the farthest (CXOU J010043.1721134, 62 . 4 ± 1 . 6 kpc) and nearest (Swift J1822.3-16.06, 1 . 6 ± 3 kpc) magnetars in the McGill catalog, placing an upper limit ∼ 125 kpc on the sample. The query was further constrained to Galactic longitudes -80 · ≤ 𝑏 ≤ 20 · based on magnetar locations. Figure 1 shows the resulting Gaia sample of 156572 potential OB star sources.  \nWeuse a Monte Carlo simulation to estimate p-values for the Gaia search. The criteria for association for a given sampled trajectory are below:  \n- (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 · .\n- (ii) The absolute 𝑉 -band magnitude of the Gaia source is consistent with that of an OB (possibly A) star ( -10 ≲ 𝑀 𝑉 ≲ 3) 14 .\n- (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.\n- (iv) The magnetar's physical travel time T 3D , mag is consistent with its age T .\n- (v) The source's parallax distance is consistent with magnetar's distance, such that Δ 𝑑 ≈ 0.  \nwhich 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:  \n𝑝 = 1 -(1 -𝑝 1 )(1 -𝑝 2 )(1 -𝑝 3 )(1 -𝑝 4 ) (1)  \nWe 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.  \nEach 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.", '3.2 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  \nTable 2. Magnetar Age Upper Limits Used for the Gaia Search  \nT decay is estimated in this work (see Appendix D1 in this work and Ferrario &Wickramasinghe 2008) using the magnetic field 𝐵 𝑑 and X-ray luminosity 𝐿 𝑋 reported in the McGill Magnetar Catalog (Olausen & Kaspi 2014, and references therin). Spindown ages are those from the McGill catalog, and are approximated as the characteristic age T rot = 𝑃 (2 / 𝑃 ) -1 . SNR ages T SNR are used only for magnetars with confident SNR associations. References for the spindown and SNR ages are in the rightmost column.  \nsample, 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 trajectories 15 . 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.", '3.2.1 PSR J0538+2817': "For PSRJ0538+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 DR33441732292729818752, was assigned a p-value 𝑝 >  \n<!-- image -->  \n<!-- image -->  \nFigure 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).  \n<!-- image -->  \n99% 16 and therefore was not identified as a candidate. Upon fur-  \n16 The minimum and maximum p-values are set by the number of trials, which varies for each magnetar depending on what constraints on its position, distance, and proper motion are available. 1000 and 100 trials were used for PSRB1951+32 and PSRJ0538+2817, respectively.  \nther inspection, we find the high p-value is the result of a high false positive rate; among the simulated sources used to estimate the H0 distribution, ∼ 95% are found with all Monte Carlo samples within the { Θ , 𝑀 𝑉 } association region ( ˆ 𝑓 = 0; see Appendix C1). We attribute this to the large association region in { Θ , 𝑀 𝑉 } space which extends to ≳ 1 ' (see Figure 2). While pulsars have larger confidence  \nregions in { Θ , 𝑀 𝑉 } space than magnetars due to the higher age limit and lower distances, the high proper motions of PSR J0538+2817, PSRJ1809-2332, and PSR B2334+61 further extend the region, resulting in unusually high false positive rates ≳ 99%, ≳ 91%, and ≳ 76%, respectively.  \nWetherefore 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 ˆ 𝑓 < ˆ 𝑓 𝑐 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.", '3.2.2 PSR B1951+32': "One Gaia companion candidate is identified for PSR B1951+32 with 𝑝 < 0 . 1%, Gaia DR3 2035109887486067456 (catalogued as TYC2673-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 kms -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.", '3.3 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 WR77F 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.738131, 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.", '3.3.1 SGR J1822.3-1606': "The primary candidate for an unbound system is SGR J1822.3-1606 and Gaia DR34097832458955594880, 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 . 1 17 , which are both reasonable for a low-mass B-type star. SGR J1822.31606 does not have a measured proper motion, and our simulation suggests 𝜇 𝛼 cos( 𝛿 ) = -29 +15 -24 masyr -1 , 𝜇 𝛿 = -25 +13 -19 masyr -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 kms -1 . Notably, the permitted stellar velocity | -→ 𝑣 | = 60 +7 -5 kms -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.", '3.3.2 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 CXOU164710.2-455216 based 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-455216 lies 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 𝜎 . Finally, the valid trajectories require velocities | -→ 𝑣 mag |≈ 390 ± 170 km s -1 and | -→ 𝑣 |≈ 80 +20 -13 kms -1 for CXOU164710.2-455216 and WR77F, respectively. The former is consistent with Hobbs et al. (2005) while the latter suggests that WR77F is a runaway star ( | -→ 𝑣 | ≳ 30 km s -1 Renzo et al. 2019). 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-  \nupporting the results of Clark et al. (2014).", '4 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.  \nWhile 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.", '4.1 Pan-STARRS and SkyMapper Data and Optical Search Method': "The Pan-STARRS Data Release 2 (PS1 DR2) catalog contains optical images of the northern sky (above 𝛿 > -30 · ) in 5 passbands ( 𝑔 -, 𝑟 -, 𝑖 -, 𝑧 -,and 𝑦 -bands, or ' 𝑔𝑟𝑖𝑧𝑦 ') with 0 . 25 '' pixels (Flewelling et al. 2020; Chambers et al. 2019). 1 ' × 1 ' image cutouts were obtained from the PS1 Image Server 20 at the position of each magnetar and candidate. 17 magnetars and all 7 magnetar candidates fall within the PS1 survey range. For the remaining 7 with 𝛿 ≤ -30 · , we use images from the SkyMapper Southern Survey Data Release 4 (SkyMapper DR4 21 ), which records in 𝑢 -, 𝑣 -, 𝑔 -, 𝑟 -, 𝑖 -, and 𝑧 -bands (' 𝑢𝑣𝑔𝑟𝑖𝑧 ') with 0 . 5 '' pixels (Wolf et al. 2018; Onken et al. 2024). Note the point source detection limit is ∼ 18 . 0 in comparison to ∼ 21 . 5 for PS1.  \nIn 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:  \n- (i) The PS1 source and magnetar are located at the same RA and declination.  \n(ii) If we assume the source absolute magnitude is that of an OB star ( -10 ≲ 𝑀 𝑉 ≲ 3), the distance is equal to the magnetar's distance.  \n7-8-ir-sensitivity  \n20  \n21  \nhttps://ps1images.stsci.edu/ps1image.html  \nhttps://skymapper.anu.edu.au/  \n,  \ndoi.org/10.25914/  \n5M47-S621  \n(iii) The difference of the Kron and point source (PSF) magnitude ( 𝑀 𝑥, Kron -𝑀 𝑥, PSF ) exceeds the Farrow et al. (2014) point source cutoff for x = ' 𝑢𝑣𝑔𝑟𝑖𝑧𝑦 ' bands as available (Kron 1980):  \n(2) Δ 𝑀 min ( 𝑀 𝑥, Kron ) = -0 . 319 + 0 . 129( 𝑀 𝑥, Kron -21) + 0 . 007( 𝑀 𝑥, Kron -21) 2  \nfor which we define p-values 𝑝 1 , 𝑝 2 , and 𝑝 3 , respectively. The three p-values are combined to get the total association p-value:  \n𝑝 = 1 -(1 -𝑝 1 )(1 -𝑝 2 )(1 -𝑝 3 ) (3)  \nAn explicit derivation of 𝑝 1 , 𝑝 2 , and 𝑝 3 is provided in Appendix B1.  \nThe 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 magnitudes 22 , 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.", '4.2 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. 2012, 2013; Smith et al. 2018). Additional sources were identified 'by-eye' in each image to include in the search.  \nThe criteria used for the IR search are below:  \n- 24 https://vizier.cfa.harvard.edu/viz-bin/VizieR?-source= II/316  \nFigure 2. ( Left ) Monte Carlo sampled trajectories for candidate Gaia DR3 sources with p-values 𝑝 < 5% possibly associated with 3 magnetars ( right , bottom ), PSRJ0538+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 CXOUJ164710.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 ˆ 𝑓 less than the critical threshold, ˆ 𝑓 𝑐 , have colorbar corresponding to ˆ 𝑓 (see Appendix B4). ˆ 𝑓 𝑐 is indicated by a red horizontal line on each colorbar. We recover the proposed companions of CXOUJ164710.2-455216 (WR77F; upper right ) and PSR J0538+2817 (HD 37424; top left ) with 𝑝 < 5%. We also identify potential runaway companions of SGR1822.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.  \n<!-- image -->  \n- (i) The source and magnetar are located at the same RA and declination.\n- (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.  \nfor 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:  \n𝑝 = 1 -(1 -𝑝 1 )(1 -𝑝 2 ) (4)  \nAn 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 .", '4.3 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 SGRJ1822.3-1606, but are most likely false positives as discussed in Appendix C3.", "4.3.1 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 companionstarCPD-29 2176,2MASSJ07554248-2933535.Wealso recover a PS1 source within the 2MASS source's 60 mas error circle, PSOJ118.9270-29.5649, and therefore conclude this is the optical counterpart of CPD-29 2176.  \nHowever, we note that the association of SGR 0755-2933 with CPD-292176 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  \n8 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.", '4.3.2 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 CXOUJ171405.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 VIRACS384090103 as a potential companion of CXOU J171405.738131.', '5.1 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 2018, 2021, 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.  \nThe 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 Portal 27 . For two magnetars (SGR0501+4516 and 1E2259+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 (CIRADA 28 ) cutout service (Gordon et al. 2020, 2021, 2023a). VLASS QuickLook images 29 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 SGR2013+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 ''  \n- 27 https://data.csiro.au/\n- 28 https://cirada.ca/vlasscatalogueql0\n- 29 Epochs 3.1 and 2.2 are used for magnetars SGR 0501+4516 and 1E2259+586, respectively.  \nFigure 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)  \n<!-- image -->  \npixels 31 . One magnetar, SGR 1627-41, is not covered by any of the radio surveys discussed 32 .  \nThe 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:  \n- (i) The SNR and magnetar are located at the same RA and declination.\n- (ii) The magnetar's proper motion direction is consistent with the direction required to intersect the SNR.\n- (iii) The magnetar's proper motion and age are consistent with the observed angular separation of the magnetar and SNR.\n- (iv) The SNR's dynamical age is consistent with an SNR shell rather than a point source.  \nfor 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 :  \nˆ 𝜒 24 = -2(ln( 𝑝 1 ) + ln(1 -(1 -𝑝 2 )(1 -𝑝 3 ))) (5)  \nThis 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:  \n𝑝 = 1 -(1 -𝑝 𝜒 )(1 -𝑝 4 ) (6)  \nNote 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.  \n- 31 VLASS data is taken 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/resolution for details.\n- 32 While 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. 33 https://simbad.u-strasbg.fr/simbad/sim-fbasic\n- 34 See https://simbad.cds.unistra.fr/guide/otypes.htx for details on SIMBAD object-types.", '5.2 Radio SNR Search Test with Known Magnetar-SNR Associations': "8 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; Corbel et al. 1999), which is likely due to the large angular offset ∼ 15 . 5 ' . No proper motion is known for SGR1627-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. (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.", '5.3 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 SNRswerecoverfive(andmarginally recover a sixth with 𝑝 = 5 . 6%), one of which is nearby the Galactic Center magnetar SGR J17452900andanotherofwhichwouldsuggest3XMMJ185246.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.', '5.3.1 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).  \nSix 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 J18082057 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 (DuarteCabral 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  \ncoincident 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'.  \nSource 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 180620, 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 180620, 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 . 2 kyr (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.  \nSNRG009.7-00.0, on the contrary, was detected by Frail et al. (1994) and confirmed through its OH maser by Hewitt & YusefZadeh (2009), which constrained its distance to ∼ 4 . 7 kpc. Yeung et al. (2016) find its distance is inconsistent with W31 and SGR 180620. 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.0 appear 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.", '5.3.2 3XMM J185246.6+003317 and CXOU J185238.6+004020 as an Unbound Double Neutron Star Binary': "For 3XMMJ185246.6+003317 we recover the Kes79 SNR (SNRG033.6+00.1) with 𝑝 < 0 . 1%, which is associated with a central X-ray pulsar CXOU J185238.6+004020 (PSR J1852+0040; e.g. Gotthelf et al. 2005; Malofeev et al. 2005; Bogdanov 2014). Bogdanov (2014) derive a Sedov-Taylor expansion age of ∼ 5 kyr for the SNRwhichcontrasts with both the decay age T decay ≈ 0 . 7 Myr 37 and spindown age T rot ≳ 1 . 3 Myr of 3XMMJ185246.6+003317. 38 Magnetars are typically expected to be unobservable beyond ∼ 10 4 years; however,the spindown age T rot is a well-constrained lower limit derived by Rea 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  \ninitial discovery that the X-ray luminosity ( ∼ 10 33 -34 erg s -1 ) far exceeds its spindown luminosity ( ∼ 10 30 erg s -1 ) making its magnetar nature robust. Therefore we take both the spindown age and the decay age to be reasonable age estimates. This suggests that Kes 79 is associated with CXOU J185238.6+004020, but allows a scenario wherein the pulsar and magnetar were gravitationally bound in the past (Zhou et al. 2014).  \nConsidering this case, let the initially bound massive star progenitors of 3XMM J185246.6+003317 and CXOU J185238.6+004020 be called stars 1 and 2, respectively. When the CCSN of star 1 disrupts the binary 39 , 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's Gaia counterpart, Gaia DR34266508881354196736. 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.  \nThe 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 | -→ 𝑣 mag |≈ 398 +208 -165 kms -1 (proper motion 𝜇 𝛼 cos( 𝛿 ) ≈ -2 . 1 +1 . 2 -5 . 8 masyr -1 , 𝜇 𝛿 ≈ -9 . 8 +4 . 9 -20 . 8 masyr -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.  \nThe pulsar requires a velocity | -→ 𝑣 psr , 1 |≈ 8 +11 -4 kms -1 following the first CCSN, and | -→ 𝑣 psr , 2 |≈ 422 +192 -172 kms -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  \nexplain 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 3XMMJ185246.6+003317 and CXOUJ185238.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.", '5.3.3 SGR J1745-2900, the Galactic Center Magnetar': "SGRJ1745-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 ∗ andhighdensity 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 J17452900 (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.  \nGiven 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 0FGLJ1746.0-2900 found 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.", '5.3.4 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 052666, 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, S24RACSJ190707+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 (SNRG10.0-0.3)wasalsorejected 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).  \nSGR0526-66islocatedwithin 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.  \nFor SGR0501+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.", '6 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 CXOUJ164710.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.  \nTable 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.  \nFigure 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 3XMMJ185246.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 the Hobbs 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.  \n<!-- image -->  \nFigure 5. Radio continuum image cutouts for four magnetars ( middle, bottom ), and two magnetar candidates ( top ) with possible SNR associations. Red pointers indicate each magnetar's location and the 1 𝜎 , 2 𝜎 , and 3 𝜎 magnetar error ellipses are shown in red (although they are too small to be visible except for SGR1808-20 and SGR1801-23. For the latter, the 3 . 8 · error ellipse is too large to be effectively shown.). Note the difference in axis extents for each image. RACS 888MHz continuum images with 25 '' × 25 '' pixels are shown for SGR 1808-20, SGR 1801-23, SGR J1745-2900, and SGR 0526-66 ( top, middle ). For SGR0501+4516, a 3 GHz 1 ' × 1 ' VLASS image with 1 '' × 1 '' pixels (smoothed to 25 '' × 25 '' with a boxcar kernel) is shown enclosed in a red box indicating the queried region, while the remainder of the image is from the 1.4 GHz NVSS survey to show SNR HBH 9 (e.g. Göˇgüş et al. 2010b). However, NVSS is only sensitive to angular sizes < 970 '' ( ≲ 16 . 7 ' ), making SNR HBH9 (120 -140 ' ) unobservable in the image shown. No radio images from RACS, VLASS, or NVSS were available at SGR 1627-41's position ( bottom left ). Circular markers indicate the locations of known radio sources while square markers indicate unknown sources identified by-eye. 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; sources outlined in red (green) have 𝑝 > 5% ( 𝑝 ≤ 5%). The most likely or previously proposed SNRs (see Section 5.3) are labelled. 'Source 1' ( top left ) is S24-RACSJ180842-202849. The angular extents of SNRs, if available, are shown as ellipses. For SGR 1745-2900 ( middle left ), the proper motion direction is marked by a purple arrow, and the purple shaded region indicates the 1 𝜎 , 2 𝜎 , and 3 𝜎 past travel direction and extent; for SGR 1627-41, SGR 0501+4516, and SGR0526-66 ( bottom, middle right ), the Hobbs et al. (2005) velocity distribution is used to show possible past travel direction and extent.  \n<!-- image -->", '6.1 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:  \n- (i) -2 . 76 ≲ 𝑀 𝐻 ≲ 4 . 34 (Chrimes et al. 2022b)\n- (ii) -2 . 61 ≲ 𝑀 𝑖 ≲ 5 . 55 (derived by converting the Chrimes et al. 2022b, constraints on 𝑀 𝐻 and 𝐽 -𝐻 using the 2MASS relationships 40 )  \nSimilarly in the Gaia search for unbound companions we expect:  \n- (i) -10 ≲ 𝑀 𝑉 ≲ 1 . 185\n- (ii) 𝑣 mag following the Hobbs et al. (2005) distribution\n- (iii) 𝑣 src following the Renzo et al. (2019) distribution  \nand for the SNR search:  \n- (i) Radio luminosity 𝐿 ≈ 10 42 erg s -1 in the Rayleigh-Jeans limit (Barker et al. 2022, 2023; Hamuy 2003)\n- (ii) Angular size 25 '' ≲ ΔΘ ≲ 1 · (from observed SNR sizes in the Green (1987) catalog and considerations of free- and SedovTaylor-expansion in Truelove & McKee (1999))  \nBy 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.  \nFor 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 CXOUJ010043.1-721134 is partially incomplete for faint stars due to its large distance. However, four magnetars lie along high-extinction sightlines with 𝐴 𝑉 ≳ 17 (CXOUJ171405.7-381031, 1E1547.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.  \nThe 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).  \nFinally, the SNR search is complete in all cases for typical radio luminosities and angular sizes in the free- or Sedov-Taylor expansion  \nphases (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.", '6.2 The Magnetar Binary and Massive Star Merger Fractions': "6.2.1 Markov-Chain Monte Carlo Analysis of Magnetar Binary Fractions  \nWith 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 MarkovChain 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 𝑓 𝑛 = 1 -𝑓 𝑢 -𝑓 𝑏 and 𝑓 𝑏 = 𝑓 𝑖 + 𝑓 𝑝 .  \nThe 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 SGRJ1822.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 07552933 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:  \nP i = 𝑓 21 𝑛 𝑓 3 𝑢 𝑓 1 𝑖 (1 -𝑓 𝑢 ) 2 ( 𝑓 𝑛 + 𝑓 𝑢 ) 3 (7)  \nwhere (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 𝑓 𝑖 .  \nIn 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:  \nP c = 𝑓 10 𝑛 𝑓 1 𝑢 𝑓 15 𝑛𝑐 ( 𝑓 𝑛 + 𝑓 𝑢 ) 2 (1 -𝑓 𝑢 ) 2 (1 -𝑓 𝑛𝑐 ) 1 (8)  \nwhere 𝑓 𝑛𝑐 is the fraction of magnetars which do not originate in CCSNe. This modifies the definition for 𝑓 𝑛 to 𝑓 𝑛 = 1 -𝑓 𝑢 -𝑓 𝑏 -𝑓 𝑛𝑐 .  \nTable 3. Search Results and Completeness Ratios  \nFor 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, 'S24survey -J hhmmss ± 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. The radio search completeness is labelled with 'R', 'V', or 'N' to indicate whether the RACS, VLASS, or NVSS survey is used. The optical/IR bound search is considered partially or fully complete if either the optical or IR search are complete. These are labelled with 'P' or 'S', to indicate whether the PS1 or SkyMapper survey is used, and 'T', 'U', or 'V' to indicate whether the 2MASS, PS1, or VVV survey is used. The Gaia search is considered partially or fully complete if both the 𝑉 -band magnitude and expected relative proper motion are detectable. Four magnetars or candidates have no distance measures and the completeness cannot be determined; these are marked with '-'.Candidates for SNR, bound, or unbound companions from previous analyses are provided in the last column with references matched to the following: (1) Gaensler & Chatterjee (2008) (2) Klose et al. (2004) (3) Gelfand & Gaensler (2007) (4) Anderson et al. (2012a) (5) Corbel et al. (1999) (6) Clark et al. (2014) (7) Halpern & Gotthelf (2010a) (8) Chrimes et al. (2022b) (9) Corbel & Eikenberry (2004) (10) Kargaltsev et al. (2012) (11) Vasisht & Gotthelf (1997) (12) Vrba et al. (2000) (13) Gaensler (2014) (14) Fahlman & Gregory (1981) (15) Gaensler et al. (1999b) (16) Gotthelf et al. (2000) (17) Zhou et al. (2014) (18) Ponti et al. (2015) (19) Yalinewich et al. (2017) (20) Cline et al. (2000) (21) H. E. S. S. Collaboration et al. (2018) (22) Ding et al. (2020b)  \nNote this case is more restrictive, as SGR 0755+2933, SGR J1822.31606, 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 3XMMJ185246.6+003317 as unbound systems, and the likelihood becomes:  \nP uc = 𝑓 10 𝑛 𝑓 3 𝑢 𝑓 13 𝑛𝑐 ( 𝑓 𝑛 + 𝑓 𝑢 ) 2 (1 -𝑓 𝑢 ) 2 (1 -𝑓 𝑛𝑐 ) 1 (9)  \nResults 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), and Renzo et al. (2019) population synthesis results, the unbound fraction significantly disagrees with their prediction  \nFigure 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).  \n<!-- image -->  \nof 𝑓 𝑢 ∼ 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.  \nIn 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), and Renzo 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.  \nA similar result is obtained in the 'SNR/unbound complete' case; notably, the larger unbound fraction 𝑓 𝑢 ≈ 8 +5 -3 %allows a lower nonCCSN 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).", '6.2.2 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.  \nThe 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. 41 While some simulations suggest that CCSNe  \nTable 4. Constraints on magnetar binaries from this work  \nNote: We define two broad cases based on the completeness of the RACS SNR search: 'SNR incomplete' assumes that all magnetars originate in CCSNe, whether one is detected or not. 'SNR complete' 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.  \n† 𝑓 𝑛 , 𝑓 𝑏 , 𝑓 𝑢 , and 𝑓 𝑛𝑐 are defined such that 𝑓 𝑛 + 𝑓 𝑏 + 𝑓 𝑢 = 1 in the 'SNR incomplete' case, and 𝑓 𝑛 + 𝑓 𝑏 + 𝑓 𝑢 + 𝑓 𝑛𝑐 = 1 in the 'SNR complete' and 'SNR/unbound complete' cases. 𝑓 𝑖 and 𝑓 𝑝 sum 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.  \n‡ 𝑓 𝑚 is estimated assuming the progenitor stars are early B-type stars (9 𝑀 ⊙ < 𝑀 ∗ < 16 𝑀 ⊙ ) with binary fraction 𝐹 0 = 84 ± 9% (Moe & Di Stefano 2017).  \ncan 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:  \n𝑓 𝑢 + 𝑓 𝑏 1 -( 𝑓 𝑏 + 𝑓 𝑢 ) 𝑓 𝑞 = 𝐹 0 (1 -𝑓 𝑚 ) (10)  \nwhere 𝑓 𝑞 = 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.  \nThe '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'  \ncase is one step removed, taking unbound companions as evidence of a disruption event (CCSN) but also allowing for 𝑓 nc ≈ 45%toform from other formation channels (e.g. AIC, neutron star merger). For this reason, 𝑓 𝑛 ≈ 40%islowered, 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 07552933, 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.  \nAlternatively, 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.", '6.2.3 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 𝑀 ⊙ ),  \nFigure 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 Table 4 in 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.  \n<!-- image -->  \nand 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 . 30 for masses 0 . 1 -0 . 5 𝑀 ⊙ and 𝛼 = -2 . 35 for 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.  \nAn 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 180820, 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.", '6.2.4 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. 2014, 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.  \nFurthermore, 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.  \n2019). 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.", '6.3 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.  \nIt 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).  \nDirect 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).", '7 CONCLUSION': '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 07552933 and marginally detect the proposed unbound companion of CXOUJ164710.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), and Renzo et al. (2019).  \nOur 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.  \nSignificant 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 (DSA2000) 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.', '8 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 .', '9 ACKNOWLEDGEMENTS': "The authors would like to thank the members of the Deep Synoptic Array (DSA-110) team for their support and insight with regards to this search effort. The authors thank staff members of the Owens Valley Radio Observatory and the Caltech radio group, including Kristen Bernasconi, Stephanie Cha-Ramos, Sarah Harnach, Tom Klinefelter, Lori McGraw, Corey Posner, Andres Rizo, Michael Virgin, Scott White, and Thomas Zentmyer. Their tireless efforts were instrumental to the success of the DSA-110. The DSA-110 is supported by the National Science Foundation Mid-Scale Innovations Program in Astronomical Sciences (MSIP) under grant AST-1836018. The authors would also like to thank Casey Law, Alexa Gordon, Alice  \nCurtin, Maxim Lyutikov, and an anonymous referee for their insightful discussions and feedback on the early draft of this manuscript. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1745301.  \nThe Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. The UKIDSS project is defined in Lawrence et al. (2007b). The national facility capability for SkyMapper has been funded through ARC LIEF grant LE130100104 from the Australian Research Council, awarded to the University of Sydney, the Australian National University, Swinburne University of Technology, the University of Queensland, the University of Western Australia, the University of Melbourne, Curtin University of Technology, Monash University and the Australian Astronomical Observatory. SkyMapper is owned and operated by The Australian National University's Research School of Astronomy and Astrophysics. The survey data were processed and provided by the SkyMapper Team at ANU. The SkyMapper node of the All-Sky Virtual Observatory (ASVO) is hosted at the National Computational Infrastructure (NCI). Development and support of the SkyMapper node of the ASVO has been funded in part by Astronomy Australia Limited (AAL) and the Australian Government through the Commonwealth's Education Investment Fund (EIF) and National Collaborative Research Infrastructure Strategy (NCRIS), particularly the National eResearch Collaboration Tools and Resources (NeCTAR) and the Australian National Data Service Projects (ANDS).  \nUKIDSS uses the UKIRT Wide Field Camera (WFCAM; Casali et al. (2007)). The photometric system is described in Hewett et al. (2006), and the calibration is described in Hodgkin et al. (2009). The pipeline processing and science archive are described in Irwin et al. (2004) and Hambly et al. (2008). The VISTA Data Flow System pipeline processing and science archive are described in Irwin et al. (2004), Hambly et al. (2008) and Cross et al. (2012). We have used data from the 5th data release of the Vista Variables in the Via Lactea (VVV) Survey, which is described in detail in Nikzat et al. (2022), and the VVV Infrared Astrometric Catalogue, described in detail in Smith et al. (2018). This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.  \nThis work also uses data obtained from Inyarrimanha Ilgari Bundara / the Murchison Radio-astronomy Observatory. We acknowledge the Wajarri Yamaji People as the Traditional Owners and native title holders of the Observatory site. CSIRO's ASKAP radio telescope is part of the Australia Telescope National Facil-  \nity (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.", 'REFERENCES': '- Greisen E. W., Calabretta M. R., 2002, Astronomy & Astrophysics, 395, 1061 Güver T., Özel F., 2009, Monthly Notices of the Royal Astronomical Society, 400, 2050  \n- Shimwell T., et al., 2022, Astronomy & Astrophysics, 659, A1', 'APPENDIX A: GAIA ADQL QUERY AND PHOTOMETRIC CONVERSIONS': 'Below is the ADQL query used to compile the initial sample of Gaia sources:', 'select * from (': 'select *,phot\\_g\\_mean\\_mag\\_absolute - (-0.02704 +  \n0.1424*bp\\_rp\\_absolute -  \n0.2156*POWER(bp\\_rp\\_absolute,2) +  \n0.01426*POWER(bp\\_rp\\_absolute,3))  \nAS phot\\_v\\_mean\\_mag\\_absolute from (  \nselect source\\_id,  \nra,  \nra\\_error,  \ndec,  \ndec\\_error,  \nparallax,  \nparallax\\_error,  \nparallax\\_over\\_error,  \npm,  \npmra,  \npmra\\_error,  \npmdec,  \npmdec\\_error,  \nruwe,  \nphot\\_g\\_n\\_obs,  \nphot\\_g\\_mean\\_flux,  \nphot\\_g\\_mean\\_flux\\_error,  \nphot\\_g\\_mean\\_flux\\_over\\_error,  \nphot\\_g\\_mean\\_mag,  \nphot\\_bp\\_n\\_obs,  \nphot\\_bp\\_mean\\_flux,  \nphot\\_bp\\_mean\\_flux\\_error,  \nphot\\_bp\\_mean\\_flux\\_over\\_error,  \nphot\\_bp\\_mean\\_mag,  \nphot\\_rp\\_n\\_obs,  \nphot\\_rp\\_mean\\_flux,  \nphot\\_rp\\_mean\\_flux\\_error,  \nphot\\_rp\\_mean\\_flux\\_over\\_error,  \nphot\\_rp\\_mean\\_mag,  \nphot\\_bp\\_rp\\_excess\\_factor,  \nbp\\_rp,  \nbp\\_g,  \ng\\_rp,  \nradial\\_velocity,  \nradial\\_velocity\\_error,  \nteff\\_gspphot,  \nteff\\_gspphot\\_lower,  \nteff\\_gspphot\\_upper,  \nlogg\\_gspphot,  \nlogg\\_gspphot\\_lower,  \nlogg\\_gspphot\\_upper,  \nmh\\_gspphot,  \ndistance\\_gspphot,  \ndistance\\_gspphot\\_lower,  \ndistance\\_gspphot\\_upper,  \nazero\\_gspphot,  \nazero\\_gspphot\\_upper,  \nag\\_gspphot,  \nag\\_gspphot\\_lower,  \nag\\_gspphot\\_upper,  \nebpminrp\\_gspphot,  \nebpminrp\\_gspphot\\_lower,  \nebpminrp\\_gspphot\\_upper,  \nl,  \nb,  \nphot\\_g\\_mean\\_mag-ag\\_gspphot-(  \n5*LOG10(1000*ABS(parallax))-5)  \nAS phot\\_g\\_mean\\_mag\\_absolute,  \nbp\\_rp - ebpminrp\\_gspphot  \nAS bp\\_rp\\_absolute  \nfrom gaiadr3.gaia\\_source  \nAS m)  \nAS mm)  \nAS mmm  \nWHERE parallax+parallax\\_error>0.008 and parallax-parallax\\_error<1 and phot\\_v\\_mean\\_mag\\_absolute<1.18 and b>-80 and b<20  \nThe magnitude limit was re-applied, along with the color limit, after the initial query to incorporate the Bayestar19 extinction estimates. From Bayestar19 we obtain the 𝐸 𝐵 -𝑉 color excess and convert it to 𝐸 45  \n𝐺 -𝑉 using the Johnson-Cousins Relations :  \n(A1) 𝐸 𝐺 -𝑉 = -0 . 04749 -0 . 0124 𝐸 𝐵 -𝑉 -0 . 2901 𝐸 2 𝐵 -𝑉 + 0 . 02008 𝐸 3 𝐵 -𝑉 + 0 . 04772 𝐸 4 𝐵 -𝑉  \nand similarly the 𝐺 𝐵𝑃 -𝐺 𝑅𝑃 excess:  \n45 See the Gaia DR3 Documentation at https://gea.esac.esa.int/ archive/documentation/GDR3/ Section 5.5.1 Table 5.9 for details.  \n𝐸 𝐺 𝐵𝑃 -𝐺 𝑅𝑃 = 0 . 06483 + 1 . 575 𝐸 𝐵 -𝑉 -0 . 7815 𝐸 2 𝐵 -𝑉 + 0 . 5707 𝐸 3 𝐵 -𝑉 -0 . 176 𝐸 4 𝐵 -𝑉 + 0 . 01916 𝐸 5 𝐵 -𝑉 (A2)  \n(e.g. Fernie 1983; Stetson 2000). The 𝑉 - and 𝐺 - band extinction are then computed using:  \n𝐴 𝑉 ≈ 3 . 2 𝐸 𝐵 -𝑉 (A3)  \n𝐴 𝐺 ≈ 𝐴 𝑉 + 𝐸 𝐺 -𝑉 (A4)  \n(e.g. Savage & Mathis 1979; Seaton 1979). The 𝐺 -band absolute magnitude is then:  \n𝑀 𝐺 = 𝑚 𝐺 -𝐴 𝐺 -(5log 10 ( 𝑑 ) -5) (A5)  \nwhile the absolute color is:  \n𝐺 𝐵𝑃 -𝐺 𝑅𝑃 = 𝑚 𝐺 𝐵𝑃 -𝑚 𝐺 𝑅𝑃 -𝐸 𝐺 𝐵𝑃 -𝐺 𝑅𝑃 (A6)  \nFinally, we derive the absolute 𝐺 -𝑉 color from the Johnson-Cousins relations:  \n(A7) 𝐺 -𝑉 = -0 . 027 + 0 . 01424( 𝐺 𝐵𝑃 -𝐺 𝑅𝑃 ) -0 . 2156( 𝐺 𝐵𝑃 -𝐺 𝑅𝑃 ) 2 + 0 . 07426( 𝐺 𝐵𝑃 -𝐺 𝑅𝑃 ) 3  \nwhich is used to get the absolute 𝑉 -band magnitude:  \n𝑀 𝑉 = 𝑀 𝐺 -( 𝐺 -𝑉 ) (A8)  \nFor source with 𝛿 ≤ -30 · , the extinction and color excess from Gaia GSP-Aeneas photometry are used as described in Section 3.', 'APPENDIX B: DEFINITION OF ASSOCIATION P-VALUES': 'In this appendix, we explicitly define p-values for association of magnetars with SNRs, bound companions, and unbound companions from the radio, optical, IR, and Gaia searches. The p-values are based on the criteria defined in Sections 4.1-3, and are estimated through Monte Carlo simulations.', 'B1 Optical Search for Bound Companions': "Recall for the bound companion search, the null hypothesis ( H0 ) is that the magnetar and source are not bound in a binary system, and we use the position, distance, and difference between PSF and Kron magnitudes (' 𝑔𝑟𝑖𝑧𝑦 ' for PS1, ' 𝑢𝑣𝑔𝑟𝑖𝑧 ' for SkyMapper) as association criteria. For each source, we draw 1000 samples of the RA and declination of the source and magnetar ( 𝛼 PS1 , 𝛿 PS1 , 𝛼 mag , 𝛿 mag), and the distance of the magnetar ( 𝑑 mag), assuming each are normally distributed. Using Δ 𝛼 = | 𝛼 PS1 -𝛼 mag | , Δ 𝛿 = | 𝛿 PS1 -𝛿 mag | , we compute the search statistic ˆ 𝑔 1 as the fraction of samples that fall below the 95 th percentiles of folded zero-mean normal distributions with 𝜎 𝛼 = ( 𝜎 2 𝛼, PS1 + 𝜎 2 𝛼, mag ) 1 / 2 and 𝜎 𝛿 = ( 𝜎 2 𝛿, PS1 + 𝜎 2 𝛿, mag ) 1 / 2 , respectively. The H0 distribution of ˆ 𝑔 1 is estimated by repeating the simulation 100 times with mean Δ 𝛼, Δ 𝛿 selected uniformly between 1 -5 𝜎 𝛼/𝛿 , and the p-value 𝑝 1 is taken as the percentile of ˆ 𝑔 1 given the resulting ˆ 𝑔 1 distribution.  \n𝑝 2 is computed by a similar process; since the errors on apparent  \nmagnitudes are small, we assume they are negligible, and use 𝑔 -and 𝑟 -band magnitudes to derive the apparent 𝑉 -band magnitude:  \n𝑚 𝑉 = 𝑚 𝑟 + 0 . 005 + 0 . 462( 𝑚 𝑔 -𝑚 𝑟 ) + 0 . 013( 𝑚 𝑔 -𝑚 𝑟 ) 2 (B1)  \n(Tonry et al. 2012). From Bayestar19 we obtain the 𝐵 -𝑉 color excess 46 along the source's sightline, 𝐸 𝐵 -𝑉 , for distances between 60 pc -60 kpc, converting to the 𝑉 -band extinction with 𝐴 𝑉 ≈ 3 . 2 𝐸 𝐵 -𝑉 . If Bayestar19 extinction estimates are not available, we use the 𝑁 𝐻 -𝐴 𝑉 relation from Predehl & Schmitt (1995) from Table 1. Then, after samples of the absolute magnitude are drawn from a uniform distribution between -10 ≲ 𝑀 𝑉 ≲ 3, linear interpolation is used to obtain a sample of the PS1 source distance, 𝑑 PS1 , from the difference 𝑚 𝑉 -𝑀 𝑉 . Finally, we use the magnetar distance samples to get samples of Δ 𝑑 = | 𝑑 PS1 -𝑑 mag | , which is again compared to the 95 𝑡 ℎ percentile of a folded zero-mean distribution and H0 samples between 1 -10 𝜎 Δ 𝑑 to get the statistic ˆ 𝑔 2 and p-value 𝑝 2 .  \nFinally, the Kron magnitude describes the flux of extended sources within a characteristic Kron radius (Kron 1980). As described by Farrow et al. (2014), the difference 𝑚 𝑥, Kron -𝑚 𝑥, PSF can be used to distinguish between point sources and extended sources by imposing a cutoff 𝑚 𝑥, cutoff ( 𝑚 𝑥, Kron ). Point sources are expected to exceed this limit, and we adopt the Farrow et al. (2014) fit equation for the cutoff:  \n(B2) 𝑚 𝑥, cutoff ( 𝑚 𝑥, Kron ) = -0 . 319 + 0 . 129( 𝑚 𝑥, Kron -21) + 0 . 007( 𝑚 𝑥, Kron -21) 2  \nSince the apparent magnitude error is negligible, we set the p-value 𝑝 3 , 𝑥 = 0 if 𝑚 𝑥, Kron -𝑚 𝑥, PSF > 𝑚 𝑥, cutoff ( 𝑚 𝑥, Kron ) and 𝑝 3 , 𝑥 = 1 if 𝑚 𝑥, Kron -𝑚 𝑥, PSF ≤ 𝑚 𝑥, cutoff ( 𝑚 𝑥, Kron ), combining results in all available passbands to get 𝑝 3 :  \n𝑝 3 = 1 -' 𝑢𝑣𝑔𝑟𝑖𝑧𝑦 ' Y 𝑥 (1 -𝑝 3 , 𝑥 ) (B3)  \nFor SkyMapper sources, multiple measurements for the Kron and PSF magnitudes are available in each band, each of which is used to compute 𝑝 3 . We additionally exclude SkyMapper sources with stellar classification probability ( class\\_star ) < 90% to avoid extended objects (Wolf et al. 2018; Onken et al. 2024).", 'B2 Infrared Search for Bound Companions': "Like the optical search, the IR search uses the position, distance, and PSF 𝐽 -and 𝐻 -band magnitudes as criteria for association, with H0 that the magnetar and source are not bound in a binary system. 𝑝 1 is computed in a similar way to the optical data for UKIDSS and VVV using their respective PSF FWHMs, 1 '' and 0 . 75 '' , as errors in RA and declination. For 2MASS sources, the RA and declination are sampled as a bi-normal distribution using the reported major axis ( 𝑎 IR ), minor axis ( 𝑏 IR ), and position angle ( 𝜙 IR ) of the error ellipse. The covariance matrix of the bi-normal distribution is given by:  \n(B4) 𝜎 pos  \n= GLYPH<20> ( 𝑎 IR sin( 𝜙 IR )) 2 + ( 𝑏 IR cos( 𝜙 IR )) 2 ( 𝑎 2 IR -𝑏 2 IR )cos( 𝜙 IR )sin( 𝜙 IR ) ( 𝑎 2 IR -𝑏 2 IR )cos( 𝜙 IR )sin( 𝜙 IR ) ( 𝑎 IR cos( 𝜙 IR )) 2 + ( 𝑏 IR sin( 𝜙 IR )) 2 GLYPH<21>  \n46 Note that Bayestar19 reports reddening in arbitrary units and is converted to 𝐸 𝐵 -𝑉 by multiplying by 0 . 995 as described at http://argonaut. skymaps.info/usage (Schlafly & Finkbeiner 2011). We proceed assuming this factor is approximately 1, which is accurate to the first decimal place.  \nWe 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.  \nFor 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 SNRis 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.  \nThe 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 ' .  \nFinally, 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  \nuniform distribution between 2 . 75 -120 𝑀 ⊙ 47 and use an explosion energy 10 51 ergs 48 (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:  \n(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  \nFor 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.  \nFor 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  \nsamples 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.  \nNext, 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:  \n-→ 𝑟 o , mag , i = -→ 𝑟 mag , i -𝜏 𝑖 -→ 𝑣 mag , i (B6)  \n-→ 𝑟 𝑜,𝑖 = -→ 𝑟 𝑖 -𝜏 𝑖 -→ 𝑣 𝑖 (B7)  \nWe 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\\_separation method, 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 𝑀 𝑉,𝑖 :  \n𝑀 𝑉,𝑖 = ( 𝑚 𝐺 -𝐴 𝐺 -(5log 10 ( 𝑑 𝑖 ) -5)) -( 𝐺 -𝑉 ) (B8)  \nΘ 𝑖 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 ˆ 𝑓 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 ˆ 𝑓 𝑐 were derived through a Receiver Operating Characteristic (ROC) analysis for each magnetar, which is summarized in Appendix B4.1. By comparing the final ˆ 𝑓 to samples in the H1 hypothesis (obtained during the ROC test), we estimate a p-value, 𝑝 1 , as the percentile of ˆ 𝑓 within the H1 sample.  \nTo further narrow down the candidate Gaia sources with ˆ 𝑓 <  \nˆ 𝑓 𝑐 , we iteratively compute the true intersection points -→ 𝑟 N ,𝑖 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:  \n𝜏 2D ,𝑖 ≈ ΔΦ N ,𝑖 √︃ 𝜇 2 𝛼,𝑖 cos 2 ( 𝛿 𝑖 ) + 𝜇 2 𝛿,𝑖 (B9)  \n𝜏 2D , mag ,𝑖 ≈ ΔΦ N , mag ,𝑖 √︃ 𝜇 2 𝛼, mag ,𝑖 cos 2 ( 𝛿 mag ,𝑖 ) + 𝜇 2 𝛿, mag ,𝑖 (B10)  \nwhere ΔΦ N 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:  \nT 3D , mag ,𝑖 ≈ | -→ 𝑟 mag ,𝑖 --→ 𝑟 N ,𝑖 | | -→ 𝑣 mag ,𝑖 | (B11)  \nWe 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; DohmPalmer & Jones 1996; Anderson et al. 2012b; Blumer et al. 2019; Borkowski & Reynolds 2017; Lyman et al. 2022; Gaensler et al. 1999a). If no SNR association is known or its age is unconstrained, we estimate the age T decay from from continuum X-ray luminosity measurements using the crustal field decay model (see Ferrario & Wickramasinghe 2008, and Appendix D1). This applies to 18 magnetars in our sample. If no SNR association is known and the X-ray luminosity or surface magnetic field strength are unknown, the spindown age T rot reported in the McGill magnetar catalog is used in its place (Olausen & Kaspi 2014). 5 magnetars in the sample require spin-down ages to be used, decreasing confidence in their p-values since T rot is predicted to over-estimate magnetar ages (e.g. Kaspi & Beloborodov 2017). For 5 magnetars, none of the three age estimates are known, and the final p-value is computed from 𝑝 1 and 𝑝 2 alone.  \nFor 𝑝 4 , we use samples of the magnetar and source distances from the initial Monte Carlo search to estimate the distribution of the difference, Δ 𝑑 𝑖 = 𝑑 𝑖 -𝑑 mag , i . The width of this distribution is dominated by the magnetar's distance error. We estimate the H1 distribution, for which the magnetar and source are at the same distance, by shifting the Δ 𝑑 samples to have median 0. We then compute the search statistic as the fraction of Δ 𝑑 samples within the 3 𝜎 width (0 . 15 -99 . 85 th percentiles) of the H1 distribution. We compare this to 100 realizations of the 𝐻 0 distribution which are shifted to have medians between 1 -6 𝜎 𝑑 , where 𝜎 𝑑 is the width of the H1 distribution, and compute the search statistic for each trial. The p-value 𝑝 4 is computed as the percentile of the observed Δ 𝑑 statistic among the 𝐻 0 distribution of statistics.", 'B4.1 Receiver Operating Characteristic Analysis': "The optimal critical threshold for the Gaia Monte Carlo search, ˆ 𝑓 𝑐 , was derived through an ROC analysis. To compute the true positive  \nTable B1. Critical Thresholds for Magnetar Association Tests  \nrate (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 ˆ 𝑓 were computed for each source. The TPR for each trial ˆ 𝑓 𝑐 was computed as the fraction of sources with ˆ 𝑓 < ˆ 𝑓 𝑐 . The false negative rate (FNR) is computed as FNR = 1-TPR.  \nTo 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 ˆ 𝑓 𝑐 was computed as the fraction of total sources (from all tiles) with ˆ 𝑓 < ˆ 𝑓 𝑐 .  \nWe 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  \n0%, the threshold is ˆ 𝑓 𝑐 = 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 ˆ 𝑓 𝑐 = 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 ˆ 𝑓 𝑐 for each magnetar in Table B1.", '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': "For 3XMMJ185246.6+003317, four candidates, Gaia DR34266311617807818624 (Source 1), Gaia DR34266313683745630080 (Source 2), Gaia DR342663155048110807488 (Source 3), and Gaia DR34258978601120897664 (Source 4), are identified with 𝑝 < 5%, all requiring the magnetar to have proper motion 𝜇 𝛼 cos( 𝛿 ) ∼ -5 and 𝜇 𝛿 ∼ -3 . 5 mas yr -1 . For Sources 1 and 2, the stellar velocities are consistent with walkaway stars with | -→ 𝑣 |≈ 28 +8 -5 kms -1 for Source 1 and | -→ 𝑣 |≈ 30 +8 -5 kms -1 for Source 2; the travel times are consistent with the magnetar's at the 1 𝜎 and 2 𝜎 levels, respectively. Higher velocities more consistent with runaway stars are required for Sources 3 and 4: | -→ 𝑣 |≈ 68 +10 -15 kms -1 and | -→ 𝑣 |≈ 143 +6 -10 kms -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. 3XMMJ185246.6+00317's distance was estimated from the nearby SNR Kes79 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 3XMMJ185246.6+003317 distinguish 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 Kes79 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 SGR0755-2933 were rejected as false positives.  \nEight 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 SGR0755-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 PSRJ1622-4950 and SGR J1822.3-1606 were rejected as false positives.', 'C3.1 PSR J1622-4950': "Through the optical search, the SkyMapper source SMSSJ162244.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.", 'C3.2 SGR J1822.3-1606': "Three distinct IR sources nearby SGR J1822.3-1606 are found to have 𝑝 < 5%: 2MASS18221839-1604241 ( 𝑝 < 0 . 1%; UKIDSSJ182218.39-160424.0, UKIDSSJ182218.39-160424.1), 2MASS18221794-1604259 ( 𝑝 < 0 . 1%; UKIDSSJ182217.94160425.9, UKIDSS J182217.94-160426.0), and 2MASS 182217971604360 ( 𝑝 ≈ 5%; UKIDSSJ182217.97-160436.0, UKIDSSJ182217.97-160435.9). The first is unlikely based on  \nits angular separation ∼ 6 . 7 '' ± 0 . 7 '' ; at the magnetar's distance 1 . 6 ± 0 . 3 kpc, this corresponds to ∼ 0 . 05 pc which is too large for a bound orbit. While this could be an unbound companion, it was not recovered in the Gaia search, and the absolute 𝐽 -𝐻 color is inconsistent with OB stars at the 5 𝜎 level when placed at the magnetar distance, with 𝐽 -𝐻 ≈ 0 . 85 ± 0 . 10 (Chrimes et al. 2022b). The latter source was identified with 𝑝 < 5% only for the UKIDSS counterparts ( 𝑝 = 1% , 2%), while the 2MASS source had a marginal 𝑝 ≈ 5%. Again, we find it unlikely this is associated with SGR J1822 based on its angular separation ∼ 9 . 2 '' ± 0 . 7 '' and color 𝐽 -𝐻 ≈ 1 . 69 ± 0 . 14. We conclude that both 2MASS18221839-1604241 and 2MASS18221797-1604360 are false positives for bound companions.  \nThe 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 SGRJ1822.3-1606 may be a Be star X-ray binary, like SGR 07552933, when initially discovered by e.g. Halpern (2011), Gogus & Gorosabel (2011), Gorosabel et al. (2011). The source's Gaia counterpart, Gaia DR34097826476059022848, 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:  \n𝐿 𝑋 ( 𝑡 ) = 𝐿 0 𝐵 𝑑 10 13 G ! 𝑒 -𝑡/𝜏 𝑑 (D1)  \nwhere 𝐵 𝑑 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:  \n𝜏 𝑑 = 𝜏 𝑑 0 𝐵 𝑑 10 13 G ! -𝛿 (D2)  \nwhere 𝜏 𝑑 0 = 500 kyr. Plugging this into Equation D2 and solving for 𝑡 = T decay , which we take to be the present-day:  \nT decay = 𝜏 𝑑 0 𝐵 𝑑 10 13 G ! -𝛿 ln 𝐿 𝑋 ( 𝑡 = T decay ) 𝐿 0 10 13 G 𝐵 𝑑 ! (D3)  \nwhere 𝐿 𝑋 ( 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.  \nThe 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:  \n𝑅 FE 1 km ≈ (7090 km / s) √︄ 𝐸 10 51 ergs 𝑀 ⊙ 𝑀 T 1 s (D4)  \nwhere 𝐸 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:  \nT ch = 𝑀 5 / 6 𝐸 1 / 2 ( 𝜇𝑛 𝑒, 0 ) 1 / 3 (D5)  \nwhere 𝜇 ≈ 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  \n𝑅 ch = 𝑀 𝜇𝑛 𝑒, 0 ! 1 / 3 (D6)  \n(Truelove & McKee 1999; Sedov 2018).  \nBeyond this timescale and radius, the system enters the SedovTaylor 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:  \n(D7) 𝑅 ST 1 km ≈ 𝑅 ch 1 km ! T 1 s ! -2 / 5 = (7090 km / s) √︄ 𝐸 10 51 ergs 𝑀 ⊙ 𝑀 T ch 1 s ! -3 / 5 T 1 s ! -2 / 5  \nwhich 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.", 'APPENDIX E: PAST SNR AND NEUTRON STAR COMPANION SEARCH THRESHOLDS': 'Here we list references for the previous surveys used for comparison in Figure 6. For past radio SNR searches (see Green 2019), we include the Canadian Galactic Plane Survey (CGPS; Gerbrandt et al. 2014; Reich et al. 1990; Kothes et al. 2006, bottom ), the Parkes Galactic Plane Survey (PGPS; Duncan et al. 1997, bottom left ), the Parkes-MIT-National Radio Astronomy Observatory (PMN; Stupar et al. 2007; Condon et al. 1993, bottom left ), the 20 cm VLA supernova remnant search (VLA92; Sramek et al. 1992, bottom middle ), the Multi-Array Galactic Plane Imaging Survey (MAGPIS; Helfand et al. 2006, bottom middle ), The HI, OH, Recombination line (THOR) and VLA Galactic Plane Survey (VGPS) combined survey (THOR+VGPS; Bihr et al. 2016; Anderson et al. 2017, bottom middle ), and the Molonglo Observatory Synthesis Telescope (MOST) Galactic Plane survey (MGPS; Whiteoak & Green 1996, bottom right ). We include limits from the Naval Observatory Merged Astrometric Dataset (NOMAD; Zacharias et al. 2004, 2005) and Tycho-2 (Høg et al. 2000b) point source and proper motion catalogs ( middle ), the former of which Kochanek (2018) use to search for stellar companions of the Crab and Cas A pulsars. We also include limits from the Asteroid Terrestrial-impact Last Alert System (ATLAS) Photometric All-Sky Survey (APASS; ) and the HST used by Kochanek (2021) (in conjunction with 2MASS, UKIDSS, PS1, and Gaia ) and Chrimeset al. (2022b), respectively, to search for pulsar and magnetar companions.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure D1. X-ray Luminosity ( 𝐿 𝑋 ) vs. Dipole Magnetic Field Strength ( 𝐵 𝑑 ) for Known Magnetars using the Decay Model of Ferrario & Wickramasinghe (2006, 2008). The magnetic field estimate is derived from the pulse period 𝑃 and spindown rate / 𝑃 using 𝐵 𝑑 ≈ (3 . 2 × 10 19 𝐺 ) √︁ 𝑃 / 𝑃/ 1 s as reported in the McGill catalog (Olausen & Kaspi 2014). The X-ray luminosity is derived from the 2 -10 keV flux and reported distance, also from the McGill catalog. The colorbar corresponds to the decay age T decay derived from Equation D3, and lines of constant age are shown in grey. This estimate is preferred to the spindown age, T rot, which likely overestimates magnetar ages through neglect of field decay and glitching (e.g. Nakano et al. 2015; Kaspi & Beloborodov 2017).  \n<!-- image -->'}