B-modes from Post-inflationary Gravitational Waves Sourced by Axionic Instabilities at Cosmic Reionization

We show that axion-like particles that only couple to invisible dark photons can generate visible B-mode signals around the reionization epoch. The axion field starts rolling shortly before reionization, resulting in a tachyonic instability for the dark photons. This generates an exponential growth of the dark photon quanta sourcing both scalar metric modes and gravitational waves that leave an imprint on the reionized baryons. The tensor modes modify the cosmic microwave background (CMB) polarization at reionization, generating visible B-mode signatures for the next generation of CMB experiments for parameter ranges that satisfy the current experimental constraints.


I. INTRODUCTION
The discovery of gravitational waves (GW) at LIGO [1] and VIRGO [2] has motivated the search for other possible GW sources beyond the mergers of astrophysical objects. Among those, GWs from cosmological sources, such as strong first order phase transitions [3] and the presence of cosmic strings [4], are of particular interest in elucidating the early history of the universe (e.g., [5,6]). The cosmological GW signals can have a wide range of possible frequencies: interferometer experiments can detect GWs with frequencies above ∼ 10 −5 Hz [7][8][9][10], and lower frequency signals down to ∼ 10 −8 Hz are relevant for pulsar timing experiments [11,12]; if GWs have frequencies lower than ∼ 10 −15 Hz, we can search for the B-mode polarization signals from GW imprints on the cosmic microwave background (CMB) [13]. Such lowfrequency signals have wavelengths comparable to the visible universe's size and must have a cosmological origin. As a result, the B-mode signal is mainly considered to come from GWs produced during cosmic inflation (see [14] and the references therein).
In this letter, we propose a new source for B-mode generating GWs produced by axion-like particles (ALPs) around the time of cosmic reionization. Axions were originally proposed to solve the strong CP problem [15,16] and realized to be a viable dark matter (DM) candidate [17][18][19][20]. ALPs generalize the cosmological phenomenology of axions without a necessary connection to strong CP. For example, an ALP can serve as the inflaton field responsible for the period of cosmic inflation [21][22][23] or as the relaxion, addressing the hierarchy problems in nature by varying the fundamental constants of nature with time [24][25][26]. On the experimental side, several new direct detection experiments have been put into action [27][28][29][30] or have been proposed [31][32][33][34] to look for ALPs. Part of the theoretically-favored axion parameter space has already been experimentally excluded.
In the particular case where ALPs couple to dark photons, the ALP field's rolling leads to a "tachyonic instability" that amplifies vacuum fluctuations of one of the dark photon helicities. The process generates exponential dark photon production, and similar phenomena have been studied under the context of inflation [23,35], production of dark photon DM [36], depletion of axion DM to avoid overclosure [37], and friction for the relaxion models [24,25]. Recently it has been shown that the stochastic GW background generated through this process in the early universe may be detectable in interferometers or pulsar timing arrays [38,39]. In [40], a similar mechanism at the recombination period is studied within the context of early dark energy solutions to the Hubble tension [41,42] and is shown to produce visible GW signals in the CMB.
In this work, we consider a similar effect of producing a GW background from tachyonic particle production late in the universe's history -after recombination and around the time of galaxy formation. As a tensor perturbation of the metric, the GW background leaves an imprint on the photon energy distribution. When the universe enters the reionization era at z rei ≈ 8, CMB photons propagating in the line-of-sight direction get polarized by the last Thomson scattering, and a combination of the tensor perturbation and the angular distribution of the photon polarization produces the B-mode signal in the large-scale CMB spectrum. In particular, we will show that for parameter ranges of our model not currently excluded by existing or past experiments [43,44], we predict a B-mode signal accesible to the next generation of B-mode detectors.
The B-mode signals sourced by the axionic instability have a power spectrum which could be distinguished from those produced by inflationary GWs. An observation of such unique B-mode signals will be a discovery of dark sector physics and will shed light on the nature of dark energy by revealing that dark energy is changing at late times. In particular, a revelation that dark energy has recently changed by an amount close to its current value is suggestive of some dynamics related to the cosmological constant (CC) problem. As we will show, next-generation experiments will be able to probe such shifts on a scale similar to the current value of the cosmological constant.
We remark that our calculation utilizes a linear semiclassical approximation and therefore our results need to be confirmed by a full lattice study. We expect this to affect the precise predictions of the spectral shapes, but not our ultimate conclusions. This paper is organised as follows: After reviewing the mechanism of tachyonic production of dark photons, we describe our setting and set up the calculation of dark photon's energy density fluctuations in Sec. II. We then discuss the metric perturbation sourced by the dark photon fluctuation in Sec. III and show the derivation of the resulting CMB spectra. Subsequently, we present our results, comparing the predicted signals within two benchmark ALP models to the sensitivity of the future B-mode experiments and to the current constraints from Planck in Sec. IV. Finally, we conclude in Sec. V.

II.1. Tachyonic production of dark photons
We consider an axion field φ coupled to a U (1) dark photon, with the Lagrangian given by where V (φ) = 1 2 m 2 φ 2 , and f is the axion constant. We assume the dark photon is massless and is produced only after inflation. The quadratic potential V (φ) can naturally arise from an axion-like potential Λ 4 cos(φ/f ), which implies m ∼ Λ 2 /f . We consider m close to the Hubble scale right before the reionization. We will see that in our setting, the CMB probes Λ ∼ O(meV), which also coincides with the order of magnitude of the cosmological constant, so that an observation of the signal we discuss may lead to new insights into dark energy 1 .
The equation of motion of the axion field is then in which a is the scale factor of the FRW metric ds 2 = a 2 (τ )(dτ 2 − δ ij dx i dx j ), and H is the Hubble parameter. The prime symbol denotes derivatives with respect to the conformal time τ . On the right hand side of the equation, the dark electromagentic field serves as friction for the rolling of the axion φ.
The rolling of the axion will cause the dark photon modes within a certain momentum range to grow exponentially, a phenomenon known as the tachyonic instability. This can be shown by examining the equation of motion of the dark photon field, which in the Coulomb gauge is written as where Dk ≡ d 3 k/(2π) 3 .
The creation and annihilation operators obey the commutation relation , and the polarization vectors obey k · ± = 0, k × ± = ∓ik ± , ± · ± = 0, ± · ∓ = 1. The dark photon field equation can then be written in terms of the mode function v as with the dispersion relation ω 2 ± (k, τ ) = k 2 ∓ kαφ /f . As long as the axion starts rolling and develops a non-zero φ , the dark photon modes of a certain helicity in the momentum band 0 < k < α|φ |/f will have ω 2 (k, τ ) < 0 and therefore grow exponentially. Specifically, the v + modes can grow when φ > 0 and the v − modes grow when φ < 0, and the growth of the two helicities are alternating as the axion field oscillates around the minimum of its potential. The helicity experiencing the tachyonic instability right after axion starts the rolling will be significantly more enhanced than the other, as it spends more time in the tachyonic band.
To solve the axion and the dark photon coupled equations of motion, we treat the dark photon mode functions v ± (k, τ ) as discretized modes of fixed k. And to the leading order, the reaction from dark photon field E·B on the right hand side of Eq. (2) is replaced by the expectation value E · B , which is calculated as

II.2. Calculation setup
In our calculation, we assume the dark photon to be non-thermal such that its abundance comes only from the tachyonic production described above. We use 200 dark photon k-modes equally spaced on a logarithmic grid in the momentum range [k min , k max ]. The value of k max is chosen such that k max α|φ | max /f , and we perform a consistency check with several choices of k max to determine the number used for each calculation. The minimum value of the momentum range is set to be k min = H 0 /4. With these choices, we make sure that the entirety of the momentum range of interest is covered.
In this work we do not include the back reaction of the gauge modes on the axion perturbations that requires a full lattice study (see [46]). This is mainly important for the axion abundance calculation which is not of interest in this setup. For the GWs, we expect the lattice results to be roughly consistent in magnitude [46] and to be mainly important for the spectral shape (see also Left panel: The evolution of the axion potential energy of the two benchmark models we use, normalized by the dark energy density today ρCC ≈ 37 meV 4 . Right panel: The evolution of the dark sector energy density ρe and its perturbation δρ 2 e 1/2 induced by the tachyonic particle production, normalized by the total energy density of the universe, [ [47][48][49]). We therefore treat the calculation here as a preliminary estimate to motivate a full lattice study which is left for a subsequent work. The produced dark photon k-modes are assumed to be in the Bunch-Davis vacuum v ± (k, τ ) = e ikτ / √ 2k before the axion rolling starts and the axion field is released at an initial misalignment of |φ 0 | = f . We choose two benchmark models for which the tachyonic instability becomes significant after recombination but before reionization, taken as z rei = 8. Note that keeping the axion mass fixed and varying the height of the initial misalignment (i.e. Λ in our parameterization, keeping f ∝ Λ 2 ) will only rescale the energy in the dark sector, and with it the energy density in the perturbations. Therefore we will think of the two benchmarks as two classes of models where the total energy in the dark sector remains a free parameter which can be constrained by current experimental data. We give the benchmark values of the parameters in Table I, where we also show bounds on the energy scale in the dark sector Λ bound which we derive in the following sections based on the uncertainty of the current power spectra measurements.
We plot the time evolution of the axion potential and the resulting inhomogeneities in the gauge modes in Fig. 1, using the two benchmark mark models and energy scale Λ bound . We see that the axions start their rolling shortly before the reionization and begin to oscillate around the minimum, producing the gauge mode inhomogeneity in the process. In Fig. 1 (right), we see that the dark photon energy is always below O(5 × 10 −4 ) of the total energy. The process therefore gives negligible corrections to the angular diameter distance that relates to the CMB spectra. However, even though the average ρ e is small comparing to ρ tot that is dominated by the matter density ρ m , the density contrast of the dark photon energy is of O(1) as can be seen in the transparent and opaque curves. The energy perturbation δρ 2 e 1/2 is thus comparable to the matter density perturbation (∼ 10 −5 ρ m ) that enters the horizon around the same time and can therefore generate visible signals as we show below.

III. CMB SPECTRA CALCULATION
Although the axion starts rolling only after recombination, remarkably it can still modify the CMB perturbation observed today. The dark photon field enhanced by the tachyonic instability generates isocurvature perturbations that also source GWs [38] affecting the CMB power spectrum through the late integrated Sachs-Wolfe (ISW) effect. The produced GWs also leave an imprint in the CMB B-mode which will serve as our target signal for the discovery of this setup. Here we present the calculation of CMB T T , EE and BB spectra, C T T , C EE and C BB .

III.1. Scalar mode contribution
Perturbations of the axion and dark photon energy density δρ e generate a gravitational potential Φ through the linear Boltzmann and Einstein equations [50] δ m + θ m = 3Φ , Here δ m = δρm ρm is the matter energy density contrast induced by the perturbation from the dark photons, and θ m is the velocity divergence of matter. For the metric perturbations we set Φ = −Ψ and ignore the stress tensor from the free streaming radiation. Once the tachyonic production starts, the dark photons dominate the energy perturbation of the dark sector, and hence: Here the energy density fluctuation is defined as an operator by subtracting the expectation value from the energy density operator [23]. Through the ISW effect, the gravity perturbation Φ, obtained by solving Eq. (6), sources a temperature perturbation today Θ 0 (n) = δT /T (−n; τ 0 ) as [51] where τ 0 and τ rec are the conformal time today and at recombination respectively. Since the dark photon perturbation from the tachyonic production is uncorrelated with the adiabatic perturbation, the cross correlator between Θ 0 (n) and the adiabatic CMB temperature perturbation is negligible. Therefore, the dark photon contribution to the temperature power spectrum is calculated as dn dn Θ 0 (n )Θ 0 (n )P l (n · n ) . (10) Using functions T r l and T i l defined in Eq. (A6) and (A7) as convolution integrals between the dark photon mode function v(k, τ ) and the spherical Bessel functions, we find that where the vector k 2 = k − k 1 . We give more details of the derivation in the Appendix A.
The scalar perturbations can also source the CMB Emode, which can be calculated as after taking the narrow width approximation of the visibility function in time. Here τ rei is the conformal time at reionization, and T rei ≈ 0.08 is the photon optical depth in the reionized universe. We find that the E-mode contribution from the scalar perturbations is subdominant to that of the tensor perturbations.

III.2. Tensor mode contribution
The tensor perturbation h(k, τ ) is obtained from the linear Einstein equation, which is written in terms of h ij = ah ij as where Π ij (k, τ ) is the anisotropic part of the energy momentum tensor T ij . The tensor perturbation then generates the B-mode power spectrum as [51] C BB where with κ = (τ 0 − τ rei )k. We take the narrow width approximation of the visibility function in time as in the E-mode calculation. In contrast to the calculation of T T , the Bmode signal relies on having the last photon scattering at the reionization. This can be seen by the presence of j 2 [(τ rei − τ ) k] that comes from expanding the photon propagation within the time interval [τ, τ rei ] into spherical harmonics and then matching the angular mode to the polarization signal. As can be seen from Eq. (13), the spectrum h (k, τ )h (k, τ ) is related to Π ij (k, τ )Π ij (k , τ ) , which again can be expressed in terms of the dark photon mode function v(k, τ ). The B-mode spectrum can therefore be rewritten as The function Θ defined in Eq. (B5) comes from the scalar products of the dark photon polarization, while B r and B i are convolutions between v(k, τ ) and the spherical Bessel function j 2 (see Eq. (B6) and Eq. (B7)). More details of the derivation appear in Appendix B.
The tensor perturbation also contributes to the EE and T T spectrum. The E-modes have a similar generation mechanism as the B-modes, and we can calculate the C EE by simply replacing J B,l (k) in Eq. (14) with  Table. I).
In the left panel we also show the 1σ error bar of the binned Planck 2018 power spectrum [43] up to l = 2000 (when l > 2000 the uncertainty increases) and in the right panel -the measurement from BICEP2/Keck Array [13] (digitized from [52]) as well as the projected instrumental noise of several future experiments, including LiteBIRD [53], CMB-S4 [52], PICO [54] and CORE [55] (digitized from [56]). Additionally, we plot on the right panel the primordial and lensing B-mode spectra for several different tensor-to-scalar ratio r (dotted gray, taken from [53]). and the changing the pre-factor from 36π to 9π. The T T spectrum can be calculated by As we can see, the power spectra contributed by the tensor and scalar perturbations are proportional Π ij (k, τ )Π ij (k , τ ) , which scales with the mode function of dark photon as v 4 . The energy density of the dark photon field is ρ X = 1 2a 4 Dk |v (k)| 2 + k 2 |v(k)| 2 − k , where the last term comes from subtracting the vacuum energy [37]. When the mode function grows due to the tachyonic production, the axion initial potential energy ρ φ = 1 2 Λ 4 quickly transfers into ρ X and generates v ∝ Λ 2 . When fixing the axion mass m and the axion-dark photon coupling α, the magnitude of the resulting spectrum is proportional to Λ 8 .

IV. RESULTS: THE T T AND BB SPECTRA
In Fig. 2, we show the C T T and C BB spectra from the two benchmark models defined in Table. I. In particular, the value of Λ is rescaled (keeping m and α fixed) so that the T T spectrum roughly saturates the error bars from the Planck 2018 data [43], as can be seen in the plot. This shows the rough bounds on Λ from the current CMB measurements. On the right of Fig. 2 we plot the corresponding C BB signals for the two benchmark models saturating the C T T constraints. This gives the upper range of the predicted B-modes within our setting.
Below we discuss the shape of the calculated spectra. As we can see in Fig. 2 (right), the two axion B-mode curves are roughly parallel to each other. This can be explained by the spherical Bessel function from the angular integral that projects the photon polarization tensor to the B-mode perturbation. The function peaks at the origin and is suppressed by x −3 when x 1, so the integral is dominated by the k-modes that minimize x. At the same time, the tachyonic production mainly produces k-modes larger than τ −1 rei . This results in D BB getting most of its contribution from perturbations at τ ∼ τ rei . This explains why the difference in the dark photon production at early times between the two benchmarks does not significantly modify the shape of the spectra even though the axions in the two models start rolling at different times (as shown in Fig. 1).
This behavior does not hold, however, for the D T T spectra which are sensitive to the starting time of the particle production. The T T spectra in Eqs. (11) and (18) are not affected by the reionization and the spherical harmonic projection has the form j [(τ 0 − τ )k], receiving contributions from a wider τ window for different -modes. Our numerical results show that the T T spectra are dominantly contributed by the early period of the dark photon production. This is why they no-longer peak at lower -modes as D BB , and the peak of the spectrum for the BM1 model, where the particle production starts earlier, is accordingly at higher compared with the peak of the BM2 spectrum.
In the D T T plot, we compare signals from the two benchmark axion models to the Planck 2018 data [43], establishing a rough bound on Λ. We find that Λ bound ≈ 15(9) meV (see also Table I) for the BM1 (BM2) that saturates the error bar of the Planck data following the same binning as in [43]. We also find a similar sensitivity from the Planck E-mode polarization data, not shown here.
In the D BB plot, we first note that the BICEP2/Keck measurement [13] does not exclude the benchmark models 2 . We have accordingly chosen the parameters to satisfy the existing constraints and find that the signal from the late time tachyonic production is well within sensitivities of next-generation CMB B-mode experiments, such as LiteBIRD [53], CMB-S4 [52], PICO [54] and CORE [55].
The B-mode signals from axions peak at low-, similarly to those from inflationary tensor modes, in both cases due to reionization. In this region, the inflationary model with r ∼ 0.01 produces B-mode signals that dominate over the gravitational lensing signal (see e.g. Fig. 1 of [55]). This suggests that the axion signals can also dominate the lensing background. The scientific goal of LiteBIRD, for example, is to achieve an uncertainty of δr ∼ 0.001 on the range 2 100 [53]. It has been shown that even with the contamination from diffuse galactic foreground, LiteBIRD can still be sensitive to D BB l ∼ 10 −4 µK 2 [57] for 10. Such sensitivity is close to the BM1 signal, and it is also comparable to the BM2 signal even with a lower Λ ≈ 7 meV, which is close to the scale of the observed dark energy ρ 1/4 CC . This signal, if observed, might have intriguing implications for the nature of dark energy.

V. CONCLUSION
We have studied the CMB power spectra generated by ALPs via a tachyonic instability and the ensuing production of dark photon quanta close to the cosmic reionization epoch. The ALP-dark photon system produces GWs that leave an imprint in the CMB, including its B-mode polarization spectrum. The signal is visible to future CMB polarization detectors while remaining compatible with the bounds from current measurements. Moreover, we find that future experiments can be sensitive to ALP potential energies similar in order of magnitude to the value of CC, which, if discovered, may lead to progress in discerning the nature of dark energy. We note that our setting may potentially also generate a signal in measures of cosmic non-Gaussianity that could be visible to future experiments. We leave this analysis for a future study.

VI. ACKNOWLEDGEMENT
We thank Pedro Schwaller, Ben Stefanek, Chen Sun for useful discussions, and especially Gustavo Marques Tavares for useful comments to the draft. MG and SL are supported in part by Israel Science Foundation under Grant No. 1302/19. MG is also supported in part by the US-Israeli BSF grant 2018236 and the GIF grant I-2524-303.7. YT is supported by the NSF grant PHY-2014165. YT was also supported in part by the National Science Foundation under grant PHY-1914731 and by the Maryland Center for Fundamental Physics. The authors also thank the KITP institute (Enervac19 program), and the Munich Institute for Astro-and Particle Physics (MI-APP) of the DFG Excellence Cluster Origins, where part of this work was conducted, for hospitality.
where we have defined g(k, τ ) ≡ −4πG N a 2 δρ e for convenience. With the dτ integration in Eq. (9), we reorganize the expression by switching the order of the conformal time integrals as Θ l (k) = 2 τ0 τosc dτ G Φ (τ, τ )g(k, τ )j l [k(τ 0 − τ )] + 2 We can then convert the correlation function Θ l (k)Θ l (k ) into δρ e (k, τ )δρ e (k , τ ) . The δρ e operator can be expressed with the dark photon fields by replacing the X i and X ij in Eq. (7) with the definitions in Eq. (3), and its spectrum is obtained as 0|δ [O] δ [O] |0 = 0|O 2 |0 − 0|O|0 2 . After a lengthy but straightforward calculation one arrives at Eq. (11). The expression of the function T r and T i are