Dipole Radiation and Beyond from Axion Stars in Electromagnetic Fields

We investigate the production of photons from coherently oscillating, spatially localized clumps of axionic fields (oscillons and axion stars) in the presence of external electromagnetic fields. We delineate different qualitative behaviour of the photon luminosity in terms of an effective dimensionless coupling parameter constructed out of the axion-photon coupling, and field amplitude, oscillation frequency and radius of the axion star. For small values of this dimensionless coupling, we provide a general analytic formula for the dipole radiation field and the photon luminosity per solid angle, including a strong dependence on the radius of the configuration. For moderate to large coupling, we report on a non-monotonic behavior of the luminosity with the coupling strength in the presence of external magnetic fields. After an initial rise in luminosity with the coupling strength, we see a suppression (by an order of magnitude or more compared to the dipole radiation approximation) at moderately large coupling. At sufficiently large coupling, we find a transition to a regime of exponential growth of the luminosity due to parametric resonance. We carry out 3+1 dimensional lattice simulations of axion electrodynamics, at small and large coupling, including non-perturbative effects of parametric resonance as well as backreaction effects when necessary. We also discuss medium (plasma) effects that lead to resonant axion to photon conversion, relevance of the coherence of the soliton, and implications of our results in astrophysical and cosmological settings.

Since the g aγ φE · B interaction between axions and electromagnetism is expected to be very weak, one might seek to compensate the tiny coupling g aγ by searching for signatures in systems with a strong electromagnetic field and/or a large axion field amplitude. Moreover, if φ is oscillating it can also have an enhanced effect through resonances [39][40][41]. With these considerations, it is natural to explore the impact of axion stars in strong electromagnetic fields. Axion stars are large amplitude, spatially localized and oscillating φ field configurations (also known as scalar solitons, oscillons, axitons etc. [42][43][44][45][46][47][48][49][50][51][52][53][54][55]). Such exploration is the main purpose of this paper. We aim to understand how the interplay between the axion-photon coupling strength, parameters defining solitons, and electromagnetic fields influences the production of electromagnetic radiation from such solitons.
The study of axions in astrophysical and cosmological electromagnetic fields has a long history [56,57]. Some of the strongest constraints on the coupling of axions to matter come from considering the production of axions in the hot and dense stellar interiors. Alternatively, if axions were produced in the early universe and survive today as dark matter, then the flux of these cold axions onto magnetized compact stars could result in a distinctive radio emission [58][59][60][61][62][63][64][65]. As much as an O(1) fraction of the axion dark matter could be in the form of axion stars, and therefore it is also important to develop strategies for detecting the encounter of axion stars with magnetized compact stars [66][67][68][69][70][71][72][73][74][75][76][77][78]. Furthermore, the collision of and collapse of axion stars can amplify even small fluctuations in the electromagnetic fields [41,53,66,79].
In the work being presented here, we consider the coherent emission of electromagnetic radiation from an axion star in an electromagnetic field. Because of the high occupation number of the axions in the solitons, it is natural to treat the axion field classically. We calculate the spectrum and luminosity of the resultant electromagnetic radiation using both analytical techniques and 3+1 dimensional lattice simulations. Based on Floquet analysis, we argue that different qualitative behaviour of the electromagnetic radiation is determined by a dimensionless effective coupling parameter C ∼ (g aγ ϕ 0 )(ωR), where ϕ 0 is the field amplitude, ω and R are the oscillation frequency and radius of the axion star respectively.
In the small effective coupling C 1 regime, the analytical analysis is based on the observation that an axion star in an external electromagnetic field develops a charge and current dipole, which oscillates in time and produces dipole radiation. In the absence of resonance effects, we find that the signal has a strong dependence on the axion star's radius and oscillation frequency, which leads to a suppression (that goes like the Fourier transform of the spatial profile) at large ωR > O(1). This understanding leads us to focus our attention on compact axion-star configurations and oscillons. Our dipole approximation is validated with numerical simulations of the axion electrodynamics on a 3 + 1 dimensional lattice.
Floquet analysis and lattice simulations also allow us to study the regime of moderate to large C ∼ 1, where perturbative analytical results are difficult to obtain. We are able to capture a non-trivial transition from a steady photon production rate to an explosive (exponential) one as we vary the coupling strength and axion field configuration. We find qualitative and quantitative differences between the photon production rate in the presence of external electric and magnetic fields (including significant suppression of the radiated power at moderate coupling). We also analyse the backreaction of photons on the axion configuration when necessary.
Most earlier work on axion stars in astrophysical magnetic fields relies on a 'resonant' axion-to-photon conversion, when the plasma frequency approximately matches the energy of the axion particles (see, for example, [7]). While our simulations do not include the effects of a plasma, we are able to incorporate this resonant conversion in our calculation analytically in the small coupling regime. We also comment on the relevance of a coherent solitonic configuration compared to an incoherent collection of dipoles, as well as the connection of our results to the well-known quantum mechanical calculation related to the axion-photon conversion probability (see, for example, [56]). The remainder of this article is organised as follows: In Sec. 2 we briefly introduce the model of interest, namely axion electrodynamics, and in Sec. 3 we introduce axion stars. In Secs. 4 and 5 we employ analytical and numerical techniques to calculate the spectrum of electromagnetic radiation that arises from an axion star in an external electromagnetic field. In Sec. 6 we comment on a few supplemental topics such as finite density and coherence effects, and in Sec. 7 we discuss several possible observational signatures. Finally, we summarize and conclude in Sec. 8. We include an Appendix A with details of the dipole radiation calculation.

Axion electrodynamics
Our system consists of a real valued, pseudo-scalar field φ coupled to the electromagnetic field. The action for our system is given by where we adopt the − + ++ signature of the metric. The electromagnetic field-strength tensor, and its dual are: where 0123 = 1. The equations of motion for the axion and the gauge fields are given by Note that ∂ µ j µ = 0, and that we have assumed that there are no free currents or charges in our system. The above four-current arises from axion-electromagnetic interactions. We define electric and magnetic fields in the usual way with ijk = ijk . The coupled Klein-Gordon and Maxwell equations are then given by [80] φ Note that the effective charge and current densities are In the above equations, we have ignored gravitational interactions. If one wishes to include weak-field gravity (gravitational potential |Ψ| 1), the substitution ∂ φ V → (1 + 2Ψ)∂ φ V in the equation of motion for the scalar field captures the most relevant gravitational contributions. Moreover, we would need to include a Poisson equation ∇ 2 Ψ = (1/2m 2 pl )ρ φ where ρ φ is the density of the axion field to close the system. This prescription allows certain gravity-supported scalar field configurations to exist, but ignores gravitational effects (such as redshifts) in the dynamics of electromagnetic fields and also ignores the contribution of electromagnetic fields in determining the gravitational potential. 1

Compact axion stars in constant electromagnetic fields
We are interested in electromagnetic radiation generated by a spatially localized, spherically symmetric, coherently oscillating axion field configuration of the approximate form φ(t, r) ≈ ϕ(r) cos(ωt) . (3.1) Such solutions of the nonlinear Klein-Gordon equation (with and without gravity), which we generically refer to as solitons, are a result of a balance between the tendency of the field configurations to disperse and (i) attractive self-interactions in the potential V (φ) and/or (ii) gravitational interactions. The detailed form of ϕ(r) depends on the potential V (φ) as well as ω. For most of our purposes, we use an ansatz of the form ϕ(r) = ϕ 0 sech (r/R) so that φ(t, r) = ϕ 0 sech (r/R) cos ωt . (3.2) The above form is motivated by the fact that it has the correct large distance behavior: Note that there is also polynomial dependence of the profile on r at large radii multiplying the exponential, which are ignoring here [81,82]. Typically, ω is not too different from m, however, ϕ 0 and R can vary significantly for small changes in ω close to m. In a typical scenario, ϕ 0 , R and ω are not independent. Usually we are free to chose only one, and even that has constraints from stability analyses.
To understand what to expect for ϕ 0 and R, we consider two relevant cases below. Fig. 1), exceptionally long lived spatially localized configurations of the above form exist, and are called oscillons (setting g aγ → 0 for the moment). Typically, for very long-lived oscillons, we have a field amplitude ϕ 0 ∼ f , a spatial width R ∼ few×m −1 and an oscillation frequency ω m [83]. In detail, there is a one parameter family of long lived configurations for a given potential V (φ). Moreover, typically the solution also includes higher frequencies and a very small radiating tail (scalar radiation [84]).

Self-interaction supported solitons
The scalar field potentials that support solitons. For the quadratic potential and cosine potential, gravity is essential for supporting long-lived solitons, whereas the "flattened" potentials can support solitons without gravity, but typically require amplitudes ∼ f . For any potential where solitons have a small amplitude compared to f , gravity is essential for long-lived stable solitons. Right: A schematic representation of a solitons. Dilute solitons have ϕ 0 f and R m −1 . Dense solitons have ϕ 0 ∼ f and R ∼ few × m −1 . The frequency is always ω m.
Because of the scalar radiation, oscillons are not perfectly stable, and they exhaust their energy on a time scale τ . This lifetime depends sensitively on the scalar potential V (φ). For example, a cosine potential leads to a lifetime τ ∼ 10 3 m −1 , whereas some other potentials shown in Fig. 1 give much longer lifetimes τ 10 12 m −1 [83,85,86]. The formation of such oscillons from cosmological initial conditions (especially in the early universe) has been explored in detail before [87][88][89][90][91], and typically happens when H ∼ m. An almost homogeneous, oscillating condensate naturally fragments into oscillons. As a result, compared to the H −1 at the time of formation, oscillons can be exceptionally long lived, and have important cosmological implications [85,[92][93][94][95]. Using similar arguments, oscillons in ultralight axions might potentially survive until today [85,96,97]. Furthermore, oscillons appear to be attractors in the space of solutions [98], and might also nucleate inside dark matter halos [79], or even near black-holes [99], although most of these analyses are done in the context of gravitationally supported solitons so far.
If we are interested in a population of oscillons in the contemporary universe which have a primordial origin, their lifetimes will likely be short compared to H −1 0 ∼ 10 33 eV −1 (unless m 10 −21 eV). While challenging, claims exist in the literature for oscillons that have lifetimes comparable to the present age of the universe [100,101]. As discussed earlier, late universe formation mechanisms can also ameliorate this problem (also see [33]). The detailed investigation of oscillon production and population is not considered in this paper; we simply take these objects to exist and examine their consequences due to encounters with external electromagnetic fields.
Less compact configurations can also exist if the amplitude of the field is not so large, and R m −1 -they are referred to as dilute axion stars [52,55]. In this regime the central field amplitude ϕ 0 ∝ 1/R 2 , the frequency ω ≈ m, and the radius R ∼ 1/ √ m 2 − ω 2 . Such dilute configurations are well described by a Schrödinger-Poisson system (with a conserved particle number), and are often the focus in fuzzy dark matter studies [85,108,109]. Dilute axion stars have the benefit of being cosmologically long-lived. However, as we will see, electromagnetic radiation from dilute axion stars in external electric and magnetic fields is heavily suppressed. As a result, most of our focus will be on the dense, smaller radius solitons. 2

Analytic calculation of electromagnetic radiation
In this section, we calculate the electromagnetic radiation generated by spatially localized, coherently oscillating axion configurations (solitons) discussed in the previous section. In the presence of external electromagnetic fields, such configurations can be effectively thought of as time-dependent charge densities and currents that produce electromagnetic radiation. We provide analytic results for the produced radiation at leading order in the coupling g aγ , and discuss difficulties with going beyond the leading order analytically. We also discuss the expected non-perturbative (in the coupling) results in general terms.
The first-order Maxwell equations (2.6) can be rearranged into the following differential equations:Ë The 4-current (ρ, J ) defined in (2.7) is spatially localized because the axion field configuration φ given by eq. (3.1) is spatially localized. Note that (ρ, J ) depend on φ as well as the E and B via eq. (2.7). Beyond the spatial extent of the axion stars, both E and B propagate like free waves.

Floquet analysis
Because the system is linear in E and B fields, and we assume φ to be periodic in time, we expect the solutions to obey Floquet's Theorem [111,112]. That is, the solutions are either bounded and periodic, or have exponential growth in time. However calculating Floquet exponents (µ), or explicit solutions is a tall order because of the large number of coupled degrees of freedom associated with each spatial point (formally infinite, and usually a rather large number in discretized three dimensions). Equivalently, the modes in Fourier space are coupled because of the spatial variation in φ. 3 While the explicit calculation of the Floquet exponents is non-trivial, we can get a physical understanding of their scaling with parameters and the parametric boundary between bounded and unbounded solutions as follows. For the homogeneous axion field with amplitude ϕ 0 and oscillating harmonically with a frequency ω, the electromagnetic fields are always unstable, with the k ≈ ω/2 electromagnetic field modes growing as e µ hom t where µ hom ≈ g aγ ϕ 0 ω/4 at least when g aγ ϕ 0 is not too large [113] (for larger amplitudes, it is model dependent [41]). In contrast, for the localized soliton configuration, we expect a threshold value of the coupling g aγ ϕ 0 for which we get exponentially growing solutions. The parameter ϕ 0 should now also be thought of as the central amplitude of the soliton. The threshold can be determined by comparing µ −1 hom to the width of the soliton R [41,113,114]. Essentially, if the produced photons can escape the system quickly enough (ie. R is small enough), they do not lead to exponential growth due to parametric resonance (equivalently, Bose-enhancement). This motivates the definition of a dimensionless effective coupling In terms of this effective coupling: C 1 −→ bounded periodic solutions, steady radiated power , C 1 −→ unbounded exponential solutions and radiated power . We remind the reader that C is independent of background electromagnetic fields. Note that for C > C crit ∼ 1, the power in radiated electromagnetic fields In Sec. 5, we will confirm this behaviour, and provide the numerical coefficient in front of this expression for µ eff based on a specific soliton profile. We remind the reader that soliton configurations do not allow us to specify ϕ 0 , ω and R independently. For example, dilute and gravitationally supported solitons have ϕ 0 ∝ R −2 . For dense, self-interaction supported axion stars/oscillons, ϕ 0 ∼ f . For the dilute case, we have ωR 1, so we can get C ∼ 1 for g aγ ϕ 0 1. For the dense case, we typically have R ∼ few × m −1 , so we can get C ∼ 1 with g aγ ϕ 0 ∼ 1. The C 1 can be achieved, for example, by simply making g aγ smaller in each case.
Before moving on to a quantitative analytical analysis, we briefly discuss the connection of C 1 and C 1 regimes with effective field theory (EFT) considerations. The action in Eq. (2.1) represents the leading operators in an EFT with cutoff Λ ∼ g −1 aγ describing axionphoton interactions. 4 The EFT also contains sub-leading operators that are suppressed by additional powers of the cutoff, e.g. L sub ⊃ c sub g 2 aγ φ 2 F 2 or c sub g 3 aγ φFF . Validity of the EFT requires the sub-leading operators to be negligible. As discussed above, it is possible to have g aγ ϕ 0 1 to get C 1. For dilute axion stars, C ∼ 1 can be obtained for g aγ ϕ 0 1 also. However, for C ∼ 1 in the dense case, we need g aγ ϕ 0 ∼ 1, which threatens to break the EFT if higher-order operators are only suppressed by additional powers of g aγ ϕ 0 . Even in this case, the EFT can remain reliable even for g aγ ϕ 0 ∼ 1 if the numerical coefficient of the higher-order operators is small, e.g. c sub 1. For some theoretical work on models with a large axion-photon coupling, see [115][116][117][118][119][120][121].
BĒ Figure 2. The effective charge and current density (dipoles) induced by the presence of a soliton in an external electromagnetic field background. The left image shows a charge dipole aligned with the external magnetic field, and the right image shows a current dipole in a plane normal to the external electric field. The charge density and current density oscillate in time, generating dipole radiation.

Perturbative analysis
With the expectation of bounded solutions for C 1, we pursue an analytic treatment in the limit of small g aγ ϕ 0 . With this small parameter in mind, we expand the fields, densities and currents as follows: Here we use the subscript (n) to denote the terms containing n-th power of g aγ ϕ 0 . At the lowest order, the E (0) and B (0) stand for the electric and magnetic backgrounds and are sourced by (ρ (0) , J (0) ) which are independent of the axion field configuration. For example such background fields could be the fields in the magnetosphere of a neutron star or in the intergalactic medium. To make the physics more transparent, we will consider spatio-temporally constant background electromagnetic fields which we denote by We are essentially assuming that the spatial extent of the axion star is much smaller than the coherence length of the background fields, and that the time variation of the background fields is slow compared to the time that configuration spends in the given volume of the fields.

Leading order in g aγ ϕ 0 : dipole radiation
At leading order in the coupling g aγ ϕ 0 , we havë At this order in g aγ ϕ 0 , the background electromagnetic fields along with the axion configuration φ(t, x) = ϕ(r) cos ωt induce an effective charge and current density: with (1)  Due to the spatial derivative acting on ϕ along the direction ofB field, the positive and the negative charges are distributed separately along theB field axis like a dipole (see left panel in Fig. 2). And with its oscillating nature of the axion configuration, such an oscillating dipole will lead to dipolar electromagnetic radiation. A constantĒ field results in an oscillating azimuthal current, which also results in dipolar radiation (see right panel in Fig. 2). It is a standard textbook problem to compute the excited electric and magnetic fields caused by the harmonic, spatially localized sources of the form (4.10), as well as the associated Poynting flux S (2) ≡ E (1) × B (1) and power emitted per unit solid angle. See for example [122,123]. We review some of the relevant details of the derivation in Appendix A. Here, we directly write down the solution for the flux below. At a position x far from the source, and at sufficiently late times, the power per unit solid angle dP γ (2) /dΩ = |x| 2x · S (2) , is given by wheref (k) is the spatial Fourier transform of f (x). Using the specific forms of the charge and current densities in (4.11), we have˜ (1) The radiation spectrum is a delta function in frequency, with the radiation emitted at the frequency ω. The spatial pattern of radiation energy density is shown in Fig. 3. The total power emitted, and its time average are given by (4.14) It is important to note that the emitted power is proportional to the squared Fourier transformφ(ω) of the oscillon radial profile evaluated at ω (which is the frequency of the oscillon, and that of the emitted electromagnetic radiation): for a spherically-symmetric oscillon. Then the ratio F (ω) =φ 2 (ω)/φ 2 (0) is a form factor for the oscillon profile. If the wavelength of the radiation is large compared to the scale radius of the oscillon, R λ = π/ω, then the form factor approaches F (ω) ≈ 1 as ω → 0, corresponding to radiation from a point-like dipole. Shorter wavelength radiation probes the structure of the oscillon and F (ω) → 0 as ω → ∞. This behavior is illustrated in Fig. 4 for a few representative oscillon profile functions.
Using the profile ϕ(r) = ϕ 0 sech (r/R), we get where the second equality assumes ωR 2. When the radius of the axion configuration R ∼ ω −1 , there is no suppression of the emitted power fromφ(ω). However, when R ω −1 we get an exponential suppression. We have checked that the exponential suppression also exists for numerically obtained spatial profiles for dilute axion stars where ωR 1. The physical origin of this suppression is destructive interference between the emitted electromagnetic waves which are emitted in phase from different locations within the oscillon. Also see discussion of coherence and interference in section 6.2. Note that this suppression is more severe than the suggested by [72], where a power law suppression is obtained because of where ω is the frequency of oscillation of the axion field (and of the emitted electromagnetic radiation), andφ(ω) is the Fourier transform of the soliton's spatial profile at k = ω. This form factor determines the amplitude of the dipole radiation, and has a very strong dependence on the radius of the soliton. The blue curve corresponds to the sech profile ∝ sech (r/R), with the correct exponential behaviour at large R. The solid gray curve is for a Gaussian profile ∝ e −r 2 /R 2 , the dashed one for an exponential profile ∝ e −r/R with a cusp at the origin, and the dotted line correspond to a top-hat profile of the soliton with radius R. While the form factor is identical at small ωR for the different profiles, it is very sensitive to the profile choice at large ωR.
to a cusp in their ϕ(r) at the origin. This can make a rather large difference in the radiated power even for ωR few. Compare the blue curve for the sech profile with the dashed gray curve for an exponential profile with a cusp at the origin. 5 Also see Sec. 3 for further discussion of the expected form of the axion star profiles.

Summary of dipole radiation
Finally, to make the dipole nature of the radiation apparent, let us set the background electric field to zero. In this case where θ is the angle with respect to theB direction. The same formula holds for the electric field also. The second equality is a good approximation for ωR 2 for the sech profile. To get significant emitted power, it is essential to have Rω not be too large, and g aγ ϕ 0 not be too small, which provides motivation for considering dense axion stars and oscillons. At the same time, it is also beneficial to have a small ω ∼ m a which pushes us towards pursuing lighter axions.

Higher orders in g aγ ϕ 0 : beyond dipole radiation
Our organization of the calculation using powers of g aγ ϕ 0 is fraught with subtleties as we go beyond the leading order in g aγ ϕ 0 , with the system best dealt with non-perturbatively using Floquet theory (with a large number of coupled degrees of freedom). However, to appreciate these subtleties, we try to follow our nose and proceed with the calculation order by order in g aγ ϕ 0 . While we will be unable to complete the calculation, the set up also provides some physical insight into how the radiated power deviates from the dipole estimate of the previous section as we increase the coupling strength.
The field equations, charge and current densities are given bÿ for n ≥ 1. Recall that (n) denotes the order in g aγ ϕ 0 , E (0) =Ē and B (0) =B are assumed to be constants, and φ = ϕ(r) cos ωt. We will continue ignoring backreaction of the produced electromagnetic fields on the axion field configuration in this subsection. 6 Order by order this represents a system of periodically forced oscillators. We can use these to understand the possible frequency structure of the fields at different orders (ignoring resonances for the moment). Since φ oscillates with a frequency ω, so do (ρ (1) , J (1) ), which in turn source (E (1) , B (1) ) which also oscillate with a frequency ω. However, because of the products of oscillating terms coming from φ and oscillating electromagnetic fields, (ρ (2) , J (2) ) will include frequency components 0ω and 2ω. Similarly, (ρ (3) , J (3) ) will contain ω and 3ω and so on.
The above arguments reveal that if we are interested in the radiated electromagnetic fields at O[(g aγ ϕ 0 ) n ], they will contain multiple frequencies. Conversely, if we want to consider fields with a fixed frequency, they will contain terms with many different orders in g aγ ϕ 0 . This latter fact does mean that there is a possibility that the generated electromagnetic fields (and power) at a given frequency, or in total, can be enhanced or decreased as we go to higher couplings. That is, the power radiated can be non-monotonic in the coupling (when the coupling is not too small), and its dominant frequency content might also change with coupling strength. We observe these effects in our numerical simulations. Finally, note that fields with higher frequencies beyond ω always come with higher powers in g aγ ϕ 0 ; this is because higher frequencies are sourced by the electromagnetic fields already sourced by the axion field configuration. Again, we confirm this behavior in the simulations.
The above discussion is incomplete because we ignored the possibility of resonances that should be present in a system with periodic forcing terms. These resonances make it notoriously difficult to carry out our perturbative scheme for long time scales. We can get a rough idea of the difficulties and subtleties by trying to solve the (4.18) equations in momentum space. First, let us consider the n = 1 case: As we see, the result is periodic and bounded, except when p equals ω. This is just the behavior of a periodically forced harmonic oscillator. As expected, for p = ω, there is a term that is linear in t (secular growth). After sufficient time, such a term will dwarf the zero order terms. This in turn can limit the reliability of the perturbative expansion we used in the first place. To maintain the validity of the perturbative expansion, besides requiring a small g aγ ϕ 0 , one should further restrict ourselves to small times. This scenario with secular terms is reminiscent of the challenge of solving the Mathieu equation via perturbative methods, and in principle, there exist mathematical tools to deal with such situations. For example, one can go beyond the naive perturbation theory (4.5), and resort to the Renormalization Group [124] or resurgent resummation [125]. But the question here is more complicated than in the Mathieu equation because of the vast number of coupled momentum degrees of freedom. Furthermore, there is another subtlety. While the individual mode for B (1) (p = ω) has a secular term, when we obtain the fields in position space via a Fourier transform, the secular term disappears. Note that the secular term above can be reached by the expression of p = ω, in the limit p → ω, and therefore this imposes no pole or singular point in the momentum integral.
Nevertheless, there are good reasons to believe the secular terms will appear at high order in g aγ ϕ 0 terms for the fields. One reason is that Floquet theory predicts (and we observe in simulations) the existence of exponentially growing solutions that can be constructed out of solutions with different power law (unbounded) in time dependencies at various orders. More generally, such terms can combine in non-trivial ways to give real and imaginary Floquet exponents corresponding to bounded and unbounded solutions.

Results of numerical lattice simulation
The axion-photon system can be simulated numerically [41,[126][127][128][129]. In general, when there exist charged matter fields, the usage of the electromagnetic potential A µ is unavoidable, as the gauge covariant derivatives require A µ explicitly. But in the simple axion-photon system where there is no charged field, the electromagnetic scalar potential φ and vector potential A are not necessary. Instead, we can directly evolve the electric field E and the magnetic field B in the axion background through the Maxwell's equations, and as a byproduct, there is no gauge fixing needed. We will present the details of our numerical scheme in a separate paper.
Simulation Parameters and Initial Conditions: Our benchmark simulation has a physical volume m 3 V = 64×64×64, with N 3 = 160 3 lattice sites, and the resolution is mdx = 0.4. We also used N 3 = 320 3 for convergence tests (for more details, see [41]). We typically run our simulations up to mt max = 50. We also used mt max = 100 when using an eight times larger simulation volume. We employ periodic boundary conditions but make sure that we do not have our results contaminated by radiation cycling through the box.
For initial conditions we start with the axion field in a solitonic configuration of the form φ(t) = ϕ 0 sech (r/R) cos(ωt + θ 0 ) with θ 0 = −π/2. That is, φ(t = 0) = 0 andφ(t = 0) = ωϕ 0 sech (r/R). We include either a constant magnetic field or a constant electric field through the box. For our fiducial values of our parameters, we use ϕ 0 = 2.6f , ω = 0.82m , R = 1.6m −1 , andĒ = 10 2 m 2 orB = 10 2 m 2 . (5.1) With these above values (see Eq. (4.2)): 2) The soliton parameters are consistent with those of dense solitons found in [83], although the precise values can differ based on the functional form used to fit the true profile. Apart from transients, our results are insensitive to the chosen values ofĒ andB apart from a trivial scaling of the radiated power in the C 1 regime, and a change in logarithmic time scale of backreaction (see below) in the C 1 regime. 7 The chosen values of the background fields are for numerical convenience. We also varied these parameters within factors of two or even an order of magnitude to delineate general statements from those that are sensitive to this particular choice of fiducial parameters.
While we focus on dense solitons, we could have carried out simulations in the dilute soliton regime as well. However, the exponential suppression expected from (4.14) would make this uninteresting (at least for dipole radiation, though not necessarily for the case with parametric resonance [113]).
Backreaction Considerations: For most of our simulations, it is unnecessary to evolve the axion field using its equation of motion numerically (although it is still necessary to solve for the electromagnetic fields numerically for C ∼ 1). That is, the soliton sources electromagnetic fields, but it is not significantly affected by them. To see this, recall that the energy extracted from a dense soliton with C 1 grows linearly with time with P γ ∼ 10B 2 (f g aγ ) 2 /m 2 (see (4.14)). Hence it will take mt br ∼ mM sol P γ ∼ 10 for backreaction on the soliton to be relevant. We have used M sol ∼ 10 2 f 2 /m above for the energy of a dense soliton. Note that since even for the strongest fields possible around neutron starsB m 2 e ∼ 10 −1 MeV 2 , and with g aγ 10 −10 GeV −1 and m 10 −7 eV, we have m/(g aγB ) 1. Hence, we can safely ignore backreaction on the axion field configuration in our simulations when C 1 and mt max 10 2 . Note that another way of thinking about backreaction is at the level of the equation of motion for the axion field (first equation in (2.6)). The correction to the axion field evolution due to the source term g aγ E · B is given δφ/ϕ 0 ∼ g 2 aγB 2 /m 2 1. This is essentially the same ratio that appears in the discussion above.
For C 1, the exponential growth in the radiated power can lead to backreaction on the soliton within mt max . At the end of this section, we provide simulation results where backreaction eventually shuts down the resonant electromagnetic field production.
Numerically Calculated Power: The main output from our simulations will be the radiated power in electromagnetic fields. We define this radiated power as the surface integral of the Poynting vector over a spherical surface whose radius is much larger than the size of our soliton. In our set-up, we compute the luminosity by a sum where the sum is over all sites of index j with a distance within (r − , r + ), given m −1 . Note that we excludeĒ orB when we compute the Poynting vector since these background fields are not part of the radiation that escapes to infinity. In our simulations, the radiated power is measured at mr = 16 with m = 0.1. Note that for mt max = 50, the radiation does not have sufficient time to cycle through our periodic box and contaminate the radiated power calculation.
In general, we found that the numerically calculated power was not very sensitive to the lattice size or the radius of the sphere where we calculated the radiated power as long as this radius R. Our finite lattice spacing, mdx, leads to a slightly smaller numerically evaluated power in comparison to the power calculated in the continuous limit (when such a calculation is possible). For mdx = 0.4, the discrepancy with the analytic expectation is ∼ 1% for C 1.

Bounded vs. unbounded radiating fields
As we discussed at the beginning of Section 4, we expect periodic solutions for C 1 and exponential growing ones for C 1 based on Floquet theory. To confirm this behaviour, we numerically solve for E and B fields sourced by the same soliton configuration (with our fiducial values of R, ϕ 0 and ω), but different f g aγ . /Ē 2 Figure 6. The radiated electromagnetic power from a dense soliton in constant background E or B field as function of time for C = 0.09. In this regime the radiation is expected to be described well by dipole radiation P γ (2) provided in (4.14). The black lines are the analytic results for the time-averaged power which matches nicely with the numerical results. Note the frequency of 2ω is expected from analytics as well.
In Fig. 5 (left panel) we plot the power radiated as a function of time for different g aγ . Note that the power radiated is constant for C 1 but increases exponentially as P γ ∝ e 2µ eff t for C 1 as expected. Also note the different scales on the vertical axes for different parts of the panel. In the right panel of Fig. 5, we plot the maximum Floquet exponent µ eff from the numerically obtained time dependence of the radiated power. This plot reveals that Note that the Floquet exponent is a property of the axion configuration and coupling g aγ (though the combination in C), but is independent of the presence or absence of background electromagnetic fields. While we expect C crit ∼ 1, its precise value and the numerical coefficient appearing in the expression for µ eff will depend on the details of the soliton solution.

Small coupling: dipole estimate
In Section 4.2.1, we provided an analytic calculation for the power radiated by the soliton configuration in the presence of background E and/or B fields. This result is expected to hold for C 1. We confirm this expectation in detail with numerical simulation for C ≈ 0.09. In Fig. 6, we present the time-dependent power radiated by the soliton in the presence of a constant background B field (left) and E field (right). The radiated power oscillates with a frequency 2ω, consistent with our analytic result in the first line of (4.14). Moreover, the magnitude of the time averaged power is also consistent with our analytic calculation in the second line of (4.14) (horizontal black lines). Finally, the spatial pattern of the radiated energy density is consistent with dipole radiation as predicted in eq. (4.17) (see Fig. 3). The dependence of the time-averaged radiated power P γ t on the effective dimensionless coupling C = g aγ ϕ 0 ωR/4. The above plot is based on a dense soliton of field amplitude ϕ 0 = 2.6f , frequency ω = 0.82m and radius R = 1.6m −1 . For changing C ≈ 0.85(f g aγ ), only f g aγ is varied. The black line is the dipole estimate P γ (2) t from eq. (4.14). The orange and green dots show the numerically evaluated P γ t for E and B field backgrounds respectively. Note that the numerics agree with the dipole estimate at C 1 as expected. The deviation becomes more and more pronounced as we move from C 1 towards C ∼ 1. In particular, note the difference in the radiated power between E and B field backgrounds. For the B background, note the significant suppression of the radiated power compared to the dipole estimate and the non-monotonic behaviour with C. Finally, as C > C crit = 1.3 (grey shaded), we have an exponentially growing (in time) power due to parametric resonance.
While we do not do so here, we can easily include both electric and magnetic field backgrounds together to confirm eq. (4.13).

Intermediate couplings
We now consider the power radiated when 0.1 C 1. As we saw in Section 4.2.2, it is difficult to extend the perturbative calculation to include terms that are higher order in g aγ ϕ 0 . Numerically, we of course have no issues probing this regime. In this regime, new phenomenon emerge, which have not been reported before to the best of our knowledge.
As we move away from the C 1 regime towards C = C crit , the radiated power is still constant in time (ie. we have periodic solutions for the radiated fields at each point in space). However, we find that the dipole estimate P γ (2) t from eq. (4.14) starts deviating significantly from our numerical results. See Fig. 7, where the black curve is P γ (2) t whereas the orange and green curves represent the numerically obtained P γ t for the case of background B and E fields respectively. In the following, when we vary C, we hold all parameters fixed apart from f g aγ .
Background B: For the case of the background B field, we find that the P γ (2) t overestimates the power in this regime (see Fig. 7). The frequency content of the radiated power  The photon particle number distribution in momentum space obtained by taking a Fourier transform of the electromagnetic fields in our simulation volume at a fixed time (for more details, see [41]). Three gray lines represent p = ω, 2ω and 3ω. Left panel is for C = 0.09 whereas the right panel is for C = 0.85. Note that at C = 0.85, the dominant radiation frequency in ω is absent, and is connected to the suppression of power in the magnetic field background at this coupling. continues to be dominated by 2ω. As discussed in Section 4.2.2, we believe that this is due to higher order contributions of O((g aγ ϕ 0 ) 4 ) to the radiated power at frequency 2ω. These next-to-leading order in g aγ ϕ 0 contributions have an opposite sign compared to the leading order result. This is confirmed by our fits to P γ t as a function of C(∝ g aγ ). As C ∝ g aγ increases further, P γ t even shows a non-monotonic behavior -first increasing with g 2 aγ at small g aγ and then turning over and decreasing as g aγ increases further. Increasingly higher order terms in g aγ ϕ 0 can no-longer be ignored as C ∼ 1/2.
For our fiducial parameters for the profile, we see that the dipole estimate can differ by more than an order of magnitude as C becomes order unity (see Fig. 8, left panel). We also find that for certain C, even the dominant frequency content of the radiated power can change from 2ω to 4ω signalling a cancellation between various higher order terms in g aγ at the frequency 2ω! In Fig. 9 we show a comparison between the frequency(=wavenumber) content of the radiation in our simulation box at for C = 0.09 and C = 0.85. Notice the vanishing of the dominant k ≈ ω as we go from C = 0.09 to C = 0.85.
Background E: For the case of the background E field, we find that the P γ (2) t underestimates the power as we move to larger C. See Fig. 7. The frequency content of the radiated power continues to be dominated by 2ω. The detailed reason for difference in behaviour of the radiated power between the background B and E fields is not entirely clear to us. However, we do note that it might be sensitive to the fiducial parameters chosen. By reducing the radius of the soliton, we were able to also find a regime where the P γ (2) t overshoots the numerical results even with an E field background.

Large coupling with backreaction
When the coupling C > C crit ≈ 1.3, we transition to exponentially growing E and B fields. This exponential growth can be rapid enough so that an order unity fraction of the energy of the soliton is extracted from the soliton within the duration of our simulations. The duration to backreaction depends on the initial energy in the soliton, as well as initial conditions. For the present case the background B and E fields generate fluctuations in the electromagnetic fields, which are then enhanced via parametric resonance.
To include this backreaction on the soliton, we evolve the fully coupled axion-photon system dynamically. The results are shown in Fig. 10. The key point to note is that backreaction naturally regulates the exponentially growing radiated power once sufficient energy has been extracted from the soliton.

Medium effects and coherence
In order to illuminate the connection between our calculation and calculations of axionphoton conversion in the literature, we address three issues here. First we discuss how our calculation should be modified in the presence of a medium, such as the dense plasma around a compact star, and the associated phenomenon of resonant conversion. Second we clarify how the axion star's coherence has affected the final radiation power. Third we compare our calculation with a perturbative calculation of the axion-photon conversion probability. Here, we will restrict our attention to leading order in g aγ ϕ 0 to simplify the discussion.

Medium effects and "resonant" conversion
Our previous calculation was done with the axion star in the presence of background electromagnetic fields, with no other medium present at the background level. However, in applications to astrophysical scenarios, such as axion-stars in the magnetosphere of a neutron star, a plasma is present that leads to the photon having an effective mass ω p . The effect of such a constant effective mass ω p < ω can be approximately taken into account by modifying eq. (4.8) and (4.9) making the replacement ∇ 2 → ∇ 2 − ω 2 p . Note that if ω < ω p the propagating mode would be exponentially suppressed. 8 Following through with the same calculation as before, but carefully keeping track of κ and ω separately, we arrive at the generalization of our equation for the radiated power (4.14): whereφ(κ) is the Fourier transform of the axion field profile at |k| = κ. Note that for ω p = 0, κ = ω and we recover our earlier result without the medium (4.14). However, for ω p → ω ("resonant conversion" domain), we have to be careful. Importantly, in the limit κ → 0, we 8 Since we assume the background E and B fields to be constant, it is natural to assume ωp is constant, although in practice it does depend on the spatially varying free charge density also. In neutron star atmospheres, an approximation to the plasma frequency is given by ωp = 4παemne/me where ne is the Goldreich-Julian charge density [130], and me is the mass of the electron. Note that ne(x) ∼ Ω · B(x) where Ω is the angular velocity of the neutron star.  Figure 11. The dependence of the time averaged power radiated by the axion star as a function of the plasma frequency ω p and radius R of the star. Note the "resonant conversion" when ω p → ω (the soliton frequency), and the strong dependence on the radius of the star. The plot should be understood under the assumption that the plasma frequency is constant within the radius of the star. In the above plot, ωR varies between 1 and 10, with the peak scaling as (ωR) 6 . Also note that while the power increases with radius for ω p ≈ ω, it decreases with radius when ω p ω for large radii.
can remove the exponential suppression inφ(κ) (see eq. (4.16)), and we obtain In the above expression we have assumed that ω p does not vary within R. Notice that the medium effects point us to move to larger axion stars to get a large amount of power emitted, whereas without the medium, we must limit ourselves to a smaller radius because of the exponential suppression (assuming ω p is constant on the scale R).
To see the detailed dependence of the radiated power on the radius R and ω p , see Fig. 11. Note the large enhancement as we increase the radius for ω p approaching ω from below. This is the resonant conversion. The ratio of the power emitted when ω p ≈ ω above, compared to that when ω p = 0 is (1/16)(πωR) 2 e πωR 1 − ω 2 p /ω 2 .

Spatio-temporal coherence
The dipole radiation power that we calculated in Sec. 4.2.1 was derived with the assumption that the value of the axion field at different points across the soliton all oscillate in phase.
That is, we have phase coherence across the entire configuration. It is worth exploring the importance of this coherence for our results. To this end, suppose that the charge density generated by our axion configuration is replaced by N idealized, equally spaced charge dipoles. Each dipole oscillates with frequency ω but a random phase θ j . The total charge density is 3) B B Figure 12. A soliton of the axion field in an external magnetic field creates a coherently oscillating dipole configuration, which leads to dipole electromagnetic radiation. If instead, we replace the solition with N oscillating dipoles with random phases, we can get less power radiated than the coherent soliton case for sufficiently large N . as shown in Fig. 12. Recall that (1) (x) = −g aγ ∇ϕ ·B. Then, under the assumption that the phases are random, we have where k = ωx, we have assumed that N is large, V is the volume in which the (1) (x) is non-zero. Note that (∆x) 3 = V /N , and we defined Since the power radiated is proportional to | N (1) (k)| 2 , we can now compare the coherent and incoherent cases P γ That is, if N > Q 2 0 |˜ (1) (k)| −2 , the radiated power will be larger from a coherent configuration. For our sech profile and a constant background B field, we have |˜ (1) (k)| 2 ≈ (4π 2B2 /ω 4 )(g aγ ϕ 0 ) 2 (πωR) 4 e −πωR and Q 2 0 ∼ 10 −1 (πωR) 4 (g aγ ϕ 0 ) 2B2 /ω 4 , which tells us that we need N 10 −3 e πωR for coherence to win. 9 For some localized configuration of radius R with a characteristic density 0 ∼ Q 0 /R 3 , we can define a coherence length: There can be purely numerical coefficients in front that depend on choice of V and the details of the profile. The scaling with ωR, is the main result we want to focus on.
If we subdivide the volume of our coherent configuration into N incoherent regions, each with a volume smaller than λ 3 C , then the power radiated from the coherent configuration will be larger. For our specific case of interest related to our soliton profile, we get λ C ∼ e −πωR/3 R. Hence for large radius configurations, incoherent emission will typically dominate over the coherent one.

Axion-photon conversion probability
We have calculated the electromagnetic radiation from an axion star in an external magnetic field by recognizing that the soliton behaves like a coherent dipole antenna. However, the conversion of axions into photons can be understood from a different point of view. The axion and photon fields mix with one another in the presence of an external magnetic field [131], and an incident axion develops a nonzero probability to be detected later as a photon. This phenomenon is the basis of many laboratory probes of axions [80].
In the presence of an external magnetic fieldB, the action from eq. (2.1) contains a mixing, L mix = −g aγ φȦ ·B, where we work in the Weyl gauge with A 0 = 0. To leading order in the coupling, the probability for an axion with momentum p µ = {E p , p} to convert into a photon with momentum k µ = {ω k , k} is [132] whereB(q z ) is the Fourier transform ofB(z). We have assumed that the magnetic field is static and only varies in the z direction. The longitudinal momentum transfer, q z = k z − p z , is restricted by energy and transverse momentum conservation to only take on two values, q ± = ± p 2 z + m 2 − p z . If the magnetic field has a top-hat shaped profile, meaning that it is only nonzero (and takes value B 0 ) for a region of longitudinal distance , then the Fourier transform isB(q z ) = (2B 0 /q z ) sin(q z /2). Moreover if the incident axion is non-relativistic, then q ± ≈ m. For a fiducial volume V containing N a axions, the power per unit area of photons being emitted in the z-direction is P γ /A = ω k N a P a→γ /V , which corresponds to P γ = ρ a P a→γ A where ρ a = mN a /V ∼ m 2 ϕ 2 0 is the axion energy density and ω k = m. It is interesting to compare this calculation with the classical dipole radiation formula from eq. (4.17). They display a similar parametric behavior, and both powers scale as P γ ∝ g 2 aγ ϕ 2 0B 2 R 2 . If we further identify the size of the B-field filling region with the radius of the axion star, = R, then we have a factor of sin 2 (mR), which also arises from an axion star with a top-hat density profileφ 2 (ω = m) ∝ sin 2 (mR). However, the two calculations are not necessarily equivalent, since eq. (4.17) is the power output from an (inhomogeneous) axion star in a homogeneous magnetic field, whereas the estimate above is for a homogeneous axion flux in an inhomogeneous (longitudinally-varying) magnetic field. We believe that these two approaches will yield consistent results when put onto the same footing (soliton structure and coherence), and we defer this investigation to future work. 10

Observational signatures
The central goal of our work is to understand the emission of electromagnetic radiation that occurs when an axion star passes through a strong electromagnetic field. We have seen that the radiation spectrum peaks at E ∼ ω ∼ m a , which corresponds to radio frequencies for typical axion masses. We have also seen that the radiation power grows as P γ ∝B 2 with the strength of the external magnetic field. In this section we will discuss how this phenomenon could lead to a variety of observational signatures in different environments with strong magnetic fields. We again restrict our attention to results at leading order in g aγ , although the large g aγ results might lead to more radiated power in some cases.
Using the dipole approximation from eq. (4.14) and (6.2), the luminosity (L ≡ P γ (2) t ) of an axion star in a background magnetic field (with strengthB) is estimated as where we have normalized the axion-photon coupling g aγ to the 95% CL upper limit from the CAST helioscope [133], and we have set ω = m. We also remind the reader that 1 W = 10 7 erg/sec. The function F holds information about the soliton shape and plasma effects: where ω p is the plasma frequency. Note the beneficial dependence on large radius and the lack of exponential suppression in the "resonant" (ω p ≈ ω) case. As long as the radius of the star is smaller than the size of the resonant region, our calculation holds, and leads to a large enhancement in the radiated power compared to the non-resonant case.
For the estimates in this section, we approximate the radiation spectrum as monochromatic, corresponding to a single spectral line. The frequency of this line is taken to be and we take the width of the line, i.e. the signal bandwidth, to be ∆ν γ ∼ ν γ . For these fiducial parameters, we also note that the mass scale and radius of a very dense axion star (soliton) are expected to be on the order of M sol ∼ 10 2 f 2 /m 2 × 10 9 kg f 10 10 GeV Note that 2 × 10 9 kg ≈ 10 −21 M .

Compact stars
The strongest magnetic fields in the universe today can be found in the magnetospheres of compact stars. The magnetic field strength at the surface of a white dwarf star is typically 10 6−8 G [134] whereas the smaller neutron stars can reach 10 12−14 G [135]. If an axion star were to encounter these extreme magnetic fields, the result would be a sudden and extreme release of electromagnetic radiation [67].
If the compact star is a distance d away, then the flux of radiation reaching Earth is F = L/(4πd 2 ), which can be measured in erg/cm 2 /sec. The corresponding spectral flux density is calculated as S = F/B where B = ∆ν γ = ω/2π is the signal bandwidth. For a nearby star, the spectral flux density evaluates to S 2 × 10 7 µJy d 100 pc whereas the flux from a star at the galactic center (d ≈ 8 kpc) would be reduced to S 3 × 10 3 µJy. For reference, an hour-long observation with a current or planned telescope (such as GBT, JVLA, or SKA) would have a flux sensitivity of δS ∼ 1 µJy; see the estimates in Refs. [61,72]. If an axion star were to pass through the magnetosphere of a compact star while it was being observed by a radio telescope, then the signal could be quite striking, even for modest couplings and field strengths.
Since the compact star is surrounded by a plasma, this must be taken into account for the signal strength estimates. In Sec. 6.1 we have argued that the finite plasma density modifies the radiation spectrum, which is captured by F in eq. (7.2). This factor depends on the plasma frequency ω p , which is grows larger at points closer to the star, and F peaks near to where the plasma frequency matches the soliton's oscillation frequency, ω p ≈ ω, as shown in Fig. 11. For example, using the fiducial parameters in Ref. [61], the width of the resonance region is estimated to be ωL ∼ O (100). If the axion star's radius is R ∼ 0.1 L then Fig. 11 implies an enhancement of F ∼ 10 4 to the spectral flux density estimate from eq. (7.5), which further increases the detectability.
Even if an axion star's encounter with a compact star could be detected, we must address the expected rate of these encounters [67,72,74]. The encounter rate between a particular compact star and the ambient population of axion stars is estimated as Γ = σ eff v rel n as where σ eff is the effective cross sectional area for the scattering, v rel is the typical relative velocity, and n as is the number density of axion stars (near the target compact star). We can also write n as = ρ as /M sol where ρ as is the local mass density in axion stars and M sol is the typical energy per axion star (soliton). The effective cross sectional area is further enhanced by the gravitational focusing factor, and we estimate σ eff = (1+v 2 esc /v 2 rel )πR 2 where v 2 esc = M /4πm 2 pl R is the escape velocity at the surface of the neutron star. Combining these factors allows us to estimate the encounter rate of axion stars with a particular white dwarf star to be whereas the rate for encountering a neutron star (with R = 10 km and other fiducial parameters unchanged) is Γ 5 × 10 −8 hr −1 . The fiducial axion star density is taken to equal the local dark matter energy density near Earth, ρ dm = 0.3 GeV/cm 3 , although axion stars are not expected to compose an O(1) fraction of the total dark matter density, which is typically dominated by a diffuse population of axion particles.
The estimate in (7.6) appears very unfavorable. For the fiducial parameters we expect a particular white dwarf star to encounter an axion star approximately once every 3 years (or once every 2000 years for a neutron star). However, there are several reasons why the rate might be enhanced over these estimates. First, the rate increases for compact stars at the galactic center or within dark matter subhalos (where ρ as is higher if it tracks the dark matter density). Second we estimated σ eff using the star's geometrical cross section, ∼ R 2 , whereas the magnetic field extends far beyond the boundary of the star and scales like B ∼ r −3 for a magnetic dipole. Third, depending on the nature of the observation, it may be necessary to integrate over a finite region of the sky, such as toward the galactic center, which could contain many neutron stars, further increasing the encounter rate [136]. Fourth, the fiducial axion star mass M sol = 10 9 kg is a free parameter, and a smaller value implies a larger encounter rate.

Direct detection in our solar system
The phenomenon of electromagnetic radiation from an axion star in an external magnetic field could be used to develop a strategy for detecting axion stars when they encounter our solar system. The strongest magnetic fields generated in laboratories on Earth can reach strengths of a few Tesla, corresponding to ∼ 10 4 G. However, the flux of axion stars at Earth is expected to be quite low, making these signals very unlikely to be observed. The flux is estimated as Φ = ρ as v rel /M sol , and the expected encounter rate with a 1 meter-scale detector is Γ = Φ(100 cm) 2 ∼ (10 −17 yr −1 ) (ρ as /0.3 GeV/cm 3 )(v rel /10 −3 )(M sol /10 9 kg) −1 . Going beyond the confines of the laboratory, the Earth sustains its own magnetic field with a strength of ∼ 1 G. The smaller field strength would lead to a weaker signal, but the larger volume implies an increased encounter rate, Γ ∼ (10 −4 yr −1 )(M sol /10 9 kg) −1 . Finally, axion star encounters with the Sun's ∼ 1 G magnetosphere could also provide a channel for detection. The encounter rate is enhanced by the Sun's much larger surface area, giving Γ ∼ (10 0 yr −1 )(M sol /10 9 kg) −1 , but the radiation power is much weaker, approximately P ∼ 10 −24 L , making this signal undetectable for the fiducial parameters.

Galactic magnetic field
The axion stars in our Milky Way galaxy are continuously exposed to its ∼ 10 −6 G magnetic field. The corresponding electromagnetic radiation power is estimated using eq. (7.1). For the fiducial parameters used above we find P ∼ 4 × 10 −10 W for a single axion star. This power output is incredibly weak. For reference, if we sum the power output from all of the axion stars in a galaxy like the Milky Way (assuming that they make up all the dark matter), then the net power output is still only 10 −3 L ! However, see also Ref. [137] for a discussion of resonant axion-photon conversion in the intergalactic magnetic field.

Early universe
We know very little about the extreme environment of the Universe during the first fractions of a second after the Big Bang. Some theories predict that a magnetic field may have arisen during the period of cosmological inflation, post-inflationary reheating, or during a subsequent cosmological phase transition [138,139]. The strength of this primordial magnetic field may have been incredibly large by our every-day standards. For instance a study of magnetogenesis from axion inflation [140] concluded that magnetic field generation could be so efficient as to transfer an O(1) fraction of the inflaton's energy into the magnetic field, leading to field strengths as large as ∼ 10 52 G at the end of inflation (for an inflaton mass of m inf ∼ 10 14 GeV).
Formation of oscillons and dense axion-star configurations has been explored in earlier works [87][88][89][90][91]141]. In well-motivated, observationally constrained models of inflation, the universe can become dominated by such solitons at the end of inflation (if the coupling to other fields is sufficiently weak) [88]. Similar phenomena are possible in moduli fields and other (pseudo-)scalars in the early universe. Typically, such configurations are long-lived compared to the age of the universe then, although they are not expected to survive until the present day.
If the early universe were to contain both a strong primordial magnetic field and a population of axion stars [142][143][144][145][146], then their interaction will induce electromagnetic radiation from the axion stars, thereby precipitating their decay (but also raise the question of whether the solitons would form in the first place). Recall from the estimates in eq. (5.3) that an axion star with mass M sol ∼ 100f 2 /m emitting with a power P γ ∼ 10B 2 (f g aγ ) 2 /m 2 would exhaust an O(1) fraction of its energy on a time scale of τ ∼ M sol /P γ ∼ 10m/g 2 aγB 2 . Sincē B can be very large in the early universe, this axion star lifetime can potentially drop below the Hubble time scale at that time, which is t H ∼ m pl /T 2 during radiation domination at temperature T . We also note that even without strong magnetic fields, the solitons might be able to decay into photons rapidly due to collisions via mechanisms similar to those discussed in [41].
The phenomenon of magnetic-induced axion star decay would be challenging to test, since we have only a few handles on early universe physics. If the decay happens to occur during primordial nucleosynthesis, then the injection of electromagnetic radiation into the primordial plasma could potentially disrupt the formation of the light nuclei [147], and measurements of the light element abundances would provide an indirect constraint on this scenario. Nucleosynthesis also provides strong constraints on the QCD axion, even in the absence of a primordial magnetic field [148]. The interplay between solitons and electromagnetic fields can affect the rate of energy transfer and equation of state during reheating, as well as gravitational wave production, and spectral distortions during the early universe [90,126,127,[149][150][151][152].

Summary and conclusion
A spatially localized, periodically oscillating axion configuration (soliton: oscillon, axion star etc.) in background electromagnetic fields, sources electromagnetic radiation. We investigated this production analytically and numerically (with 3+1 dimensional lattice simulations when necessary), focusing in particular on the dependence of the emitted radiation on the characteristics of the axion field configuration as well as the strength of the coupling to the electromagnetic field. We also pointed out how the coherence of the soliton configuration, as well as the plasma effects, change the radiated energy in electromagnetic fields.
Our key results regarding the radiated power (luminosity) in electromagnetic waves are as follows: • We delineated and verified the boundary between bounded, constant luminosity solutions and exponentially growing ones based on axion-photon coupling and soliton properties. For a soliton with central amplitude ϕ 0 , oscillation frequency ω and radius R, this boundary lies at C ≡ g aγ ϕ 0 ωR/4 ∼ 1. This boundary is independent of the background electromagnetic fields.
• For C 1, we get dipole radiation with a constant time-averaged luminosity. We derived an explicit formula for this dipole radiation, including an understanding of the strong (exponential) dependence on the radius of the solitons.
• For the dense solitons (which we explore in detail in the numerics), we see a rich behavior of the radiated power as g aγ is varied to explore all scenarios from C 1 to C crit ∼ 1. Although the time-averaged radiated power remains constant in time, the details of the magnitude of the radiated power differ between background E and B field cases, and they are also sensitive to the details of the soliton configuration. For the B field case, for all cases we have considered, we see a suppression compared to the dipole estimate and a non-monotonic behavior with C. The same is not true for the E field background. 11 • For C 1, parametric resonance leads to an exponentially (in time) growing luminosity based on Floquet Theory. In the parametric resonance regime, background electromagnetic fields are unnecessary, small fluctuations can be sufficient. The exponential transfer of energy can be significant enough to cause backreaction on the soliton, and regulate photon production.
• We explained the relevance of the coherently oscillating axion field configuration compared to an incoherent collection of dipoles and defined a critical coherence length which allows us to determine whether the coherent or incoherent configuration would radiate more efficiently.
• We explored how the presence of a plasma affects the radiation from a soliton. In particular, we find that when the plasma frequency is approximately equal to the oscillon frequency we get an enhanced resonant conversion to photons ('resonant' conversion, which is different from parametric resonance).
There are a number of avenues for future work to extend our results. Our formalism and code includes background electric and magnetic fields together, however, we presented detailed numerical results for each separately. Considering them both together would introduce additional rich phenomenology which might be necessary, for example, when axion stars are boosted through static fields or when the astrophysical background fields themselves are time-dependent as is the case with neutron stars. In future work, we also plan to numerically include gravitational effects, such as tidal disruption, and take time-dependent medium effects around compact stars into account.

Acknowledgements
The numerical simulations were carried out on the NOTS cluster supported by the Center for Research Computing at Rice University. MA is supported by a NASA ATP theory grant NASA-ATP Grant No. 80NSSC20K0518, and PMS acknowledges support from STFC grant ST/P000703/1. We thank Yang Bai for helpful comments on the draft. We also thank Kun Hu, Mudit Jain, Siyang Ling and Hongyi Zhang for useful discussions regarding the radiated power as well as neutron star atmospheres. 11 We have checked that the radiated power is constant in this regime by doubling the linear size of the box and the duration of the simulation. While this constancy is expected from Floquet theory, it is not quite a proof since the possibility of band structure in C with complicated boundaries also exists which we might have missed out on numerically. Furthermore, the boundary at C = Ccrit might be richer than just going from a constant-in-time radiated power to an exponentially growing one. We cannot exclude the possibility of a power-law behavior with time for the radiated power at this boundary. We leave this investigation to future work.

A Dipole radiation Green's function
In this appendix we solve the field equations, (4.8) and (4.9), using the method of Green's functions, and we derive the time-averaged Poynting vector S t . Consider the retarded Green's function Solutions of eqs. (4.8) and (4.9) are written as where the integration contour M is understood as the upper half plane (time t > 0), while the boundary ∂M is located at time t = 0. In each expression above, the second and third terms enforce the initial conditions. However, for the purposes of calculating the late-time radiation at a point far away from the localized charge distribution, we can safely ignore these terms. Notice that these terms depend on the Green's function through G(t, x; 0, x ) ∝ δ(t − |x − x |). When x is restricted in the oscillon region and x is fixed for the observation, the delta function in the Green's function can not be satisfied given t is big enough.
The first two terms are independent of time t, while the third term oscillates with period π/ω and it vanishes upon taking the time average (over many oscillations cycles). Thus, the time-averaged Poynting vector is where k = ωx, for a spatially-localized, spherically-symmetric charge distribution that oscillates with period 2π/ω. A similar calculation goes through for a photon with nonzero mass 0 < m γ < ω. In a medium m γ = ω p is the plasma frequency. The field equations, (4.8) and (4.9), are extended to include the mass term. The Greens function in eq. (A.1) involves a massive propagator that enforces the dispersion relation, k 2 0 − |k| 2 = ω 2 p . Ultimately the Poynting vector is expressed as in eq. (A.11) but with the wavevector k = ωx replaced by k = κx with κ = [ω 2 − ω 2 p ] 1/2 . For larger photon masses, ω < ω p , the fields are exponentially damped.