Time-domain properties of electromagnetic signals in a dynamical axion background

Electromagnetic waves in a dynamical axion background exhibit superluminal group velocities at high frequencies and instabilities at low frequencies, altering how photons propagate through space. Local disturbances propagate causally, but unlike in ordinary Maxwell theory, propagation occurs inside as well as on the lightcone. For the unstable modes, the energy density in the electromagnetic field grows exponentially along timelike displacements. In this paper we derive retarded Green functions in axion electrodynamics in various limits and study the time-domain properties of propagating signals.


Introduction
Axion models provide some of the most well-motivated extensions to the Standard Model, providing a mechanism to resolve the strong CP problem and a class of dark matter candidates. Through its coupling to gluons, the vacuum expectation value of the axion field cancels theθ parameter of quantum chromodynamics (QCD) and restores CP symmetry [1,2], explaining the surprisingly small experimentally measured value, θ < 6 × 10 −11 [3,4]. Through the misalignment mechanism [5][6][7], axions can also be produced in the early universe in sufficient abundance to comprise most or all of the dark matter.
Many of the most stringent constraints on axion models utilize a coupling between the axion and the Standard Model electric and magnetic fields, This coupling enables axion production through the Primakoff process [8,9]; axion decay to photons; and axion-photon interconversion in the presence of electromagnetic fields. As such it can be used either to detect or to produce axions in the laboratory. Stellar observations [10][11][12] constrain g aγγ < 10 −10 GeV −1 for a wide range of axion masses, and additional constraints set by the power spectra of bright X-ray point sources provide a more stringent limit of g aγγ 10 −12 GeV −1 for light axions of mass m a 10 −12 eV [13][14][15]. At large occupation numbers and de Broglie wavelengths, axion dark matter behaves as a classical, oscillating background field that induces small time-dependent perturbations to electrodynamics, which can be probed with a variety of different sensitive experimental techniques [9,[16][17][18][19][20][21] .
The coupling in Eq. (1) also affects the propagation of classical electromagnetic radiation [22][23][24]. Electromagnetic plane waves traveling through an axion background acquire modified phase velocities for left-and right-polarizations, an effect which may be observable in interferometers [25,26], atomic clocks [27] or astrophysical sources [28], for some ranges of axion masses and couplings. Furthermore, low frequency modes exhibit tachyonic instabilities, while at high frequencies, group velocities for both polarizations are superluminal.
Despite the presence of plane wave solutions with superluminal group velocities, axion electrodynamics is a causal theory: local disturbances do not propagate outside the lightcone. This was first shown long ago in the case of tachyonic scalar field theory by Aharonov, Komar, and Susskind [29]. Here we show that the same mechanism is at work in axion electrodynamics, with consequences that include the exponential growth of local disturbances.
In this paper we calculate the classical electromagnetic retarded Green function in a coherent, dynamical axion background in several disparate regimes of axion parameter space. Our results are organized based on the hierarchical ordering of three different scales: • m a , the axion mass; • µ 0 ≡ 1 2 g aγγ √ ρ a , a mass scale that determines the rate of exponential growth of the electromagnetic fields, based on the axion density ρ a ; and • T −1 and L −1 , the inverses of the characteristic propagation time and distance T ∼ L of a signal.
Rather than focusing only on those axion models that provide a natural solution to the strong CP problem, we consider the broader realm of axion-like particles (ALPs), where the axion mass m a and decay constant f a are not required to satisfy m a f a ∼ m π f π , and the value of g aγγ is not determined by m a . These axions can still provide a wide range of dark matter candidates (see, e.g., [30][31][32]). For the lowest-mass "fuzzy dark matter" candidates, generic constraints on ultra-light scalars from Lyman-α data [33][34][35][36] impose a lower bound on the axion mass of m a 2 × 10 −21 eV, though in the "large misalignment" regime of Ref. [37] the Lyman-α bound is altered by the effect of ALP self-interactions. Our analysis in this paper encompasses the nearly twenty decades of ALP parameter space above this bound, where the axion is still light enough that it can be treated as a coherently oscillating background field. In terms of the axion virial velocity v, the mass sets an upper bound on the characteristic L and T , L (m a v) −1 and T (m a v 2 ) −1 , after which any analysis must incorporate the effects of decoherence. In the example of fuzzy dark matter with m a ≈ 2 × 10 −21 eV, and taking v ∼ 10 −3 as the ALP virial velocity, the coherence length and time are respectively L c 10 17 m ≈ 3 pc and T c 3 × 10 11 s ≈ 10 4 yr. On the other extreme, for m a much larger than 10 −4 eV, even a tabletop experiment will encounter significant decoherence. Our analysis is focused on the nonrelativistic limit, neglecting these decoherence effects and ignoring spatial gradients in the axion field.
We begin in Section 2 with the simplest analysis, the m a µ 0 , 1/T limit. In this case the value of the axion field changes at an approximately constant rate, ∂ t a(x, t) ≈ const, and we find an analytic solution to the Green function valid for all values of µ 0 T . This Green function exhibits exponential growth inside the lightcone of the disturbance when µ 0 T 1. Despite the potentially catastrophic consequences of this unbounded growth, the dilute density of dark matter and experimental constraints on the axion-photon coupling ensure that the timescales for the genuinely exponential phase of the growth are outside the reach of all but the lightest ALP candidates, unless the local ALP density ρ a is enhanced by several orders of magnitude above 0.4 GeV/cm 3 . In Section 2.3 we highlight some of the curious and potentially detectable perturbations to classical electrodynamics induced by the axion background.
For almost all allowed values of g aγγ and ρ a , the hierarchy µ 0 m a is more realistic, and we explore this limit in Section 3. In the case of the oscillating background axion field it is no longer possible to derive an exact analytic expression for the Green function using the methods of Section 2. Instead, we construct perturbative expansions for the µ 0 T 1 and µ 0 T 1 limits by expressing the Green function as a continued fraction. When the frequency support of the radiation includes ω ≈ 1 2 m a , a narrow resonance induces exponential growth for large T . In Section 3.1 we calculate the dominant part of the Green function in the µ 0 T 1 limit. In this late-time limit the resonant enhancement dwarfs the contribution from frequencies ω = 1 2 m a . This resonant emission has been previously studied in [38,39], although the resonant band is so narrow that dispersion effects and gravitational redshifting may completely prevent the exponential growth [40]. For the non-resonant limit µ 0 T 1, and for electromagnetic signals which do not include support near the resonant frequency ω ≈ 1 2 m a , Section 3.2 provides a continued fraction expression for the Green function that is valid to arbitrary order in g aγγ . In Sections 2.3 and 3.3, we provide numeric examples to illustrate the behavior of signals in various corners of ALP parameter space, and to verify our analytic expressions.
The primary results of this paper are collected in Eqs. (30), (65) and (73) in Section 4. Despite the significant differences between the two limits, the Green functions of Sections 2 and 3 both exhibit the novel inside-the-lightcone propagation and exponential growth in certain modes.

Axion Electrodynamics
In terms of θ(x µ ), the local value of the effective CP violation induced by the axion background, the Lagrangian for electrodynamics includes the interactions where F µν is the electromagnetic field strength tensor, A µ and J µ are the vector potential and 4-current, and θ(x µ ) is related to the value of the axion field via In Lorenz ∂ α A α = 0 gauge, the equations of motion for A µ reduce to which depends explicitly on the derivatives of θ(x µ ) rather than θ itself. Taking the external source to be neutral and transverse, J 0 = 0 and ∇ · J = 0, and neglecting any spatial gradients in the background axion field, |∇a| |ȧ|, the equations of motion for the scalar and vector potentials decouple, 1 For a typical model of ALP dark matter,θ(t) is given by where the value of θ 0 is set by the local axion dark matter density, ρ a ∼ (0.042 eV) 4 , and wherė For future reference, we note that 1.75 · 10 −23 eV 2.66 · 10 −8 Hz 0.84 yr −1 .

Green Functions for the Steady-State Background
For timescales that are short compared to the period of the axion oscillation, m a t 1 and m a θ , it is appropriate and instructive to consider the approximation ∂ 2 t θ ≈ 0, whereθ(t) assumes a nearly constant value −θ 0 ≤θ(t) ≤θ 0 . In this case the differential equation for A can be solved using a Fourier transform, so that the differential equation for A becomes a set of algebraic equations for A i . For planar waves propagating in theẑ direction, k = kẑ, the polarization basis where j ± (k) is the Fourier transform of the transverse source (ˆ ± · J) in the polarization basis. Circularly polarized plane waves propagate with the dispersion relations 1 The equations of motion with ∇a = 0 are discussed in e.g. [23,41].
producing subluminal or superluminal phase velocities depending on the sign ofθ and the polarization of the radiation. The group velocities for both modes are superluminal, for both positive and negativeθ, with dω ± /dk = 1 only forθ = 0. The effect appears at quadratic order inθ/k, so that forθ k the modification of the group velocity is subdominant to the O(θ/k) change in the phase velocity, We demonstrate below that the retarded Green function vanishes outside the light cone, preserving causality despite the presence of superluminal group velocities. As a necessary consequence, disturbances in the field induced by local sources grow exponentially in timelike directions. To leading order inθ, the phase velocities alternate about a central value ω ± /k = 1 based on the polarization of the light and the sign ofθ. After multiple periods of the axion oscillation, the perturbations to the phase velocity tend to cancel each other. On the other hand, the group velocity is superluminal for both positive and negativeθ, so the exponential growth is not ameliorated by any periods of exponential decay when the sign ofθ changes. The effects from the modified group velocities should grow over time.
With the approximation thatθ(t) is nearly constant, the two-dimensional Green function can be obtained analytically to all orders inθ. Experimental constraints on g aγγ < 10 −10 GeV −1 and m a exclude the m a θ possibility unless the local axion density is significantly enhanced, ρ a O(GeV/cm 3 ), so the results in this section are directly applicable primarily to situations involving to dense clumps of ultralight axions. When we calculate the Green function for the more broadly relevant m a θ hierarchy of scales in Section 3, the steady-state case with constantθ also provides a helpful consistency check in the limit where the exponential growth becomes important.

Green Function Solution in Two Dimensions
The Green function can be calculated analytically for the simplified case of plane waves k = kẑ with a spatially homogenous (∇θ = 0), steady-state (∂ 2 t θ = 0) axion background. Imposing translational symmetry in x and y effectively reduces the system from (3 + 1) dimensions to (1 + 1). The A ± equation of motion in Eq. (9) admits a Green function solution g ± (z, t) of the form where In this section it is convenient to fold a factor of 1/2 into the definition ofθ, where µ determines the rate of exponential growth, as we show below. Defining a related Green function G 0 , Eq. (15) can be simplified to so that µ 2 acts as an effective tachyonic mass for the scalar-like Green function G 0 . Applying the Fourier transform and integrating both sides of Eq. (18) produces the integral form of the Green function, where ε > 0 indicates that the contour in the complex ω plane should correspond to the retarded Green function, which vanishes for t < 0. For k 2 > µ 2 the ε → 0 + limit can be recovered easily. However, for k 2 < µ 2 one of the poles in ω is located above the real axis, at To recover the retarded Green function, the contour in ω should pass above both poles, with ε → µ + on the imaginary axis.
The dω integral can be completed using the residue theorem, where ω 0 = k 2 − µ 2 = i µ 2 − k 2 and Θ(t) is the step function. The integral is simplified by a coordinate substitution k → ϕ, where the contour L in the complex ϕ plane is shown in Figure 1, and allows cosh ϕ to vary smoothly from −∞ to ∞ with Im(cosh ϕ) = 0. By mapping the coordinates (t, z) to (η, λ) via for λ ≥ 0, the integral can be written in terms of just λ and ϕ ± η: Figure 1: Above, we show the contour L that allows Eq. (??) to be written in terms of Hankel functions. The shorthand notation for this contour is R 1 1+i⇡ d', to indicate that the contour approaches the imaginary axis from negative real 1 with a constant imaginary component of +i⇡.
2 Figure 1: The contour L that allows Eq. (24) to be written in terms of Hankel functions. The shorthand notation for this contour is ∞ −∞+iπ dϕ, indicating that the contour approaches the imaginary axis from negative real ∞ with a constant imaginary component of +iπ.
Here the notation ∞ −∞+iπ dϕ indicates the imaginary offset for Re(ϕ) < 0 shown in Figure 1. In order to make this last simplification, removing the η dependence completely, note that the integrand has no poles for finite ϕ , so that the contours L ± from the coordinate substitutions ϕ ± = η ± ϕ can be shifted horizontally to compensate for η.
Outside the lightcone, for spacelike displacements z 2 > t 2 , the coordinate transformation Eq. (23) is replaced by the alternative µz =z cosh η, µt =z sinh η, In this case, with λ < 0, the two contributions to the integral cancel each other, and so the retarded Green function vanishes outside the lightcone. Using this contour notation, and the integral definition of the Hankel functions G 0 can be written as Eq. (28) simplifies for integer values of ν to recover an expression in terms of the modified Bessel function of the first kind, I ν (z): Finally, we obtain As anticipated by the discussion of tachyonic scalar theories in Ref. [29], the Bessel function I 0 ( √ λ) grows exponentially with timelike λ ∝ t 2 − z 2 . In the λ 1 limit it approaches The retarded Green function vanishes outside the light cone, t 2 < z 2 , and it explicitly satisfies the equations of motion From the solution for A ± one can derive the corresponding Green functions for the E and B fields, with E z = B z = 0.
In many of the simplest cases of interest, including idealized axion interferometers and astrophysical sources, the photon source is localized in space and produces a signal s(τ, z) = δ(z)s(τ ) or alternatively s(τ, z) = s(τ − z) that varies in time and propagates in the forward direction, z ≥ 0. With this choice, the E and B fields can be found from the Θ(z) components of the derivatives of A ± , and the familiar Standard Model limit µ → 0 is recovered by

Green Function Solution in Four Dimensions
The Green function Eq. (30) is valid for sources that are spatially uniform in the x and y directions. In many situations, including laser pulses and the light from distant stars, this approximation is sufficient. However, in other cases the fully four-dimensional Green function may be relevant. As we show in this section, most of the 4D solution can be written in terms of the 2D Green function derived in Section 2.1. One of the new integrals cannot be so easily solved analytically, but with some effort it can be put in the form of a rapidly converging infinite series for easier numerical evaluation.
In the Lorenz gauge with neutral sources, ∇ · A = 0, A 0 = 0, the Green function for A satisfies where the cross product term in Eq. (5) forces the Green function G kj to have a nontrivial tensor structure. Its Fourier transform G kj satisfies in terms of µ from Eq. (16), frequency ω, and k 2 = k i k i for i = 1, 2, 3. Inverting the operator that acts on G kj , the Green function can be written as where where the four roots of β are Since the current sources are transverse, the Green function can be simplified by the transverse projection where A, B and A ij are the Fourier transforms of A, B, and − Ak i k j /k 2 , respectively. Both A and B can be written in terms of the scalar function G 0 (t, z) from the two-dimensional case, Eq. (30), where Attempting the same technique for A ij leads to the incomplete expression an integral that does not have a simple expression in terms of Bessel functions or other hypergeometric functions. In Appendix B we show how the integral form of A ij can be replaced with an infinite series over a product of hypergeometric functions, with the result The series in converges rapidly for µr ≤ 2. For large µr 1 and µt 1 it converges for > max , for an max (µt, µr) given in Appendix B.

Application to Monochromatic Signals
The Green function for the vector potential A (Eq. (30)) and its derivatives exhibit novel insidethe-lightcone components, which induce exponentially growing, semi-static residual fields in the wake of a signal. The on-the-lightcone contribution to the signal is modified as well, as a result of the perturbed phase velocity. Both of these effects provide signatures of the axion background in the path of an electromagnetic wave. In this section, we provide a few examples to show how simple monochromatic signals can be distorted on timescales T that are shorter than the period of axion oscillation, T m −1 a .

Phase Velocity
One distinctive feature of the modified electrodynamics, the helicity-dependent phase velocities, can be quantified directly from the equations of motion as in Refs. [25,26]. From the dispersion relations for right-and left-polarized light, Eq. (10), the phase velocities can be expanded in powers This result can be also be recovered from the Green function in Eq. (30). If the signal is driven by a monochromatic source, the solutions for A ± (t, z) are proportional to trigonometric factors of exp(iΩ(t − z) ± µz), reproducing the linear order term in Eq. (47). Extremely sensitive measurements of the phase shift induced by gravitational waves are the bedrock for the remarkable recent detections of black hole and neutron star mergers. The sensitivity of Advanced LIGO [42] to the gravitational strain h = ∆L/L exceeds 10 −23 / √ Hz for gravitational waves with frequency f ∼ 10 2 Hz. An axion interferometer comparing the phases of left-and rightpolarized laser beams of angular frequency Ω observes a phase difference after the laser propagates a length L. Compared to the equivalent phase shift corresponding to a change in the path length ∆L, ∆φ = Ω∆L, the technology capable of detecting an h min ∼ 10 −23 would also be able to set a limit on µ/Ω of order in the context of an axion interferometer. To use the 1064 nm laser of LIGO [43] as an example (Ω = 1.77 × 10 15 rad/s), an interferometer capable of similar precision has a potential sensitivity to anẏ θ 2 × 10 −8 Hz, which overlaps with the parameter space indicated by Eq. (7) for 10 −11 GeV −1 g aγγ 10 −10 GeV −1 .
To achieve the sensitivity indicated by Eq. (50), the approximation thatθ(t) = 2µ is constant needs to hold only for as long as it takes the photon beam to traverse the interferometer. As long as m a L 1 and m a T 1 for the characteristic length L and time-of-flight T for the measurement, it is not necessary to insist that m a θ . A detector with cT ∼ nL ∼ O(10 3 km) (where n is the effective number of reflections of the light within the chamber) would thus be sensitive to m a 10 −12 eV, while more massive axions would be more easily detected with shorter interferometers. In the case of LIGO, the beam cavity storage time T is long compared to the length of each arm, with T ∼ nL/c for some n ≈ 70 [43].

After-Pulse Residual Fields
A substantively new effect appears when we consider the Green function to all orders inθ, for signals of finite duration T . In theθ = 0 vacuum of standard electrodynamics, such a signal propagates away from the source at speed c, maintaining the same duration T . In theθ = 0 background, this is no longer true. Rather than returning to zero after the signal passes, the E and B fields retain a residual nonzero value that grows with time. Interpreting − 1 4θ 2 as a tachyonic mass term for the photon, the growth of the E and B fields is a consequence of the tachyonic instability-a transfer of energy from the axion background into long-wavelength photons triggered by the original signal. This growth takes place inside the lightcone, rather than strictly on it.
Again specializing to a monochromatic signal j ± (τ ) = 2e iΩτ f ± (τ ) for simplicity, and taking the source at z = 0 to satisfy j ± (τ ) = 0 for τ < 0 and τ > T , the "residual" field A (r) ± at z > 0 refers to the nonzero field value after t − z > T . We have In the high frequency limit, the rapid oscillations of e iΩτ tend to cancel out the contributions from both integrals, so that the strengths of the residual E and B fields are proportional to µ/Ω. However, especially for µT ∼ O(1), the residual fields can become substantial: for λ O(1), the exponential growth of the Bessel functions becomes apparent, and eventually compensates for the µ/Ω suppression.
To demonstrate the distortion to a signal as it propagates through space, Figure 2 shows the Poynting vector S = E × B as a function of z at four snapshots in time. At t = T 0 , when the signal is newly produced, it exhibits relatively mild modifications to the original S z (t, z) = S 0 square wave, in this example with T 0 = µ −1 . By the t = 5T 0 snapshot, not only has the signal become notably distorted, but it also has developed nonzero values inside the lightcone, of magnitude S 0 , and continues to grow exponentially for t > 7T 0 .
In Figure 3, we show E x (t) for two square pulses with different durations T , to demonstrate the relationship between the pulse length and the strength of the residual fields. In this example we inflate the typical value of µ to µ = 1 kHz so that the strength of the residual pulse approaches the amplitude of the original signal for t ∼ O(few) × T in an example with 10 kHz radiation. However, the same plots can be reused for any value of µ by rescaling Ω, t and z so as to keep µΩ, µt and µz constant, as in the right hand panel of Figure 3. For another example, using µ = 10 −8 Hz instead of 10 3 Hz, the plot in Figure 3 would show t in units of 10 5 s rather than µs, and z = 1.5 × 10 13 km = 0.49 pc.

Exponential Growth
At sufficiently late times, µt 1, the growth in the E and B fields highlighted in Figure 3 becomes exponential, driven by the Bessel functions I 0 ( √ λ) and λ −1/2 I 1 ( √ λ) in the Green functions where the solid and dashed lines correspond respectively to µ = ±10 −3 Hz. For 0 ≤ t − z ≤ T 0 ("on the lightcone"), the difference in sign affects the phase velocities and spectrum of the pulses, but for t − z > T 0 ("inside the lightcone") the sign of µ is irrelevant for S z (t, z).
Eqs. (51-52). This is true not only of the semi-static residual fields in the wake of the signal, but also for the signal itself, if its duration T is long compared to µ −1 . For the phenomenologically relevant values of µ given by Eq. (7), probing the strongly exponential behavior requires extremely long coherence times T 0 > µ −1 10 8 s for both the axion field and the radiation source, unless the axion density is significantly enhanced beyond the expected ρ a ∼ 0.4 GeV/cm 3 . Nevertheless, it is a useful exercise to explore the behavior of the fields in this extreme limit.
Even in the late-time limit µ(t − z) > 1, there is a clear distinction between the received "signal", when 0 ≤ t − z ≤ T , and the "residual" fields, t − z > T . In contrast to the relatively mild modifications to the signal in Figure 3, the signal eventually becomes significantly amplified and distorted for large µt 1 and µz 1. The semi-static residual fields that appear in the wake of the signal continue to grow exponentially, until the point where the energy density in the electromagnetic fields becomes comparable to ρ a , and the back-reaction on the axion field can no longer be neglected. Figure 3 (right panel) shows the modified signal for an example with µ = 10 −3 Hz, where the pulse duration (T = 5 × 10 3 s) and propagation distance (z = 6 × 10 8 km) are both larger than µ −1 . In the µ → 0 limit of standard electrodynamics, the power density |S z | of this circularly-polarized square wave would remain constant, S m=0 z (0 ≤ t − z ≤ T ) = S 0 , returning to zero for t − z > T . Instead, the received power in the µ = 0 case varies as a function of time, oscillating with everlarger fluctuations and increasing exponentially. Once µ(t − z) 1, the fluctuations in the power ∆S z become larger than the magnitude of the power at the source, S 0 , and the exponential growth soon ensures that S z (t) S 0 for all t − z µ −1 . After the signal has passed by, the residual fields are well described by Eqs. (51-52), which in To show the dependence of the strength of the residual fields on the pulse length ΩT 0 , the left panel of Figure 4 shows S z (t, z) at t = 10 ms, z = 1 ms, as a function of pulse duration T 0 , with fixed frequency f 0 = Ω/2π = 5 kHz. Once T 0 4µ −1 , the strength of the field at (t, z) approaches a constant value. Evidently, S z (t − z > T ) is driven primarily by the first part of the signal, which has induced the longest-lived instabilities: the subsequent exponential growth of the initial instability makes this part of the signal the most consequential.
By considering the high-frequency limit of Eqs. (51-52), where f ± (τ ) and I j ( √ λ ) both vary slowly compared to e iΩτ , it is easy to see that the inside-the-lightcone electromagnetic fields scale as |E| ∝ Ω −1 and |B| ∝ Ω −1 , so that the radiation density is proportional to S z ∝ Ω −2 . In the right panel of Figure 4, the power S z (t, z) is shown for the same values of t and z, this time as a function of the frequency of the square pulse (f 0 = Ω/2π), with fixed pulse duration T 0 = 2 · 10 −3 s. At low frequencies, f 0 T −1 0 , the power approaches a constant. In this limit the phases of the integrands in Eqs. (51-52) are essentially independent of Ω, e iΩτ ±iµz ≈ e ±iµz . In the opposite limit, f 0 T −1 0 , the S z ∝ µ 2 /Ω 2 scaling suggested by Eqs. (51-52) becomes manifest. Together, the T 0 µ −1 and f 0 µ limits indicate that the power scales roughly as S z (t, z) ∼ e with t − z T 0 . In many applications involving radio, visible or X-ray radiation, the hierarchy between µ and Ω imposes an extreme penalty on the magnitude of the surplus power, so that only after multiple e-foldings would it be possible to detect the signal.
Our derivation of the Green function Eq. (30) assumes thatθ vanishes, so that µ(t) can be treated as constant. For a standard axion oscillating in a quadratic potential, µ(t) varies according to At times T comparable to the period of axion oscillation, T ∼ m −1 a , the steady-state approximation µ(t) µ is no longer valid, and it is necessary to use the methods described in Section 3. In principle, because µ 0 is set directly by ρ a and g aγγ , it is independent of m a , and the exponential growth in the late-time limit µT 1 can be explored without violating the steady-state condition m a T 1. Achieving µ 0 m a does, however, require a non-standard ALP model with an extended field range. Eq. (54) implies so µ 0 m a necessarily implies that θ 0 2.
In practice, this region of parameter space with µ 0 > m a is mostly ruled out by experiment. Based on the constraints for g aγγ [10][11][12][13][14][15] and m a [33][34][35][36], for dark matter densities of ρ a ≈ 0.4 GeV/cm 3 , Eq. (7) indicates thatθ m a for all allowed m a and g aγγ . To realize the late-time exponential growth within the steady-state axion background, three conditions must be satisfied: m a must be near the low-mass "fuzzy dark matter" extreme, m a ∼ 10 −21 eV; the coupling g aγγ must be relatively strong, g aγγ 10 −11 eV; and the value of √ ρ a in the path of the photon must be enhanced by a few orders of magnitude, for example by concentrating some fraction of the axions into dense clumps or "axion stars." If these conditions are not satisfied, then m a T 1 implies µ 0 T 1, and the Green function is well approximated by its series expansion in λ.

Oscillating Axion Background
In an oscillating axion background, the Green function for the vector potential A ± in (1 + 1) dimensional spacetime satisfies As in Section 2.1 we restrict our analysis to propagating plane waves. Unlike theθ ≈ const. limit of Eq. (19), Eq. (57) cannot be inverted to find an algebraic expression for the Fourier transform of g ± . Instead, couples g ± (ω) to g ± (ω ± m a ). This complication stems from the fact that the solutions of the homogeneous equations of motion are Mathieu functions. In both of the limitsθ 0 T 1 andθ 0 T 1, the leading forms of the Green functions can be extracted from Eq. (58) without invoking Mathieu functions or their Fourier transforms. Section 3.1 focuses on the former limit, in which frequencies ω = 1 2 m a ± O(θ 0 ) are resonantly enhanced, inducing exponential growth for these unstable frequencies. In the alternate limitθ 0 T 1 where the exponential growth is not realized, and forθ 0 T 1 for signals that do not include the resonantly enhanced frequencies, Section 3.2 provides a continued fraction expression for the Green function that converges quickly for smallθ 0 m a .

Resonantly Enhanced Propagation
For smallθ 0 m a , Eq. (58) suggests that the solution for g ± (ω) could be found as an expansion in kθ 0 /(k 2 − ω 2 ). However, for ω ≈ k ± O(θ 0 ) the simple perturbative expansion is disrupted, especially near frequencies ω = ± 1 2 m a where ω 2 = (ω ± m a ) 2 . At late times, whenθ 0 T 1, the exponential growth of the unstable modes with ω = ±m a /2 + O(θ 0 ) dominates the propagation of a signal. In this late-time limit the Green function can be approximated by integrating over the resonantly enhanced modes, 2 where we expand ω and k about ±m a /2, defining The factor of ±θ 0 in Eq. (57) corresponding to right-and left-polarized light is absorbed into the definition of the parameter . Even in the corner of ALP parameter space with small masses and relatively large couplings, Eq. (7) indicates that 10 −1 remains perturbatively small. For ALPs more closely resembling a QCD axion with m a > 10 −12 eV, g aγγ < 10 −11 GeV −1 , and fixed ρ a ≈ 0.4 GeV/cm 3 , the value of drops to < 10 −11 .
The function G(ω ± , k ± ) introduced in Eq. (59) is defined to include only the O( −1 ) part of g ± in the neighborhood of the resonant frequencies, where α and β are O(m a ). By dropping the O( 0 ) portion of the Green function, Eq. (58) can be disentangled to solve for G(ω ± , k ± ), In terms of α and β, Eq. (59) reduces to As always, the retarded Green function is defined to satisfy g(t < t 0 ) = 0, so the contour in ω(β) passes above both poles at β = ± α 2 − m 2 a , even when α 2 < m 2 a . Fortuitously, this integral is nearly identical to the one encountered in Section 2. Following the example of Eqs. (22)(23), we introduce the coordinate transformations α ≡ m a cosh ϕ, As in Figure 1, ϕ is such that cosh ϕ runs smoothly from −∞ to +∞ with Im cosh(ϕ) = 0. Extending the limits of integration in Eq. (62) to −∞ < β < ∞ and −∞ < α < ∞, we find which is again causal and exhibits propagation inside the lightcone. The 0 term in the expansion is calculated separately, by considering the Maxwell theory → 0 limit. In the late-time limit, when T = (t − t 0 ) satisfies 2 m a T + log log m a T, the I 0 ( √ λ) and I 1 ( √ λ) contributions become larger than O(1), and the resonantly enhanced modes dominate the Green function.
Compared to the m a → 0 limit from Section 2, the approximate retarded Green function in the resonantly enhanced regime of m a θ is remarkably similar. The exponential growth scale µ(t) has been replaced by which is in line with what we naively expect from Section 2. For example, if we used the steadystate result to approximate the late-time exponential growth by replacing |µ(t)| with its average value, |µ(t)| = 1 πθ 0 , the resulting estimate for the growth factor is off by only 27%. Recall from Eq. (7) that for fixed axion density ρ a = 0.4 GeV/cm 3 , even in the corner of parameter space saturating g aγγ 10 −10 GeV −1 and m a 2 · 10 −21 eV, the value of is still perturbatively small, | | 2 · 10 −3 . For this roughly-maximal value of , Eq. (66) is satisfied by m a T 3600. Elsewhere in the (m a , g aγγ ) parameter space, can assume significantly smaller values, requiring larger m a T 10 3 to satisfy Eq. (66). Our treatment of the axion background as a coherently oscillating field requires the field to remain coherent throughout the signal propagation, or T < T c with T c ∼ (m a v 2 ) −1 , where v ∼ 10 −3 is the virial velocity of the axions. In the case of photons traveling freely through space, as opposed to reflecting within some cavity, the propagation distance L must also be smaller than some L c ∼ (m a v) −1 . For 10 −5 , the onset of exponential growth indicated by Eq. (66) requires m a T > 10 6 ∼ 1/v 2 , meaning that decoherence effects become important on the timescales associated with the exponential growth, and must be accounted for. This result is consistent with Ref. [40], which found that for axion models with 10 −8 eV < m a , decoherence completely obscures the exponential growth. Additionally, for these models the width is narrow enough that gravitational redshift by the dark matter halo is sufficient to detune the resonance, even if the velocity for the axion cloud is taken to be v 10 −3 in order to satisfy m a T > 1/v 2 , leading the authors of Ref. [40] to conclude that for m a > 10 −8 eV in the observable range the parametric resonance at ω = 1 2 m a never develops into exponential growth. Nevertheless, there is a window where 10 −5 occupying a couple decades of the (m a , g aγγ ) parameter space in which the resonance can develop. It applies to extremely low-frequency radiation of ω 10 −21 eV, or equivalently f 10 −7 Hz ∼ 10 yr −1 , which does not propagate through the interstellar medium.

Propagation Without Resonance
In the opposite limit to Eq. (66), where m a T ∼θT 1 and the resonance is not given time to grow, an alternate approach provides a perturbative expansion of the Green function in powers of . This approach is also valid in the m a T > 1 limit for signals that do not include support within the instability band ω = m a /2 ± O(θ), and generically for any signal where the resonant component from Eq. (65) can be approximated by its series expansion.
For values of ω = m a /2, Eq. (58) can be rearranged into a continued fraction solution for g ± (ω), as a generalization of the relation in terms of from Eq. (60). By iterating the replacement of g ± (ω±m a ) with g ± (ω) and g ± (ω±2m a ), we derive an approximation of the form g ± (ω, k) ≈ a (ω, k) as follows: where the k dependence of each function g ± (ω ± 2 m a , k) has been left implicit, and where The recursion is provided by .
The convergence of the continued fraction expression effectively depends on a small-b expansion, and for generic values of ω ∼ k O( m a ) counting powers of is relatively easy: all of the a i are O( 0 ), while b ± ∼ O( 2 ). Near the k ≈ ω poles, where k = ±ω + O( m a ), a factor of (k 2 − ω 2 ) ∝ −1 modifies the power counting to b ± ∼ O( 2 −1 ).
However, when ω ≈ n 2 m a for integer n, the power counting is overturned near the k 2 = ω 2 poles, prompting the special treatment in Section 3.1. For example, for ω = 1 2 m a ± O( m a ), the poles in a 0 (ω) and a 0 (ω − m a ) are encountered simultaneously, and rather than finding a 1 (ω) = a 0 (ω)(1 + O( )), the difference between a 1 (ω) and a 0 (ω) becomes O( 0 ); similarly, b −1 (ω ≈ m a /2) ∼ O( 0 ). Excepting this ω = ± 1 2 m a resonance and the family of higher-order, narrower resonances, the Green function can otherwise be approximated to arbitrary order in by The = 0 case is simply the Maxwell theory → 0 result. At = 2, the expression for the retarded Green function at O( 2 ) is t + t 0 ))] cos m a z 2 + m a z sin m a z 2 which agrees at O( ) and even at O( 2 m a T ) and O( 2 m a z) with the late-time expression Eq. (65). Note that the series expansions of the Bessel I ν (z) functions are To express g ± in terms ofθ 0 , recall that a ± sign is incorporated into the definition of , where As a result, the replacement → − is equivalent to switching g + ↔ g − ; taking the complex conjugate, g ± = g ∓ ; or applying the parity transformation z → −z.
In Appendix A we provide the O( 3 ) form of g(t|t 0 , z). We also demonstrate that as m a T 1 and m a z 1 approach the late-time limit, the leading ( m a T ) n and ( m a z) n terms reconstruct the series expansions of I 0 ( √ λ) and I 1 ( √ λ)/ √ λ. In Eq. (78) we verify this explicitly as far as the 1/192 coefficient of the ( m a T ) 5 term. This indicates that the continued fraction expression for the Green function provides a smooth interpolation between the m a T 1 and m a T 1 limits.
By including the O( 2 ) terms, Eq. (73) is more precise than Eq. (65) in the m a t 1 limit, and it includes a novel effect: the −4 2 term in the expansion, which is not proportional to any sinusoidal factors. When the Green function is convolved with a signal of some duration T and some spectrum of frequencies Ω, in the T m −1 a limit the sinusoidal terms act as approximate Dirac δ functions to enhance the modes with Ω ≈ n 2 m a for integers n ≥ 1. If the axion mass is heavy enough that m a coincides with observable frequencies of light, then this resonant enhancement may be the most easily visible effect. However, for very light ALP dark matter where m a ∼ ω < 2π · O(kHz) corresponds to difficult-to-detect radio waves, the frequency-independent perturbation to the Green function becomes much more significant. As 2 is proportional to ρ a , changes in the axion density can modify the strength of visible light passing through it, causing "nongravitational microlensing." Depending on the cosmological history, some fraction of the axions can clump together to form minihalos with density perturbations δρ a /ρ a potentially much larger than O(1) [45]. (For recent work, see [46][47][48][49][50].) In addition to the gravitational microlensing, the direct effect from axion electrodynamics on starlight passing through an axion cluster may be detectable if the average value of 2 = 1 16 g 2 aγγ ρ a /m 2 a is not too small.

Numeric Results
In Sections 3.1 and 3.2 we used different methods to approximate the Green function in the homogeneous, oscillating axion background. In this section we verify Eq. (65) and Eq. (73) by comparing them to the numeric solution of the differential equation (76) (Recall that ≡ ±θ 0 4ma .) For the numeric calculations, we use a small value of σ < m −1 a to approximate the delta function source.
In Figure 5, we compare the series expansion of the continued fraction g ± (ω, k) = a (ω, k) at = 2 to the numeric result calculated from Eq. (76) with σ = 0.05m −1 a . Rather than truncating the Re(g ± (t, z)) with ✏ = +1/16 Im(g ± (t, z)) with ✏ = +1/16 Im(g ± ) m a z Figure 5: Cubic approximation for Re(g ± ) and Im(g ± ). Im(g ± (t, z)) with ✏ = +1/16 Im(g ± ) m a z Figure 5: Cubic approximation for Re(g ± ) and Im(g ± ). Im(g ± (t, z)) with ✏ = +1/16 Im(g ± ) m a z Figure 6: Cubic approximation for Re(g ± ) and Im(g ± ). 4θ 0 /m a → 0 result at t = 9m −1 a and t = 45m −1 a , in the upper and lower panels (respectively). Right column: the series expansion (dashed) and numeric result (solid) for Im(g ± (t|t 0 , z)) are shown as functions of z for the same fixed values of t as in the left column, at early times t = {5m −1 a , 9m −1 a } (top) and late times t = {25m −1 a , 45m −1 a } (bottom). In the → 0 limit (ordinary Maxwell theory), Im(g ± ) → 0. Changing the sign of is equivalent to replacing g ± with its complex conjugate, g ± → g ± , or performing the parity transformation z → −z. In this example we take = +1/16 and t 0 = 0, so that the top and bottom rows represent "early" and "late" times, with m a t 1 and m a t 1, respectively. The t = 45m −1 a curve in red shows the breakdown of the series expansion for m a t 1. Re(g ± (t, z)) with ✏ = +1/16 Im(g ± (t, z)) with ✏ = +1/16 Im(g ± ) m a z Figure 7: Cubic approximation for Re(g ± ) and Im(g ± ). Im(g ± (t, z)) with ✏ = +1/16 Im(g ± ) m a z Figure 7: Cubic approximation for Re(g ± ) and Im(g ± ). For simplicity, we use t 0 = 0 in this example, with = 1/16. The agreement at early times, m a t 1, is quite good, but the approximation begins to fail by t 2m −1 a / . Even though the O( 3 ) expression includes terms that grow as ( m a T ) and i ( m a T ) 2 , by t ≥ 45m −1 a the exponential growth has begun to invalidate the series expansion in ( m a T ) for both the real and imaginary parts of the Green function. At extremely late times, when m a T −2 , even our "late-time" approximation from Section 3.1 fails to capture the dominant behavior of the Green function. In addition to the ω = ±m a /2 poles, the contribution from the next-to-leading resonance at ω = ±m a +O( 2 m a ) becomes significant. To calculate the extremely-late Green function, the methods of Section 3.1 can be repeated with some new k ± = ±m a + 2 γ and ω ± = ±m a + 2 δ. Together with the subleading terms from the ω = ±m a /2 poles, the resulting expression would include all terms of O( 2 (m a T ) n ), and should match the n = 0 and n = 1 terms that appear already in Eq. (77). For the phenomenologically viable values of with ρ a ∼ 0.4 GeV/cm 3 , the decoherence at T c ∼ (m a v 2 ) −1 ensures that 2 m a T O(1) for all T T c , and the axion field loses coherence before the behavior at these "extremely late" times can be explored.

Conclusions
In this paper we have computed a set of Green functions appropriate to various limits of axion electrodynamics. Between Section 2.1, Section 3.1, and Section 3.2, our analysis covers the phenomenologically viable parameter space for coherent axion backgrounds, as well as more extreme ALP models with stronger couplings or enhanced energy densities.
Our study of photon propagation at early times T m −1 a in Section 2 shows that the QCD axion induces gentle modifications to the propagation of local disturbances in standard electrodynamics. ALP models in the more extreme corners of parameter space can instigate more dramatic growth in low-frequency modes as signal pulses pass through space. In Section 3, extending our analysis to account for the oscillation of the axion field, we derive a Green function that preferentially enhances radiation with frequencies close to ω = n ma 2 for positive integers n. In all cases, perhaps the sharpest qualitative distinction between radiation in Maxwell theory and in axion electrodynamics is the presence of inside-the-lightcone propagation in the latter.
We collect our main results below. For plane waves symmetric in the x and y directions, the Green functions g ± satisfy the equations of motion ( for right-and left-polarized light, respectively. • Semi-static Limit at times t m −1 a , with µ = 1 2θ 0 cos(m a t) const: The phase factor e ±iµz provides polarization-dependent phase velocities, as anticipated by the dispersion relation Eq. (13). The exponentially growing I 0 function corresponds to superluminal group velocities for both polarizations of light. For nonstandard ALP models with m a θ 0 , the exponential growth at t µ −1 can begin while the semi-static condition t m −1 a is still satisfied. Otherwise, when µt 1, the Bessel function is well approximated by its series expansion in powers of µt and µz.
• Resonant, Late-Time Limit for | |m a t 1, where m a = ± 1 4θ 0 for g ± , respectively: where g ± is expressed as a function of λ = 2 m 2 a ((t − t 0 ) 2 − z 2 ). This Green function describes the behavior at late times, λ 1, when the resonant enhancement of frequencies ω m a /2 is the dominant effect. In contrast to Eq. (30), Eq. (65) is appropriate for the standard ALP models withθ 0 m a , and is valid even for t m −1 a . Both Eq. (30) and Eq. (65) exhibit exponential growth at late times t θ −1 , though in the latter case the enhancement is specific to frequency modes within the narrow resonance ω = m a /2 ± O(θ).
At extremely late times, when 2 m a t 1, the contributions from the higher resonances such as ω = m a ± O(θ 2 /m a ) also become significant. Given the small values of for most regions of ALP parameter space, and the fact that the axion background usually exhibits decoherence before 2 m a t ∼ O(1) is satisfied, we do not provide the Green function in this limit. However, the methods of Section 3.1 can be extended in a straightforward manner to describe the dominant behavior in the 2 m a t 1 limit.
• Non-Resonant Propagation for | |m a t 1, and for signals that do not include the resonant frequencies ω = m a /2 + O(θ 0 ): g ± (t|t 0 , z) = Θ(t − t 0 )Θ((t − t 0 ) 2 − z 2 ) 2 1 + 4i cos m a (t + t 0 ) 2 sin m a z 2 − 4 2 + 2 4 cos(m a (t + t 0 )) cos(m a z) + 2m a (t − t 0 ) cos m a z 2 sin m a (t − t 0 ) 2 − 2 cos m a (t − t 0 ) 2 [−2 + 2 cos (m a (t + t 0 ))] cos m a z 2 + m a z sin m a z 2 where the O( 3 ) term of the non-resonant series expansion is provided in Eq. (77) in Appendix A. As long as | |m a t < 1 and 1, the Green function can be calculated to arbitrary precision in O( n ) using this approach. For | |m a t 1, this series expansion interpolates smoothly onto the Bessel functions in Eq. (65), which include all terms of order ( m a t) n and ( m a z) n for n = 0, 1, 2 . . ..
In addition to the enhancement of specific frequencies ω ≈ n 2 m a , Eq. (73) and Eq. (77) include terms that are not sinusoidal, and which modify the propagation of all frequencies of light. For ALP models with small axion masses, where the resonant frequencies themselves are too small to detect, the frequency independent (1 − 4 2 ) part of the Green function continues to affect visible wavelengths of light, and may provide a new signal of ALP dark matter in regions where the axion density ρ a ∝ 2 varies significantly.
Both analyses in Section 3 are predicated upon = ± 1 4θ 0 /m a being a small parameter, and so the > O(1) case is generally the most difficult to address. If > 1 then the continued fraction expression for the Green function does not converge. In terms of the resonant analysis in Section 3.1, once n m a t 1, the Green function receives contributions from an infinite set of poles at ω n = n 2 m a . This unusual m a θ limit is handled by the treatment in Section 2.1 for times t m −1 a , where the semi-static approximation is valid. However, once t m −1 a , the > 1 Green function must be written in terms of Mathieu functions or calculated numerically.
If axions make up some component of the dark matter, then their presence may be discerned through their influence on propagating photons. Sensitive terrestrial experiments designed to measure the inside-the-lightcone propagation or the polarization-dependent perturbations from axion electrodynamics might complement existing detection strategies for ALP dark matter. Astrophysical observations may also be sensitive to the effects of axions. Depending on the axion mass, the resonant enhancement and frequency-independent modifications could provide additional opportunities to detect axion dark matter, especially if some fraction of the dark matter has collapsed into minihalos.
The work presented here can be developed in several directions. To fully understand the propagation of light through a galaxy containing axions, especially in the neighborhood of axion minihalos, it is important to consider the effects from spatial gradients in the axion background. For light propagating through multiple coherent patches, over distances L (m a v) −1 or for times T (m a v 2 ) −1 , decoherence effects are similarly important. Finally, Green function techniques may provide useful tools for exploring the sensitivity of terrestrial axion detection experiments. We hope to return to these issues in future work.
with the result One of the two infinite series can be replaced by a hypergeometric function, after replacing the index k with ≡ k − j, so that both j and run from zero to infinity: Even for relatively large values of µt and µr, the series expression for I converges relatively quickly. In the asymptotic v 0 1 limit, the 1 F 2 function is approximately while to leading order the Bessel functions I j (2 √ σ) approach Consequently, the number of terms in jmax j=0 required for convergence is driven primarily by the value of µt: in the limit µt 1, the series converges quickly for j > j max once 3j max 2 log j max e + j max log µt √ 2 − 5j max 2 log µr 2 µt.