# The Galactic Population of Magnetars: A Simulation-based Inference Study

M. Sautron, A E Mcewen, G. Younes, J. Pétri, P. Beniamini, D.

Huppenkothen

► To cite this version:

M. Sautron, A E Mcewen, G. Younes, J. Pétri, P. Beniamini, et al.. The Galactic Population of Magnetars: A Simulation-based Inference Study. The Astrophysical Journal, 2025, 986 (1), pp.88. 10.3847/1538-4357/add0aa. insu-05115028

HAL Id: insu-05115028

https://insu.hal.science/insu-05115028v1

Submitted on 16 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# The Galactic Population of Magnetars: A Simulation-based Inference Study

M. Sautron<sup>1,8</sup> ID, A. E. McEwen<sup>2,8</sup> ID, G. Younes<sup>2,3</sup> ID, J. Pétri<sup>1</sup> ID, P. Beniamini<sup>4,5,6</sup> ID, and D. Huppenkothen<sup>7</sup> ID

<sup>1</sup>Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France; matteo.sautron@gmail.com

<sup>2</sup>Center for Space Sciences and Technology, University of Maryland, Baltimore County, Baltimore, MD 21250, USA; almcewen0@hotmail.com

<sup>3</sup> Astrophysics Science Division, NASA/GSFC, Greenbelt, MD 20771, USA

<sup>4</sup>Department of Physics, The George Washington University, Washington, DC 20052, USA

<sup>5</sup>Department of Natural Sciences, The Open University of Israel, P.O. Box 808, Ra'anana 4353701, Israel

<sup>6</sup> Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, P.O. Box 808, Ra'anana 4353701, Israel

<sup>7</sup>Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098XH Amsterdam, The Netherlands

Received 2025 March 14; revised 2025 April 18; accepted 2025 April 23; published 2025 June 9

## Abstract

Population synthesis modeling of the observed dynamical and physical properties of a population is a highly effective method for constraining the underlying birth parameters and evolutionary tracks. In this work, we apply a population synthesis model to the canonical magnetar population to gain insight into the parent population. We utilize simulation-based inference to reproduce the observed magnetar population with a model that takes into account the secular evolution of the force-free magnetosphere and magnetic field decay simultaneously and self-consistently. Our observational constraints are such that no magnetar is detected through their persistent emission when convolving the simulated populations with the XMM-Newton EPIC-pn Galactic plane observations, and that all of the ~30 known magnetars are discovered through their bursting activity in the last ~50 yr. Under these constraints, we find, within 95% credible intervals, the birth rate of magnetars to be $1.8_{-0.6}^{+2.6}$ kyr$^{-1}$, leading to having $10.7_{-4.4}^{+18.8}\%$ of neutron stars born as magnetars. We also find a mean magnetic field at birth ($\mu_b$ is in T) $\log(\mu_b) = 10.2_{-0.2}^{+0.1}$, a magnetic field decay slope $\alpha_d = 1.9_{-1.3}^{+0.9}$, and timescale $\tau_d = 17.9_{-14.5}^{+24.1}$ kyr, in broad agreement with previous estimates. We conclude this study by exploring detection prospects: an all-sky survey with XMM-Newton would potentially allow around seven periodic detections of magnetars to be obtained, with approximately 150 magnetars exceeding XMM-Newton's flux threshold, and the upcoming AXIS experiment should allow these detections to be doubled.

Unified Astronomy Thesaurus concepts: Magnetars (992); X-ray bursts (1814); X-ray sources (1822)

## 1. Introduction

Magnetars are a peculiar subgroup of neutron stars, with magnetic fields much stronger than typical neutron stars ($10^9–10^{11}$ T, while normal pulsars have magnetic fields closer to $10^8$ T). Given that such B-fields are unattainable in terrestrial laboratories, magnetars are the only means for studying the effects of ultra-high magnetic fields on matter, radiation, and fundamental physical processes. Naturally, questions about the formation of such neutron stars emerge. One explanation utilizes a dynamo process during the birth of the neutron star, which dramatically increases the magnetic field after core collapse (R. C. Duncan & C. Thompson 1992; P. Barrère et al. 2025). Another explanation (known as the fossil field hypothesis) suggests that the magnetic field of the progenitor star was already large and subsequently amplified during the core-collapse supernova event by magnetic flux conservation.

Known magnetars are typically young objects, as evidenced by their rapid rotational decay compared to normal pulsars, as well as the placement of some within X-ray-bright, young supernova remnants (J. D. Gelfand & B. M. Gaensler 2007; G. E. Anderson et al. 2012; K. J. Borkowski &

S. P. Reynolds 2017; R. Kothes et al. 2018). There are also magnetars associated with young stellar clusters like CXOU J164710.2-455216 detected in Westerlund 1 (K. Belczynski & R. E. Taam 2008), as well as SGR 1806-20, which is surrounded by a cluster of massive stars (D. F. Figer et al. 2005). They have long spin periods ($P \approx 1–15$ s) and large spin period derivatives ($\dot{P} \approx 10^{-15} – 10^{-11}$ s s$^{-1}$). Under the assumption of a rotating dipole in a vacuum, these rotational characteristics indeed translate to dipole fields of the order of $10^{11}$ T measured at the equator. Most known magnetars are bright X-ray emitters, with X-ray luminosities in the range of $10^{25}–10^{28}$ W, which at times exceed the corresponding rotational energy losses $\dot{E}$ by a few orders of magnitude. Hence, the emission process from these objects has been demonstrated to be driven by the dissipation of their extremely strong magnetic fields (C. Kouveliotou et al. 1994), rather than being rotationally powered. Additionally, a unique trait of magnetars is their bursting properties, rarely shared with any other member of the neutron star population. These come in differing flavors. The short bursts are characterized by a 0.2 s duration and hard X-ray energies reaching $10^{34}$ J (e.g., A. C. Collazzi et al. 2015). On the other end, there are the most extreme events known as giant flares, where up to $10^{40}$ J are released by a magnetar in about 0.5 s in the hard X-ray, soft gamma-ray band (e.g., D. M. Palmer et al. 2005). At times of bursting activity, magnetars may enter weeks- to years-long episodes during which their persistent emission increases (F. Coti Zelati et al. 2018). These episodes are known as outbursts.

<figure><img src="image_3.png" /></figure>

<footnote n="8"><sup>8</sup> The two first authors contributed equally.</footnote>
Currently, there are about 27 confirmed magnetars, with 3 potential candidates (S. A. Olausen & V. M. Kaspi 2014). These are overwhelmingly discovered during bursting periods (unlike canonical pulsars, which are mostly discovered via their persistent radio emission in large-scale surveys; e.g., R. N. Manchester et al. (2001), K. Stovall et al. (2014), J. L. Han et al. (2021), A. E. McEwen et al. (2024), etc.). This is largely due to the lack of large field-of-view, sensitive X-ray surveys capable of identifying magnetar-like spin periods without the need for deep follow-up observations.

Studying magnetars can also deepen our understanding of magnetic field evolution in neutron stars, providing crucial insights into the life cycles of these stars and the role of magnetic fields in their rotational evolution and outburst activities (V. M. Kaspi & A. M. Beloborodov 2017). Additionally, magnetars are key to understanding the connection between neutron star populations, such as the possible evolutionary links between high-B pulsars, X-ray-dim isolated neutron stars (XDINSs), and central compact objects (CCOs), which could help to solve the birth rate problem since the core-collapse supernova (CCSN) rate measured in the Galaxy today, $\beta_{CCSN}$, is estimated to be 11.76–21.28 kyr$^{-1}$ within a 68.3% confidence interval (K. Rozwadowska et al. 2021). This is smaller than the sum of all the estimates of birth rates of isolated neutron stars (A. A. Gençali & Ü. Ertan 2024). Furthermore, recent detections of ultralong-period magnetars (ULPMs) such as GLEAM-X118J162759.5-523504.3 might allow us to increase the number of detected magnetars and even broaden the definition of what is a magnetar (P. Beniamini et al. 2023). Moreover, it is also plausible that a small fraction of millisecond magnetars produce gamma-ray bursts (GRBs; V. V. Usov 1992; Z. G. Dai & T. Lu 1998; D. Kasen & L. Bildsten 2010; B. D. Metzger et al. 2011, 2018; A. C. Trigg et al. 2024), although their association as central engines of GRB jets is challenged on theoretical and observational grounds (P. Beniamini & W. Lu 2021; R. Ricci et al. 2021; Z.-F. Wu et al. 2025). They can also be responsible for fast radio bursts (FRBs), the brightest events in the Galaxy, as evidenced by the association of FRB 20200428 with the known Galactic magnetar SGR 1935+2154 (C. D. Bochenek et al. 2020; W. Lu et al. 2020; B. Margalit et al. 2020). While there is not yet conclusive proof, various lines of evidence indicate that at least some, and perhaps all, extragalactic FRBs are associated with magnetars (S. B. Popov & K. A. Postnov 2010; P. Kumar et al. 2017; J. M. Cordes & S. Chatterjee 2019; A. C. Gordon et al. 2023; T. Totani & Y. Tsuzuki 2023; P. Beniamini & P. Kumar 2023, 2025). Furthermore, some FRBs originate in old stellar populations (M. Bhardwaj et al. 2021; F. Kirsten et al. 2022; F. H. Panther et al. 2023), which may challenge the magnetar scenario or suggest alternative formation channels, such as accretion-induced collapse of a white dwarf or compact object mergers (B. Margalit et al. 2019; W. Lu et al. 2022). Recently, Ó. Rodríguez et al. (2024) analyzed a sample of 54 stripped-envelope supernovae spanning all subtypes (SNe IIb/Ib/Ic) and found evidence for the presence of a central engine in most cases. If these engines are indeed magnetars, their inferred initial magnetic fields and spin periods are consistent with those expected at birth. This suggests that stripped-envelope supernovae could represent the formation events of objects known as magnetars. Their violent outbursts and persistent X-ray emission also provide constraints on magnetothermal evolution models,

shedding light on the interplay between magnetic field decay, thermal conductivity, and neutron star cooling. In this sense, magnetars serve as natural laboratories for probing the complex interactions between magnetism, gravitation, and particle physics, offering valuable clues about the behavior of matter under extreme conditions. Also, understanding the magnetar population may lend clues to the links between them and events in the Universe, as well as newly discovered populations resembling magnetars.

Neutron star population synthesis is an efficient tool that allows us to extract fundamental information about this population that would be unattainable otherwise, such as their emission processes or the parameters characterizing the population (e.g., birth rate, magnetic decay timescale, and mean magnetic field at birth, among others). In this type of study, neutron stars are generated from birth and evolved up to the present time. Once the sample is generated, detection criteria are applied to determine whether they can be observed with current and future instruments. Recently, V. Graber et al. (2024) and J. Berteaud et al. (2024) applied simulation-based inference (SBI; K. Cranmer et al. 2020) to pulsar population synthesis, demonstrating its effectiveness for such studies. By relying on simulations rather than analytical likelihoods (which may not be feasible or even possible to construct), SBI enables robust parameter estimation and model comparison.

Studies constraining the characteristic parameters of the population of magnetars have been done; for instance, see M. Gullón et al. (2015), P. Beniamini et al. (2019), and J. A. Jawor & T. M. Tauris (2022), where they constrained the mean magnetic field at birth, the birth rate, and the magnetic decay timescale of magnetars. This paper combines magnetar evolution with the detection statistics while using simulation-based inference in a complete way that has not been included in previous studies. As such, finding these characteristic parameters will both improve models for magnetar evolution and provide insight for future searches.

This paper is organized as follow: in Section 2, we describe how magnetar populations are generated, followed by an explanation of the magnetar evolution model in Section 3. Section 4 details the periodicity search conducted for this study, and Section 5 discusses the two types of detections considered in this work. Section 6 provides an overview of the simulation-based inference method, while Section 7 focuses on validating the entire pipeline. The results are presented in Section 8, discussed in Section 9, and summarized in the conclusions in Section 10.

## 2. Generation of Magnetars

We generate a population of magnetars whose ages are taken from a random uniform distribution between 0 and the maximum age of $5 \times 10^6$ yr (this is a conservative limit, as no observed magnetars have a characteristic age older than $3.6 \times 10^6$ yr and most are likely significantly younger). In each simulation, the number of sources is characterized by a birth spacing $X$, where $X$ is chosen as an integer number of years between the births of subsequent magnetars (i.e., $\frac{5 \times 10^6}{X}$ sources are drawn for each simulation). In this paper, we report this quantity as either the birth spacing $X$ in years or the birth rate $BR = (1/X)$ yr$^{-1}$. This approach can be used to verify whether a constant birth rate is realistic, by comparing the results from the simulation to observation. Studies by
P. Beniamini et al. (2019) and J. A. Jawor & T. M. Tauris (2022) found the birth rate of magnetars to be in the range 2.3–20 kyr⁻¹, corresponding to the birth spacing of 50–434 yr; the value of X is drawn from this range.

To describe the position of the magnetars in the Galaxy, we use the right-handed Galactocentric coordinate system (x, y, z), with the Galactic center at its origin, y increasing in the disk plane toward the location of the Sun, and z increasing toward the direction of the north Galactic pole. The initial spatial distribution of the magnetars is given by the radial spread ρ(R) from I. Yusifov & I. Küçük (2004), and the spread in altitude ρ_z(z) is given by S. A. Olausen & V. M. Kaspi (2014):

$$ \rho(R) = A \left( \frac{R + R_1}{R_{\odot} + R_1} \right)^a \exp \left( -b \left( \frac{R - R_{\odot}}{R_{\odot} + R_1} \right) \right) \quad (1a) $$

$$ \rho_z(z) = \frac{e^{-|z+z_{\odot}|/h_e}}{h_e}, \quad (1b) $$

where R is the axial distance from the z-axis and z is the distance from the Galactic disk. The numerical values for the constants are A = 37.6 kpc⁻², a = 1.64, b = 4.0, R₁ = 0.55 kpc, h_e = 30.7 pc, z_☉ = 15 pc, and R_☉ = 8.5 kpc. The purpose of choosing Equation (1a) to describe the radial spread is to put magnetars' birth positions within the spiral arms of the Galaxy. More precisely, the Galactic spiral structure contains four arms with a logarithmic shape function, and the azimuthal coordinate φ is calculated as a function of the distance from the Galactic center,

$$ \phi(R) = k \ln \left( \frac{R}{r_0} \right) + \phi_0. \quad (2) $$

The values of the models describing each arm are given in Table 1. These values were used by M. Ronchi et al. (2021) and taken from J. M. Yao et al. (2017) to match the shape of the arms of the Galaxy. The Local Arm is not modeled, as it is far less dense than the four other arms. Each star has an equal probability to be in one of the four arms with angular coordinate φ calculated using (2) for a given R. Furthermore, the Galaxy is not static, and its arms are moving with an approximated period T = 250 Myr (D. M. Skowron et al. 2019). Because the Galaxy is rotating in the clockwise direction, the age of a magnetar will then dictate its angular position at birth. Following the same procedures described in M. Ronchi et al. (2021), noise is added to both coordinates R and φ in order to avoid artificial features near the Galactic center. This correction factor is calculated via

$$ \phi_{\text{corr}} = \phi_{\text{rand}} \exp(-0.35 R), \quad (3) $$

where φcorr is randomly drawn from a uniform distribution between 0 and 2π. Similarly, a correction factor for the radial distance rcorr is drawn from a normal distribution with a mean of 1 and a standard deviation σcorr = 0.07. These corrections are taken into account to compute the values φ and R, respectively:

$$ R_{\text{birth}} = R \times r_{\text{corr}} \quad (4) $$

and

$$ \phi_{\text{birth}} = \phi + \phi_{\text{corr}} + \frac{2 \pi t_{\text{age}}}{T}, \quad (5) $$

where tage is the age of the magnetar. These coordinates are then converted to Galactocentric coordinates (x, y).

Next, the rotational properties of each magnetar are produced, including the angle between rotation and magnetic axes (inclination angle) α₀, the birth spin period P₀, and birth magnetic field B₀ (throughout this study, the term “magnetic field” always denotes the surface magnetic field). The birth inclination angle α₀ is assumed to follow an isotropic distribution generated from a uniform distribution U ∈ [0, 1] and given by α₀ = arccos(2 U − 1). In this work, P₀ and B₀ both follow log-normal distributions, as suggested from the results of a study of 56 young neutron stars by A. P. Igoshev et al. (2022). These distributions are used in other magnetar population synthesis studies, including in J. A. Jawor & T. M. Tauris (2022), where a log-normal distribution is used for the spin period at birth, as well as in N. Rea et al. (2015), where a log-normal distribution is used for the magnetic field at birth.

Explicitly, the probability distributions for B₀ and P₀ are given by

$$ p(\log(B_0)) = \frac{1}{\sigma_b \sqrt{2\pi}} e^{-(\log B_0 - \log \mu_b)^2 / (2\sigma_b^2)} \quad (6) $$

$$ p(\log(P_0)) = \frac{1}{\sigma_p \sqrt{2\pi}} e^{-(\log P_0 - \log \mu_p)^2 / (2\sigma_p^2)}, \quad (7) $$

where μb,p and σb,p are the mean and standard deviation of the log-normal distributions for magnetic field and period.

Studies have posited several mechanisms (asymmetry in the emission of neutrinos, hydrodynamic asymmetry during the collapse, etc.) to be at the origin of the asymmetry in the supernova explosion (E. Tademaru & E. R. Harrison 1975; R. Cowsik 1998; H. Spruit & E. S. Phinney 1998; T. Foglizzo et al. 2015). These mechanisms seem to result in an alignment between the kick velocities at birth and the rotation axis of the neutron star, as suggested by J. M. Rankin (2007) and A. Noutsos et al. (2013). Therefore, this alignment is taken into account in the simulation. Moreover, according to G. Hobbs et al. (2005), a Maxwellian distribution replicates well the observations of pulsars velocities, and more recently, H. Ding et al. (2024) showed that the velocity distribution of magnetars is broadly similar to that of normal pulsars, albeit slightly lower. A description based on supernova (SN) simulations is provided in V. Kapil et al. (2023), but was not implemented here, as it would have required progenitor parameters for the simulated neutron stars, which are not considered in this work. Additionally, here, we assume that magnetars originate from single-star progenitors, neglecting potential binary interactions that could affect their velocity distribution (P. Beniamini & T. Piran 2024). Thus, in this work, magnetar velocities at birth are given by

$$ p(v) = \sqrt{\frac{2}{\pi} \frac{v^2}{\sigma_v^3}} \exp \left( -\frac{v^2}{2\sigma_v^2} \right). \quad (8) $$

The mean velocity of the distribution v̄ is related to the standard deviation by v̄ = σv√(8/π) with a standard deviation σv = 265 km s⁻¹ for magnetars. The velocity vector is then
distributed along the unit vector $\mathbf{n}_{\Omega}$ of the rotation axis, which is drawn from an isotropic distribution. The Cartesian coordinates of the unit rotation vector $\mathbf{n}_{\Omega}$ are $(\sin \theta_{n_{\Omega}} \cos \phi_{n_{\Omega}}, \sin \theta_{n_{\Omega}} \sin \phi_{n_{\Omega}}, \cos \theta_{n_{\Omega}})$.

### 3. Evolution of Magnetars

#### 3.1. Evolution of the Spin Period, the Magnetic Field, and the Inclination Angle

The magnetar's initial period $P_0$, inclination angle $\alpha_0$, and magnetic field $B_0$ are evolved in time in a fully self-consistent way that takes into account spin-down losses and internal magnetic field dissipation within the neutron star. The neutron star period evolves according to the force-free magnetosphere model, which takes into account the electric current and charge flowing inside the magnetosphere. The spin-down luminosity $\dot{E}$ is

$$ \dot{E} = \frac{dE_{\text{rot}}}{dt} = -I \Omega \dot{\Omega} = L_{\text{ffe}} \quad (9a) $$

$$ L_{\text{ffe}} = \frac{4\pi R^6 B^2 \Omega^4 (1 + \sin^2 \alpha)}{\mu_0 c^3}, \quad (9b) $$

where $I \approx 10^{38}$ kg m$^{-2}$ is the neutron star moment of inertia, $\Omega = 2\pi/P$ the rotation frequency of the pulsar ($P$ is the spin period), $\dot{\Omega}$ is its time derivative, $R = 12$ km is the typical radius of a neutron star as found by recent NICER observations (T. E. Riley et al. 2019; S. Bogdanov et al. 2019), $B$ is the magnetic field at the equator, $\alpha$ is the inclination angle, $\mu_0 = 4\pi \times 10^{-7} \text{H m}^{-1}$ is the vacuum permeability constant, and $c$ is the speed of light. The term $L_{ffe}$ for the spin-down luminosity was given by A. Spitkovsky (2006) and J. Pétri (2012). Equation (9a) clearly shows the correlation between the obliquity, the magnetic field, and the rotation frequency. Combining Equations (9a) and (9b) leads to

$$ \dot{\Omega} = -K_{\text{ffe}} \Omega^n \quad (10a) $$

$$ K_{\text{ffe}} = \frac{4\pi R^6 B^2 (1 + \sin^2 \alpha)}{\mu_0 c^3 I}, \quad (10b) $$

where $n$ is the braking index ($n=3$ for magnetic dipole radiation).

The most speculative part of our work regards the evolution of the inclination angle, as it is still a hotly debated issue. From a theoretical point of view, the study of A. Philippov et al. (2014) suggests that the force-free model allows for an evolution of $\alpha$ different from the exponential decay, that applies when considering a neutron star in vacuum (as applied to observed data by, e.g., J. A. Jawor & T. M. Tauris 2022). Another theory points out that the inclination angle should tend to $\pi/2$, to minimize energy loss (V. S. Beskin et al. 1993; V. S. Beskin & E. E. Nokhrina 2007; V. S. Beskin et al. 2013). However, from an observational point of view, there is no clear evidence in favor of one of these approaches; even when there is a tendency toward alignment with increasing age, there is a large spread (T. M. Tauris & R. N. Manchester 1998; P. Weltevrede & S. Johnston 2008). Finally, the free precession of neutron stars studied in A. Melatos (2000), L. Arzamasskiy et al. (2015), and A. P. Igoshev & S. B. Popov (2020) shows that there should be oscillations but long-term decay in the evolution of $\alpha$, which could explain why we

observe neutron stars that seem to have an increasing inclination angle (A. G. Lyne et al. 2015). Since we are considering the force-free model in this work, we will use the integral of motion between $\Omega$ and $\alpha$ found in A. Philippov et al. (2014), which reads

$$ \Omega \frac{\cos^2 \alpha}{\sin \alpha} = \Omega_0 \frac{\cos^2 \alpha_0}{\sin \alpha_0}, \quad (11) $$

where quantities with subscript 0 indicate their initial value and those without a subscript are the value at present time.

The evolution of neutron stars' magnetic fields is also debated. Few pulsar population syntheses have claimed that the magnetic field of neutron stars is not evolving (D. Bhattacharya et al. 1992; C.-A. Faucher-Giguère & V. M. Kaspi 2006). However, as discussed earlier, magnetars cannot be powered by spin-down. They are therefore highly likely to be powered by magnetic field decay. Several phenomena strengthen the idea that the magnetic field should decay: ohmic dissipation in the crust, ambipolar diffusion in the core, vortex movements in the crust, internal turbulence, and the loss of energy via radiation emission (A. Cumming et al. 2004; S. Dall'Osso et al. 2012). In the specific case of magnetars, the X-rays could not be powered by spin-down and therefore be powered magnetically, implying the decay of the magnetic energy reservoir (V. M. Kaspi & A. M. Beloborodov 2017). In addition, a recent study, which takes into account ohmic heating (A. P. Igoshev & S. B. Popov 2024), period derivative measurements, and spectral information about the neutron star RX J0720.4-3125, has suggested that its magnetic field is decaying rapidly. Furthermore, it has been shown in other and more recent pulsar population studies that a decaying magnetic field better reproduces the observed pulsar/magnetar population (P. Beniamini et al. 2019; J. A. Jawor & T. M. Tauris 2022; M. Sautron et al. 2024). Therefore, we prescribe a magnetic field decay according to a power law, i.e.,

$$ B(t) = B_0 \left( 1 + \frac{\alpha_d t}{\tau_d} \right)^{-1/\alpha_d}, \quad (12) $$

where $\alpha_d$ is a constant parameter controlling the speed of the magnetic field decay and $\tau_d$ the typical decay timescale.

In line with the spin-down luminosity, the inclination angle satisfies another evolution equation that after integration was found by A. Philippov et al. (2014) for a spherically symmetric neutron star with a constant magnetic field. With our decaying prescription, the inclination angle $\alpha$ is found by solving for the root of

$$ \ln(\sin \alpha_0) + \frac{1}{2 \sin^2 \alpha_0} + \frac{K \Omega_0^2 \tau_d B_0^2 \cos^4 \alpha_0}{\alpha_d - 2 \sin^2 \alpha_0} \times \left[ \left( 1 + \frac{\alpha_d t}{\tau_d} \right)^{1-2/\alpha_d} - 1 \right] = \ln(\sin \alpha) + \frac{1}{2 \sin^2 \alpha}, \quad (13) $$

where $t$ is the age of the pulsar and $K = 4\pi R^6 / I\mu_0 c^3$. The typical decay timescale for a mean magnetic field of $10^{10}$ T is $3 \times 10^4$ yr in D. Viganò (2013), for instance. However, as we will show in Section 6, we do not assume that $\tau_d$ depends on magnetic field; instead, we use statistical inference to determine $\tau_d$.
## 3.2. Description of the Galactic Potential

Finally, the motion of the magnetar within the Galactic potential is calculated. Each magnetar evolves in the gravitational potential $\Phi$ subject to an acceleration $\ddot{x}$ according to

$$ \ddot{x} = -\nabla\Phi. \quad (14) $$

Equation (14) is integrated numerically thanks to a Position Extended Forest Ruth-Like (PEFRL) algorithm (I. P. Omelyan et al. 2002), a fourth-order integration scheme (see Appendix A of M. Sautron et al. 2024 for a detailed analysis of the PEFRL scheme and its accuracy). We now discuss the Galactic potential model used.

The Galaxy is divided in four distinct regions with different mass distributions and associated gravitational potentials. The four potentials are: the bulge $\Phi_b$, the disk $\Phi_d$, the dark matter halo $\Phi_h$, and the nucleus $\Phi_n$. The total potential of the Milky Way $\Phi_{tot}$ is the sum of these potentials:

$$ \Phi_{tot} = \Phi_b + \Phi_d + \Phi_h + \Phi_n. \quad (15) $$

The expressions for these potentials are taken from A. T. Bajkova & V. V. Bobylev (2021) with parameters given in Table 2. The nucleus mass was found in J. Bovy (2015).

The potential for the bulge $\Phi_b$ and for the disk $\Phi_d$ both have the forms proposed by M. Miyamoto & R. Nagai (1975), which are typically used for models of the gravitational potential of the Milky Way. They read

$$ \Phi_i(R, z) = - \frac{GM_i}{[R^2 + (a_i + \sqrt{z^2 + b_i^2})^2]^{1/2}}, \quad (16) $$

where $R^2 = x^2 + y^2$, $i = b$ is for the bulge, and $i = d$ is for the disk. In those formulae, $R$ depends on the coordinates $x$ and $y$, $a_i$ and $b_i$ are the scale parameters of the components in kpc, $M_i$ is the mass (for the disk or the bulge), and $G$ is the gravitational constant. The values of these constants can be found in Table 2.

The potential for the dark matter halo $\Phi_h$, according to the frequently used potential of J. F. Navarro et al. (1997), is expressed as

$$ \Phi_h(r) = -\frac{GM_h}{r} \ln\left(1 + \frac{r}{a_h}\right), \quad (17) $$

where $M_h$ is the mass of the halo, $r = x^2 + y^2 + z^2$, and $a_h$ is the length scale whose values are found in Table 2.

Finally, for the last part of the Galaxy, the nucleus of the Milky Way is simply represented by a Keplerian potential $\Phi_n$, such as

$$ \Phi_n(r) = -\frac{GM_n}{r}, \quad (18) $$

where $M_n$ is the mass of the nucleus.

## 3.3. X-Ray Emission Mechanism

X-ray emission from magnetars is usually associated with resonant Compton scattering (RCS) or thermal emission (N. Rea et al. 2008; S. Zane et al. 2009; M. Gullón et al. 2015; N. Rea et al. 2015). For the former, electrons in the magnetosphere absorb nearby photons and re-emit them as X-rays afterward. In this work, we focus exclusively on the dipolar magnetic energy losses ($\dot{E}_B$), which account for a

fraction of the total magnetic energy losses ($0.1 \le f_{dip} \le 0.3$). Most of the magnetic energy is stored in the toroidal component of the magnetic field, which is not considered here. The persistent X-ray luminosity of the star appears to represent a fraction $f_{per} < 1$ of the total magnetic energy losses (P. Beniamini et al. 2025). Based on current understanding, it is reasonable to assume that $f_{dip} = f_{per} = 1/3$, leading to the following estimation for the persistent X-ray luminosity of a simulated magnetar: $L_X \approx \dot{E}_B$. Furthermore, the dipole magnetic energy ($E_B$) and the dipole magnetic energy losses can be written as

$$ E_B = \frac{2\pi R^3 B^2}{3\mu_0}, \quad (19a) $$

$$ \dot{E}_B = \frac{4\pi R^3 B \dot{B}}{3\mu_0}, \quad (19b) $$

where $R$ is the radius of the neutron star (12 km in our simulations), $\mu_0$ is the permeability constant, and $B$ is the magnetic field. This luminosity is then converted into an unabsorbed flux,

$$ S = \frac{L_X}{4\pi d^2}, \quad (20) $$

where $d$ is the distance of the magnetar to us. By assuming blackbody emission and under an assumption of the magnetar's emitting surface fraction in X-rays, the surface temperature can be calculated as such in the simulation:

$$ T = \left( \frac{L_X}{\xi \times 4\pi R^2 \sigma} \right)^{1/4}, \quad (21) $$

where $\xi$ is the emitting surface fraction in X-rays, $\sigma \approx 5.67 \times 10^{-8} \text{ W m}^{-2} \text{ K}^{-4}$ is the Stefan-Boltzmann constant, and $T$ is the temperature. This computation is very complicated observationally, because of the deviation of the emission from a pure blackbody spectrum, the complex emission geometry on the surface, and the poorly known distances to the sources. Nonetheless, we explored varying estimates of the emitting surface fraction for our simulated magnetars, which were considered detected in the process. This is discussed further in Section 8.

# 4. Observations and Data Reduction

A crucial component of the simulation process is the comparison of the summary statistics measured from each simulated population to equivalent measurements from the “true” population of magnetars. One of these statistics comes from a search for new magnetars within the observed XMM-Newton sky; in particular, the number of magnetars discovered by an untargeted search. As there is no complete survey of the Galactic plane in X-rays, archival XMM data were utilized as an ad hoc survey; a periodicity search was conducted on point sources identified in the 4XMM-DR13 catalog (N. A. Webb et al. 2020). Unlike modern surveys for pulsars, where large regions of the Galaxy are uniformly searched to identify new sources and place constraints on the entire population, the XMM observations cover a far smaller region of the sky. Additionally, these observations are primarily targeted observations of previously identified sources, as opposed to “blind” searches. For these reasons, the effective “survey” used is very
uneven and biased. Nevertheless, searching for sources with well-understood sky coverage can begin to inform population and evolution models for the Galactic magnetar population.

Our search is based on the 4XMM-DR13 catalog. We exclude all known bright sources (e.g., the Crab pulsar) and only consider sources that are within ±5° of the Galactic plane; magnetars age rapidly and are not expected to be detected beyond such a high Galactic latitude. Also, we only consider sources with a total number of counts ≥200. Finally, we exclude observations in small window and timing mode.

These restrictions reduced the number of sources to search to 26,073 sources contained in 1584 observations. For each candidate source, we performed data reduction using tools in the XSPEC (K. A. Arnaud 1996) and XMM-Newton Science Analysis System (XMM-SAS) version 21.0.0 tools. First, the raw observation data in which the candidate was observed are downloaded. Using the full 4XMM catalog, all known sources (aside from the candidate of interest) are extracted from the data set using the “EXTENT” parameter as the extraction radius to create a “Swiss cheese” image; this is the background for each observation. High-energy (≥10 keV) background sources associated with proton flaring are also removed. The 4XMM catalog provides observed count rates in several energy bands: 0.2–0.5 keV (band 1), 0.5–1 keV (band 2), 1–2 keV (band 3), 2–4.5 keV (band 4), 4.5–12 keV (band 5), and the full band (band 8). To determine the best band for the search, the signal-to-noise (S/N) for each band and all contiguous combinations were calculated using the following expression:

$$S/N = \frac{\sum_{i} r_{i}}{\sum_{i} \sigma_{i}} \sqrt{\tau_{\text{exp}}}, \qquad (22)$$

where $r_i$ is the rate in band i, $\sigma_i$ is the uncertainty on the rate in band i, and $\tau_{\text{exp}}$ is the duration of the exposure. The band with the highest S/N was used for the search, and an image was created using this energy band. Next, a King function was fit to the source to find an optimal radius of extraction:

$$K(r) = A \left( 1 + \left( \frac{r}{r_0} \right)^2 \right)^{-\beta} + DC, \qquad (23)$$

where A is a scaling factor, $r_0$ is a characteristic angular radius, $\beta$ is the index of radial decay, and DC is a constant offset from zero. These four parameters were fit to the radial intensity curve produced by the eradial routine in XMM-SAS. Using this fit, the optimal extraction radius was chosen to include 99% of the area under the curve, and a time-series data set was created from the background-subtracted source region. This created our final “cleaned event list” for each source.

With the time series in hand, a search for periodicity was conducted using an H-test (O. C. de Jager & I. Büsching 2010). The search spans frequencies between 50 mHz and 2 Hz; step sizes are chosen to equal the inverse of the observation duration or 5e-5, whichever is larger. This limit is chosen to prevent excessive computation time for observations lasting over 20 ks. All calculated H-statistics are saved, and significant detections are highlighted. For these candidates, spectral files are also created for additional confirmation. To choose a significance threshold, the same H-test was conducted on simulated time-series data containing only noise, and all resulting H-statistic measurements were binned in a

histogram. From this, a CDF was measured, and the value of the H-statistic below which 99.999% of samples were included was determined to be 47. Using the PDF as measured by O. C. de Jager & I. Büsching (2010), 47 sits even closer to 100%, and so it was accepted as a very conservative limit. The search turned up no new sources, though 14 previously known magnetars were detected.

## 5. Simulated Detection

### 5.1. Magnetar Outbursts

Given that all the known magnetars were discovered through outbursts, modeling this phenomenon in our simulation was very important in order to match the results to the observed population. The process used here is very similar to the one in P. Beniamini et al. (2019), where it is assumed that the burst energy distribution per magnetar is dominated by the most energetic bursts ($dN/dE \propto E^{-\gamma}$ with $1.4 \le \gamma \le 1.8$ (E. Göğüş et al. 1999, 2000), and the total energy released is $\propto E^2 dN/dE$). In addition, the energy released during a GF is on the order of $10^{37}$ J, and the burst rate scales with the normalized magnetic energy loss $\dot{E}_B/E_{B,0}$. We assume

$$\frac{\dot{E}_B(t)}{E_{B,0}} = \frac{2B\dot{B}}{B_0^2}. \qquad (24)$$

For a magnetar with age $5 \times 10^3$ yr and initial magnetic field of $3 \times 10^{10}$ T, giant flares are expected to occur at a rate of approximately $10 \text{ kyr}^{-1}$ for a typical energy of $10^{37}$ J. Given that hard X-ray burst-monitors have been active for approximately 50 yr, for each magnetar in the simulation, we determine whether it has experienced an outburst within the last 50 yr using its burst rate. For these simulations, we assume that all bursts in this timeframe are detected, although doing so might slightly overestimate the number of detections; no burst detection threshold is applied. We explore the implications of these assumptions in the discussion.

### 5.2. Convolution with XMM-Newton Survey and Periodicity Search

After producing a simulated population of Galactic magnetars, we used the XMM coverage of the Galactic plane to determine which (if any) of the sources would be detected by periodicity search. The source positions were matched to the catalog of XMM observations, and sources that fall within the observation footprint were retained. In cases where multiple observations overlap the same position, we used the longest available. Vignetting in XMM observations reduces the effective exposure away from the boresight, and some pixels are corrupted. This information is encoded in the exposure maps, which we produced for each XMM observation using the expmap routine in XMM-SAS. From these, a precise exposure was found for each source. Using the distances to each simulated source and its sky position, we estimated the expected column density of hydrogen ($N_H$) using the 3D-NH tool<fnref n="9" /> (V. Doroshenko 2024). This tool utilizes X-ray, visible, and radio data to estimate absorption along a given line of sight for a provided distance. At this point, we determined the optimal temperature band for a search by simulating source and background count rates in each 4XMM-

<footnote n="9">⁹ http://astro.uni-tuebingen.de/nh3d/nhtool</footnote>
DR13 band by modeling an absorbed blackbody in XSPEC and choosing the band with the highest S/N, per Equation (22).

With the exposure, optimal band, and count rates in hand, we simulated a time series for each candidate source. We drew pulsed fractions from an empirical PDF derived from a double-Gaussian fit to known magnetar pulse fractions (C.-P. Hu et al. 2019). A sinusoidal pulse profile was produced with individual phase bins that are required to be greater or equal to the XMM sampling time. Source counts were distributed across each rotation of the magnetar according to its pulse fraction, and background counts are distributed according to Poisson noise. As was done in the “real” search, we conducted an H-test on the lightcurve with bin sizes corresponding to the inverse of the exposure or 20 µHz, whichever is larger. Because the source spin frequency is known and the simulation is intended to determine simply whether or not the source is detected, only the frequency bins within 5% of the source frequency are used in the search, to save processing time. Any sources detected with H-statistics above 47 (the search threshold) at their spin frequency are considered detections via periodicity search.

# 6. Simulation-based Inference

## 6.1. Concept

To constrain the parameters that characterize the Galactic magnetar population, SBI (K. Cranmer et al. 2020) was utilized. Consider a model M with a vector of parameters $\theta$ that produces data D. In traditional statistical inference, we assume that the data-generating process $p(D|\theta)$, also called the likelihood, can be evaluated analytically. We can then find the probability of the model parameters, also called the posterior probability, via the Bayes theorem:

$$p(\theta|D) = \frac{p(D|\theta)p(\theta)}{p(D)}, \qquad (25)$$

where $p(\theta)$ denotes the prior probability distribution of the parameters, and $p(D)$ is called the marginal likelihood. In parameter estimation problems, where we are interested primarily in the shape of the distribution, we generally only need to consider $p(\theta|D) \propto p(D|\theta)p(\theta)$, since $p(D)$ is a normalization constant and does not depend on the parameters.

Evaluating the posterior requires an analytical form for the likelihood. This may not be available in practice, for example, when the measurement process introduces complex effects that are not easily parameterized, or when the output of the data-generating process is inherently stochastic. In these cases, it may be possible to simulate data sets that are somewhat faithful representations of the real data, but we cannot easily relate the model to the observed data via the likelihood.

In this context, SBI provides principled methods for relating simulated data to the underlying model parameters such that (unbiased) inference becomes possible. In SBI, we use a large number of pairs of parameters $\theta$ and simulated data D to learn a mapping between them to approximate the likelihood (neural likelihood estimation; G. Papamakarios et al. 2018; J. Alsing et al. 2019), the posterior (neural posterior estimation (NPE); G. Papamakarios & I. Murray 2016; S. Mishra-Sharma & K. Cranmer 2022; M. Vasist et al. 2023; V. Graber et al. 2024), or the ratio of likelihood to marginal likelihood (neural ratio estimation; J. Hermans et al. 2019; U. Bhardwaj et al. 2023). Flexible neural network-based density estimators called

normalizing flows, designed to approximate complex probability distributions, are often used to learn the mapping between parameters and data.

Because the data are often high-dimensional and noisy, it can be much easier to first compute summary statistics $x$ of the raw (simulated) data: these summaries can be engineered to explicitly encode properties of the data that are informative with respect to the parameters, or they can be learned automatically, often again using neural networks (V. Graber et al. 2024).

Here, we use NPE, since we are particularly interested in the posterior distributions from our observed data. We generate training data from flat priors for all parameters (see Section 6.2 and Table 3 for details), compute summary statistics (explained in more detail in Section 6.2), and train a normalizing flow to approximate the posterior probability of the population model parameters given the observed sample of magnetars. We use the Python package sbi<fnref n="10" /> (A. Tejero-Cantero et al. 2020) for neural network setup and training.

## 6.2. Simulation Setup

The characterizing model parameters of the Galactic magnetars population are $\vec{\theta} = \{\sigma_b, \mu_b, \sigma_p, \mu_p, \alpha_d, \tau_d, BR\}$ (see Sections 2 and 3 for descriptions). The ranges for the flat priors on the model parameters are given in Table 3. Note that $\alpha_d$ can be negative, and by the definition of the magnetic field decay in Equation (12), this allows for a simulated magnetar to reach a negative magnetic field. Because this is nonphysical, a magnetar that reaches a negative magnetic field is ruled out for detection. These ranges were chosen to explore a parameter space similar to—and even larger than—those of P. Beniamini et al. (2019) and J. A. Jawor & T. M. Tauris (2022).

Since the parameter space is relatively large, the Latin hypercube sampling (LHS) method (M. D. McKay et al. 1979) is used to explore it most efficiently. This method divides the parameter space into a specified number of intervals (depending on how many simulations are made) and samples the model parameters such that only one value is drawn per interval. Therefore, this method ensures the values chosen for the model parameters in different simulations are representative of the whole space.

The simulations each produce eight summary statistics (five vectors and three scalars) to characterize the population. In order to compare the simulated population with the observed population (which has been discovered almost entirely through bursts), the first six statistics are obtained from the magnetars that experienced at least one outburst episode in the last 50 yr when burst monitoring has been active. Five of these statistics are feature vectors from the following density maps (which are 2D histograms of size 32 × 32):

1. $P - \dot{P}$,

2. $P - L_x$,

3. $P - g_b$,

4. $\dot{P} - L_x$, and

5. $\dot{P} - g_b$.

The feature vectors from the density maps of the simulated magnetars detected through outbursts are obtained after the density maps were passed to a CNN, as detailed later in this

<footnote n="10">¹⁰ https://github.com/sbi-dev/sbi</footnote>
<caption>**Table 1**</caption>

Parameters of the Milky Way Spiral Arm Structure, Adapted from J. M. Yao et al. (2017)

<table><thead><tr><th>Arm Number</th><th>Name</th><th>k<br>(rad)</th><th>r<sub>0</sub><br>(kpc)</th><th>φ<sub>0</sub><br>(rad)</th></tr></thead><tbody><tr><td>1</td><td>Norma</td><td>4.95</td><td>3.35</td><td>0.77</td></tr><tr><td>2</td><td>Carina-Sagittarius</td><td>5.46</td><td>3.56</td><td>3.82</td></tr><tr><td>3</td><td>Perseus</td><td>5.77</td><td>3.71</td><td>2.09</td></tr><tr><td>4</td><td>Crux-Scutum</td><td>5.37</td><td>3.67</td><td>5.76</td></tr></tbody></table>

subsection. The sixth population statistic is the slope of the log N-log S curve, dubbed *K*.

The final two summary statistics are related to detection and are therefore compared with the results of the search conducted in this project (as well as the observational history of magnetars):

1. $n_{\text{detect}}$ is the number of magnetars detected by a periodicity search in XMM-Newton data, and

2. $n_{\text{outburst}}$ is the total number of outburst episodes of magnetars in the last 50 yr.

These eight outputs of the simulations are then compared to the same characteristics of the observed data in order to make the best possible match.

The physical properties of magnetars (e.g., *P* and *Ṗ*) are often found to be correlated, mirroring the underlying physical relationships (e.g., of the aging process of magnetars). Consequently, simple summaries of the data such as the mean and the standard deviation of the simulated population are insufficient to capture important relationships with respect to the input parameters. In order to characterize the quantities and their correlation, we generate density maps of *P* − *Ṗ*, *P* − *L*<sub>X</sub>, *P* − *g*<sub>b</sub>, *Ṗ* − *L*<sub>X</sub>, and *Ṗ* − *g*<sub>b</sub>, and then a convolutional neural network (CNN) is used to generate informative features from these maps (each type of map will have its own dedicated CNN for training). These features are 1D vectors that represent the latent features extracted by the CNNs that can be parsed easily by the density estimator.

The CNNs used in this work have the following architecture (see Figure 1; V. Graber et al. 2024 proposed a similar CNN architecture):

1. Two-dimensional convolution layer with kernel size 3 × 3, 3 input channels, 32 output channels, stride 1, padding “same.”

2. Two-dimensional Max pooling layer with size 2 × 2, stride 2, no padding.

3. Two-dimensional convolution layer with kernel size 3 × 3, 32 input channels, 64 output channels, stride 1, padding “same.”

4. Two-dimensional Max pooling layer with size 2 × 2, stride 2, no padding.

5. Fully connected linear layer with the flattened output from the second pooling layer as input and 32 output neurons encoding the latent representation.

After each convolution and the fully connected layer, a rectified linear unit (ReLU) activation function is applied.

The summary statistics (the five vectors from the density maps and the three scalars) of the simulation were then passed to the default density estimator of SBI, a masked autoregressive flow (MAF; G. Papamakarios et al. 2017), a type of

<table><caption>**Table 2**

Constant Values Used for the Different Potentials, with the Solar Mass $M_{\odot} = 1.99 \times 10^{30}$ kg</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>$M_h$</td><td>2.9 × 10<sup>11</sup> ± 7.6 × 10<sup>10</sup>$M<sub>☉</sub></td></tr><tr><td>$a_h$</td><td>7.7 ± 2.1 kpc</td></tr><tr><td>$a_d$</td><td>4.4 ± 0.73 kpc</td></tr><tr><td>$b_d$</td><td>0.308 ± 0.005 kpc</td></tr><tr><td>$M_d$</td><td>6.50 × 10<sup>10</sup> ± 1.9 × 10<sup>9</sup>$M<sub>☉</sub></td></tr><tr><td>$a_b$</td><td>0.0 kpc</td></tr><tr><td>$b_b$</td><td>0.267 ± 0.009 kpc</td></tr><tr><td>$M_b$</td><td>1.02 × 10<sup>10</sup> ± 6.3 × 10<sup>8</sup>$M<sub>☉</sub></td></tr><tr><td>$M_n$</td><td>4.10<sup>6</sup> ± 0.42 × 10<sup>6</sup>$M<sub>☉</sub></td></tr></tbody></table>

<table><caption>**Table 3**

Ranges of Input Parameters for SBI Procedure</caption><thead><tr><th>Quantity</th><th>Parameter</th><th>Range</th><th>Unit</th></tr></thead><tbody><tr><td>Standard deviation of birth magnetic field</td><td>$\sigma_b$</td><td>[0.1,0.9]</td><td></td></tr><tr><td>Standard deviation of birth period</td><td>$\sigma_p$</td><td>[0.1,0.9]</td><td></td></tr><tr><td>Mean birth period</td><td>$\mu_p$</td><td>[1,100]</td><td>ms</td></tr><tr><td>Mean birth magnetic field</td><td>$\mu_b$</td><td>[1 × 10<sup>10</sup>,2 × 10<sup>11</sup>]</td><td>T</td></tr><tr><td>Field decay scaling</td><td>$\alpha_d$</td><td>[-7.4]</td><td></td></tr><tr><td>Field decay timescale</td><td>$\tau_d$</td><td>[2,100]</td><td>kyr</td></tr><tr><td>Birth rate</td><td>$\text{BR}$</td><td>[1/1000,1/50]</td><td>yr<sup>-1</sup></td></tr></tbody></table>

normalizing flow that models complex distributions through a sequence of autoregressive transformations. The model was trained by minimizing the Kullback-Leibler (KL) divergence (S. Kullback & R. A. Leibler 1951), ensuring that the learned distribution closely matches the true posterior. This approach is particularly well-suited for NPE, as it enables flexible and efficient density estimation.

In this work, 10,490 simulations were employed to train the density estimator of SBI for parameter inference. This number was chosen based on a specific criterion: following the recommendation from other SBI users, who suggest using at least 1000 × *n* simulations, where *n* is the number of parameters. Additionally, due to the complexity of simulating magnetars, the number of simulations for inference was extended to 10,000 in order to ensure confidence in the results.

## 7. Model Validation

We performed validation checks to test the robustness of our SBI pipeline. First, a population of magnetars was simulated following the technique described in Sections 2 and 3, with a known set of model parameters $\theta_{\text{in}}$. This data set replaced the “true” magnetar population, and the SBI pipeline was run under this assumption to check whether the model parameters obtained with SBI $\theta_{\text{out}}$ match the known set of model parameters $\theta_{\text{in}}$. We repeated this process many times. The results from these tests are given in Figure 2, which shows that the ground-truth values were largely recovered by the SBI process, except for the birth-period parameters ($\mu_p$ and $\sigma_p$). The birth rate (1/X) and the field decay parameters ($\alpha_d$ and $\tau_d$) were the parameters that were better recovered, with the majority of results matching parameters within the 95% credible interval. To further validate the results, a second
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Convolutional neural network architecture used for this work.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Validation plot for the SBI pipeline. In the upper panels, input values are plotted directly against their recovered counterparts with 68% credible interval error bars. The lower panels directly compare the CDFs of both the input and output parameters. If the pipeline has successfully recovered parameters, these lower plots should track the dotted line (unity); deviations indicate regions where the parameters are not recovered accurately.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** *P − *P diagrams following pipeline validation checks. The main panel includes contour maps for two collections of 5000 populations: the black contours correspond to the simulation constructed using a known set of parameters drawn from the prior (θ_in), and the red contours are the populations generated from parameters drawn from the posteriors following a full run of the SBI pipeline (θ_out). Along the top and right side are the CDFs for *P* and *P* for both populations.</figcaption></figure>

comparison was performed using 5000 *P* − *P̂* diagrams from populations where contours were created for bursting sources. These populations were generated using both the known set of model parameters θ_in and the model parameters inferred from the posteriors obtained with SBI. These are overlaid in

Figure 3, and the CDFs of both distributions are given as well. It can be seen here that the contour of the reference population (generated with θ_in) has been quite well reconstructed from the distributions of parameters found by SBI, validating the pipeline’s ability to find the parameters. It shows the pipeline’s insensitivity to birth period, suggesting that the current magnetar population is insensitive to the parameters μ_p and σ_p (as reported in other works, e.g., M. Gullón et al. 2015, P. Beniamini et al. 2019, and J. A. Jawor & T. M. Tauris 2022). This is most likely due to the fast decay of the magnetar spins and magnetic fields, as a magnetar will be detected with ages comparable to τ_d but the timescale of modification of *P* is much shorter than τ_d. Because of this, the birth spins are very hard to determine with a small sample size; with only 30 known magnetars, we were unable to sufficiently probe the initial spin period. We further detail our validation steps in the Appendix.

## 8. Results

We have applied the SBI pipeline to the observed magnetar population. The posterior distributions for the model parameters are shown in Figure 4. The best estimates are presented in Table 4 with their 95% credible intervals. We do not report the values found for μ_p, as the validation checks suggest that this parameter is not recovered well with the available data. However, while we do not show the tests here, using the best estimates from Table 4 with different μ_p and σ_p made clear that 1 s for μ_p and 0.9 for σ_p are appropriate upper limits, as
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Corner plot of results obtained with SBI. Vertical yellow lines indicate the best estimate.</figcaption></figure>

<table><caption>**Table 4**
Results of the SBI Pipeline</caption><thead><tr><th>Parameter [Unit]</th><th>Value</th></tr></thead><tbody><tr><td>σ<sub>b</sub></td><td>0.2<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td>log (μ<sub>b</sub>) (μ<sub>b</sub> [T])</td><td>10.2<sup>+0.1</sup><sub>-0.2</sub></td></tr><tr><td>α<sub>d</sub></td><td>1.9<sup>+0.9</sup><sub>-1.3</sub></td></tr><tr><td>τ<sub>d</sub> [kyr]</td><td>17.9<sup>+24.1</sup><sub>-14.5</sub></td></tr><tr><td>Birth spacing [yr]</td><td>549<sup>+292</sup><sub>-323</sub></td></tr></tbody></table>

**Note.** Uncertainties indicate the 95% credible interval.

models far from these values lead to *P* − *Ṗ* diagrams where simulations and observations no longer align. Note that the measured birth spacing reported in Table 4 corresponds to birth rates BR = 1.8<sup>+2.6</sup><sub>-0.6</sub> kyr<sup>-1</sup>; our result provides a significantly tighter constraint compared to the bounds of 2.3–20 kyr<sup>-1</sup> reported in the previous study of P. Beniamini et al. (2019).

The contour plot presented in Figure 5 shows that the posterior distributions on the parameters obtained with the SBI

pipeline produce magnetars that broadly align with observations. Note that, to compare with known magnetars, the only sources included in these contours are those that have undergone an outburst in the last 50 yr; we have assumed that all of these would be detected, so we refer to this subpopulation as the “simulated” group. The known magnetars (from the McGill catalog; S. A. Olausen & V. M. Kaspi 2014) are shown as red stars in the *P* − *Ṗ* diagram. Furthermore, even magnetars with low *Ṗ* are generated and have undergone outbursts in the simulation with a low occurrence rate, matching the observations. When conducting a KS test on the *P* and *Ṗ* distributions of each simulation and the observations, 91% of the measured *p*-values are greater than 0.05 for *Ṗ* and 87% are greater than 0.05 for *P*, suggesting a significant similarity between these simulations and the observations. This is highlighted in the upper and right panels of Figure 5, where we show the CDFs of the observed population and those of 5000 simulations for *P* and *Ṗ*. The CDFs start to diverge from each other slightly at high *P* and low *Ṗ*, where we are largely hampered by low-number statistics.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Contour P − ˙ P diagram. To produce the contours, 5000 simulations were obtained using the posterior distributions of the parameters obtained at the end of the SBI. Only sources that have undergone an outburst in the past 50 yr are included, as this is the subpopulation we assume to have been detected. Plotted in red are the observed Galactic magnetars. Along the side and top of the plot are the CDFs of each distribution for comparison.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** LogN–logS plots of the observations compared with the mean logN–logS of the mean of 5000 simulations.</figcaption></figure>

As another layer of comparison, we generated logN–logS distributions for the observed and simulated populations in Figure 6. The solid blue line indicates the observed logN–logS curve (unabsorbed fluxes were used), and the solid green line indicates the mean logN–logS curve of the 5000 post-SBI simulations. The shaded area represents the 1σ uncertainty. We can accurately predict the shape of the logN–logS within 1σ for most of the flux ranges. There is a small flux range (10−14–10−13.5 W m−2) within which the observed values deviate by slightly more than 1σ compared to our simulations.

The temperature of the blackbody model for the X-ray emission serves as another point of comparison, as illustrated in Figure 7. The emitting surface fraction of observed magnetars, defined as $\xi = (R_{eff}/R)^2$, with R being the neutron

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Contour of the temperature versus characteristic age of the simulated magnetars in comparison with the observed magnetars, while choosing $\xi = 0.01$ (in blue) and $\xi = 0.08$ (in orange) as the emitting surface fractions (normalized by the total number of magnetars over 1000 simulations in each bin).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Probability distribution of the inclination angles for simulated magnetars: initial (red), final (green), and those younger than 17.9 kyr (blue), normalized by the total number of magnetars across 5000 simulations in each bin.</figcaption></figure>

star radius and $R_{\text{eff}} = \sqrt{\frac{d^2S}{\sigma T^4}}$ (following the blackbody emission definition in Equation (21)), yields an average value of $\xi = 0.06$ with a standard deviation of 0.08. The observed emitting surface fraction rarely exceeds 8%, which motivates tests with $\xi$ ranging from 1% to 8%. However, for clarity, the plots presented here focus only on these two extreme values. The observations indeed fall within the simulated contours in the temperature versus characteristic age plot ($\tau_c = P/2\dot{P}$) when $\xi$ is within this range. These contours suggest that the emitting surface fraction might be relatively low for most magnetars, as lower values of $\xi$ tend to place observed magnetars in more plausible regions of the simulated parameter space. However, the precision of this estimate is primarily constrained by the limited number of known magnetars. Additionally, the lack of detailed knowledge regarding the X-ray emission region of magnetars further hinders the ability to draw robust conclusions on this matter.
<figure><img src="image_12.png" /><figcaption>Figure 9. Ratio of real and characteristic age as a function of characteristic age, with the parameter $\alpha_d$ used to simulate each magnetar in 120 simulations.</figcaption></figure>

The distribution of inclination angles in the simulated population is shown in Figure 8. Around 80% of the population exhibits either alignment at 0° or anti-alignment at 180°. The analytical spin-alignment timescale, following L. Dirson et al. (2022) and A. Philippov et al. (2014), is given by

$$ \tau_{\text{align}}^{\text{MHD}} = \frac{I\mu_0 c^3 P_0^2 \sin^2 \alpha_0}{16\pi^3 R^6 B_0^2 \cos^4 \alpha_0}, \quad (26) $$

where subscript 0 denotes initial values (except for the vacuum permeability). For a magnetar with $B_0 = 1.7 \times 10^{10}$ T (the best estimate of $\mu_b$ with SBI), $P_0 = 50$ ms, and $\alpha_0 = 89^\circ$, we obtain $\tau_{\text{align}}^{\text{MHD}} \approx 7 \times 10^6$ yr. Conversely, an initial $\alpha_0 = 91^\circ$ results in anti-alignment over the same period. A slight decrease to $\alpha_0 = 85^\circ$ leads to $\tau_{\text{align}}^{\text{MHD}} \approx 11$ kyr, and smaller $\alpha_0$ values further shorten the alignment timescale. Given the best estimate for the magnetic decay timescale $\tau_d = 17.9$ kyr, most magnetars align before significant magnetic field decay, as shown by the blue distribution representing the magnetars that are younger than $\tau_d = 17.9$ kyr. Figure 8 also confirms an initially uniform $\alpha_0$ distribution among magnetars, implying rapid alignment regardless of initial inclination. This is consistent with the computed alignment timescales, which become significantly shorter as $\alpha_0$ deviates from 90°. Consequently, only magnetars with $\alpha_0$ near 90° retain a final inclination different from alignment or anti-alignment. As a result, only 20% of magnetars exhibit nonaligned angles, with fewer than 1.5% having $\alpha_0$ between 60° and 120°. This behavior contrasts with canonical pulsars, where the distribution is more uniform and peaks at 90°, making initial inclination a critical factor in their detection (see Figure 13 of L. Dirson et al. 2022).

Figure 9 presents the ratio of real to characteristic age for magnetars across 120 simulations, highlighting the impact of $\alpha_d$ on their evolution. Negative $\alpha_d$ leads to a rapid divergence between real and characteristic age once the magnetar surpasses the characteristic magnetic decay timescale. In contrast, for positive $\alpha_d$, the larger the value, the longer it takes for this deviation to occur. The posteriors favor $\alpha_d \approx 2$, implying that only magnetars older than $10^4$ yr—around 20%

<figure><img src="image_13.png" /><figcaption>Figure 10. Positions of the simulated magnetars in the Galactic plane (normalized by the total number of magnetars over 5000 simulations in each bin) in comparison with the observed magnetars.</figcaption></figure>

of the detected population—experience a significant discrepancy between real and characteristic age. A few young magnetars with $\alpha_d \approx 2$ also show such differences, due to rapid alignment of their rotation and magnetic axes, which quickly results in high $P$ and low $\dot{P}$.

Recovering the spatial distribution of magnetars presents a challenge in population synthesis. Figure 10 compares the fraction of magnetars occurring in a given region of the Galactic plane across 5000 simulations with the observed magnetar positions. While the simulated and observed positions do not align, this is expected, as the outburst detection model used here is based solely on energy (which is dependent on magnetic field strength and age) without accounting for distance. The simulations reproduce the shape of the Galactic spiral arms, reflecting the youth of the magnetars that have not migrated significantly from their birth sites. Furthermore, even if distances were incorporated, discrepancies between simulations and observations would persist due to the considerable uncertainties in distance estimates. Including distance in the detection model could nonetheless provide valuable insights into the true spatial distribution of observed magnetars.

# 9. Discussion

## 9.1. Comparison with Other Works

In this study, we estimate the birth rate of magnetars to be $BR = 1.8^{+2.6}_{-0.6}$ kyr$^{-1}$, given a CCSN rate of $11.8-21.3$ kyr$^{-1}$ (K. Rozwadowska et al. 2021), within a 68.3% confidence interval—and given the fact that $30.4 \pm 5\%$ of CCSNe are classified as SNe IIb/Ib/Ic (I. Shivvers et al. 2017), the estimated rate of SNe IIb/Ib/Ic is approximately $5.0 \pm 1.6$ kyr$^{-1}$ (1 $\sigma$ uncertainty). This estimate is consistent with the birth rate of magnetars within the $1.6 \sigma$ uncertainty found in this work. Moreover, if not all SNe IIb/Ib/Ic lead to magnetar formation, then the $5.0 \pm 1.6$ kyr$^{-1}$ rate should be interpreted as an upper limit on the magnetar birth rate, which may indeed be the case. Given the CCSN rate of $11.8-21.3$ kyr$^{-1}$ (as mentioned earlier), with approximately 80% of these CCSNe leading to the formation of isolated neutron stars (C.-A. Faucher-Giguère & V. M. Kaspi 2006),
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Median period and source counts for varied decay parameters. Holding all other population parameters fixed, $\alpha_d$ and $\tau_d$ were varied and many populations were produced. The median simulated period was measured from each combination when enough sources had been produced for it to stabilize. In the right-most panel, the average number of sources per population is shown instead. The left three panels are split by detection criteria: the leftmost shows all sources, whether they are detectable or not; the second includes sources detected with S/N > 6 in a 10 ks XMM observation; the third includes those detected with S/N > 6 in a 10 ks AXIS observation; and the fourth shows sources that have undergone an outburst in the past 50 yr.</figcaption></figure>

and with the estimate from M. Sautron et al. (2024) of a 24 kyr⁻¹ birth rate for isolated canonical pulsars, it places the overall birth rate of isolated neutron stars between 9.4 and 24 kyr⁻¹. Meanwhile, the current star formation rate in the Milky Way is approximately $1.65 \pm 0.19 M_\odot \text{yr}^{-1}$ (T. C. Licquia & J. A. Newman 2015). Assuming a Kroupa initial mass function (P. Kroupa 2001), this corresponds to the birth rate of stars with $M > 8M_\odot$, most of which evolve into neutron stars, yielding an estimated neutron star birth rate of $17 \pm 2 \text{ kyr}^{-1}$ within a 95% confidence interval. By comparing the last approximation, which is consistent with the previous estimates, with our magnetar birth rate, we find that $10.7_{-4.4}^{+18.7}\%$ of neutron stars become magnetars. This is close to the 10% estimate from S. B. Popov et al. (2010) and provides a well-defined uncertainty on this value. It also provides more constraints compared to the estimate of P. Beniamini et al. (2019), who reported that more than 12% of neutron stars become magnetars. This result is on the lower end of the previous estimates, slightly relaxing the birth rate problem of neutron stars (A. A. Gençali & Ü. Ertan 2024).

To obtain a similar simulated population, the magnetic field decay timescale must be quite short; in this study, we find $\tau_d = 17.9_{-14.5}^{+24.1} \text{ kyr}^{-1}$. This estimation is in agreement with previous studies of the magnetar population that found $\tau_d \sim 10^4 \text{ yr}$ (C. Kouveliotou et al. 1998; M. Colpi et al. 2000; S. Dall’Osso et al. 2012; D. Viganò et al. 2013; P. Beniamini et al. 2019).

The mean magnetic field at birth, $\mu_b$, is estimated between $1.15 \times 10^{10} \text{ T}$ and $2.40 \times 10^{10} \text{ T}$ in our work. This is much greater than the estimate of M. Gullón et al. (2015) and lower than the estimate of P. Beniamini et al. (2019). Interestingly, this range of values is consistent with the magnetic field strengths required for millisecond magnetars to power superluminous supernovae (SLSNe; B. D. Metzger et al. 2015, 2018) and GRBs (P. Beniamini et al. 2017; B. D. Metzger et al. 2018). The other key factor is the initial spin period; however, the only constraint on its distribution suggests upper limits of 1 s for $\mu_p$ and 0.9 for $\sigma_p$.

P. Beniamini et al. (2019) suggested a relationship between the magnetic field decay parameters and the period distribution of the magnetar population that can be seen in simulations from this study. To demonstrate this, we generated populations of sources while varying $\tau_d$ and $\alpha_d$ and fixing the remaining population parameters. For each combination of decay parameters, many populations were produced and the median period and total number of sources per population were measured. The number of populations used varied, as

increasing alpha dramatically increased the number of magnetars in the population. For larger $\alpha_d$, fewer populations could be used without sacrificing the stability in median/ source count measurements. Additionally, the populations were filtered by the detection criteria used above to show the difference in detectable and “full” populations. These results are shown in Figure 11, and they illustrate the clear dependence on both maximum period and number of Galactic magnetars.

## 9.2. Model Improvements and Future Plans

While the magnetar evolution model proposed in this work features most of the known aspects of the evolution of these objects, it can be improved. For instance, similarly to canonical pulsars, magnetars have been observed to glitch (A. G. Lyne et al. 2009; D. Eichler & R. Shaisultanov 2010; G. Younes et al. 2023) in many studies. Furthermore, J. R. Fuentes et al. (2017) found a relationship for all neutron stars between the glitch activity (a quantity relating the increase in rotation frequency over the time of observations for all the glitches of the concerned neutron star) and the spin-down rate. In the future, we could incorporate glitches in our model with this relationship and test if we are still able to recover a population similar to the observations. Adding this phenomenon would form a more realistic model and may tighten constraints on the parameters that characterize the population of magnetars.

Another improvement for the model would be to refine the bursting model used in this paper. Outburst events are very energetic, but we are certainly missing some of these (especially less-energetic bursts). For example, magnetars XTE J1810−197 and 3XMM J185246.6+003317 were both discovered serendipitously in observations of other targets (A. I. Ibrahim et al. 2004; P. Zhou et al. 2014), with detections showing that these older sources emit weaker bursts with longer tails compared to their younger counterparts. This may suggest that our technique for assigning bursts should weight the burst detectability with the magnetars’ ages. Running simulations with a parameter that informs us about the fraction of detected outbursts (this topic is addressed in P. Beniamini et al. 2025) could enable to estimate how many we miss, further improving the accuracy of our results. To achieve a more realistic simulation, it would not be sufficient to model only the rate of outbursts. Instead, each outburst should be characterized (in addition to the rate of outburst) by an energy release, $E_{\text{out}}$, and a duration, $t_{\text{out}}$. Both $E_{\text{out}}$ and $t_{\text{out}}$ should be drawn from appropriate distributions that have to be
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Detectable magnetar point sources in 1000 simulated populations. The sky region is limited to $|b| \le 5$, and detection is limited to the 0.4–7 keV band. The black band indicates the mean number of sources in all realizations for each limiting flux. In red, known Galactic magnetars are included (only those in the plane with published fluxes are shown). The green band corresponds to sources detected in the periodicity search in XMM data, and blue indicates sources that are detected above 6 S/N in those data. Also included for comparison are comparable results using simulated AXIS response files; these are plotted in orange and magenta. All shaded regions correspond to $1\sigma$ uncertainties.</figcaption></figure>

determined, guided by studies such as R. Perna & J. A. Pons (2011), C. Dehman et al. (2020), S. K. Lander (2023), and P. Beniamini et al. (2025), which have investigated the bursting phenomenology of magnetars. The corresponding flux, $F_{\text{out}}$, could then be computed based on the source distance and compared to the detection threshold of a specific instrument (e.g., Swift-BAT). To refine this comparison, one could also account for the instrument's attitude and on-times, ensuring that the source was within the field of view and the instrument was active during potential detection periods. Detection would be accepted only if the flux exceeds this threshold and the observational conditions were met. Importantly, the distributions of $E_{\text{out}}$ and $t_{\text{out}}$ should also be correlated with the age of the magnetars, such that older magnetars tend to produce weaker and longer events as mentioned previously. As a consequence, this more refined modeling will likely lead to higher estimates of the birth rate and $\tau_d$, since a fraction of distant or intrinsically faint magnetars should now go undetected. In order to still detect a few magnetars with high P and low $\dot{P}$ in the simulations, $\tau_d$ might need to increase to allow magnetic fields to decay more slowly, thus maintaining a higher probability of bursts—even if they are weaker. Ultimately, this would improve the match between the spatial distribution of simulated and observed magnetars, reducing the excess of simulated sources found across the whole Galaxy, an effect not seen in the observations, since the observed magnetars seem to be near to us. Consequently, $F_{\text{out}}$ should also incorporate the column density along the line of sight, as it plays a key role in shaping detection biases. However, as noted in Section 8, the observed positions still suffer from significant uncertainties. Therefore, even in this refined model describing the bursts, the similarity between observations and the simulation is imperfect and might inform us better with regard to the real positions of the population of magnetars. This analysis needs to be studied in more detail, and it will be explored in a future work.

Finally, another perspective to improve the simulations would be to consider the possible precession of magnetars (A. Melatos 2000; L. Arzamasskiy et al. 2015; S. K. Lander & D. I. Jones 2018; A. P. Igoshev & S. B. Popov 2020).

### 9.3. Detection Prospects

With a suitable model for the population of magnetars, we can assess the detectability of new magnetar sources. To that end, using realizations of magnetar populations generated with the optimal parameters, we examined all sources above a nominal flux threshold of $10^{-20} \text{ W m}^{-2}$ (i.e., $10^{-17} \text{ erg/(s cm}^2$)) in the plane ($|b| > 5^\circ$) using the same detection procedure described in Section 5. This formed a mock survey of the sky in the plane, which informs future searching techniques and detection probabilities. For each source, absorption was modeled and a simulated XMM lightcurve was generated for a 10 ks observation. We applied an H-test and detectable sources were recorded. Here, detectability is broken into three subgroups: the first is simply whether the source has undergone an outburst as predicted using the model described in Section 3; the second is defined as “point-source” detection, with S/N > 6 in a simulated lightcurve; and finally, whether the source is detectable in a periodicity search. Bursting sources are expected to most closely match the true, Galactic magnetar population that has been detected. Sources with S/N > 6 would not necessarily be identified as magnetars. But periodic detections should include several of the Galactic discoveries as well as a possibly undetected-but-detectable population. This subsample indicates how effective a full-plane search would be, and it can be tailored to various instrumentation. For consistency, we simulate XMM observations. Figure 12 shows the number of sources above a given limiting flux $S_{\text{lim}}$ as a function of $S_{\text{lim}}$ for 1000 realizations. This number was chosen based on the stability of the mean value in each bin; above 500 simulations, this mean did not vary substantially, so 1000 ensures accuracy. This search resulted in seven periodic detections per simulation on average, and the average number of magnetars above XMM’s flux threshold ($\approx 10^{-14} \text{ erg/(s cm}^2$) $= 10^{-17} \text{ W m}^{-2}$) is approximately 150. More sensitive instruments will probe more deeply, and experiments like the Advanced X-ray Imaging Satellite (AXIS; R. Arcodia et al. 2024) are expected to detect objects fainter by 2 orders of magnitude (S. Safi-Harb et al. 2023). To test this, we repeated the procedure described above with data simulated using AXIS predictions; the results are shown in Figure 12. It is apparent that AXIS will detect many more sources, with the potential of doubling periodic detections and more than doubling the number of detectable magnetar point sources.

It is apparent from Figure 12 that a significant fraction of sources are detectable as point sources (i.e., S/N > 6) but are not detected in the periodicity search, including those on the brightest end of the plot. Examining those sources more closely revealed that the evolution model was producing many sources outside of the periodic search range of 0.5–20 s. Sources in the $10^{-14}-10^{-13} \text{ W m}^{-2}$ range correspond to young, bright magnetars with periods below 0.5 s; at the low-flux end, the gap between the two detection populations corresponds to a large number of older magnetars that have evolved to periods longer than 20 s. Notably, the former period range includes the periods of both the “magnetar-like” high-B pulsars J1119–6127 (0.4 s period) and J1846–0258 (0.3 s
<figure><img src="image_16.png" /><figcaption>**Figure 13.** CDFs of simulated sources, split by detection criteria. The dashed lines indicate the period range that was used in the search pipeline, where the left side of the plot indicates periods below 0.5 s, those in the center are between 5 and 20 s, and those on the right are above 20 s. The green, black, red, and cyan curves are cumulative distributions of the sources from 1000 simulated populations. Error bars indicate $1\sigma$ uncertainties on the distributions. Shown in each period range are the fraction of each detection subpopulation that are within that range and the mean characteristic age of those sources.</figcaption></figure>

period). Neither of these sources were initially discovered during outburst (F. Camilo et al. 2000; E. V. Gotthelf et al. 2000), but both have undergone bursts since (R. F. Archibald et al. 2016; H. Blumer et al. 2021). In fact, this model suggests that as many as 80% of magnetars that could be detected as point sources (and possibly as periodic sources) would effectively be omitted from the search simply because of the search bounds. Figure 13 shows the CDFs of the spin periods of those sources that are “detected” according to the three detection criteria. In order to assess whether the sources outside of our search range would be detectable in a periodicity search, we have included counts for periodic detections outside of that range. This result is corroborated by the recent discoveries of slowly rotating radio emitters that may be old magnetars (P. Beniamini et al. 2020; M. Caleb et al. 2022; N. Hurley-Walker et al. 2022; P. Beniamini et al. 2023; Z. Wang et al. 2024), and it may link magnetars to high-B pulsars (F. Camilo et al. 2000; B. Zhang & A. K. Harding 2000). Granted, the search region does cover the majority of sources (≥68%, assuming all bursts are detected) that have undergone bursting in the last 50 yr, reflecting the detection bias toward younger sources with shorter periods.

## 10. Conclusions

Simulation-based inference was used in this study to constrain the parameters characterizing the Galactic population of magnetars for the first time. By changing only the parameters for the birth distributions and the magnetic field decay compared to normal pulsars, we used an evolution model that is otherwise identical to normal pulsars and showed it was possible to produce a population of magnetars. This study considered both persistent X-ray detection and outburst

detection, taking into account biases from sky coverage and observation information. In the future, with telescopes sensitive enough to detect magnetars through persistent emission alone, this pipeline can be used as-is with only minor changes to the reference population statistics.

Furthermore, the use of CNNs was essential for extracting features from the density maps, preserving the correlations between the different quantities, and adequately informing the density estimator of SBI. First, the pipeline was validated by a population that was obtained via a known set of parameters. After validation, the following population parameters were constrained: the mean magnetic field at birth $\log(\mu_b) = 10.2_{-0.2}^{+0.1}$ and standard deviation $\sigma_B = 0.2_{-0.1}^{+0.2}$, the parameter controlling the speed of the magnetic field decay $\alpha_d = 1.9_{-1.3}^{+0.9}$, the typical decay timescale $\tau_d = 17.9_{-14.5}^{+24.1}$ kyr, and the birth rate $\text{BR} = 1.8_{-0.6}^{+2.6} \text{ kyr}^{-1}$. Validation checks showed that the current population is globally insensitive to the initial spin period distribution parameters, $\mu_p$ and $\sigma_p$. However, upper limits of 1 s for $\mu_p$ and 0.9 for $\sigma_p$ seem to be appropriate. The recovered blackbody temperatures of the simulated magnetars suggest that the effective emission radius of magnetars may be quite low for most magnetars, and will require more study to match the temperature versus characteristic age plot with the observations. An analysis of the inclination angle showed that around 80% of the detected population have aligned (or anti-aligned) rotation and magnetic axes, which implies that seeing pulsation from this population could be quite difficult compared to pulsars.

Nonetheless, a limitation of the simulation is that the positions of the detected magnetars do not correspond to those of the real magnetars. When accepting a magnetar as detected via outburst, no discrimination through the energy nor the
distance is made. Therefore, in order to improve this work in the future, these factors should be taken into account as discussed in the previous section.

Predictions from these results suggest that a full-scale Galactic survey for magnetars using XMM would be able to uncover several dozen magnetars with a modest exposure duration. From this relatively shallow mock survey, the population could feasibly double; using upcoming improved instruments like AXIS could potentially quadruple it. Given that many of the results of this study (and other magnetar population studies) are challenged by the small sample size, large-scale searches like those simulated here would be crucially important to fully understanding this population of neutron stars.

## Acknowledgments

This work has been supported by the grant ANR-20-CE31-0010 and is based upon work supported by NASA under award number 80GSFC24M0006. We thank Ósmar Rodríguez for his insightful comments regarding the birth rate of magnetars and stripped-envelope supernova. We also thank the anonymous referee for the helpful comments and suggestions, which helped improve the paper.

**Facilities:** XMM.

**Software:** Astropy (A. M. Price-Whelan et al. 2018), SciPy (P. Virtanen et al. 2020), NumPy (T. E. Oliphant 2006), PyGEDM (D. C. Price et al. 2021), PyXSPEC (C. Gordon & K. Arnaud 2021), SBI (A. Tejero-Cantero et al. 2020), TensorFlow (M. Abadi et al. 2016).

## Appendix

### Extended Model Validation

To complete the model validation, the Kolmogorov–Smirnov (KS) test (N. Smirnov 1948) is applied to the *P* and

<figure><img src="image_17.png" /><figcaption>**Figure 14.** *p*-values from two-sample KS tests of *P* and *P̂* between populations produced using $\theta_{in}$ and $\theta_{out}$. In total, 20,000 population comparisons were made. The vertical dashed line indicates a threshold *p*-value of 0.01; for the two distributions, 2.68% and 6.08% lie below this threshold.</figcaption></figure>

$P̂$ distributions between the reference population and each validation simulation. The process was repeated 20,000 times, and the distribution of the *p*-values is plotted in Figure 14. Those results suggest that, in the vast majority of cases, the null hypothesis cannot be rejected.

The posterior distributions obtained in Figure 15 are an example of the results from one of the validation runs mentioned in Section 7. In this figure, it is apparent that, in each distribution, the ground truth (in red, representing $\theta_{in}$) is close to the median of the posterior found (shown in yellow, $\theta_{out}$), except for the distribution of $\mu_p$. Even though in this run $\sigma_p$ is well-recovered, the distribution for this parameter is quite spread, indicating that it also seems insensitive to this parameter, as evidenced by Figure 2.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Corner plot of the posterior distributions of the parameters corresponding to the reference simulation obtained after one validation run. Posterior medians are in yellow, and the ground truths are in red.</figcaption></figure>
ORCID iDs

M. Sautron https://orcid.org/0009-0006-7913-1186

A. E. McEwen https://orcid.org/0000-0001-5481-7559

G. Younes https://orcid.org/0000-0002-7991-028X

J. Pétri https://orcid.org/0000-0003-3790-8066

P. Beniamini https://orcid.org/0000-0001-7833-1043

D. Huppenkothen https://orcid.org/0000-0002-1169-7486

References

Abadi, M., Agarwal, A., Barham, P., et al. 2016, arXiv:1603.04467

Alsing, J., Charnock, T., Feeney, S., & Wandelt, B. 2019, MNRAS, 488, 4440

Anderson, G. E., Gaensler, B. M., Slane, P. O., et al. 2012, ApJ, 751, 53

Archibald, R. F., Kaspi, V. M., Tendulkar, S. P., & Scholz, P. 2016, ApJL, 829, L21

Arcodia, R., et al. 2024, Univ, 10, 316

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

Arzamasskiy, L., Philippov, A., & Tchekhovskoy, A. 2015, MNRAS, 453, 3540

Bajkova, A. T., & Bobylev, V. V. 2021, A&AT, 32, 177

Barrère, P., Guilet, J., Raynaud, R., & Reboul-Salze, A. 2025, A&A, 695, A183

Belczynski, K., & Taam, R. E. 2008, ApJ, 685, 400

Beniamini, P., Giannios, D., & Metzger, B. D. 2017, MNRAS, 472, 3058

Beniamini, P., Hotokezaka, K., van der Horst, A., & Kouveliotou, C. 2019, MNRAS, 487, 1426

Beniamini, P., & Kumar, P. 2023, MNRAS, 519, 5345

Beniamini, P., & Kumar, P. 2025, ApJ, 982, 45

Beniamini, P., & Lu, W. 2021, ApJ, 920, 109

Beniamini, P., & Piran, T. 2024, ApJ, 966, 17

Beniamini, P., Wadiasingh, Z., Hare, J., et al. 2023, MNRAS, 520, 1872

Beniamini, P., Wadiasingh, Z., & Metzger, B. D. 2020, MNRAS, 496, 3390

Beniamini, P., Wadiasingh, Z., Trigg, A., et al. 2025, ApJ, 980, 211

Berteaud, J., Eckner, C., Calore, F., Clavel, M., & Haggard, D. 2024, ApJ, 974, 144

Beskin, V. S., Gurevich, A. V., & Istomin, Y. N. 1993, Physics of the Pulsar Magnetosphere (Cambridge: Cambridge Univ. Press)

Beskin, V. S., Istomin, Y. N., & Philippov, A. A. 2013, PhyU, 56, 164

Beskin, V. S., & Nokhrina, E. E. 2007, Ap&SS, 308, 569

Bhardwaj, M., Gaensler, B. M., Kaspi, V. M., et al. 2021, ApJL, 910, L18

Bhardwaj, U., Alvey, J., Miller, B. K., Nissanke, S., & Weniger, C. 2023, PhRvD, 108, 042004

Bhattacharya, D., Wijers, R. A. M. J., Hartman, J. W., & Verbunt, F. 1992, A&A, 254, 198

Blumer, H., Safi-Harb, S., McLaughlin, M. A., & Fiore, W. 2021, ApJL, 911, L6

Bochenek, C. D., Ravi, V., Belov, K. V., et al. 2020, Natur, 587, 59

Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJL, 887, L26

Borkowski, K. J., & Reynolds, S. P. 2017, ApJ, 846, 13

Bovy, J. 2015, ApJS, 216, 29

Caleb, M., Heywood, I., Rajwade, K., et al. 2022, NatAs, 6, 828

Camilo, F., Kaspi, V. M., Lyne, A. G., et al. 2000, ApJ, 541, 367

Collazzi, A. C., Kouveliotou, C., van der Horst, A. J., et al. 2015, ApJS, 218, 11

Colpi, M., Geppert, U., & Page, D. 2000, ApJL, 529, L29

Cordes, J. M., & Chatterjee, S. 2019, ARA&A, 57, 417

Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. 2018, MNRAS, 474, 961

Cowsik, R. 1998, A&A, 340, L65

Cranmer, K., Brehmer, J., & Louppe, G. 2020, PNAS, 117, 30055

Cumming, A., Arras, P., & Zweibel, E. 2004, ApJ, 609, 999

Dai, Z. G., & Lu, T. 1998, A&A, 333, L87

Dall'Osso, S., Granot, J., & Piran, T. 2012, MNRAS, 422, 2878

de Jager, O. C., & Büsching, I. 2010, A&A, 517, L9

Dehman, C., Viganò, D., Rea, N., et al. 2020, ApJL, 902, L32

Ding, H., Lower, M. E., Deller, A. T., et al. 2024, ApJL, 971, L13

Dirson, L., Pétri, J., & Mitra, D. 2022, A&A, 667, A82

Doroshenko, V. 2024, arXiv:2403.03127

Duncan, R. C., & Thompson, C. 1992, ApJL, 392, L9

Eichler, D., & Shaisultanov, R. 2010, ApJL, 715, L142

Faucher-Giguère, C.-A., & Kaspi, V. M. 2006, ApJ, 643, 332

Figer, D. F., Najarro, F., Geballe, T. R., Blum, R. D., & Kudritzki, R. P. 2005, ApJL, 622, L49

Foglizzo, T., Kazeroni, R., Guilet, J., et al. 2015, PASA, 32, e009

Fuentes, J. R., Espinoza, C. M., Reisenegger, A., et al. 2017, A&A, 608, A131

Gelfand, J. D., & Gaensler, B. M. 2007, ApJ, 667, 1111

Gençali, A. A., & Ertan, Ü. 2024, MNRAS, 534, 1481

Gordon, A. C., Fong, W., Kilpatrick, C. D., et al. 2023, ApJ, 954, 80

Gordon, C., & Arnaud, K., 2021 PyXspec: Python Interface to XSPEC Spectral-fitting Program, Astrophysics Source Code Library, ascl:2101.014

Gotthelf, E. V., Vasisht, G., Boylan-Kolchin, M., & Torii, K. 2000, ApJL, 542, L37

Göğüş, E., Woods, P. M., Kouveliotou, C., et al. 1999, ApJL, 526, L93

Göğüş, E., Woods, P. M., Kouveliotou, C., et al. 2000, ApJL, 532, L121

Graber, V., Ronchi, M., Pardo-Araujo, C., & Rea, N. 2024, ApJ, 968, 16

Gullón, M., Pons, J. A., Miralles, J. A., et al. 2015, MNRAS, 454, 615

Han, J. L., Wang, C., Wang, P. F., et al. 2021, RAA, 21, 107

Hermans, J., Begy, V., & Louppe, G. 2019, arXiv:1903.04057

Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974

Hu, C.-P., Ng, C. Y., & Ho, W. C. G. 2019, MNRAS, 485, 4274

Hurley-Walker, N., Zhang, X., Bahramian, A., et al. 2022, Natur, 601, 526

Ibrahim, A. I., Markwardt, C. B., Swank, J. H., et al. 2004, ApJL, 609, L21

Igoshev, A. P., Frantsuzova, A., Gourgouliatos, K. N., et al. 2022, MNRAS, 514, 4606

Igoshev, A. P., & Popov, S. B. 2020, MNRAS, 499, 2826

Igoshev, A. P., & Popov, S. B. 2024, MNRAS, 535, L54

Jawor, J. A., & Tauris, T. M. 2022, MNRAS, 509, 634

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Kapil, V., Mandel, I., Berti, E., & Müller, B. 2023, MNRAS, 519, 5893

Kasen, D., & Bildsten, L. 2010, ApJ, 717, 245

Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261

Kirsten, F., Marcote, B., Nimmo, K., et al. 2022, Natur, 602, 585

Kothes, R., Sun, X., Gaensler, B., & Reich, W. 2018, ApJ, 852, 54

Kouveliotou, C., Fishman, G. J., Meegan, C. A., et al. 1994, Natur, 368, 125

Kouveliotou, C., Dieters, S., Strohmayer, T., et al. 1998, Natur, 393, 235

Kroupa, P. 2001, MNRAS, 322, 231

Kullback, S., & Leibler, R. A. 1951, The Annals of Mathematical Statistics, 22, 79

Kumar, P., Lu, W., & Bhattacharya, M. 2017, MNRAS, 468, 2726

Lander, S. K. 2023, ApJL, 947, L16

Lander, S. K., & Jones, D. I. 2018, MNRAS, 481, 4169

Licquia, T. C., & Newman, J. A. 2015, ApJ, 806, 96

Lu, W., Beniamini, P., & Kumar, P. 2022, MNRAS, 510, 1867

Lu, W., Kumar, P., & Zhang, B. 2020, MNRAS, 498, 1397

Lyne, A. G., Jordan, C. A., Graham-Smith, F., et al. 2015, MNRAS, 446, 857

Lyne, A. G., McLaughlin, M. A., Keane, E. F., et al. 2009, MNRAS, 400, 1439

Manchester, R. N., Lyne, A. G., Camilo, F., et al. 2001, MNRAS, 328, 17

Margalit, B., Beniamini, P., Sridhar, N., & Metzger, B. D. 2020, ApJL, 899, L27

Margalit, B., Berger, E., & Metzger, B. D. 2019, ApJ, 886, 110

McEwen, A. E., Lynch, R. S., Kaplan, D. L., et al. 2024, ApJ, 969, 118

McKay, M. D., Conover, W. J., & Becker, R. L. 1979, Technometrics, 21, 239

Melatos, A. 2000, MNRAS, 313, 217

Metzger, B. D., Beniamini, P., & Giannios, D. 2018, ApJ, 857, 95

Metzger, B. D., Giannios, D., Thompson, T. A., Bucciantini, N., & Quataert, E. 2011, MNRAS, 413, 2031

Metzger, B. D., Margalit, B., Kasen, D., & Quataert, E. 2015, MNRAS, 454, 3311

Mishra-Sharma, S., & Cranmer, K. 2022, PhRvD, 105, 063017

Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Noutsos, A., Schnitzeler, D. H. F. M., Keane, E. F., Kramer, M., & Johnston, S. 2013, MNRAS, 430, 2281

Olausen, S. A., & Kaspi, V. M. 2014, ApJS, 212, 6

Oliphant, T. E. 2006, A Guide to NumPy, 1 (USA: Trelgol Publishing)

Omelyan, I. P., Mryglod, I. M., & Folk, R. 2002, PhRvE, 66, 026701

Palmer, D. M., Barthelmy, S., Gehrels, N., et al. 2005, Natur, 434, 1107

Panther, F. H., Anderson, G. E., Bhandari, S., et al. 2023, MNRAS, 519, 2235

Papamakarios, G., & Murray, I. 2016, arXiv:1605.06376

Papamakarios, G., Pavlakou, T., & Murray, I. 2017, arXiv:1705.07057

Papamakarios, G., Sterratt, D. C., & Murray, I. 2018, arXiv:1805.07226

Perna, R., & Pons, J. A. 2011, ApJL, 727, L51

Pétri, J. 2012, MNRAS, 424, 605

Philippov, A., Tchekhovskoy, A., & Li, J. G. 2014, MNRAS, 441, 1879
Popov, S. B., Pons, J. A., Miralles, J. A., Boldin, P. A., & Posselt, B. 2010, MNRAS, 401, 2675

Popov, S. B., & Postnov, K. A. 2010, in Evolution of Cosmic Objects through their Physical Activity, ed. H. A. Harutyunian, A. M. Mickaelian, & Y. Terzian (Yerevan: Gitutyun), 129

Price, D. C., Flynn, C., & Deller, A. 2021, PASA, 38, e038

Price-Whelan, A. M., Hogg, D. W., Rix, H., et al. 2018, AJ, 156, 123

Rankin, J. M. 2007, ApJ, 664, 443

Rea, N., Gullón, M., Pons, J. A., et al. 2015, ApJ, 813, 92

Rea, N., Zane, S., Turolla, R., Lyutikov, M., & Götz, D. 2008, ApJ, 686, 1245

Ricci, R., Troja, E., Bruni, G., et al. 2021, MNRAS, 500, 1708

Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL, 887, L21

Rodríguez, Ó., Nakar, E., & Maoz, D. 2024, Natur, 628, 733

Ronchi, M., Graber, V., Garcia-Garcia, A., Rea, N., & Pons, J. A. 2021, ApJ, 916, 100

Rozwadowska, K., Vissani, F., & Cappellaro, E. 2021, NewA, 83, 101498

Safi-Harb, S., Burdge, K. B., Bodaghee, A., et al. 2023, arXiv:2311.07673

Sautron, M., Pétri, J., Mitra, D., & Dirson, L. 2024, A&A, 691, A349

Shivvers, I., Modjaz, M., Zheng, W., et al. 2017, PASP, 129, 054201

Skowron, D. M., Skowron, J., Mróz, P., et al. 2019, Sci, 365, 478

Smirnov, N. 1948, The Annals of Mathematical Statistics, 19, 279

Spitkovsky, A. 2006, ApJL, 648, L51

Spruit, H., & Phinney, E. S. 1998, Natur, 393, 139

Stovall, K., Lynch, R. S., Ransom, S. M., et al. 2014, ApJ, 791, 67

Tademaru, E., & Harrison, E. R. 1975, Natur, 254, 676

Tauris, T. M., & Manchester, R. N. 1998, MNRAS, 298, 625

Tejero-Cantero, A., Boelts, J., Deistler, M., et al. 2020, JOSS, 5, 2505

Totani, T., & Tsuzuki, Y. 2023, MNRAS, 526, 2795

Trigg, A. C., Burns, E., Roberts, O. J., et al. 2024, A&A, 687, A173

Usov, V. V. 1992, Natur, 357, 472

Vasist, M., Rozet, F., Absil, O., et al. 2023, A&A, 672, A147

Viganò, D. 2013, PhD thesis, University of Alacant

Viganò, D., Rea, N., Pons, J. A., et al. 2013, MNRAS, 434, 123

Wang, Z., Rea, N., Bao, T., et al. 2024, arXiv:2411.16606

Webb, N. A., Coriat, M., Traulsen, I., et al. 2020, A&A, 641, A136

Weltevrede, P., & Johnston, S. 2008, MNRAS, 387, 1755

Wu, Z.-F., Damoulakis, M., Beniamini, P., & Giannios, D. 2025, ApJL, 980, L28

Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29

Younes, G., Baring, M. G., Harding, A. K., et al. 2023, NatAs, 7, 339

Yusifov, I., & Küçük, I. 2004, A&A, 422, 545

Zane, S., Rea, N., Turolla, R., & Nobili, L. 2009, MNRAS, 398, 1403

Zhang, B., & Harding, A. K. 2000, ApJL, 535, L51

Zhou, P., Chen, Y., Li, X.-D., et al. 2014, ApJL, 781, L16