{'fraction_non_alphanumeric': 0.03684477667035806, 'fraction_numerical': 0.013819675156884459, 'mean_word_length': 4.426264396594893, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Axion configurations around pulsars

Axion-like particles lead to a plethora of new phenomena relating to compact astrophysical objects including stellar and black hole superradiance, axion stars and axion clusters. In this work, we investigate a new scenario in which macroscopic axion configurations are sourced by the electromagnetic fields of pulsars via the axion-photon coupling. We solve the inhomogeneous axion field equation with an explicit source term given by the electromagnetic fields of the pulsar in a rotating magnetic dipole approximation. We find that the axion profile either forms a localised boundstate or radiates as outgoing waves, depending on whether the pulsar frequency is smaller or greater than the axion mass, respectively. We derive the total mass of the scalar cloud generated around the pulsar, and the power loss for radiative solutions. We point out that it will be necessary to incorporate pulsar magnetosphere effects and their partial screening of the axion-photon interaction for a more accurate quantitative prediction. Finally, we suggest some observational signatures which should be investigated in future work.


Introduction
Axion physics currently represents one of the most active areas at the interface of particle theory, astrophysics and cosmology. Motivated by the axion solution to the strong CPproblem [1][2][3] and axion-like particles arising from string theory [4,5] we consider a generic scalar field which interacts with electromagnetism (in Lorentz-Heaviside units µ 0 = ε 0 = 1) via the following Lagrangian minimally coupled to gravity: whereF µν = ε ρσµν /2F ρσ . Using F µνF µν = −4E · B, the axion equation of motion becomes This paper seeks to make the first steps towards understanding a simple question: what is the nature of scalar clouds sourced by an E · B component of the strong electromagnetic fields of pulsars, and what might its observational signatures be?
In comparison to many other astrophysical sources, the electromagnetic fields of pulsars [16] can reach some of the strongest known values in the Universe, with field strengths of up to 10 15 G in the case of magnetars [17]. This makes pulsars an excellent playground in which to study the axion-photon coupling, whose strong magnetic fields can compensate for the large supression of g γφ . It should be noted that we consider macroscopic axion field configurations, rather than individual scattering processes in the cores of stars [18] used to place stellar cooling bounds on axion parameters. Thus our solutions are more in the vein of other classical axion configurations around astrophysical objects which we now discuss.
There already exists a wide array of novel astrophysical predications should the axion turn out to exist in nature. These are of merit not only as tools to further constrain the physics of ultralight scalar degrees of freedom but also as astrophysical phenomena in their own right. Perhaps the most widely studied effect involving axions and compact astrophysical objects is superradiance, in which scalar fields extract rotational energy from a spinning black hole [19][20][21][22][23][24] or star [25][26][27] via an instability for massive perturbations. Other notable signatures relate to pulsar timing, in which pressure oscillations from a background axion field distort spacetime, causing an apparent frequency change in pulses [28]. Similarly axion stars [29][30][31][32] present an entirely new class of astrophysical objects whose detection has recently been considered in [33,34]. Clusters [35] also represent yet another kind of axion-based astrophysical object whose radio signatures have been examined in [36]. It is remarkable to think that such a rich array of phenomena arise from something as simple as a scalar field and it is in this spirit that we explore axion clouds around pulsars.

Contrasting with superradiance
It is worth remarking that the nature of axion solutions explored in this paper is qualitatively different from discussions of superradiance in the context of black holes and stars mentioned in the preceding paragraph. The superradiant scenario essentially describes a homogeneous wave equation for a particle "trapped in a box" which exhibits a fundamental set of eigenfrequencies ω = ω R + iω I with an imaginary part ω I whose sign is such that solutions grow exponentially in time. Crucially, the scalar field does not corotate and has ω = Ω, where Ω is the frequency of the spinning object. In fact the superradiant constraint is more stringent, requiring ω R < mΩ, where m is the azimuthal quantum number of the mode in question. The physical picture underlying this mechanism is that rotating objects can dissipate energy onto an incoming wave, leaving the scattered wave with more energy so that |A out |/|A in | > 1 where A out and A in are the amplitudes of the outgoing and incoming waves, respectively. If one repeats this process ad infinitum by placing a "mirror" near the black hole or star, the incoming wave comes back with more energy each time. The analogue of such a mirror is to give the field a non-vanishing mass, which has the effect of confining it to a region near the black hole.
By contrast, in our setup, the scalar field continually extracts energy from the E · B source term in the inhomogeneous wave equation (1.2). This generates a configuration which, under certain circumstances (specifically when M φ < Ω), consists of a purely outgoing wave (|A in | = 0) at infinity as the axion solution radiates away energy from the system. Furthermore, by linearity, the scalar field contribution sourced by the electromagnetic fields inherits the same frequency as the pulsar so that by contrast with the superradiant case, we have instead a corotating solution ω = Ω. From a calculational point of view then, we must go further than simply solving a homogeneous wave equation with particular boundary conditions. We shall actually derive Green functions and integrate them with a source to construct axion profiles. It is also worth noting that the scenario considered here has the advantage that it probes the axion-photon coupling directly.

Summary of results
• We compute axion solutions to the inhomogeneous wave equation (1.2) by a method of Green functions, with the source given by the electromagnetic fields of a pulsar in an oblique rotating magnetic dipole approximation. Solutions are shown in fig. 1.
• The maximum density of solutions grows with decreasing axion mass. For M φ > Ω, solutions are localised near the pulsar with an amplitude falling off exponentially at large distances. For M φ < Ω they radiate as outgoing axion waves. For very low masses M φ Ω, densities tend to a maximum value ( fig. 3) with the largest axion densities (see eq. (5.3)) reaching ∼ 10 −1 g cm −3 for typical neutron star and axion parameters. We find that in the extremal case, for the largest known pulsar magnetic fields of B * ∼ 10 15 G and the upper bound for the axion photon copuling g γφ = 10 −13 GeV −1 , the total mass of the cloud can reach ∼ 10 −13 M , a comparable size to certain axion stars and clusters.
• In the radiative regime M φ < Ω, we compare the axion power output to the magnitude of electromagnetic dipole radiation -see eq. (6.8) and fig. 5.
The structure of this paper is as follows. In the following two sections we introduce our basic model, beginning with the gravitational metric in sec. 2. We then outline the structure of electromagnetic fields around the pulsar in sec. 3 in a simple rotating magnetic dipole model. There we explain how the axion field satisfies an inhomogeneous wave equation with a source term derived from the dot product of the pulsar's electric and magnetic fields. In sec. 4 we describe axion solutions to this wave equation, deriving the radial component of solutions and offering an approximate analytic solution in sec. 4.4. We then compute the total mass of the cloud in sec. 5 which for reference we compare to the typical masses of axion stars. Figure 1: Axion contours sourced by the electromagnetic fields of the pulsar in the vacuum dipole model (VDM), as seen in the M -Ω plane. We took M φ R = 0.1 and r s /R = 0.3. The left column shows a boundstate with ΩR = 0.01 whilst the right shows radiating solutions with ΩR = 0.5. The upper two rows show the axion radial function (solid) together with its analytic approximation (dashed) given by neglecting gravity and far field corrections to EM fields described in sec. 4.4. In the bottom plots, the magnetic field lines are shown in grey.
We also describe the total power radiated by the axion cloud in sec. 6 and compare this to the power of electromagnetic dipole radiation. In sec. 7 we suggest how our stellar modelling could be improved to incorporate plasma around the pulsar and how this will localise the axion cloud to regions in which E · B is poorly screened by the pulsar magnetosphere. Finally in sec. 8 we offer our conclusions and give some possible observational effects which should be pursued in future work.

Gravitational Metric
In our choice of gravitational metric we shall neglect rotational effects (i.e. frame dragging etc.) that would be encountered in e.g. a Kerr metric [37]. This has the benefit of greatly simplifying the form of the wave equation of the scalar field. We shall therefore consider a Schwarzschild geometry matched to an interior solution corresponding to uniform stellar density. The Schwarzschild metric is a good approximation provided the ratio of the star's total angular momentum to mass J/M , is much smaller than its radius R, in which case Kerr gives the limiting case of the Schwarzschild metric. As a rough approximation, treating the pulsar as a sphere of uniform density one has J = 2 5 ΩM R 2 , so that J/M = 2 5 R(ΩR). Note that ΩR is precisely the tangential speed of the pulsar at the equator and therefore must be less than the speed of light. For typical pulsars, R ∼ 10km and Ω 10 2 Hz so that ΩR 1. In this case, J/M R making the Schwarzschild metric an excellent approximation for most pulsars.
Explicitly then, our metric is of a Schwarzschild form matched to a perfect fluid of uniform density in the stellar interior [38,39] so that the full metric takes the form where r s = 2GM is the Schwarzschild radius of the star. Using the relation it is easy to see the axion wave equation then satisfies where (r, θ, ϕ) are spherical polar coordinates. Consider a monochromatic source for the wave equation of the form where Y m (θ, ϕ) are spherical harmonics. In accordance with eq. (2.4) we shall consider axion solutions of the form Of course, since the wave equation (2.4) is linear, a monochromatic source will only excite axion modes of the same frequency and quantum number as the source, so that, up to the addition of a homogeneous solution arising from some other unspecified physics, the axion cloud should have the same time dependence as the source. However, in a regime where axion self-interactions are important, one could have further axion modes excited. One could envisage a situation in which the source first fills some quantum numbers, which then begin to fill other harmonics at different frequencies due to axion self-interactions, as happens in superradiance and bosanova effects [40,41]. For instance, in the case of a cosine potential non-linear effects would become important if the axion cloud reached densities such that φ ∼ f a . However, for the sake of progress we shall consider the simplest scenario in which we focus only on a linear setup with corotating axion solutions of the form (2.6), whose substitution into (2.4) gives the following radial wave equation for the ψ m (r): where primes denote derivatives with respect to r and This can be thought of as a Schödinger equation with an extra source term with effective potential given by V (r) = e 2Ψ l(l + 1) The homogeneous part of the radial equation has been well-studied in the Schwarzschild case [42][43][44][45][46] we consider here and also in its extension to Kerr backgrounds [47][48][49][50], the main interest being in superradiance and quansinormal modes. There are many approaches to solving the homogeneous equation [51], ranging from series solutions to numerical shooting methods. In general though, with the exception of a few special cases, no simple analytic solutions are known and one must resort to some form of numerics to get the full solutions.
At first glance, one might also worry that by neglecting Kerr effects, one is missing important instabilities of the scalar field arising from its extraction of rotational energy from the star -so called superradiance [23] encountered in the context of rotating black holes. However, this mechanism arises from the special form of boundary conditions which are imposed on the wave equation at the black hole horizon. It is this boundary condition which gives rise to a quantised set of eigenmodes exhibiting an imaginary frequency part which leads to instabilities which grow exponentially in time.
Since stars do not possess a horizon, the wave equation in its simplest form does not give rise to superradiant effects of this nature. However, it is worth noting that if there is some source of dissipation for the scalar field then superradiance can occur for stars. For instance, if a friction term term αφ is added to the left hand side of (1.2) (where α is some dissipation constant), then superradiant solutions are possible, as shown in [25]. In principle, such a dissipative term could arise from quantum corrections to the classical equations of motion. For instance, it could be generated by the imaginary part of the one-loop photon correction to the axion propagator, which by the optical theorem corresponds to φ → γγ decays, in analogy with the mechanism described in ref. [52]. Similarly, for dark photons around a star, subject to some version of Ohm's law, a finite conductivity provides a mechanism for dissipation [26,27], in which case dark photon solutions exhibit superradiance provided they have a non-vanishing mass.
In the present case, we shall not explore these possibilities and consider the simplest setup, consisting only of a wave equation with a source E · B. Our calculation will proceed as follows. For the exterior solution, we shall generate the axion profile by first constructing a Green function for the radial equation in the region r ≥ R, which we then integrate with the source terms to get the full exterior solutions, subject to appropriate boundary conditions at infinity. In these exterior solutions it remains to fix a boundary condition at the stellar surface, which we do by imposing continuity of the solution at r = R. For black hole spacetimes, one imposes a stringent boundary condition by demanding that modes at the horizon should be purely ingoing, so that the horizon acts as a one-way membrane. This gives rise to a set of discrete eigenfrequencies -the quasinormal modes of the black hole [49]. Since stars do not possess a horizon, the imposition of such boundary conditions at the stellar surface would be artificial. Instead, we should demand only that the solutions satisfy appropriate continuity conditions at the stellar surface and match smoothly onto an interior solution.

Pulsar electromagnetic field structure
Next we must choose a model for the electromagnetic fields inside and outside the pulsar. In order to solve the axion equation of motion, we shall make two simplifying assumptions. Firstly we assume that the interior of the star is a perfect conductor, in the sense that the electric field E in the rest frame of the star vanishes: Using the transformation laws for electromagnetic fields in rotating frames [53,54], the fields E and B as seen in the inertial frame are related to E by where r is the position vector with respect to an origin at the centre of the star. Eq. (3.2) then implies that the source of the axion wave equation (2.4) vanishes in the stellar interior: Of course in reality, the axion experiences in medium effects due to the nuclear matter inside the star which could modify its potential, as has recently been examined in [55]. This could change the profile of the interior axion solution. However, for the sake of making progress we leave such considerations for future work. Next we discuss the structure of the electromagnetic fields outside the star. We shall choose perhaps the simplest, and one of the earliest models of pulsar magnetic fields. This is the classic vacuum dipole model (VDM) of Deutsch [56] which, in accordance with the assumption above of a perfectly conducting stellar interior, models the pulsar by a rotating dipole, inclined by an angle α relative to the axis of rotation. In this setup, one neglects any plasma surrounding the pulsar, so that the external electromagnetic fields are in vacuo. The effects of plasma surrounding a pulsar and its implications for E · B are discussed briefly in sec. 7.
Working in Coulomb gauge with ∇ · A = 0 and A 0 = 0, the vacuum equations for the gauge potential become with the electric and magnetic fields are given by The VDM is characterised by the rotating dipole equations where M(t) is the rotating, inclined dipole moment. Explicitly one has (see e.g. [57]) it follows that in the "near zone", characterised by rΩ 1, the second term in (3.7) and the last two terms in the expression for (3.8) are sub-leading, giving the following near-zone approximations where hats denote unit vectors and is the value of the magnetic field at the dipole on the stellar surface. The solutions (3.10) can also be derived by taking retarded time t − r → t in (3.6), so that they correspond to the neglect of retarded effects and electromagnetic radiation. One can find an explicit expression for M(t) by solving the second constraint in (3.6). This can be achieved by taking a dipole aligned with Ω (which we choose as the z-axis) and performing two successive SO(3) rotations, the second of which depends on time. One first takes a vector of constant magnitude M , pointing in the z-direction and rotates it by an angle α relative to the z axis -this makes the dipole oblique. Then one sets it spinning by performing a rotation about the z-axis through an angle Ωt. This can be written as Inserting the expression (3.12) into (3.7) and (3.8), one finds, after a lengthy but straightforward calculation that which we recognise as the azimuthal coordinate in the rest frame of the star. There is therefore no explicit t dependence other than that contained in λ, from which we see that the electromagnetic fields corotate with the pulsar. Notice that the terms O(rΩ) in the arguments of the trigonometric functions correspond to electromagnetic radiation, giving outgoing spherical waves ∼ e iΩ(r−t) , such effects are therefore neglected in the near field limit. Note also that E · B consists of a stationary axisymmetric part and time-dependent part with non-trivial azimuthal dependence. We notice immediately that the expression (3.13) can be rewritten in terms of the harmonics so that where we have absorbed the explicit e −iΩt dependence into the spherical harmonics via λ.
The profile of E · B is shown in figure 2.

Axion Solutions -Solving the radial equation
In accordance with our discussion in sec. 2, a source of the form (3.16) corresponds to an axion solution The prefactor depending on axion and stellar parameters is inserted to simplify the final form of the radial equation. It also means that the radial functions ψ m are dimensionless. Explicitly one finds the radial equations (with = 1) where V is the effective potential given by eq. (2.9) and θ(r−R) is the heaviside step function. The m = 0 mode corresponds to a stationary axissymmetric solution, whilst m = 1 gives a corotating configuration which violates axisymmetry via a non-trivial ϕ dependence. Notice that the m = 0 mode therefore depends only on the axion mass (and not Ω) whilst for m = 1, the radial function depends depends both on Ω and M φ . Notice that the radial equation (4.2) can be written purely in terms of dimensionless variables by taking (r, M, Ω) → (r/R, MR, ΩR), so that the shape of the axion profile is determined purely by the two dimensionless parameters ΩR and MR.

Exterior axion solution
The exterior wave equation is inhomogeneous and can be solved by method of Green functions satisfying We shall construct the Green function in the standard way [44] from two linearly independent mode functions f and g which satisfy the homogeneous version of (4.3). This leads to a Green function of the form (see appendix A) We have also defined Notice that the m = 1 mode has distinct behaviour depending on relative size of Ω and M φ . For M > Ω, the quantity k is real and positive, whilst for M < Ω we see k is pure imaginary. The boundary conditions (4.5) and the sign choice in (4.6) therefore ensures that the full solutions given by (4.4) either decay or give purely outgoing waves at infinity in the cases Ω < M φ and Ω > M φ , respectively. In the latter case, |k| should be thought of as the 3-momentum of the outgoing waves at infinity, whilst for the Ω < M φ case it sets the scale over which solutions decay. Thus as expected, axion radiation can be produced by an m = 1 mode.
To compute the mode functions f and g we solve the homogeneous version of (4.3) and impose the boundary conditions (4.5) at sufficiently large r = r max . We perform a consistency check of these numerical solutions by checking their Wronskian matches the analytic result W = 2k/A(r). We also test the mode functions against their flat space limit r s = 0 for which there are exact analytic expressions given in sec. 4.4.

Interior axion solution
The interior solution satisfies a homogeneous wave equation where the metric components take the interior form given by eq. (2.2). Since r = 0 is a regular singular point, one can use Frobenius' method (see Appendix B) to write solutions in the form where the prefactor R −s ensures both ψ andψ are dimensioneless.ψ then satisfies r 2ψ + r (P + 2s)ψ + [s(s − 1) + sP + Q]ψ = 0,ψ(0) = 0, (4.9) with We have obtained the interior solution by solving equation (4.9) numerically with appropriate boundary conditions and arbitrary overall normalisation by choosingψ(0) = 1. The remaining normalisation of ψ int , is fixed by imposing continuity at the stellar surface, as we now describe.

Full solution: matching at the stellar surface
The full solution is then given by where β gives the overall normalisation of the solution in the stellar interior and α = γ ∞ R drf (r)j(r)/A(r). The source j corresponds exactly to the right hand side of eq. (4.2). The constant α and β are uniquely determined by continuity of ψ and its first derivative at the stellar surface. This leads immediately to where With the equations (4.11) and (4.12), we have now fully determined the solution for the axion field both inside and outside the pulsar in the VDM, together with the assumption of an interior Schwarzschild metric corresponding to a stellar interior of constant mass density. The steps involved in computing the full solution thus go as follows: for each value of M and ω, compute the homogeneous mode functions f and g (in our case numerically) and compute the interior homogeneous solution viaψ − . Then insert f and g into the integrals with the source in (4.11) as well as equation (4.13). Then one can use these to determine the matching coefficients via (4.12). This should then be repeated for each mode ( , m) after which one will have completely determined the axion solution for a given Ω and M φ .
Notice that in the radiative regime (Ω > M φ ) the solution at infinity consists of purely outgoing waves, unlike in superradiance which involves continual reflection off the central object. The solution takes the form ψ = Ae i|k|r r → ∞, (4.14) with Notice that the amplitude of the outgoing wave is essentially determined by what is going on near the stellar surface, since every quantity except the integral is evaluated at r = R and even that is dominated by this region, since the source j decays as 1/r 4 .

Approximate analytic solutions: near field, flat space limit
We can also derive approximate analytic solutions for the axion profile by neglecting gravitational effects (r s = 0) and restricting to the near field limit (3.10) corresponding to rΩ 1. where

The axion equation then becomes
where γ E is Euler's constant and κ 1,2 are constants determined by demanding continuity of u and u at the stellar surface given by Plots of the different solution types are shown in fig. 1 together with a comparison with the approximate analytic result (4.17).

Energy density and total mass
The energy density of the scalar field is given by where T µν is the energy momentum tensor for φ. The largest values of the energy density correspond to the scenarios in which ψ is greatest. This occurs for RM φ RΩ 1 as shown in fig 3. In this limit,φ 2 ∼ Ω 2 φ 2 and M 2 φ 2 are much less than (∇φ) 2 ∼ R −2 φ 2 , in which case the energy density is dominated by the gradient term, so that, from (4.1) Thus near the surface of the star r ≈ R, O(ψ 2 ) ≈ 10 −2 (see fig 3), so that the largest axion densities are characterised by ρ axion = 9.5 × 10 −2 g cm −3 B * 10 13 G 4 · g γφ 10 −13 GeV It is easy to see then that for typical pulsars, this can be considerably more than the ambient halo density of 10 −24 gcm −3 [6], even when g γφ is considerably below the exclusion boundary of roughly 10 −13 GeV −1 .
The energy density, δρ(r), contained in a spherical shell lying between r and r + dr is given by δρ(r) = drr 2 dΩT 00 . (5.7) The total mass of the axion cloud is plotted in figure 4 for the maximal case of a magnetar and the upper bound on the axion photon coupling.
It is interesting to compare this to other astrophysical objects composed of axions, for instance, dilute axion stars examined in [30] have a range M/M ≤ 10 −13 , and axion clusters have similar masses at around M/M ∼ 10 −13 [35]. What is of particular interest is the density of our axion configurations in comparison to axion stars in [30]. Our densities are of order ρ ∼ 10 −1 gcm −3 , for typical pulsars, slightly more dense than dilute axion stars shown in fig. 1 of ref. [30] which have ρ 10 −2 gcm −3 . The main advantage we have is that whilst existence of axion stars remains open, configurations around pulsars are a simple consequence of the axion-photon coupling. In this sense, they will inevitably be generated in regions around the pulsar in which E · B is sizeable. Furthermore, one has an axion cloud in the vicinity of strong magnetic fields, which could have interesting observational consequences [6] in terms of decays and stimulated decays via the axion-photon vertex.
This comparatively high axion density could have important implications in terms of observational signatures. However it is worth bearing in mind that a quantitative analysis will require a better understanding of the E · B profile under the influence of plasma screening in the magnetosphere of the pulsar.

Axion power output
A further interesting question is how much power the axion configuration is extracting from the pulsar. This can be determined by computing the flux of the momentum density Π i ≡ T 0i through a sphere at radius r. The momentum density of the axion field is given by Π =φ∇φ. (6.1) The energy flux per unit area across a sphere of radius r is therefore Note that f has dimensions [f ] = Energy · Area −1 · time −1 and is different from the axion number flux. The energy flux per solid angle is therefore given by Inserting the solution (4.1), the energy flux (6.3) averaged over the pulsar period is therefore given by * Ω 3 128 sin 2 2α sin 2 θ Im ψ ψ * . (6.4) Notice therefore that there is only radiation if ψ 11 is complex, which requires Ω > M. In this case, from eq. (4.14) one has which is constant, so that the radiation given by (6.3) falls off as 1/r 2 , as required by energy conservation. The power radiated by the axion cloud is therefore given by integrating this quantity over a sphere at infinity, which gives We now contrast this with the electromagnetic dipole radiation given by integrating the Poynting vector S = E × B (recalling we set µ 0 = 1) over the sphere at infinity, which gives the power in dipole radiation From this, it follows that This ratio is plotted in fig. 5. In principle, since the derivativeΩ is determined by dissipation processes, as the pulsar loses rotational kinetic energy, one could use pulsar spindown measurements to place bounds on axion parameters. However in order to make reliable quantitative statements in this direction a better modelling of the screening of E · B in the pulsar magnetosphere will be needed.

Magnetospheric screening of E · B
The vacuum dipole model [56] considered in previous sections assumes there is no plasma surrounding the neutron star, and that the electromagnetic fields satisfy the vacuum Maxwell equations (neglecting axion backreaction). However, this represents only a first approximation to a realistic pulsar model. Firstly, the VDM has only a trivial emission spectrum -there is only monochromatic dipole radiation with frequency identical to that of the pulsar and this ultra-low frequency radiation cannot be observed directly. That said, dissipation of rotational kinetic energy via dipolar radiation does play a role pulsar spin down [58]. In reality, observations show pulsars have a rich emission spectrum in observable frequency ranges arising from the acceleration of charged particles in their magnetic fields. Thus a realistic model of pulsars must take into account the collection of surrounding charges, currents and their corresponding electromagnetic fields, which together form the pulsar magnetosphere. Recent reviews of pulsar magnetospheres can be found in [59][60][61], and an excellent and condensed history of the key developmental stages is given by [62].
Even in the early days of pulsar theory it was pointed out by Goldreich and Julian [63] that the electromagnetic forces on particles at the stellar surface exceed that of gravity by many orders of magnitude. As a result, particles will be inevitably lifted from the stellar surface and into the surrounding atmosphere. Early work [63] treated the magnetosphere as perfect conductor in the following sense. As explained at the beginning of sec. 3 a perfectly conducting plasma must satisfy E + v × B = 0. (7.1) As a result, particles with charge q, corotating with the pulsar have velocity v and experience no Lorentz force Furthermore, one sees immediately that E · B = 0 so that there is no component of force in the direction of B, and particles slide freely along magnetic field lines. The electric field satisfying (7.1) has a corresponding charge density given by Gauss' law [57] Notice that we have not specified anything about the form of the magnetic field, so that eq. (7.3) remains entirely general. If one specialises to a magnetic dipole, one obtains the classic Goldreich Julian density [63]. The interpretation is that this distribution corresponds to charges arranging themselves in the plasma in such a way as to exactly screen the component E · B/B of the induced electric field arising from the time-dependent B field. If it were true that everywhere in the pulsar magnetosphere E · B vanishes, two things would follow. Firstly, there would be no acceleration of particles along magnetic field lines and secondly, the source term for the axion field (3.4) would vanish. However, in order to explain observed emission spectra of pulsars, one must have regions in which particles do experience such acceleration, producing radiation in the process. Regions in which E · B is non-vanishing are known variously as acceleration regions, or acceleration gaps, with the term "gap" arising from the fact the plasma density may deviate from (7.3).
Understanding the exact location of acceleration gaps, the size of E · B within them, and the acceleration of charges in these regions has been a major effort in pulsar magnetosphere modelling in recent decades. One approach is to relax the constraint of perfect conductivity embodied in the condition (7.1). This can be achieved by incorporating a finite conductivity σ into the electrodynamic equations. This is known as resistive pulsar magnetosphere modelling, see e.g. refs. [64,65]. Such a setup allows one to formulate a semi-analytic approach (the equations of motion must still be solved numerically). Although, resistive pulsar magnetosphere models are able to reproduce observed properties of pulsar emission spectra [66,67], their drawback is that the exact form of the conductivity has essentially to be added by hand.
It therefore becomes necessary to treat the pulsar magnetosphere and its constituent plasma in a fully self-consistent way, by solving the full set of electrodynamic equations for a general current and EM fields. In addition one should allow for the injection of particles from the stellar surface and pair creation in the magnetosphere. Recent success in this direction utilises the "particle-in-cell" approach [68][69][70][71][72][73][74][75] which enables one to determine in a more precise way the location of the acceleration gaps and regions in which E · B = 0.
In light of this, our paper presents some of the key methodology needed to model axion clouds around pulsars and identifies some of the main qualitative features, e.g. the radiative and bound state solutions, and their scaling with axion mass and pulsar frequency. Generally speaking, we would expect the size and spatial extent of E · B (and therefore our effect) to be strongest for pulsars with low plasma densities, but higher magnetic fields, such as young pulsars -see figures 9 and 13 and table 2 of ref. [73].
However, in moving towards a more realistic pulsar model and its corresponding axion clouds, one should precisely identify the acceleration gaps with E · B = 0. One possibility, as a first step, would be to produce the E · B profile in one of the resistive pulsar magnetosphere models, decompose into spherical harmonics as in eq. (2.5) and extract the corresponding axion solutions. One could then (at least in the linear regime) determine the axion profile on a mode-by-mode basis. Such an approach has the advantage that it is a straightforward generalisation of the present calculation and also allows for a simple semi-analytic treatment of the electrodynamics with finite conductivity. Eventually this could in principle be extended to a full particle-in-cell simulation, though we shall leave such considerations for future work.

Summary and discussion and possible observable effects
In this paper we examined the size of axion solutions sourced by the axion-photon coupling in the electromagnetic fields of neutron stars, using a simple vacuum dipole model. To determine the axion profile, we solved the inhomogeneous wave equation with a source term given by E · B. We found that in the extremal case, for the largest known magnetic fields, achieved by magnetars, B * = 10 15 G and an axion-photon coupling of g γφ ∼ 10 −13 GeV −1 , that the total mass of the cloud can be of order 10 −13 M , and density reaching ρ ∼ 10 −1 g cm −3 , which compares favourably to certain classes of axion stars. We computed the total power output in the radiative case and speculated on how this might be used to constrain axion physics if sufficiently accurate measurements of pulsar spin down can be made.
This paper therefore makes the first important steps towards mapping out the structure of axion clouds around pulsars sourced by the axion coupling to electromagnetism. As alluded to in sec. 7 the next step is to incorporate magnetospheric screening of E · B to understand more precisely in what regions around the pulsar, and for what classes of pulsar E · B, is most sizeable. This will allow a better quantitative understanding of the effects presented here allowing for a more accurate discussion of observational signatures. There are a number of possible effects which should be investigated which we mention below.
• Polarisation and birefringence. Since the axion-photon coupling violates parity, photons passing through a medium of axions experience different dispersions relations depending on whether they are left or right polarised. Thus axion clouds behave as a birefringent medium. This has already been examined for light passing through axion clusters [36] and also through a superradiant cloud around a black hole [76]. There are a number of possibilities which should be pursued in the present case. Firstly, it is in precisely the "acceleration regions" where E · B = 0 that pulsar emission spectra are generated due to accelerated charges. Therefore the presence of non-vanishing axion configurations in these regions may have polarisation effects on pulsar emissions. Similarly, in the same vein as [76], one could look at light from a known linearly polarised source passing close by the pulsar and experiencing birefringence. One final possibility is that in the Ω > M regime where outgoing axion wave fronts are produced, one could examine how photons emitted from the pulsar "ride" the outgoing axion waves, experiencing polarisation effects on their journey from the pulsar. This would be an integrated effect over the light's journey from the pulsar to a detector.
• Luminosity and spontaneous decay of axions. It was noted some time ago that axion clouds can glow via the decays of axions into two photons with frequency M φ /2. This was exploited in the case of axion clusters [77,78] and more recently for superradiant axion clouds around black holes [24]. However, one would need an axion cloud of significant mass to produce an observable luminosity and also a high axion mass to give a large decay rate which goes as Γ φ→γγ = M 3 φ g 2 γφ /64π. In our case, large axion masses give highly suppressed cloud sizes, whilst smaller masses give a tiny decay width, thus it seems unlikely that the luminosity will be significant. One also would need to understand from a theoretical point of view, how dissipation into photons can be incorporated into the field equations as quantum corrections. Presumably one would expect the imaginary part of the one-loop photon correction to the axion propagator to give rise to a dissipative term as in e.g. [52].
• Stimulated Decays. One major advantage of studying axion configurations in the vicinity of pulsars is that they possess tremendously strong electromagnetic fields. This could have interesting consequences in terms of stimulated axion to photon conversions in the presence of the classical electromagnetic fields. The signatures arising from such processes in the vicinity of neutron stars have been examined in [79][80][81][82]. This exploits the same physics as light shining through a wall and haloscopes experiments [6,15].
• Generation of secondary electromagnetic fields. The axion solution can backreact onto electromagnetic fields via their equations of motion ∂ µ F µν = −g γφ ∂ µ φF µν . Thus the primary magnetic fields of the neutron star and the axion cloud form an effective current j µ eff ≡ −∂ ν φF νµ which can source secondary electromagnetic fields. For instance, in the present VDM model for the pulsar, both φ and A µ contain a Y 11 harmonic, so the effective charge ρ eff = −g γφ (B · ∇φ) will contain a corotating quadrupole moment Y 22 (θ, λ), with λ = ϕ − Ωt. This could give rise to a secondary sub-pulse from the pulsar in between the main beam. Firstly though, one would need a better pulsar model incorporating magnetospheric effects and, from an observational point of view, one should also take into account any higher multipole moments of neutron stars [83] which may obscure any signal coming from axion back-reaction.
Further effects which should be investigated are bosanova phenomena for high axion densities when non-linear self-interactions become important [40,41]. This can lead to supernova style explosions and bouncing solutions. Another interesting possibility is to study how our clouds modify the forces between neutron stars as in [55] and how our discussion might modify the nature of the solutions present therein.

Appendices A Green Functions
The Green functions for eq. (4.3) can be expanded in terms of the two linearly independent mode functions f (r) and g(r) which we choose to have boundary conditions given by eq. (4.5). The Green function therefore takes the form for some functions a,b,c and d depending on r . The Green function must also satisfy the usual continuity and jump conditions given by The source term and the Wronskian decay only as an inverse power of r, however g(r) grows exponentially. It follows that when integrating over the region r > r, the g(r ) contribution gives a divergent integral. We therefore must cancel this behaviour from g by taking a(r ) = γ(r ) + g(r ) A(r ) 2 W (r ) , where γ(r ) is some function decaying suitably fast so as to make all integrals over r convergent. The Green function therefore takes the form G(r , r ) = γ(r )f (r) + g(r )f (r) A(r ) 2 W (r ) r > r γ(r )f (r) + g(r)f (r ) A(r ) 2 W (r ) r < r . (A.8) Next we recall that the differential operator appearing on the left hand side of (4.3), is of the Sturm-Liouville type, so that it is self-adjoint with respect to following inner product with weight A(r) −1 , given by ψ 1 , ψ 2 = dr A(r) ψ 1 (r)ψ 2 (r). (A. 10) In other words, ψ 1 , Lψ 2 = Lψ 1 , ψ 2 . It therefore follows that since G(r, r ) is the inverse of a self-adjoint operator that G is also self-adjoint, since Gψ 1 , ψ 2 = Gψ 1 , LGψ 2 = ψ 1 , Gψ 2 . Since this holds for any ψ 1 and ψ 2 , it therefore follows that G(r, r )/A(r) is symmetric with respect to r ↔ r . By this symmetry, it follows that γ(r ) = γf (r )/A(r) for some constant γ. Furthermore, since W (r)A(r) = 2k, one sees that G(r, r ) now takes the form G(r , r ) = γ A(r ) f (r)f (r ) + 1 2kA(r ) θ(r − r )g(r )f (r) + θ(r − r)g(r)f (r ) , (A.11) which manifestly has the desired symmetry properties, where the constant γ is fixed by matching the exterior and interior solutions.

B Interior Solution
The interior wave equation (4.7) can be recast as r 2 u + rP (r)u + Q(r)u = 0, (B.1) where P and Q are given by eqs. (4.10). Notice that P and Q are analytic about r = 0. It therefore follows that r = 0 is a regular singular point since The general solution in the stellar interior can then be expressed as a linear combination of solutionsψ ± corresponding to these two roots. Note that for the specific case = 1 we consider here, the roots are integers with s + = 2 and s − = −1 so that the full solution takes the form ψ int (r) = r 2ψ + (r) + r −1ψ − (r). (B.7) Expressing eachψ as a power seriesψ A short calculation shows that P (0) = Q (0) = 0 and soψ (0) = 0. If we impose the further constraint that u is regular at the origin, the s = −1 contribution ψ ∼ r −1ψ − , implies that ψ − (0) = 0, which means thatψ − therefore vanishes everywhere in the stellar interior. This means that the interior solution must take the form ψ int (r) = r 2ψ + (r), 0 ≤ r ≤ R (B.11) whereψ + (r) satisfies eq. (B.6) with s = 2 subject to the boundary conditionψ + (0) = 0.