Novel Constraints on Axions Produced in Pulsar Polar-Cap Cascades

Axions can be copiously produced in localized regions of neutron star magnetospheres where the ambient plasma is unable to efficiently screen the induced electric field. As these axions stream away from the neutron star they can resonantly transition into photons, generating a large broadband contribution to the neutron star's intrinsic radio flux. In this work, we develop a comprehensive end-to-end framework to model this process from the initial production of axions to the final detection of radio photons, and derive constraints on the axion-photon coupling, $g_{a\gamma\gamma}$, using observations of 27 nearby pulsars. We study the modeling uncertainty in the sourced axion spectrum by comparing predictions from 2.5 dimensional particle-in-cell simulations with those derived using a semi-analytic model; these results show remarkable agreement, leading to constraints on the axion-photon coupling that typically differ by a factor of no more than $\sim 2$. The limits presented here are the strongest to date for axion masses $10^{-8} \, {\rm eV} \lesssim m_a \lesssim 10^{-5} \, {\rm eV}$, and crucially do not rely on the assumption that axions are dark matter.


I. INTRODUCTION
Axions are among the best-motivated candidates for physics beyond the Standard Model; these particles are a fundamental prediction of the leading solution to the strong CP problem [1][2][3][4], are an ideal candidate to explain the 'missing matter' in the Universe (i.e.dark matter) [5][6][7][8], and are expected to arise naturally in string theory from the compactification of gauge fields on topologically nontrivial manifolds [9,10].
One environment in which axion-photon mixing is particularly strong is the magnetosphere of a neutron star, where the large coherent magnetic field and the ambient plasma allow for highly efficient resonant transitions 2 .If axions contribute to the dark matter, one expects this conversion to generate narrow spectral lines that can be observed using radio telescopes [50][51][52][58][59][60][61][62][63][64][65][66]; this idea has ignited numerous observational efforts [50][51][52], with the most recent study setting world-leading limits for axion masses near ∼ 30 µeV [52].Despite the success and future promise of these spectral line searches, they are limited by the assumption that axions contribute significantly to dark matter, that the dark matter is smoothly distributed throughout the Galaxy, and that the Galactic dark matter halo is cuspy.Furthermore, such searches are confined to masses 10 −7 eV m a 10 −4 eV, as the mass must be sufficiently high to produce observable radio emission and sufficiently low so that resonances may be encountered.
Recently, Ref. [67] proposed an alternative way to detect axions in neutron star magnetospheres that overcomes the aforementioned challenges.The idea is based on relativistic axions sourced locally in the magnetosphere from the spacetime oscillations of E • B 3 .Ref. [67] showed that the electromagnetic fields are strong enough in the polar caps of neutron stars (situated above the magnetic poles, and spanning distances of order r pc ∼ O(10 − 100) meters -see Fig. 1) to produce an enor-mous flux of axions.As these axions traverse the magnetosphere they may encounter resonances, generating a broadband radio flux in the MHz − GHz regime.This axion-induced radio flux provides an alternative observable in the search for these evasive particles.
In this work, we construct the first pipeline to compute the intrinsic spectrum of axions produced in neutron star polar caps, their resonant conversion to photons, and the nonlinear evolution of these radio photons as they escape the magnetosphere.Our analysis uses both state-of-the-art numerical simulations as well as a newly developed semi-analytic model to predict the axion production rate; the overall agreement of these approaches illustrates that our procedure is not strongly sensitive to reasonable modeling uncertainties in the gap dynamics.We use this pipeline to constrain the axion-photon coupling by comparing the predicted radio flux with measurements of 27 nearby pulsars.The constraints derived here are the strongest to date for axion masses spanning 10 −8 eV m a 10 −5 eV.

II. AXION PRODUCTION FROM VACUUM GAP DISCHARGES
The e ± pair plasma populating the magnetospheres of neutron stars is expected to efficiently screen the component of the background electric field along the magnetic field lines (E ), except in small localized regions referred to as vacuum gaps which are responsible for particle acceleration and production of the pair plasma itself [70,71].
Vacuum gaps are expected to arise in a variety of locations, including in the polar caps [72], the slot gap (located near the neutron star along the last closed field lines) [73], and the outer gap (located near the light cylinder along the last closed field lines) [74].Recent progress in global particle-in-cell (PIC) simulations of the magnetosphere has shown that e ± pairs can be produced in the current sheet near the light cylinder, and efficiently screen gaps situated along the last closed field lines [75,76].We therefore choose to focus on the dynamics of the vacuum gaps in the polar caps, which also have the highest local values of E • B.
Recently, local PIC simulations of the neutron star polar-cap pair cascade have been performed [69,77,78].These simulations show that the discharge is an oscillatory process: e ± pairs accelerated in the gap produce gamma-rays that convert to more e ± pairs, which proceed to screen the gap, shutting down pair production [79].It was proposed that this process can produce coherent electromagnetic radiation, potentially answering the long-standing puzzle of the origin of pulsar radio emission [78].
The oscillatory pair discharge process will in general lead to inductively driven oscillations of E at a frequency set by the local plasma frequency, which is expected to evolve from the Goldreich-Julian (GJ) os-cillation frequency ω GJ = 4πα n GJ /m e4 to a value O(10 − 100) times greater as the density of the pair plasma increases [79][80][81][82].Here we have have introduced the GJ charge density, given by n GJ ≡ 2 B 0 • Ω NS /e, with B 0 and Ω NS being the surface magnetic field strength and rotational frequency of the neutron star, e the elementary charge, m e the electron mass, and α the fine structure constant.The presence of an oscillating E directly enters as a source term in the axion's equation of motion, The differential production rate of axions per momenta k can subsequently be expressed as [83] d Ṅ where ω( k) is the axion energy, T is the quasiperiodic timescale of the gap collapse, and ( k) is the Fourier transform (FT) of the source term, The rate and spectrum of axions produced during the collapse of the vacuum gap is fully determined by the spacetime evolution of E .The condition that axions be on shell inherently connects the spatial and temporal evolution of E , meaning a given k-mode can only be produced if the spatial component of the FT contains support on scales ∼ k −1 and the temporal component contains support on scales ∼ (k 2 + m 2 a ) −1/2 .In order to provide a quantitative understanding of axion production, we briefly estimate the production rate for a pulsar with B 0 = 10 12 G and Ω NS = 2π Hz.The polar-cap radius for such a pulsar is r pc ∼ 150 m, and the maximum value of the unscreened electric field is roughly E ∼ 6 × 10 −6 B 0 , see the Supplementary Material (SM).We expect production to be most efficient when E is largest, which occurs prior to the screening phase when the characteristic scale is k c ∼ 2π/r pc ∼ 10 −8 eV.Neglecting the phase in Eq. ( 3), and taking d 3 k/ω ∼ k 2 c , T ∼ 10 −7 s, and g aγγ ∼ 10 −11 GeV −1 , we find Ṅ ∼ O(10 50 ) axions per second, which is comparable to the values obtained from more careful calculations in the SM.

III. MODELING GAP COLLAPSE
We adopt two approaches to estimate the axion spectrum produced from the gap dynamics, one based on a semi-analytic model and the other on a numerical simulation.

A. Semi-analytic Model
We begin by highlighting the key physical features entering our 2+1D semi-analytic model; a more detailed description of each step is deferred to the SM.
We model two stages of gap evolution: pair cascade and gap collapse.The pair cascade phase begins with an initially unscreened E and low charge density.As the pair cascade progresses, the number density increases exponentially until it reaches n GJ .At this point the gap locally collapses, marking the end of the pair cascade.
A single seed particle (electron) is assumed to initiate the pair cascade.Because of the presence of E , the electron is accelerated to the radiation-reaction-limited Lorentz factor, γ max , in time t acc .Since particles move along curved field lines, they then emit curvature radiation (CR) photons with characteristic energy ε CR ∝ γ 3 max /ρ c , where ρ c is the field line curvature [80,84].CR photons may be absorbed by the magnetic field to produce new pairs and synchrotron photons.Newly produced electrons (positrons) are accelerated away from (towards) the neutron star surface, producing their own CR photons, which may again produce new pairs.Synchrotron photons can also produce new pairs if their mean free path is less than the size of the gap (see Fig. 1).
In our model, we compute the creation positions, times, and energies of all photons and pairs produced by the single seed particle.Our model iterates over 'generations' of particles, where the first generation is the seed particle, and generation n particles are sourced by generation n − 1.We run five generations of the cascade for different seed particle locations within the gap, identifying points where the plasma density locally reaches n GJ as the starting points for gap collapse -these points, which we refer to as 'pair production seeds', define the initial conditions of the semi-analytic model (see SM).
The initial stages of the gap collapse have been studied using analytic toy models and numerical simulations; in both cases, one expects the initial burst of particle production to induce exponentially damped oscillations in E .These oscillations initially occur locally, and subsequently propagate outwards along the magnetic axis [69,78,82,[85][86][87][88][89][90][91].The frequency of the damped oscillations is set by the local plasma frequency (which itself is set by the charge density and typical Lorentz factor), and will evolve from ω GJ (i.e. the minimum value necessary to collapse the gap) to a value 10 − 100 times greater as more pairs are produced.The characteristic screening timescale t c is set by the time required for pair production processes to yield a GJ charge density, which is roughly equivalent to the time required to accelerate charges to γ max [82]; for realistic pulsars this timescale lies between 10 −9 − 10 −6 s.
In order to capture the general features of the gap collapse process, we model the unscreened electric field E with a static profile, and describe the screening of E as the combination of outward propagating 2D plane waves.Each plane wave is exponentially damped on a timescale t c , and has a time-evolving oscillation frequency growing from an initial value ω GJ .In total, we source four propa-gating plane waves5 evenly spaced across the width of the gap, and with initial conditions determined by the onedimensional pair-production process discussed above.

B. PIC Simulation
Our second model used to compute the axion spectrum relies on a 2.5 dimensional6 PIC simulation developed in [81].Working in axisymmetric cylindrical coordinates (r, z), the authors impose a dipolar magnetic field.A rotating disk of radius r pc is established at the stellar surface to produce a potential drop in the open field line zone of the neutron star.Outside r = r pc , E is forced to zero to model the plasma-filled closed field line zone.Particles are extracted from the surface at a rate that depends on the local value of E ; these particles are accelerated to the radiation-reaction limit and emit gamma-rays through synchrocurvature radiation.Those gamma-rays subsequently produce e ± pairs in the ultrastrong magnetic field through one-photon magnetic pair production [92].This process is modeled using the stateof-the-art QED module in the PIC code OSIRIS [93].Videos showing the dynamical screening of E are available online; a snapshot of the simulation is shown in the left inset of Fig. 1.
Performing simulations from first principles of the gap collapse process is extremely challenging due to the large separation of scales between the size of the polar cap and the kinetic scale of the plasma (typically differing by 4 − 5 orders of magnitude).The simulations performed in [81] overcame this difficulty by re-scaling the quantum parameters χ ±,γ = (p µ F µν ) 2 /(B c m e ) for both photon emission and pair production, multiplying them by a numerical constant; here, p µ is the 4-momentum of the corresponding e ± or γ, and B c 4.4 × 10 13 G is the Schwinger field strength.This re-scaling effectively allows pair production to occur at a much lower voltage drop.However, it also significantly reduces the inherent scale separation in the problem: the kinetic plasma length scale becomes only a few hundredths of the polarcap size.As a result of this compression of scales, the FT necessarily compresses the power to an artificially narrower range of k-modes (the largest scales are expected to be unaffected, however power in small scales, i.e. large k-modes, will have been shifted to intermediate scales).We describe in Section I of the SM a procedure for rescaling the power of the FT, along with additional details on the PIC simulation.
Once the initial axion spectrum has been determined, we employ an updated version of the ray-tracing algorithm developed in [65] to compute the radio spectrum.We describe the general features of this procedure below, and defer further details to the SM.
We begin by calculating the rate of axion production for momenta spanning from the escape momentum to ∼ O(10 GHz) (the phenomenology of gravitationally bound axions differs markedly, and thus we leave a detailed study of bound states to a companion paper [94]).For each momentum state the corresponding axion trajectory is propagated to the light cylinder, and all resonances encountered during propagation are identified.Resonances occur when [65,95] where θ is the angle between the axion momentum and the magnetic field.At every level-crossing we compute the conversion probability, P a→γ , and axion survival probability7 , P a→a = 1 − P a→γ , and subsequently weight each photon sourced at crossing i by P (i) a→a .The conversion probability at an individual level crossing is given by the Landau-Zener formula [62] with Here, we have introduced ∆ ≡ sin θ/(1 − ω 2 p cos 2 θ/ω 2 ) and defined [95] ∂ with k ( k⊥ ) representing the parallel (perpendicular) direction to the axion momentum.Photons are propagated to the light cylinder using the dispersion relation for a Langmuir-O mode (see [65]).The observed spectrum is obtained by summing over the final distribution of binned and weighted photons -example spectra are included in the SM.

V. RESULTS AND DISCUSSION
The radio emission mechanism of active pulsars is not well understood, making it difficult to identify signatures arising from this process.Nevertheless, one can constrain the existence of axions without knowing the intrinsic pulsar flux by requiring that their contribution not exceed the observed flux.In this section we derive limits on g aγγ using observations of 27 nearby pulsars; our sample includes representative nearby pulsars, whose surface magnetic field and rotational period roughly span 10 12 G B 0 10 13 G and 10 −1 s P NS 2 s.For computational ease, we choose to focus on pulsars whose radio emission geometry is constrained by observations, thus evading the need to marginalize over the misalignment and viewing angles.Details of all pulsars used in this analysis are presented in the SM.
We present our fiducial 95% confidence level upper limits in Fig. 2, which are obtained by computing the con-straints for both models (a comparison between models is left to the SM) and taking the weaker limit at each mass.The bands around the fiducial limit represent a conservative estimate of the systematic uncertainties (see SM).For comparison, we plot constraints from radio line searches using neutron stars (blue), axion haloscopes (gray), the CAST experiment (teal) and X-ray and gamma-ray telescopes (light green).The limits derived in this work significantly improve upon existing bounds, and unlike axion haloscope experiments (and radio line searches), do not assume axions contribute to the dark matter.In addition, since the radio flux scales ∝ g 4 aγγ , the constraint is largely insensitive to minor mismodeling errors.The mass range covered by our constraints is limited by the frequency of radio observations (higher frequencies could probe higher masses), and the computational expense (computing time increases at both lower and higher masses).
In the SM we show that the derived bound is controlled by observations of a few strong pulsars, with high frequency observations providing the most constraining power.A comprehensive analysis of all pulsars in the ATNF catalog, as well as more dedicated pulsar observations at high frequencies, could significantly improve upon these results; we reserve this broader analysis for future work.

ACKNOWLEDGMENTS
We would like to thank Georg Raffelt, Andrea Caputo, Ben Safdi, and Anatoly Spitkovsky for useful discussions.The authors would also like to thank Andrea Caputo  , where the PIC simulations used in this work were performed.We acknowledge the use of [99] in creating the figures containing axion constraints.
In this section we illustrate specifics concerning the computation of the initial axion spectra.We begin by discussing in detail the relevant PIC simulation developed in [81], and then describe the ingredients used in the semi-analytic model.

A. PIC Simulation
We use the results from one of the 2.5 dimensional QED particle-in-cell simulations developed in [81] and performed using the code OSIRIS [93,100].Here, we outline the details of this simulation, introduce the re-scaling methodology, and explain how we extrapolate the simulation results back to more realistic regimes.
The simulation is performed in axisymmetric cylindrical coordinates (r, z), where r = 0 is the rotation axis and z = 0 is the stellar surface (ignoring the curvature of the surface) modelled as a conductor.The stellar magnetic field is taken to be a dipole with magnetic axis aligned with the rotation axis, and magnitude given by B(r = 0, z = 0) = 0.1 × B c with B c the critical magnetic field strength (conventionally, B c 4.4×10 13 G, however in the simulation B c is reduced to 10 7 G to make the computations feasible).To mimic the potential drop induced by the rotation of the neutron star, we rotate a disk of radius r pc on the stellar surface; r pc is known as the polar-cap radius and given by where R NS is the neutron star radius and Ω NS its rotational frequency.The polar-cap radius is chosen to be 0.1R NS in the simulation, which corresponds to Ω NS of a few hundred Hz.Rotation of the disk is achieved by imposing a boundary condition on the radial electric field E r (r, z = 0) = E 0 × (r/r pc ) × g(r/r pc ), with E 0 = Ω NS B 0 R NS and g(x) = 0.5 × (1 − tanh((x − 1)/0.2)).The latter function has been chosen to allow for a smooth transition from the open to the closed field line zone.The upper boundary of the computational domain is taken to be open for both particles and fields [69,93].The plasma-filled closed field line zone is modelled as a perfect conductor where E = 0.The simulation is performed on a grid that spans L r × L z = 1.5 r pc × 2.5 r pc , with resolution 6000 × 10000.This resolution, however, does not allow us to resolve the realistic scale separation between the polar-cap size r pc and the plasma skin depth λ p , especially when pair production creates a large number of e ± particles.In order to resolve the physics correctly, the plasma skin depth corresponding to the GJ charge density, λ GJ = ω −1 GJ , is increased to a value of ∼ 5 × 10 −3 r pc , roughly O(100) times larger than expected in realistic neutron stars.This re-scaling decreases the overall potential drop induced in the open field line bundle.Nevertheless, in order to model the pair production physics the electrons must still be allowed to create pairs at the lower energies.This is achieved by also scaling the quantum parameters χ ± and χ γ , which govern the emission of gamma-ray photons and the rate of conversion to e ± pairs.These quantum parameters are defined by χ ±,γ = (p µ F µν ) 2 /(B c m e ) -here p µ is the 4-momentum of the corresponding particle, F µν the electromagnetic field tensor, and m e the electron mass -and are multiplied by a factor of ζ ±,γ = 10 3 in the simulation to enhance pair production at lower energies.The simulation proceeds as follows.At each time step in the simulation charges are injected at the surface with an effective charge density proportional to the local parallel electric field, n inj = κ(E /er pc ), with κ = 0.1 and e the unit of elementary charge.As the particles are accelerated from the surface, their quantum parameters χ ± are individually calculated, multiplied by the scaling factor ζ ± , and then used to emit synchrocurvature photons with energies sampled via Monte-Carlo methods.Upon emission of a synchrocurvature photon, the energy of the photon is subtracted from that of the emitting particle.This is a discrete implementation of the 'radiation-reaction force' (see the following subsection for additional discussion).The emitted photons are also propagated in the simulation, their quantum parameters χ γ are computed, multiplied by the scaling factor ζ γ , and then used to compute the quantum Breit-Wheeler pair production cross-section.The photon is afterwards conditionally converted to an e ± pair based on the resulting pair production probability.The initial conditions of the simulation are set such that there are no particles along the open field lines.The conducting disk is gradually spun up, releasing particles which subsequently initiate quasiperiodic particle cascades.When performing the FT, we remove the initial part of the simulation, and focus only on the quasiperiodic solution achieved in the later stages.
The re-scaling of fundamental parameters mentioned above introduces a complication -the plasma scale is artificially inflated with respect to the size of the system (instead of being separated by ∼ 4 − 5 orders of magnitude, they are only separated by ∼ 3 orders of magnitude).When calculating the FT, this re-scaling shifts power from k-modes near k ∼ λ −1 GJ to a value O(100) times smaller.In order to correct for this effect we compute the FT  at a log-spaced series of momentum modes ( kmin , • • • , kmax ), and subsequently re-stretch the momentum vector k i = (ξ) i/n ki before computing the differential rate in Eq. ( 2).Here k min and k max are determined by the largest and smallest scales resolvable by the semi-analytic model, and n is the total length of the series.The extrapolation factor ξ is found by computing the ratio between the GJ scale predicted using the neutron star parameters and the value inferred from the ratio of λ GJ /r pc used in the simulation.
Finally, note that, due to the parameter re-scalings highlighted above, one must re-scale the value of E || /B 0 from the simulation to each pulsar under consideration.In order to do this, we use the fact that NS .In the simulation we analyzed, the parallel electric field, E s || , corresponds to an effective rotational frequency of Ω s,eff NS ∼ 200 Hz.Keeping the above in mind, we use the following re-scaling relation to determine the parallel electric field, E p || , for any real pulsar

B. Semi-analytic Model
There are two primary reasons to jointly develop a semi-analytic model.Firstly, the simulation discussed above was developed for a particular neutron star.We have applied it to the more generic population by re-normalizing the polar-cap radius and the maximum amplitude of E to their theoretical values, but it is not guaranteed that this re-normalization is sufficient to capture all relevant features.The second concern is that the re-scaling procedure introduced in the simulation to reduce scale separation may generate unexpected features in the FT.In order to address both of these concerns, we have built a semi-analytic model that can: (i) be applied to the broader class of pulsars and (ii) be run on a higher resolution grid, significantly reducing (although not entirely eliminating) the need for post-analysis stretching of the FT.We describe the details of this semi-analytic model below.
Pair cascade: An initially unscreened vacuum gap becomes unstable to runaway pair production.Due to variation in the magnetic field strength and line curvature across the gap pair production occurs non-uniformly.Once the local charge density becomes comparable to the GJ value dynamical screening of E takes place, leading to exponential damping of the parallel electric field.The dynamics of E in 1+1 dimensions are described in [82,91] and summarized in the following section.In order to extend these results to 2+1 dimensions, we construct a 2+1 semianalytic model of pair cascades in the gap.The model provides the initial conditions for sourcing the aforementioned damped oscillations.
We use axisymmetric cylindrical coordinates with the lower boundary being a disk of radius r pc rotating at frequency Ω NS .The initial parallel electric field is constant throughout the gap with magnitude E = 2Ω NS B 0 r pc , where B 0 is the surface magnetic field.The domain of the gap is also permeated with a constant magnetic field, B 0 ẑ8 .We further approximate the pair cascade initiated by a single seed particle to develop along a single magnetic field line, labeled by its radial coordinate s.The characteristic radius of a flux tube containing the particles produced in a cascade is of the order ψ ∼ 1/ε γ , where ψ is the angle between the photon momentum and the local magnetic field (i.e.pitch angle) and γ is the photon energy (in units of m e ).Low energy photons that can have larger pitch angles generally have mean free paths greater than the size of the gap, preventing them from contributing to the cascade.
We inject macroparticles of charge κn GJ at the neutron star surface, where κ 1. Positrons that are accelerated away from the surface reach a maximum Lorentz factor determined by the balance between acceleration and radiative losses to curvature radiation according to the equation where + (−) applies to positrons (electrons), ρ c (s) = (4/3)R 2 NS /s is the magnetic dipole curvature at radius s, and γ(t) is the particle Lorentz factor.The first term on the RHS corresponds to acceleration by the electric field and the second term to radiative losses due to emission of curvature photons.From Eq. (S3), the radiation-reaction-limited Lorentz factor is γ max = 3E 0 ρ c (s) 2 /2e 1/4 .A key quantity is the amount of time taken for a particle with initial energy γ 0 (in units of the electron mass) to reach energy γ max .This acceleration time is largely independent of the initial Lorentz factor (see Fig. S1).We define it as the amount of time it takes a particle to reach γ = (1 − )γ max , which can be determined by integrating Eq. (S3) to yield where 1. Going forward we set = 0.1.Once a particle reaches γ max , it emits CR photons with characteristic energy ε CR = 3γ 3 max /2ρ c and initial pitch angle ψ 0 = 0.As the photon propagates its pitch angle, and hence opacity, increases.The condition for absorption by the magnetic field is [80] where a defines the mean free path and α B is the opacity for pair production by a single photon in a strong magnetic field.Once the CR photon is absorbed it produces an e ± pair, each particle having energy ε γ /2.Since ψ a ≡ ψ( a ) is non-zero at absorption the pair will be produced at a high Landau level and quickly decay to the ground state, emitting a number of synchrotron photons in the process.The final momentum (and energy assuming they are ultrarelativistic) of each fermion is given by the component of the parent CR photon energy parallel to the background magnetic field prior to absorption with χ ±,a ≡ χ ± ( a ) and b ≡ B/B c .We assume the remaining energy is split equally among synchrotron photons of characteristic energy ε syn = 3bψ a ε γ /8.Once produced, the mean free path of synchrotron photons is calculated in the same way as for the CR photons, with the exception that they inherit the pitch angle of the parent CR photon.The secondary electron (positron) is emitted upward (downward) along the magnetic field line and reaches γ max in time t acc 9 .The secondary electrons (positrons) can now emit upward (downward) propagating CR photons and the process continues.We study the above model using an iterative algorithm that works as follows.We define four categories of particles: electrons, positrons, upward-propagating (γ ↑ ), and downward propagating (γ ↓ ) photons.For each e ± we keep track of its creation position (z ± create ) and time (t ± create ) and for each photon its creation position (z γ create ), time (t γ create ), energy (ε γ ), and initial pitch angle (ψ 0 ).We assume that particles reach γ max at time t create + t acc , where t acc is given in Eq. (S4), at which point acceleration is balanced by emission of CR.We assume CR photons are emitted at regular intervals δt CR = ε CR /P CR , with P CR = 2e 2 γ 4 max /3ρ 2 c being the CR power.Thus for each electron (positron) in generation n we create new upward (downward) propagating CR photons according to The characteristic number of CR photons emitted by each electron/positron is given by the transit time across the gap divided by the time interval between emission of CR photons, For each photon (indexed by i, with i = 1, 2, . . . ) in generation n we compute the mean free path to pair production as described above and produce generation n + 1 pairs and synchrotron photons at the point of absorption.We run five generations of particle production as described in Fig. S2.The output of the algorithm gives information about the particle density as a function of position and time within the gap.We use this output to determine locations at which the charge density equals n GJ , which serve as source points for gap collapse as described in the following subsection.
Screening of E : The 2+1 dimensional semi-analytic model is constructed by fixing the spatial profile for the unscreened gap E ( r), and mimicking the screening behavior induced by pair production processes using damped plane waves -notice that this is motivated by the qualitative features observed in the simulations of [69].The physical scales driving the damping and oscillations of the plane waves are inferred using the simplified 1+1 dimensional analytic solutions, as well as the insight gained from the 2.5 dimensional PIC simulations.The procedure followed is outlined below.
We start by defining the physical dimensions of the vacuum gap.The radius of the polar cap r pc is set by the typical footprint produced from a dipolar magnetic field (see Eq. (S11)) and the height of the gap h can be estimated by following the trajectories of e ± particles (which travel along magnetic field lines), and determining the point at which pair-production generated a GJ charge density (as described above); this height is approximately given by [72] h 50 m ρ c 10 6 cm 2/7 where ρ c is the curvature of magnetic field lines near the neutron star surface.For a dipolar magnetic field, the radius of curvature is given by [101] where θ is the colatitude of the footprint of the magnetic field line, θ pc the colatitude of the polar-cap boundary, and P NS the rotational period of the neutron star.In determining the gap height we set θ to a constant characteristic value of θ pc /2 [101].
The size of E may vary within the gap itself, growing linearly along the magnetic axis and decaying off at radial distances near r pc .In order to incorporate smooth spatial boundaries we define a gap profile ψ(r, z) such that the maximum unscreened electric field is given by E ,max (r, z) = E ,max × ψ(r, z), with E ,max ≡ e n GJ r pc .We take the gap profile to be with δr = r pc /10.Notice that the radial dependence of ψ is somewhat ad hoc (although similar functional forms have been adopted e.g. in [69]).In reality, we simply require a quick and smooth fall-off of the gap near r pc , and we expect other types of smooth functional cut-offs to provide comparable results.The linear growth of ψ on the other hand is intended to mimic the inability of the extracted surface current to screen E || , the efficiency of which should decrease as one moves away from the surface (see e.g.[77,78]).We emphasize that the precise details of this gap profile do not have a large impact on the axion spectra.
The initial conditions of the semi-analytic model are that of a fully grown gap (i.e.E = E ,max (r, z)).At this point, pair production is expected to collapse the gap, sending E → 0; the nonlinear interplay between the motion of the plasma and the response of the electric field, however, induces damped oscillations in E , with each phase of the oscillation triggering the production of higher generation e ± pairs.The damping is expected to be exponential, with a typical timescale roughly given by [82] where λ is the charge multiplicity.The frequency of oscillation is set by the local plasma frequency.We assume here that the plasma frequency during the initial response is set by the GJ value (i.e. the minimal value required in order to suppress the electric field), and increases by a factor of 100 as subsequent generations of particles are produced (see e.g.[82]).The growth of the oscillation frequency should be slightly less than exponential, since only a subset of each generation's e ± pairs contribute in producing the subsequent generation (at this point, pair production is driven only by the high energy tail of the distribution).As such, we adopt a time-dependent evolution of the oscillation frequency given by Here β ≡ (ω f − 1) × e −2α , and we set α = 0.5 and ωf = 100.In the following section we illustrate that this choice does not have a major impact on the derived constraints.
As an informative illustration, we show in Fig. S3 the temporal evolution of the 1+1 dimensional solution identified in [82] (computed using the same parameters as used in their Fig. 2) as a function of (re-normalized) time t, and using the frequency evolution in Eq. (S14).In the figure, we extend the result beyond the nonlinear damping regime, including the linear oscillation and decay regime identified in [82].We note that the latter two regimes are not significant for this work, as the solution at these late times is only expected to be representative of small localized plasma over-densities located well outside of the polar cap.
The fundamental difficulty in modeling the gap collapse processes arises from the nonlinear response that appears in two dimensions.Particle production itself is inherently inhomogeneous since pair production depends heavily on the radius of curvature of the magnetic field lines, which varies across the gap.The response of the electromagnetic fields to the inhomogeneity in the particle production further induces nonlinear dependencies, which significantly complicate the evolution of E .Here, we attempt to capture the initial inhomogeneity in the pair formation front by accelerating test particles and tracing the radiation and initial pair production events (see section on pair cascade above); this procedure neglects the back-reaction on the electromagnetic fields, and thus is only valid when E is large.
We model the subsequent nonlinear screening of E by drawing well-isolated samples from the pair formation front, and inducing two-dimensional damped plane waves that propagate outward -we refer to these samples as 'pair production seeds', with coordinates denoted by (r seed,i , z seed,i ).Note that similar behavior can be seen in the simulations of [81].We approximate the evolution of E using the Ansatz where the product runs over each of the sampled seeds.Our intent with this model is to capture the typical temporal and spatial variations of E .We will validate it below by comparison with PIC simulation results.The response of each seed is given by We take k = ω (i.e.assume it travels at the speed of light), and the functional dependencies on spacetime are understood to be implicit.We have introduced in Eq. (S16) the notion of an 'effective time' ti which captures the relative spacetime dependent response of E (in other words, this captures the relative delay that would arise in the event that gap collapse occurs locally and expands outwards, as seen in [81]); this is given by ti (r, z, t) where t is the absolute time, t pp,i is the initial time at which particle production event i was initiated, and d pp,i is the relative distance from the center of particle production event i to the position (r, z).Since the electric field will push the newly produced plasma out of the gap, we adopt a drift in the initial seed location such that The factor t delay is set by the minimum time required for the plane waves to collapse the gap at all radii, which is also roughly the time over which the pair formation front forms, i.e. max(t pp,i ) − min(t pp,i ) (note that the inclusion  S1), an axion mass of 10 −8 eV, and an axion-photon coupling of gaγγ = 10 −11 GeV −1 .
of t delay is necessary to ensure that the gap collapses).We run each semi-analytic model for a total time set by t max = 2 × max(t pp,i ), which is sufficiently long to observe the gap collapse and the subsequent outward drift of the seed points; moreover, this number agrees quite well with the expected gap collapse period of typical neutron stars, which is around 10 −9 − 10 −6 s [91].
In order to be able to perform the full FT across the momentum range of interest one must have sufficient spatial and temporal resolution.Specifically, we require a spatial resolution at the level of the GJ oscillation wavelength (which is well-beyond the resolution needed to resolve momentum modes k ∼ O(10 GHz)), and temporal resolution given by where ω max is the maximum axion energy computed.For most of the pulsars and masses of interest, storing the gridded evolution of the gap at this resolution requires O(100 Gbs) to O(Tbs) of memory.
In order to make these computations more feasible, we perform the semi-analytic calculations with a slightly enhanced value of the GJ scale λGJ = ξ × λ GJ , with ξ taken to be 6 (roughly the smallest value which we can achieve given current resources).Similarly to the case of the PIC simulation, this modification effectively inflates small scales and thus mis-places the power of the FT in smaller k-modes (here, we resolve around 4.5 orders of magnitude in scale, but require around 5 -thus the effect is rather minimal).As before, we attempt to correct for this by first computing  at a series of log-spaced momentum modes ( kmin , • • • , kmax ), and subsequently extrapolating the momentum vector ki before calculating the differential rate in Eq. (2).To illustrate the sensitivity to this extrapolation factor, we plot in Fig. S4 the axion spectrum predicted using ξ = 6, 10, and 20 for pulsar B1737-30 and an axion mass of 10 −8 eV (the axion-photon coupling is set to g aγγ = 10 −11 GeV −1 ).In the left panel of Fig. S4 we also show for comparison the corresponding spectrum produced by the simulation.One can see that the relative difference between the simulation and semi-analytic model greatly exceeds the small differences that appear as one changes the re-scaling factor.All runs are computed at 35 log-spaced k-values in both k ρ and k z .In SM Section V, we return to the effect of this re-scaling factor and show that taking a value of ξ = 10 makes only minimal difference in the derived limits, thus justifying the adopted re-scaling factor used for the fiducial analysis.

II. COMPUTING THE RADIO FLUX
From the spectra computed using the PIC simulation and the semi-analytic model, we can find the initial spectrum of axions produced in the gap.With the initial axion spectrum in hand, one can subsequently calculate the radio flux generated via resonant axion-photon transitions through: (i) carefully following the trajectories of all axions created in the polar caps, (ii) identifying the resonant level crossings, (iii) computing the conversion probabilities from axions into photons, and (iv) propagating the produced photons through the dispersive magnetosphere to find their direction and energy at asymptotic infinity.The ray-tracing code premiered in [65] was developed to treat the resonant conversion and propagation of radio flux arising from axion dark matter around neutron stars; we have modified this code to include axion propagation in a Schwarzschild background and incorporated a generalization of the mixing which is valid for (ultra-)relativistic axions.Both of these new features, as well as other relevant steps in computing the axion-induced radio flux, are outlined below.

A. Axion Production Rate
In this section we review the derivation of Eq. ( 2) from Eq. ( 1).We follow the discussion in [83].In the particular case of electromagnetic production of axions the classical source is given by j(x) = −g aγγ ( E • B)(x).The average production rate can be calculated assuming the source is non-vanishing only over a time period T , which in our scenario corresponds to the gap collapse time.In the asymptotic past the axion field can then be expanded as where a † k a k is the creation (annihilation) operator for the axion and ω( k This can be thought of as the homogeneous solution to Eq. ( 1) or, equivalently, the so-called 'in' state.The effect of the source, after it has been turned off, is to shift the creation and annihilation operators where ( k) is the Fourier transform of j(x).It follows from Eq. (S21) that the average rate of particle production is given by where |0 is the asymptotic vacuum state.

B. Axion propagation
Since axions are produced close to the neutron star (where gravity is strong), we propagate axions assuming they follow Schwarzschild geodesics 10 .The relevant equations describing the axion trajectories are Here we have chosen to work in polar coordinates, with τ specifying the proper time and G being the gravitational constant.It is important to note that both the angular momentum and total energy are constant along the trajectories and given by where γ is the standard Lorentz factor, v ⊥ is the perpendicular component of the local particle velocity and r s is the Schwarzschild radius of the neutron star.
The axions are propagated in just two dimensions, as their trajectories are inherently symmetric about the magnetic axis of the neutron star.The locations of the conversion points, however, are not.Consequently, we project each calculated axion trajectory in 70 different randomly sampled azimuthal directions when computing the conversion points and probabilities.This amount of samples was chosen to maximize stability and accuracy in the final results while maintaining computational tractability.After ray-tracing, we average over the rotational period to compute the expected radio flux.

C. Axion-photon conversion
In its most general form, the axion-photon transition probability is given by the Landau-Zener equation (S27) This properly describes conversion in the adiabatic (i.e.large coupling and/or large magnetic field) regime, assuming that the length scale over which the mixing takes place is sufficiently small with respect to variations in the magnetic field and plasma density 11 .Details in calculating the form of the adiabaticity parameter Γ from the electrodynamic equations of motion can be found in [65].Here, we generalize the aforementioned analysis to include: (i) relativistic mixing and (ii) the latest approximation of the axion-photon conversion probability from [95].
Resonant axion-photon transitions occur when k a k γ [65], where the axion momentum is given by As in the main text, ω is the axion/photon energy, ω p is the plasma frequency, and θ is the angle between the axion momentum and the external magnetic field B. In the non-relativistic limit, this resonance happens at m a ∼ ω p ; the generalized condition (which is also applicable in the relativistic regime), however, has the form Recently, Ref. [95] pointed out that off-diagonal derivatives (e.g.∂ k , k⊥ E ⊥ ) in the axion-photon mixing equations modify the conversion probability.In particular, this contribution shows that the group velocity of the sourced electromagnetic mode travels in the direction ŝ, which is defined through its differential Here k and k⊥ are the parallel and perpendicular directions to the axion momentum respectively (note that we have defined k⊥ to lie in direction of B, which generates a sign difference with respect to [95]), and Since the axion-photon conversion probability depends on change of the photon momentum along the direction of the group velocity, one must derive the correction to the adiabaticity parameter arising from these off-diagonal derivatives; the result is given by [95] 11 Recent work in [102] has shown that Eq.S27 breaks down when the conversion is adiabatic and when these scales are comparable.For small axion masses, we do find that some photons are sourced from resonances in which Eq. S27 does not hold; however, removing these photons entirely leads to a modest correction to the constraints on the order of 10% for axion mass ma = 10 −9 eV and 1% for ma ≥ 10 −8 eV, implying our constraints are robust to this uncertainty.We caution that lower mass axions are likely to be more severely affected by this problem, and thus a careful treatment would be required in this regime.

D. A brief comment on adiabatic axion-photon conversion
It was recently pointed out in [52] that axion-photon conversion in the magnetospheres of highly magnetized neutron stars can become extremely efficient, generating order one conversion probabilities.In the context of axion dark matter (in which axions are assumed to fall onto the neutron star), this can lead to an exponential suppression of the radio flux -this is a consequence of the fact that each axion is expected to pass an even number of resonances, and thus the generation of an observable radio signal requires a produced photon to survive a level-crossing during escape from the magnetosphere (and P γ→γ = 1 − P a→γ 1 when the conversion probability is large).We emphasize here that the localized axion production in the polar caps naturally evades this suppression as these axions typically undergo an odd number of level crossings, making this a highly complementary probe of axions in the large coupling limit.

E. Example photon spectra
Before providing an illustration of the photon spectra produced through our models, we reiterate the full procedure used to obtain these spectra.This procedure includes: • We first compute the axion production rate across a range of momentum bins using either the re-scaled PIC simulation or the semi-analytic model.
• We then use these sampled momenta as initial conditions for our ray-tracing algorithm.Axions are traced away from the gap, and their resonances along the trajectory are identified.
• For each resonance, we compute the probability of sourcing a radio photon (which depends on the conversion probability at the resonance, as well as the survival probability of the axion at previous resonances) and ray-trace the photon to large distances.
• After all photons have been traced, we determine the relative viewing angle of Earth, θ ⊕ , and compute the average radio flux by summing over the weighted contributions of photons with final positions θ ∼ θ ⊕ ± (with taken to be very small and with photon frequencies falling in the desired bin [ν 1 , ν 2 ]. We provide an illustrate of the photon spectra produced by pulsar B1737-30 in Fig. S5, where the left panel shows a comparison between the predictions of the simulation and the semi-analytic model (at three viewing angles), while the right panel displays the observed flux density at five viewing angles (computed using only the semi-analytic model).In both cases, the axion mass is set to m a = 10 −8 eV and the axion-photon coupling to g aγγ = 10 −11 GeV −1 .Fig. S5 highlights two important points: there is reasonable agreement between the spectrum predicted by the PIC simulation and that of the semi-analytic model, and the observed flux density depends rather strongly on the observed viewing angle.S1.List of pulsars used in this work.For each pulsar the columns denote: name, rotational period (in seconds), inferred dipolar surface magnetic field B0 (in Gauss), the angle θm between the magnetic-and rotational axis (in degrees), the angle β between the magnetic axis and the observer at closest approach (in degrees), the flux density between 25 − 80 MHz [103], 50 − 80 MHz [103], 408 ± 4 MHz [104], 1.4 ± 0.032 GHz [104], and 8.6 ± 0.8 GHz [105] (all in mJy), and the distance (in kpc).The inferred geometric angles are obtained from [106][107][108].Distance measurements and errors are taken from [109][110][111][112]; following the ATNF procedure, we prioritize measurements not inferred using the dispersion measure and Galactic electron models.In some cases (marked with a †), the Galactic coordinates and dispersion measure are used to infer the distance utilizing the online tools from [112].Should the various Galactic models agree, we take the mean and standard deviation of these inferred distances; if they instead differ by more than a factor of 2, we defer to the more recent model of [110] and place a 100% error on the distance.

III. PULSAR DATA
In order to constrain the existence of axions, we compare the predicted radio flux with observations from 27 nearby pulsars.We have limited ourselves to this small carefully selected sub-sample in order to make the analysis computationally feasible.Our selection is largely based on two fundamental criteria: we want pulsars with large magnetic fields and a well-measured orientation (meaning that both the misalignment angle and the pulsar's orientation relative to Earth are known).The former is particularly important as the production of axions scales with ∝ B 4 0 and the resonant conversion process scales with ∝ B( r c ) 2 ( r c being the point of resonant conversion).We have chosen to apply the geometrical constraint as the axion-induced radio flux varies rather strongly with viewing angle (see e.g.Fig. S5), and knowledge of the geometry evades the need for marginalization.It is worth mentioning that a broader study over the entire pulsar population would naturally marginalize over the geometry, and could yield stronger constraints than those obtained here; currently, such a study is computationally prohibitive, however we hope that future improvements will make this manageable.
The radio observations were obtained using the LOFAR telescope [103], the Lovell telescope at Jordell Bank [104], and the Shanghai Tian Ma Radio Telescope [105], and span a wide range of frequencies from 25 MHz to 8.6 GHz12 .The properties of each of the 27 pulsars and the observed flux densities are listed for reference in Table S1.Below we outline the analysis used to derive upper limits on the axion-photon coupling, and highlight the impact of various systematic uncertainties in the modeling of the axion signal.

IV. PROFILE LIKELIHOOD ANALYSIS
We utilize the profile likelihood ratio to derive a one-sided upper limit on g aγγ .In general, this is accomplished by defining a test statistic [114][115][116] Here µ is the parameter of interest (in our case, g aγγ ), θ denotes the nuisance parameters (which here include the pulsar distances and the intrinsic radio flux from each pulsar at each frequency), (μ, θ) are the maximum likelihood estimators, and θ is the conditional maximum likelihood estimator for a given value of µ.
We adopt Gaussian likelihoods for the flux density and pulsar distances.The full likelihood is given by where the product is taken over each pulsar and each frequency bin for which observations are available (see Table S1).Some important comments regarding this equation: • The predicted flux density has two components, one due to the predicted axion signal and another that is intrinsic to the respective pulsar, i.e. S pred,ij = S axion,ij + S pulsar,ij .Owing to the fact that we do not currently have a firm understanding of the intrinsic radio flux produced by pulsars, we treat the parameters S pulsar,ij as nuisance parameters.This implies that our test statistic is zero, unless the predicted axion flux in at least one observing bin of one pulsar exceeds the observed flux.This procedure ensures that our constraints remain conservative.
• All observed flux densities and mean distances, as well as the variances in both of these quantities, are reported in Table S1.Note that the variance in the distance is generally asymmetric, meaning that the second exponential in Eq. (S34) is actually an asymmetric Gaussian.
• The pulsar distances show up explicitly in the likelihood, but also enter through the predicted flux densities.
With the likelihood in hand, we compute our test statistic at set couplings and derive an upper limit on g aγγ by identifying the value at which TS = 2.71.This corresponds to the 95% upper limit [116].
The derived limits for the simulation and semi-analytic model are shown in Fig. S6, where the lines represent the fiducial constraints and the bands reflect the uncertainty due to particular modeling choices -the bands are also shown in Fig. 2 and discussed thoroughly in the next section.We find that the final limits obtained through our pipeline are largely driven by observations of a few pulsars at higher frequencies.The reason for the latter is that the pulsar flux tends to drop much faster than the predicted axion flux with increasing frequency.In order to illustrate the origin of the constraining power, we plot in Fig. S7 the constraints derived from one of the strongest pulsars (in the semianalytic model) using exclusively observations of low (∼ 408 MHz), mid (∼ 1.4 GHz), or high (∼ 8.6 GHz) frequencies (left panel).We also plot in the right panel of Fig. S7 the individual constraints obtained from both models using the three strongest pulsars (utilizing data at all frequencies).We find that the relative differences between the two models is typically no larger than a factor of ∼ 2 in coupling (although larger differences appear for a small amount of outliers, like pulsar B0656+14 at intermediate masses), providing confidence that reasonable changes to our axion production model have a minimal impact on the final constraints.

V. UNCERTAINTIES
This section reports the various uncertainties that enter our calculations, including uncertainties associated with the inferred pulsar properties, uncertainties in treating the PIC simulation, and uncertainties associated with the free parameters entering the semi-analytic model.We also provide quantitative estimates of the effects due to these uncertainties on the final constraints.

A. Observational uncertainties
Let us begin by discussing uncertainties associated to the inferred properties of each pulsar.In general, the properties that are measured to high precision for every pulsar are the rotational period P NS and the spin-down rate ṖNS .The value of the surface dipolar magnetic field B 0 can be inferred from these quantities.This can be seen by comparing the net Poynting flux through the light cylinder of a misaligned rotator, given by L = µ 2 Ω 4 NS (1 + sin 2 θ m ) (where µ = B 0 R 3 NS , Ω NS = 2π/P NS , and θ m is the misalignment angle [117]), to the spin-down luminosity, yielding There are two points worth noting here.First, the ATNF adopts a fiducial moment of inertia I NS = 10 45 g cm 2 , which is actually 30-50% lower than current estimates [118].Second, the ATNF inferred field strengths assume θ m = 0.One can see that these effects somewhat offset, however likely lead to a slight underestimation of B 0 at the level of ∼ O(10%).
Another important consideration is that Eq.S35 only reflects the magnitude of the dipolar magnetic field -near the surface of the neutron star additional higher order multipoles could, and in some cases are expected to, contribute.This likely implies that the values of B 0 listed in Table S1 underestimate the true surface magnetic field strength (and thus we are likely underestimating our signal, potentially by a significant amount).Incorporating higher multipoles, however, must be done self-consistently, and thus in what follows we adopt the assumption that the magnetic fields are purely dipolar.We believe that this is likely a conservative assumption, however we plan to address this more rigorously in future work.
Since the derived limit is expected to scale like g lim aγγ ∝ B α 0 with 1 α 1.5 (depending on the efficiency of the resonance), we adopt in what follows a characteristic uncertainty on g lim aγγ of ±20%.We expect this to be conservative, as it underestimates the impact of the corrected moment of inertia and neglects the contribution of higher order multipoles, both of which are expected to enhance the signal.
Notice that in estimating the uncertainty in the dipolar field strength we have neglected the contribution coming from the unknown equation of state, which allows for ∼ 10% fluctuations in R NS [119].This is because the neutron star radius enters at multiple points in the analysis, and thus we treat this as an independent uncertainty.In particular, we expect the neutron star radius to enter with the following scalings: NS , and B 0 ∝ R −3 NS .Roughly speaking, this translates into a scaling on the limits as NS , where as before 1 α 1.5.Thus, a 10% variation in the neutron star radius amounts to a 20% shift in the limits.
There exist a variety of techniques for inferring the distances to observed pulsars.By far the most reliable of these are parallax measurements, but unfortunately such measurements are only available for a small number of known pulsars.One can also analyze the HI spectrum on-and off-pulse, and use this information to infer the characteristic rotational velocities of the HI regions behind-and in front of the pulsar.In many cases, however, this information is not available.Without additional information, one must rely on inferring the distance using the observed dispersion measure (DM = d n e ), which requires a model for the Galactic electron density.In Table S1 we compile inferred distances and uncertainties using various combinations of the techniques listed above.Priority is given to geometric measurements, and we only rely on dispersion measure distances when no other information is available (as done in the ATNF) [109][110][111][112].In some cases, no published distance could be found; here, we used a variety of Galactic electron density models to find the distances using the tools available in [112].In a few cases, uncertainties on the distance were not readily available.In these cases, we compared the inferred distance from two different Galactic electron density models -if the models showed reasonable agreement (differing by no more than a factor of 2), we inferred the standard deviation from these measurements.In all other cases we defaulted to the more recent model [110] and placed a 100% uncertainty on the distance.It is worth emphasizing that none of these pulsars (denoted in Table S1 with a †) contribute meaningfully to the derived constraints.
Finally, let us mention that the uncertainties on the pulsar geometry (the misalignment angle and orientation) as interpreted within the rotating vector model (RVM) [120] for the pulsars used in this analysis are expected to be small FIG.S8.Upper limits derived using observations of B1737-30 and the semi-analytic model for the axion production rate.We asses the impact of varying several unconstrained parameters (namely Npts, the final value of ω, α, and ξ) on the derived limits, and compare with the fiducial model (solid, gold).
(see e.g.[107]).Note that this is not a generic feature of all pulsars, however we have chosen our pulsars largely to minimize this uncertainty.There is of course a larger systematic uncertainty that pertains to the validity of the model itself, but it is unclear how this can be quantified.Given that the RVM is the state-of-the-art, we take these values to be fixed in our analysis -future analyses which are applied over large populations of pulsars may be able to evade the need for specific assumptions on pulsar geometry.We leave such an endeavour for future work.

B. PIC Simulation
There are two clear identifiable uncertainties that enter our analysis of the PIC simulation.First, the physical dimensions of the simulation are expressed in terms of the polar-cap radius.In general, the fractional surface area which leads to pair cascades depends on a variety of factors including the magnetic field geometry, the misalignment angle, the properties of the return current, the compactness of the neutron star, and various nonlinearities in the electrodynamics.While it is difficult to quantify the effects of each of these separately, it is expected that together they can lead to an O(1) modification of the polar-cap radius [77,121,122] 13 .Axion production at low masses scales with the square of the volume (with a weakening of the dependence on the volume at higher masses), meaning we expect our limit on the axion photon coupling to scale proportional to this uncertainty.In what follows we therefore adopt a characteristic uncertainty on r pc at the level of ±30%14 , which roughly corresponds to the one-sigma range inferred from a flat distribution varying between ±50% of the fiducial value.
The second source of uncertainty comes from the re-scaling factor ξ, which is intended to undo the compression of scales that is performed in the simulation.Unfortunately, without more simulations at various resolutions it is impossible to truly understand the impact of this re-scaling factor.Given, however, that the semi-analytic model also contains a slight scale compression (and subsequent re-scaling), we adopt the uncertainty inferred from the semianalytic model in the simulation as well (see following subsection).We do note that this procedure is not exactly valid given that the scale compression is not done in an equivalent manner -it is largely for this reason that the semi-analytic model was constructed in the first place.
Finally, it is worth mentioning that there also exists an unquantifiable uncertainty that enters the analysis of the PIC simulation, which arises from the fact that the physical parameters used in the simulation do not directly correspond to those of true pulsars (recall that this is because the parameter re-scalings performed in the simulation modify the fundamental scales of the problem).In the future, one would ideally like to perform a number of simulations varying the re-scaling parameters in order to understand the sensitivity of axion production to this procedure, however at the moment this is computationally unfeasible.Thus, we instead use the semi-analytic model as a cross-check, allowing us to ensure the robustness of the results.

C. Semi-analytic Model
As with the simulation, the semi-analytic model is subject to uncertainties in the characteristic polar-cap size.Additionally, the semi-analytic model has a number of free parameters which have been fixed to fiducial values in derivation of the upper limits.Here, we outline the impact of reasonable variations in each of these.The primary free parameters include: N pts (the number of seed points used in the plane wave damping), ωf (the final plasma frequency at the end of the particle cascade), α (the power law index controlling how quickly the plasma frequency evolves after the initial collapse), ξ (the re-scaling factor), (the parameter controlling the acceleration time of particles in the gap), and κ (the parameter describing the initial plasma density present in the gap prior to the pair cascade, note that this parameter also enters the simulation).
The latter two parameters (namely, and κ) have no perceptible effect on the final results.In the semi-analytic model, the acceleration time (see Eq. (S4)) is the smallest length-scale in the problem.This means that for the fields considered in this work, the acceleration time is often smaller than the time resolution of the model.As a result, we generally accelerate each particle to its final Lorentz factor within a single time step.The value of , which describes at what fraction of the radiation-reaction limited Lorentz factor each particle begins to emit CR photons, shows up logarithmically in the acceleration time and thus has at most an O(1) effect on t acc , which is already small compared to the time resolution.Therefore, even orders of magnitude changes in have negligible effect on the final result.The parameter κ is introduced to avoid having to simulate ∼ n GJ r 3 pc particles, which would be computationally infeasible.Since the number density grows exponentially during a cascade, κ can be understood as a logarithmic re-scaling of the time in the pair cascade process.However, we note that the pair cascade process is introduced to simply set the initial conditions in the gap prior to dynamical screening.The initial conditions and the dynamics of the screening phase are independent of κ, and thus the effect of κ on the final result is also negligible.
In order to assess the sensitivity of the presented limits to the choice of the other parameters we run an independent analysis on one of the most-constraining pulsars, B1737-30.In particular, we perform an extra run with N pts = 3, one with N pts = 5 (rather than N pts = 4), a run in which ωf = 10 (rather than 100), two runs in which α is varied from 0.1 to 1.0, and various runs using a larger stretching parameter ξ.The limits obtained in each of these runs are displayed alongside our fiducial limit from this pulsar in Fig. S8.The characteristic uncertainty of each parameter variation on the limit is quantified and included in Table S2.Note that in the special case of ξ, we infer the uncertainty by looking at the limit in which ξ → 1 rather than by comparing variations across runs.
Most of the uncertainties listed in Table S2 are extremely sub-dominant to that coming from the uncertainty associated to the polar-cap size, with the variation in N pts otherwise producing the largest effect.The total uncertainty for each model (estimated at the level of 36% and 45% for the simulation and semi-analytic models, respectively) is obtained by summing the errors in quadrature, and is reflected in a band about the fiducial limits in Fig. 2 and Fig. S6.

FIG. 1 .
FIG. 1. Schematic figure showing axion production in neutron star vacuum gaps.The vacuum gap is depicted by a truncated cone on the neutron star surface.The left inset shows a time snapshot of E (from the simulations of [69]), with the brown/green coloring reflecting negative/positive values of E .The right inset depicts the microphysical processes responsible for the pair cascade, with green arrows indicating the direction the cascade flows with time.Axions (red) are emitted from the gap and convert to photons (purple) in the presence of the neutron star's magnetic field, B (gray).
FIG. S1.Electron (solid) and positron (dashed) energies as a function of distance propagated along a curved magnetic field line with E0 = ΩNSB0RNS.The neutron star parameters are taken to be b ≡ B/Bc = 0.1, ΩNS = 2π Hz, and ρc = 100 km.
FIG. S2.Schematic depiction of the generations of particle creation in our model.Photons are labeled as either upward propagating (↑) or downward propagating (↓).Each electron/positron produces many CR photons (see Eq. (S9) for an estimate).Each CR photon, when absorbed, produces a single electron/positron pair and several synchrotron photons.

2 ||
FIG. S3.Illustrative evolution of various phases (excluding growth phase) of the evolution of E , following the 1+1 dimensional formalism of[82].The inset shows a zoom-in of E 2 during the the decay phase.

FIG. S4 .
FIG. S4.Number of axions produced per unit time predicted by the simulation (left panel) and semi-analytic model (right three panels) for various extrapolation factors ξ, ranging from 6 to 20.The spectra are produced for pulsar B1737-30 (see TableS1), an axion mass of 10 −8 eV, and an axion-photon coupling of gaγγ = 10 −11 GeV −1 .
FIG. S5.Left: Comparison of photon spectra produced by pulsar B1737-30 using the PIC simulation (solid) and the semianalytic model (dashed).Results are shown for three different viewing angles.Right: Photon spectra produced by pulsar B1737-30 using the semi-analytic model over a wider range of viewing angles.All results assume an axion mass of 10 −8 eV and an axion photon coupling gaγγ = 10 −11 GeV −1 .
FIG. S6.Fiducial upper limits derived for the simulation (solid) and semi-analytic model (dotted), shown together with corresponding bands depicting systematic model uncertainty.

TABLE S2 .
Parameter Characteristic uncertainty on g lim aγγ Simulation Semi-analytic Estimated breakdown of the characteristic uncertainty on the derived limits arising from systematic uncertainties in the semi-analytic model and simulation (we denote in each column whether the listed uncertainty is included and relevant for each model).The total uncertainty is obtained by combining the individual contributions in quadrature, and is reflected in the final row of each model.