Audible axions

Conventional approaches to probing axions and axion-like particles (ALPs) typically rely on a coupling to photons. However, if this coupling is extremely weak, ALPs become invisible and are effectively decoupled from the Standard Model. Here we show that such invisible axions, which are viable candidates for dark matter, can produce a stochastic gravitational wave background in the early universe. This signal is generated in models where the invisible axion couples to a dark gauge boson that experiences a tachyonic instability when the axion begins to oscillate. Incidentally, the same mechanism also widens the viable parameter space for axion dark matter. Quantum fluctuations amplified by the exponentially growing gauge boson modes source chiral gravitational waves. For axion decay constants f ≳ 1017 GeV, this signal is detectable by either pulsar timing arrays or space/ground-based gravitational wave detectors for a broad range of axion masses, thus providing a new window to probe invisible axion models.


Introduction
Axions, or more generally, axion like particles (ALPs) appear in many extensions of the Standard Model (SM). This includes solutions to the strong CP problem via the Peccei-Quinn mechanism [1,2], string theory [3,4], natural models of inflation [5], dark matter (DM) [6][7][8], the relaxion mechanism for solving the hierarchy problem [9], or just in general models where an approximate global symmetry is spontaneously broken at a high scale, resulting in a very light pseudo Nambu-Goldstone boson. The viable parameter space for ALP masses and couplings spans many orders of magnitude, which makes searching for them challenging, but also motivates new ideas and approaches to probe previously JHEP01(2019)053 inaccessible regions. In particular, if the ALP is effectively decoupled from the SM, the only bound comes from black-hole superradiance [10].
Here we show that axions and ALPs may leave a trace of their presence in the early universe in the form of a stochastic gravitational wave (GW) background. This is possible if the axion couples to a dark photon. Initially, the axion is displaced from its minimum and Hubble friction prevents it from rolling until the Hubble rate drops below the axion mass. Once it begins to roll, the axion induces a tachyonic instability for one of the dark photon helicities, causing vacuum fluctuations to grow exponentially. This induces time-dependent anisotropic stress in the energy-momentum tensor, which ultimately sources gravitational waves. In the process, a large fraction of the energy density stored in the axion field is converted into radiation in the form of dark photons and gravitational waves.
Subsequently, a period of oscillation occurs where the axion undergoes parametric resonance, further suppressing the amount of energy stored in the axion field. Observable gravitational wave signals require that the energy density stored in the axion field at the time of GW emission is large, so the combination of the tachyonic instability and the subsequent parametric resonance is necessary to prevent overabundant axion dark matter unless the axion is heavy enough to decay. This mechanism for efficient depletion of axion DM abundance via exponential production of dark photons was first pointed out in ref. [11], where it was used to increase the QCD axion decay constant without tuning the initial conditions. It was also noted that the SM photon cannot play the role of the gauge boson because its fast thermalization rate would destroy the conditions required for exponential particle production.
The tachyonic instability of rolling ALPs has been previously exploited in a variety of contexts, e.g. inflationary models [12][13][14][15], the seeding of cosmological magnetic fields [16][17][18][19][20][21], reducing the relic abundance of the QCD axion [11,22], populating vector dark matter [23][24][25][26], friction for the relaxion mechanism [27,28], and GW from the string axiverse [29,30]. Here, we assume that the axion starts rolling sometime after the end of inflation when the universe is radiation dominated, as is the case in relaxion [27,28] or axion-curvaton models. The resulting GW signal will therefore be peaked, with the peak frequency set by the axion mass and the amplitude by the Hubble rate at the time when the backreaction from the tachyonic instability becomes large.
As shown in figure 2, the tail of the GW spectrum from the QCD axion may be detectable by pulsar timing arrays while generic ALP models are in reach of LISA and other future GW detectors. At the same time, the relic abundance of the axion can be a fraction or all of the observed DM in the universe. This paper is organized as follows: after describing the model and initial conditions in section 2, we review and discuss the particle production mechanism in section 3. Section 4 contains an outline of the GW spectrum computation and gives an analytic estimate for the peak frequency and amplitude. Plots of the GW spectra from our numerical simulation can be found in section 5.

JHEP01(2019)053 2 Model
Our model consists of an axion φ and a dark photon X µ of an unbroken U(1) X gauge symmetry under which the Standard Model fields are uncharged 1,2 where f is the scale of the global symmetry breaking that gives rise to the Nambu-Goldstone field φ. We assume a cosine-like potential for the axion with mass m however our results do not depend crucially on the precise form of the potential, allowing the mechanism to be extended to other types of rolling pseudoscalars. Adopting a metric convention of ds 2 = a(τ ) 2 (dτ 2 − δ ij dx i dx j ) and assuming that φ is spatially homogeneous, the equation of motion for the axion field is where primes indicate derivatives with respect to the conformal time τ , E and B are the dark electric and magnetic fields, and the Hubble parameter is defined as H = a /a 2 . As the axion rolls towards its minimum, the coupling φX µν X µν can lead to exponential production of X µ quanta. Different from most of the existing literature, here we assume that the dynamics takes place after inflation, in a radiation dominated epoch. More precisely, we assume the following initial conditions, which can naturally arise at the end of inflation: • φ is displaced from its minimum by ∆φ = θf with the initial misalignment angle θ = O(1).
• The energy density in φ is smaller than that of the radiation bath, such that its backreaction on the geometry can be ignored.
• The gauge field X µ is not thermalized and thus has zero initial abundance.
• The initial velocity φ is negligible.
The last two assumptions may be relaxed without spoiling our mechanism. While H m, the axion is pinned by Hubble friction and no gauge bosons are produced. At the temperature T osc defined by H osc ≈ m, the axion becomes free to roll toward the minimum of its potential. In a radiation dominated universe, the Hubble expansion rate is approximately H ≈ T 2 /M P , so the temperature at which the axion begins to oscillate is T osc ≈ √ mM P , where M P is the reduced Planck scale. 1 We assume there are no light degrees of freedom which carry U(1)X charge, otherwise exponential production of the dark vector may be impeded by the resulting Debye mass. 2 We define X µν = µναβ X αβ /2 with 0123 = 1/ √ −g.

JHEP01(2019)053
The rolling axion induces a tachyonic instability in the gauge boson equation of motion. This allows dark photon modes in a specific frequency band to grow exponentially, amplifying quantum fluctuations of X µ into classical modes and transferring a large fraction of the axion energy density into dark radiation. It is this process of exponential particle production amplifying quantum fluctuations within a characteristic frequency band (or set of length scales) that we later identify as the source of gravitational radiation.

Dark photon production
We now review the particle production mechanism more concretely, following refs. [12][13][14][15]. Working in the Coulomb gauge defined by ∇ · X = 0, the equation of motion for X is To study the production of dark photons, we quantize the dark gauge field aŝ and the circular polarization vectors satisfy k · ε ± = 0, k × ε ± = ∓ikε ± , ε ± · ε ± = 0, and ε ± · ε ∓ = 1. It follows that the dark photon with a time-dependent frequency As φ begins to roll, one of the helicities will have negative values of ω 2 for modes in the range 0 < k < α|φ |/f . This corresponds to a tachyonic instability which causes the corresponding dark photon helicity to grow like v ± ∼ e |ω ± |τ . The fastest growing mode is k = α|φ |/(2f ), whereω 2 (k) = −k 2 is the most negative tachyonic frequency. Since the mode with momentumk grows the fastest and has the most energy, it will set the peak of the gauge and gravitational wave power spectra. We numerically solve the equations of motion and calculate the GW spectrum in section 5. However, let us first offer some analytic understanding of the dynamics which determine the shape and amplitude of the GW spectrum.
The dependence ofk on φ means that the scale where most of the energy is being deposited becomes larger as the system evolves. To understand the evolution ofk, we approximate the solution of eq. (2.3) in physical time t as

JHEP01(2019)053
which holds while the friction from production of dark photons is small. With this, we can approximate the envelope of φ as

Tachyonic production band
It is important to point out that since the helicity which experiences tachyonic instability flips between "+ and "− when φ changes sign every half period, no modes experience significant growth unless their growth timescale |ω| −1 is less than the conformal oscillation time (am) −1 . We define the tachyonic production band k − < k < k + as the range of modes which are both tachyonic and have growth times less than (am) −1 , where k − and k + are given by solving ω 2 = −(am) 2 . The result is from which we see that αθ > 2 is required to have the tachyonic production band open initially. As the universe expands, the band contracts while simultaneously shifting toward lower k likek ∝ a −1/2 . The band closes when (a/a osc ) = (αθ/2) 2/3 wherek = a osc m (αθ/2) 2/3 . As an example, to keep the band open until the scale factor has grown by a factor of ∼ 5, one requires αθ 20. 3 We will later use the value of the scale factor when the tachyonic band closes to estimate the scalek at the time of gravitational wave production.

Parity violation and dark photon chirality
As was pointed out in ref. [33], the gauge field helicities are not produced in equal amounts because the operator αφX X/f violates parity when φ = 0. This parity violating effect manifests itself as friction due to production of dark photons of a single helicity. From misalignment arguments we expect φ = θf initially, so the amount of parity violation is controlled by αθ. It also follows that the initial sign of φ and thus the first helicity to become tachyonic (without loss of generality we take this to be "+") are randomly selected. As φ rolls, the amplitude of φ decays due to friction from the expansion of the universe and particle production. As a result, when φ changes sign and the opposite helicity becomes tachyonic, it receives dramatically less enhancement since the growth of the mode functions depends exponentially on φ .
In figure 1 we show the dark photon spectral energy density after the tachyonic band has closed, where the exponential suppression of one helicity with respect to the other can be clearly seen. The broad peak of the spectrum is generated ask evolves toward larger scales. Thus, the value ofk when the axion begins to oscillate gives the right edge of the peak and its value when the tachyonic band closes roughly gives the left edge. As we discuss later and in appendix A, the final spectrum also features additional enhancement from parametric resonance after the tachyonic band has closed.

Gravitational waves
To study gravitational waves, we consider the perturbed metric The linearized Einstein equations for h ij ≡ ah ij in Fourier space are Since gravitational waves are sourced by the highly occupied modes of the dark photon, the quantity Π ij is an operator

JHEP01(2019)053
Neglecting the a term in the equation of motion for h ij (which vanishes in a radiation dominated universe where a ∝ τ ), the solution for h ij iŝ where G = sin [k(τ − τ )] /k is the causal Green's function for the d'Alembert operator. The gravitational wave spectral energy density is given by where we have averaged over one period ∆τ = 2π/k which gives a factor of 1/2. The function Π 2 (k, τ , τ ) is the Unequal Time Correlator and is defined as . A full calculation of this object and the final formula for the gravitational wave spectrum in terms of the dark photon mode functions can be found in appendix B. Useful for comparison to experiment is the fractional gravitational wave spectral energy density defined by the value of which today is usually plotted as h 2 Ω 0 GW with h = H 0 /100 and H 0 = 67.8 km s −1 Mpc −1 .

Estimating the gravitational wave spectrum
The energy density stored in the axion field at T osc , sets an upper bound for the amount of energy that can be transformed into gravitational radiation. Notice that large decay constants are required to have an observable GW signal in planned future experiments. Normally, for f 10 17 GeV and θ ∼ O(1), the axion mass should be less than about 10 −23 eV so that the relic abundance does not overshoot the observed dark matter density. However, the efficient transfer of axion energy density into dark radiation provided by the tachyonic instability and subsequent phase of parametric resonance reopens this parameter space. The energy taken from φ during the tachyonic phase of particle production is transferred to dark photon modes in the range k − < k < k + with the peak energy deposition JHEP01(2019)053 occuring at the scalek. Therefore, we expectk at the time of gravitational wave emission (which we denote ask * ) to set the location of the peak of the gravitational wave spectrum where the factor of 2 approximates the addition of dark photon momenta. We use the value of the scale factor when the tachyonic band closes a/a osc = (αθ/2) 2/3 to estimate the scale factor a * at the time of gravitational wave emission. With this, we estimate the location of the peak of the GW spectrum as Following refs. [34,35], we postulate a simple scaling relation for the peak amplitude of the gravitational wave spectrum at the time of emission where (a * H * ) −1 is the comoving horizon at the time of emission and we have defined the energy density fraction of the gravitational wave source as Ω s = c eff Ω * φ . Here, c eff is a model dependent factor characterizing the efficiency of converting energy in the source to gravitational waves. Causality gives an upper bound on the gravitational wave amplitude, since k peak < a * H * ≤ a * m corresponds to all modes having longer growth timescales than (ma * ) −1 , so there is no tachyonic particle production. In the radiation dominated era, we have Ω * φ /Ω osc φ = a * /a osc and H * /H osc = (a osc /a * ) 2 ≈ H * /m. Using eq. (4.8) we write the scaling relation eq. (4.11) in terms of the model parameters at the time when the axion begins to oscillate as (4.12) We note that our estimate here gives the contribution to the peak amplitude coming only from the tachyonic phase of particle production. The peak amplitude is further enhanced as the physical momemtum k peak /a redshifts and enters the narrow parametric resonance band (see appendix A). Thus, there is reason to expect this scaling relation to underestimate the actual peak amplitude.

Present time gravitational wave spectrum
To obtain the amplitude and frequency of the gravitational wave spectrum today, we need to account for redshifting. The emitted amplitude Ω * GW is redshifted by a factor with g s,eq = 2+2N eff (7/8)(4/11) = 3.938 and T 0 = 2.73 K. Assuming radiation domination at the time of emission, the amplitude today can also be written as 14)

JHEP01(2019)053
where g ρ is the number of effective degrees of freedom in the thermal bath associated with the energy density and we have g γ ρ,0 = 2 and Ω 0 γ = 5.38 × 10 −5 [36]. Additionally, for the last step we have made the approximation g s, * = g ρ, * , which is very good for emission temperatures T * > m e .
The physical peak frequency redshifts as f peak = k peak /a, so its value today is given by Inserting k peak from eq. (4.10), we see that the peak frequency is related to the axion mass via (4.16)

Results
We numerically solve the coupled axion and dark photon equations of motion, treating the backreaction from particle production by assuming the axion responds to the expectation Since we assume that the axion field is homogeneous, the equation of motion for the gauge modes only depends on | k| = k. The mode functions v λ (k, τ ) are included in eq. (5.1) by discretizing the momenta k and approximating the integral as a sum over simulated modes. We discretize using 5000 equally spaced modes with momenta ranging from 0 to θαma osc , where the mode functions were initially taken to be in the Bunch-Davies vacuum v λ (k, τ τ osc ) = e −ikτ / √ 2k. We start the simulation at the temperature defined by H = m and integrate until energy transfer has ended. All of our benchmark points have Ω φ , Ω X < 1, so the total energy density of the universe is dominated by radiation and we neglect the energy density in the axion and dark photon when calculating the background evolution. All changes in the number of relativistic degrees of freedom are fully taken into account following ref. [37]. We assume no temperature dependence of the axion mass, which is a good approximation when f 10 17 GeV in case of the QCD axion.

Benchmark gravitational wave spectra
To compute the gravitational wave spectra, we express eq. (4.6) in terms of the simulated mode functions (see appendix B for details). Discretizing the resulting double integral over momenta results in the computation time growing as O(N 2 ), so these integrals were computed using a subset of N = 100 of the total 5000 simulated modes. We checked that increasing the number of modes produced no significant changes in our results. We Both peak and amplitude of the numerically obtained GW spectra agree reasonably well with the analytic estimates. The signals are detectable in LISA and current PTA experiments if the peak falls into the most sensitive regions of the experiments. Future experiments with sensitivities significantly below h 2 Ω 0 GW ∼ 10 −13 could even detect the tails of the GW signals and thus probe larger bands of axion masses. In particular, SKA could observe a GW signal from the QCD axion. Since the axion dark matter abundance is very sensitive to small variations of the initial conditions, we only demand that our benchmark points do not grossly overproduce dark matter. Therefore all the points listed in table 1 should be considered consistent with cosmology.
It is intriguing that some of the parameter space for axion dark matter might first be probed by GW detectors. The low mass region 10 −19 eV m 10 −13 eV will be probed indirectly by the black hole superradiance with data from LISA [10], showing some unexpected complementarity of GW measurements by LISA and PTAs.

Chirality of the gravitational wave spectrum
As we discussed in section 3.2, the dark photon population is completely dominated by a single helicity and has a relatively narrow range of momenta corresponding to the modes that experienced significant tachyonic growth. Since gravitational waves are sourced by exponentially amplified dark photon quantum fluctuations, they inherit the parity violation in the dark photon population. The peak of the gravitational wave spectrum comes from the addition of two approximately parallel "+" polarized dark photons of similar momenta k, such that a "+ circularly polarized gravitational wave is produced with momentum ≈ 2k. In contrast, the low-k tail of the gravitational wave spectrum comes from two ap-  proximately anti-parallel "+" polarized dark photons of similar momenta k. This results in an approximate cancellation of the polarizations and momenta, leading to the production of unpolarized, low momentum gravitational waves. These features can be seen in figure 3, where the peak of the gravitational wave spectrum is dominated by "+ polarized gravitational waves while the tail has equal components of both helicities such that the net spectrum is unpolarized.

Relic abundance and N eff
The dark photon modes that become highly occupied in the tachyonic instability phase also allow for further efficient transfer of energy from the axion to the gauge fields via parametric resonance. A detailed discussion of the latter can be found in appendix A. In combination, the two mechanisms can suppress the axion DM abundance by up to 14 JHEP01(2019)053  orders of magnitude, significantly widening the range of axion masses and decay constants that are consistent with cosmology [11]. As an example, the red curve in figure 4 shows the suppression of the axion abundance compared to the case with no particle production for the benchmark model ALP 2.
The energy density of the dark photons dilutes as radiation and changes the number of effective relativistic degrees of freedom (N eff ). At the recombination time, the dark radiation contribution to N eff is given by where the energy density of photons is ρ γ and of dark photons is ρ X . The Planck 2018 TT,TE,EE,lowE+lensing+BAO dataset constrains ∆N eff < 0.3 at 95% confidence level [38]. The changes in N eff generated by our benchmark points, shown in the right most column of table 1, are largely consistent with this bound, with some points having a mild ∼ 1σ tension. The next generation of ground-based telescope (CMB Stage-4) experiments conservatively expect to achieve a sensitivity of ∆N eff < 0.03 [39], which is sufficient to probe much of the parameter space which gives a large gravitational wave signal. If the axion is heavy, decays to dark photons will deplete the relic abundance. The axion decay width to dark photons is given by so for f ∼ 10 17 GeV and α ∼ 50−100 the axion decays before Big Bang Nucleosynthesis for m 10 3 GeV and before recombination for m 10 GeV. The GW peak frequency scales with the axion mass as in eq. (4.16), which for m 10 3 (10) GeV translates into f 0 10 4 (10 3 ) Hz. We note that in this work we ignore the back-scattering of the gauge fields which can induce inhomogeneities in the axion field. A recent study suggests that including this effect JHEP01(2019)053 could make the suppression of the axion relic abundance less efficient [22]. This will not have a dramatic impact on the GW signal, since it is dominantly produced during the large initial drop in the axion energy density. However, the parameter space which is consistent with cosmology would be reduced unless additional mechanisms were invoked to suppress the axion relic abundance.

5.4
Gravitational waves from relaxation with particle production?
Tachyonic particle production can also be used as an alternative to Hubble friction in the relaxion mechanism [27]. As explored in detail in ref. [28], when relaxation with particle production takes place after inflation, the maximum allowed cutoff is around Λ ∼ 10 5 − 10 6 GeV. Assuming the relaxion dominates the energy density of the universe at the time of particle production such that the energy budget available for gravitational waves is maximized, the Hubble parameter at the time of emission is H * ∼ Λ 2 /M P . The most negative tachyonic frequency in the case of the relaxion is k peak ∼ v EW by construction, where v EW = 246 GeV is the electroweak scale. We then estimate the peak of the gravitational wave spectrum at the present time using eqs. (4.11) and (4.14) and find where we have used g ρ, * = g s, * = 106.8 as the number of relativistic degrees of freedom at the time of emission and the range in the result corresponds to considering a maximum cutoff of Λ = 10 5 − 10 6 GeV. The peak frequency at the time of emission is given by v EW and redshifting to the present time yields This result, while simply an order of magnitude estimation, illustrates why gravitational radiation from relaxation with particle production would be very challenging to detect. The present time peak frequency is far above the reach of present or planned future detectors, and while the spectrum has a tail which extends to lower frequencies, even the peak amplitude is many orders of magnitude below the sensitivity of any planned future experiment.

Conclusions
We propose a novel method to search for axions or ALPs that couple to dark photons, without relying on SM couplings. The essential dynamics occur in a radiation dominated era, which is distinct from the inflation and preheating scenarios explored previously in the literature. If the axion-dark photon coupling is sizeable, a tachyonic instability followed by a phase of parametric resonance occurs as the axion rolls, allowing the energy of the axion to be efficiently transferred to dark photons. This suppression of the axion relic abundance allows for larger decay constants without tuning the initial misalignment angle. In this context, it would be important to understand the effects of the gauge field back-scattering, neglected in our simulation.

JHEP01(2019)053
The same mechanism causes vacuum fluctuations of the dark photon to experience exponential growth, resulting in time-dependent anisotropies that source gravitational waves. The GW amplitude is controlled mainly by f /M P and the frequency by the ALP mass, allowing for a wide region of parameter space to be explored by future GW experiments as we show with our numerical results for the GW spectra in figure 2. In addition to the distinct shape of the spectrum, its chiral nature could help distinguish it from other cosmological sources of stochastic GW backgrounds.
Complementary to GW searches, the region of parameter space which gives an observable GW signal will also be probed with future CMB experiments which will constrain the production of dark radiation. It is exciting that these astrophysical observations provide a new window to probe invisible axion models.

A Parametric resonance
Modes which leave the tachyonic production band can still experience significant growth due to parametric resonance with the coherently oscillating axion field, which allows for additional suppression of the relic abundance of the axion. In order to get a better understanding of the interplay between the tachyonic and the parametric resonance bands, we write eq. (3.3) in the form of the Mathieu equation using eq. (3.5) (e.g. [40][41][42]) where z = mt/2, F (z) is a harmonic function with unit amplitude, and we have While the backreaction from dark photon production is small, the tachyonic instabilities and parametric resonance can be studied using the stability and instability regions in the (A k , q) plane. The tachyonic regime is set by A k − 2q < 0, or for comoving momenta k < θαam(a osc /a) 3/2 , in agreement with the discussion in section 3. Roughly speaking, the broad resonance regime is given by q 1 while the narrow resonance regime occurs for q 1. The broad resonance regime is given initially by k ma osc /(2αθ) and k ma osc (αθ/16) 2/3 when the tachyonic band closes. One can show that for αθ 1, the broad resonance regime includes entire tachyonic production band throughout its evolution. Thus, we expect that the initial tachyonic growth phase also includes effects from broad parametric resonance which cannot be easily decoupled.
Once friction from dark photon production becomes important, we expect that the corresponding rapid drop in the axion amplitude (which is not captured by eq. (A.2)) will cause a transition from broad to narrow parametric resonance. The narrow resonance bands are given by A k ≈ n 2 (n = 1, 2, . . .), so the lowest k band corresponds to n = 1 or k = ma/2. The study in ref. [42] finds that once αφ/f 1, then the narrow resonance band around k = ma/2 is the most important and that the enhancement from narrow parametric resonance falls off sharply for k < ma/2 since there are no more narrow bands for lower k. This implies that the energy transfer from the axion to dark photons should end when the largest comoving scale which experienced significant tachyonic growth becomes less than ma/2. Estimating this scale as the initial value ofk leads to the condition where a f is the scale factor when energy transfer from the axion to dark photons ends. We can also give a qualitative explanation of the features observed in the evolution of the axion comoving energy density as being due to the narrow parametric resonance band at k = ma/2 entering and then later exiting the most highly pumped part of the gauge power spectrum. This is shown more explicitly in figure 5. While the parametric resonance phase is essential for suppression of the axion relic abundance, the GW spectrum is dominantly determined by the initial tachyonic phase, where a large fraction of the axion's initial energy is transferred to radiation. This can be seen in figure 6, where we show the GW spectrum at different times, including the large initial drop in the axion's energy, the closing of the tachyonic band, and the total spectrum when energy transfer has ended which includes the contribution from the narrow parametric resonance.

B Computation of the gravitational wave spectrum
The energy momentum tensor and perturbed metric are where M P = (8πG) −1/2 = 2.44 × 10 18 GeV, Π ij (x, τ ) = Λ kl ij T kl (x, τ ) is the anisotropic stress energy, and Λ kl ij = Λ k i Λ l j − 1 2 Λ ij Λ kl is the transverse traceless projector with Λ ij = δ ij − ∂ i ∂ j /∇ 2 . One can easily see that the part of T ij which is proportional to g ij ∝ δ ij will not source gravitational waves. Thus, the relevant part of the stress energy tensor is 3) where we have used the result X ik X kj = B i B j − B 2 δ ij and thrown away the isotropic term. Defining h ij = ah ij and taking the Fourier transform of the GW equation, we find

JHEP01(2019)053
and if we throw away the a term which vanishes in the radiation dominated era where a ∝ τ , the solution for h ij is (B.5) where G is the retarded (or causal) Green's function for the d'Alembert operator. The gravitational wave power spectrum is given by . Inserting the solution for h ij , we find where we have averaged over one period ∆τ = 2π/k which gives a factor of 1/2. The function Π 2 (k, τ , τ ) is called the Unequal Time Correlator (UTC) and is defined via . We now turn to computing this object.

B.1 Unequal time correlator
The Fourier transform of the anisotropic stress requires a convolution Working in the Coulomb gauge defined by ∇ · X = 0, we have X 0 = 0. Taking the Fourier transform of E and B yields E i (q) = X i (q) = v λ (q, τ )ε λ i (q). The transform for B gives B i (q) = −i ijk q j X k (q) = λ q v λ (q, τ )ε λ i (q). Putting it all together and leaving the sum over helicities implied, we have with an angular function defined as and a source function We now promote Π to an operator via v λ (q) →â λ (q)v λ (q), where the operators satisfy the commutation relation

JHEP01(2019)053
The object we need to compute is (B.13) Using the commutation relation, one can show that (B.14) The first term matches the helicities and q to q whereas the second term sends q → k − q and exchanges the helicities. It is easy to see that S is invariant under a transformation of the second type, and because Λ kl ij = Λ lk ij , so is Θ. Put explicitly, we have therefore we can just take twice the first term in eq. (B.14) to arrive at the result (B.18)

B.2 Polarization vectors and angular function
Using Λ abcd = Λ ac Λ bd − 1 2 Λ ab Λ cd where the individual projectors can be written in terms of the polarization vectors as and keeping in mind that the projection operator is only acting on symmetric tensors we can write We notice that the term ε λ a (k)ε λ b (k) picks up a phase exp(i2λφ) under a rotation about k by an angle φ, so we identify the two summands in the equation above as helicity projectors that leave us with the part of the anisotropic stress sourcing one particular GW circular JHEP01(2019)053 polarization. Since the wave equation for h ij is linear, the different polarizations do not interfere and it is sufficient to introduce ) in order to study the polarization of the gravitational wave spectrum. We use the relation to arrive at with the functions Φ 1 ≡ Φ 1 (k, q, θ) and Φ 2 ≡ Φ 2 (k, q, θ) defined as Φ 1 (k, q, θ) ≡ k · q |k||q| = cos θ , Φ 2 (k, q, θ) ≡ k · (k − q) |k||k − q| = k − q cos θ k 2 + q 2 − 2qk cos θ , (B.24) where θ is the angle between k and q. The function |Θ| 2 has a few nice properties which are worth pointing out. The first is an exchange symmetry (1 ↔ 2) under which it is invariant. Perhaps unsurprisingly, the transformation Φ 1 → Φ 2 is equivalent to q → k − q, so the symmetry we identified in eq. (B.15) for Θ has been preserved. Additionally, we see that |Θ λ λ 1 λ 2 | 2 is invariant under λ 1 → −λ 1 , λ 2 → −λ 2 and λ → −λ. The domain of both Φ 1 and Φ 2 is [−1, 1], which one can use to check that the range of |Θ| 2 is [0, 1], thus the function is positive definite and unitary.

B.4 Numerics
In the case where one vector helicity dominates (taken to be "+" ), we have  where q max = θαma osc and L = L(k, q) is a list of simulated modes between |k − q| and |k + q|.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.