Analytic study of dark photon and gravitational wave production from axion

Axion-like fields heavier than about $10^{-27}$eV are expected to oscillate in the radiation dominated epoch when the Hubble parameter drops below their mass. Considering the Chern-Simons coupling with a dark gauge boson, large amount of dark photons are produced during a short time interval through tachyonic resonance instability. The produced dark photons then source gravitational tensor modes leading to chiral gravitational waves. Through this process, one can indirectly probe a large parameter space of coupled axion-dark photon models. In this work we first find an analytic expression for the number density of the dark photons produced during the tachyonic resonance regime. Second, by using the saddle point approximation we find an analytic expression for the gravitational wave spectrum in terms of the mass, coupling and misalignment angle. Our analytic results can be used for the observational analysis of these types of scenarios.


Introduction
Most of the energy content of the universe is in the dark sector and, nowadays, the so-called dark energy and dark matter problems are the big challenges for the standard model of cosmology. It is then reasonable to assume that the dark sector, analogous to the ordinary visible matter sector, has a rich structure of its own particles and forces such as (pseudo-)scalars and gauge bosons [1].
Axion is one of the well-motivated examples of pseudo-scalar fields in the dark sector. While from the bottom-up point of view it can be a solution to the strong CP problem in the standard model of particle physics, from the top-down viewpoint it can arise in UV complete theories, such as string theory, in which a global Peccei-Quinn symmetry is spontaneously broken [2,3,4,5]. The properties of an axion are characterized by two independent parameters: the scale f a of symmetry breaking and the mass m which arises due to nonperturbative effects. There are different bounds on the coupling and mass of axions. For example, avoiding the cosmological overabundance of the QCD axion puts an upper bound f a 10 12 GeV. However, bounds on the coupling can be relaxed by considering a coupling of the axion to hidden photons or monopoles [6,7,8] as such a coupling is predicted in some models based on string theory [9], late time entropy production before big bang nucleosynthesis [10,11,12], and dynamical axion misalignment [13,14,15,16,17]. Having larger values of f a makes it possible to consider scenarios in which axion-like fields play significant roles in the earlier times in the expansion history of the universe which was the subject of many studies [18,19,20,21,22,23,24]. The important assumption we make in the following is that the symmetry breaking happens before the end of inflation so that the axion acquires a homogeneous background produced by the misalignment mechanism [25].
While for the case of the QCD axion the mass is determined by the symmetry breaking scale (see Ref. [26]) it is generally expected to be an independent parameter within a wide range of scales. When the Hubble expansion rate drops and becomes comparable to the axion mass during the expansion of the universe, the homogeneous background of axion starts to oscillate. In our scenario, for concreteness, we assume that m 10 −27 eV which corresponds to the Hubble scale at matter radiation equality [27] such that the axion oscillations take place in the radiation dominated (RD) universe. However, our analysis will be more or less the same for axions with smaller masses that start oscillation later, i.e. in the matter dominated era. Furthermore, we ignore the temperature dependence of the mass as well as the self-interactions in the axion field which does not make any qualitative differences in what follows.
Oscillation of the axion results in an energy transfer to fields that are coupled to it via the mechanism of parametric resonance very similar to the mechanism of preheating after inflation [28,29]. The amplification of the fields coupled to the oscillating axion-like fields leads to some interesting features which were the subject of many recent investigations [30,31,32,33,34,35,36,37,38]. In this work, we consider a U (1) gauge boson in the dark sector coupled to the axion field. A dark gauge boson may arise in string compactifications [39,40,41] as well as in vector dark matter models [42,43]. The natural interaction between the axion and dark photon is the so-called Chern-Simons coupling. Such a coupling is very well studied in different contexts like magnetogenesis [44,22,45,46], inflation [47,19,48], preheating [49] and in axion dark matter models to solve the problem of the overabundance of the axions [6,7]. The Lagrangian density for the dark sector takes the form where φ is the axion field, F µν is the field strength tensor of the dark photon with the U (1) gauge symmetry, andF µν = 1 2 µναβ F αβ is the dual field strength tensor with µναβ being the Levi-Civita tensor associated to the metric with the convention 0123 = 1/ √ g. The coupling to the dark photon is characterized by f a and a dimensionless and model-dependent numerical factor α a . The latter can be designed to take large values (even O(100)) for example in clockwork models [50,51,52]. In this work, we are especially interested in large values α a where, as we will see in more detail below, the dark photons become tachyonic and the particle production is very efficient. However, since the amplitude of axion oscillations decreases due to both expansion of the universe and energy transfer to photons, particle creation will soon stop being efficient. The whole process starts and ends very fast so that it takes place entirely in the RD universe. While the axion field interacts weakly with standard model particles, most of the attempts towards detection of the axion is based on its coupling to the standard model (see Refs. [53,54,55,56,57,58]). However, the natural probe of the structure of the dark sector is detection of the gravitational wave (GW) signal. The reason is that according to the equivalence principle, any types of energy or matter would at least minimally couple to gravity. As a result, in the coupled axion-dark photon model described above, it is important to investigate the GW signal from huge amount of dark photons produced via tachyonic instability. This is the main idea investigated in the pioneering works of Refs. [59,60] in which the authors explored the dark photon production and the corresponding distinct signal in the GW spectrum mostly based on numerical analysis. According to their results, there is a good chance to detect the GW signal by the future experiments for a wide range of parameters [60].
Note that the coupled system of axion and dark photon in an expanding background is very complicated. Fig. 1 shows the interplay between different degrees of freedom in the system under consideration. Here is a brief qualitative description of the system. Let us for the moment ignore the metric fluctuations. The initial energy stored in the homogeneous part of the axion field is transferred to its inhomogeneous part as well as the dark photons due to the parametric resonance. The former is due to the self-coupling 1 of the axion potential while the latter is because of the Chern-Simons coupling. As the universe expands the energy transfer becomes less efficient. In the reverse order, the presence of dark photons and axion fluctuations affect the dynamics of the background axion. The former is usually called the backreaction effect. There is also energy transfer between dark photon and axion particles which is sometimes called backscattering. Furthermore, one must consider the fluctuations in the metric and their backreaction on the other fields' dynamics. The complicated and nonlinear dynamics described above is the main reason why most of the analyses performed in the literature are based on the numerical methods.
The aim of this paper is to study the model mostly with analytical techniques. Indeed, one should not expect to be able to obtain analytic understanding of the model in the nonlinear regime due to the complicated dynamics. Instead, we focus on the two legs in Because of the axion oscillations, its own fluctuations and dark photons are produced due to parametric resonance. However, the dynamics is nonlinear as A and δφ affect each other's dynamics (backscattering) as well asφ (backreaction). Besides, we must include metric perturbations. Here we have shown only the tensor modes (h) which can be produced from photons and axion particles. The analytic investigations of this work is focused on the orange arrows.
and also the tensor modes from the dark photons. For the case of dark photon production we use the method of successive scattering matrices used in preheating scenarios [29] to find an analytic expression for the number of dark photons. We will try to find an estimate of the time when backreaction and backscattering effects become important and the analytic results cannot be trusted anymore. As for the GW production, we use our understanding from dark photon production to find an analytic expression for its spectrum. We obtain an approximate closed form expression for the spectrum of GWs by using the method of saddle point approximation as well as by using the fact that dark photons are mostly produced at a specific momentum. These analytic results, though being approximate, not only are useful in giving insights but also improve the template of the GW spectrum for the observational analysis [60]. The rest of the paper is organized as follows. In section 2, based on the method previously used in the investigation of preheating scenarios after inflation, we find analytic solutions for the mode functions of the dark photons during the tachyonic and semi-tachyonic regimes and then we find an analytical expression for the number density of the produced dark photons. We also find an analytic expression for the correction to the effective number of relativistic degrees of freedom induced by the produced dark photons. In section 3, we analytically compute the GW spectrum. Section 4 is devoted to a summary of the results and conclusions. In appendices A, B, and C, we present some useful formulas which are used in the paper while in appendix D we present an analytic study for the possibility of the perturbative decay of the axion.

Production of dark photons
In this section we explore non-perturbative production of the dark photons caused by the oscillations of the axion field φ. As we explained in the Introduction, the model is defined by the Lagrangian density (1.1) for the dark sector that is minimally coupled to the standard Einstein-Hilbert action. Varying the corresponding total action with respect to the axion field φ we find where the field strength associated to the gauge field A µ is F µν = ∂ µ A ν − ∂ ν A µ . Taking the variation of the action with respect to the gauge field, we find the corresponding equation of where we have used the identity 1 √ g ∂ µ ( √ gF µν ) = 0. With the assumed range of the mass of the axion m 10 −27 eV, oscillation starts in the RD era when the Hubble expansion rate drops below the axion mass scale. Moreover, the process of non-perturbative dark photon production and the emission of GW are very short and fits entirely into the RD epoch. The energy density of the universe is then dominated by its radiation content and any contribution from the dark sector can be ignored. As a result, from the Friedmann equation we find that the background geometry is where t is the cosmic time, η = dt/a(t) is the conformal time and a ∼ t 1/2 ∼ η is the scale factor. The electric and magnetic fields associated to the gauge field are defined from the field strength tensor as F 0i = −aE i and F ij = a 2 ijk B k such that they coincide with those measured by an inertial observer 2 . The equation of motion for the axion field φ (2.1) in the background (2.3) takes the formφ with the solution Note that if we neglect spatial variation of the field φ then we get A 0 = 0 which means that the Coulomb gauge condition also implies the temporal gauge condition at leading order [49]. The 2 This definition can be deduced by looking at the components of the strength tensor in the tetrad basis where the metric takes the form gµν = η ab e a µe b ν with η ab being the local inertial metric and the spacetime curvature is encoded in tetrads e a µ. From Eq. (2.3), the nonzero components of the tetrads are e t t = 1 and e a i = a(t)δ a i . Expanding the strength tensor field in the tetrad basis Fµν = f ab e a µe b ν , we find that Fti = aftaδ a i and Fij = a 2 f ab δ a i δ b j . We then define the electric and magnetic fields with respect to the flat space counterpart of the strength tensor f ab as Ei ≡ −ftaδ a i and Bi ≡ 1 equation for the spatial part of the gauge field can be obtained from Eq. (2.2) which turns out to beÄ where A 0 is given by Eq. (2.6). As mentioned before, we assume that the Peccei-Quinn symmetry breaking happened before the end of inflation such that φ has a homogeneous expectation value denoted byφ and in the following we ignore perturbations of the axion field. Note that if we consider self-interactions of the axion field or backscattering effects, perturbations of φ will be produced (see Ref. [28]). We will ignore these effects in our analysis. The equation for the homogeneous vacuum expectation value of the axionφ then reads where the subscript 0 in the right hand side shows the zero mode of the electromagnetic source. Initially the energy density of the dark photons is assumed to be small and we can ignore the source. When the Hubble friction term becomes subdominant due to the expansion, the homogeneous background of the axion starts to oscillate. In this regime, the approximate solution is of the formφ (t) = φ os a 3/2 cos(mt) , where φ os is the initial amplitude at the time of the beginning of the oscillation t os and for simplicity we have set a(t os ) = 1 by rescaling of spatial coordinates. The displacement of axion from its minimum before the time of the beginning of the oscillation then determines the misalignment angle θ ≡ φ os /f a . We define the time of the beginning of the oscillation t os as the time when m = 2H os [61], where m is the mass of the axion. We then choose the origin of the time coordinate so that t os = 0 and that the scale factor can be written as a(t) = (mt + 1) 1/2 . Note that since the oscillation of the axion is semi-harmonic in cosmic time, the analysis is more transparent in terms of the cosmic time compared to the conformal time coordinate. However, the equations of motion for the tensor modes take simpler forms in conformal time which is the topic of the next section.
The quantity of interest is the number density of the dark photons produced by an oscillating axion in an expanding background. To define the number of particles we need to quantize the gauge field. In this regard, first we define the canonical field X i ≡ √ aA i which is then expanded in terms of the creation and annihilation operators in momentum space as follows where ε λ i (k) are the polarization vectors properties of which are presented in appendix A. The effects of time-dependent background are completely encoded in the mode functions χ k,λ (t). Note that due to isotropy the mode functions χ k,λ are functions of the magnitude of the wave vector k = |k|. Substituting Eq. (2.10) in Eq. (2.7), the mode functions satisfy the equations of motion of the form of a harmonic oscillator with a time-dependent frequencÿ where we have omitted the indices of the mode function(s) for brevity while we keep in mind that the behaviour of the two polarizations are different as the frequency is an explicit function of λ which can be ±1 corresponding to right or left-handed photon. The mode function must also satisfy the normalization condition χχ * − χ * χ = i to ensure the standard commutation relations for creation and annihilation operators. Furthermore, we assume that photons live in the Bunch-Davies vacuum state annihilated by allâ k,λ .
If ω 2 ph > 0 then we can always write the solution of Eq. (2.11) for the mode function and its time derivative as [62,63,64] which imply a set of differential equations for time-dependent coefficients α(t) and β(t) as followṡ 14) The normalization condition of the mode function implies |α| 2 − |β| 2 = 1. It can be shown that by defining the instantaneous annihilation operator the Hamiltonian can be diagonalized and the number density of photons with definite momentum and helicity at each time is We assume that the initial state has no particle which means β(0) = 0 and α(0) = 1. As we shall see shortly later, ω 2 ph can also become negative for some time intervals. In that case, we define Ω 2 ph = −ω 2 ph > 0 by means of which we can express the mode function as where the time-dependent coefficients a(t) and b(t) satisfẏ with the constraint ab * − a * b = i. We do not extend the notion of particle number to the regime where ω 2 ph < 0 which does not seem to be well defined [65] and only measure the number density in the regime when ω 2 ph > 0.

Solving the mode function
In order to obtain the number density of the produced dark photons, we need to solve the equations of motion for the mode function Eq. (2.11) with the initial conditions of zero particle number density. A solution in the form of Eq. (2.12) for ω 2 ph > 0 or Eq. (2.17) for ω 2 ph < 0, despite being exact, may not be very useful since solving equations for the coefficients α, β and a, b are as difficult as the equation for the mode function itself. In this subsection we explore the general properties of Eq. (2.11) which is characterized by the effective photon frequency ω 2 ph and then find approximate solutions based on the analytic method employed in preheating scenarios after inflation in the next subsection [28,29].
Let us rewrite Eq. (2.11) in a more appropriate form by substituting from (2.9) and defining the dimensionless variable x ≡ mt as where we have defined the dimensionless frequency in terms of the dimensionless parameters κ ≡ k/m and ϕ ≡ α a φ os /f a . In deriving the above result, we have neglected the term suppressed by H/m in computingφ from Eq. (2.9). Further, in the last step in Eq. (2.21) we have also neglected the mass term induced by the expansion of the universe. The justification is as following: in RD universe we haveḢ ∼ H 2 and therefore the mass term induced by the expansion is proportional to H 2 /m 2 which is less than unity after the axion field started its oscillation and evolves like a −4 which falls off faster than the other terms. Moreover, we will show that the interesting range of parameters are κ 1 and ϕ 1. As a result, we can safely neglect the mass term induced by the expansion of the universe in our upcoming analytic investigations. Since the Hubble parameter is smaller than the mass scale it is reasonable as a first approximation to set a = 1 and neglect the effects of expansion. We then add the effects of the expansion in an appropriate manner. In this approximation we can write Eq. (2.20) in the standard form of the so-called Mathieu equation [66] where we have defined z ≡ 1 2 x − π 4 , A m ≡ 4κ 2 and q m ≡ 2κϕ. In deriving the above equation we have set λ = −1. The equation for the other polarization is obtained by the transformation q m → −q m , but this does not change characteristics of Mathieu equation since it can be realized by a shift in time z → z + π/2. Hence, if we neglect the expansion and also the backreaction of photons to the dynamics ofφ (right hand side of Eq. (2.8)), then there is no parity violation. We will see that this is no longer true if we take into account the expansion of the universe [19,49] and the backreaction.
The interesting feature of Eq. (2.22) is that although its coefficient is periodic with period equal to π, the solution is not necessarily periodic but in general has the form e µz p(z) where p(z) is a periodic function and µ is a complex number called the characteristic exponent (see appendix B for more details). If µ happens to have a positive real part then the amplitude of χ grows exponentially and the solution is called unstable while if µ is purely imaginary then χ is bounded and stable. The boundary, in the parameter space, between stable and unstable solutions can be found by seeking periodic solutions of Eq. (2.22) with period equal to π or 2π (see appendix B). Fig. 2 shows stable regions as shaded area for different values of ϕ and κ while the white area corresponds to unstable solutions. In the case of instability, the number of particles grows exponentially like e 2µz . This phenomena is usually called parametric resonance. The process of exponential particle production continues until backreaction of photons becomes important in the dynamics of the oscillating fieldφ(t). For very small ϕ, unstable solutions are obtained for a narrow band around κ = n/2 for n = 1, 2, . . . which is called narrow resonance while for larger ϕ the instability bands include a wider range of κ thus called broad resonance. The nature of instabilities are different in narrow and broad resonance regimes [29,67]. In the narrow resonance band particles are produced gradually all the time while in the latter case they are produced explosively in a short interval of time in each period of oscillation. The thick red line in Fig. 2 corresponds to κ = ϕ. Below this line ω 2 becomes negative for a time interval in each period of oscillation while above this line ω 2 is always positive in each period. The general expectation is that when ω 2 < 0 the solutions grow exponentially and particle production can be very efficient. This is consistent with Fig. 2 where most of the instability bands are below the line κ = ϕ. However, even below this line there are narrow stable regions which occur due to the destructive interference between positive and negative frequency solutions during the time ω 2 > 0 [68]. Parametric resonance becomes more and more inefficient in the course of time if we take into account the expansion of the universe. This can be studied adiabatically by considering a time dependence in the parameters of the form (2.23) which are related to each other as κ(t) ∝ ϕ(t) 2/3 where the constant of proportionality is fixed by initial conditions. The dashed blue lines in Fig. 2 show this time evolution. Even if a point (ϕ(t), κ(t)) is initially located in an unstable region, the background expansion flows this point along the curves κ(t) ∝ ϕ(t) 2/3 through many stability/instability bands and finally to the region where there is no particle production. In particular, if ϕ > κ so that the point is initially (x = 0) located below the thick black line in Fig. 2, after the time x p given by it crosses the red line in Fig. 2 which means that ω 2 becomes always positive. For a fixed ϕ, different momenta corresponding to different κ have different fates. If κ is very large (but still less than ϕ) then x p is very small and the mode crosses the line κ = ϕ very soon and spends most of its time in the stable shaded area of Fig. 2. Its last chance for particle production is a short period in the first instability band, i.e. the one opening up at around κ = 1/2. As a result, we expect that κ with the most efficient particle production is the one that leaves the first band at a time around x p which means that it has spent all its particle production history below the line κ = ϕ. Since this usually happens in the last stages of the parametric resonance, we expect that the resonance band is narrow where Eq. (B.16) is valid for the width of the first band. Then particle production terminates at x = x en when the lower bound of the inequality in Eq.
If we demand a en = a p then the corresponding momentum κ * is given by Since the amplitude of the unstable mode function is of exponential form we expect that the number of photons corresponding to other momenta be exponentially suppressed compared to κ * . As a result, the spectrum of produced photons has a peak at κ * . Another way to obtain a similar result is an intuitive reasoning due to [59] which we review here for the completeness. At each time the range of momenta which are below the thick red line of Fig. 2, is the range The growth rate is intuitively controlled by the amplitude of |ω| which is maximized in the middle of the above range at ϕ/(2 √ a) with the maximum value |ω max | = ϕ/(2a 3/2 ). On the other hand, the time scale for one oscillation is m −1 which means that in order to have efficient particle production we need |ω| 1 (or equivalently |ω ph | m). As a result, we find that the efficient particle production ends when the scale factor approaches (ϕ/2) 2/3 . The specific momentum which maximizes |ω| at this time is (ϕ/2) 2/3 ∼ κ * which is consistent with Eq. (2.26). As we will see in more detail in the following section, for the tachyonic regime we will have κ * 1. This seems to be counter-intuitive since from perturbation theory we expect the axion particle with mass m can decay into two photons with momentum m/2 corresponding to κ = 1/2. However, parametric resonance is a nonperturbative effect and the naive intuition from the perturbative interactions no longer works here. Even in the narrow resonance regime where parameters might fit into the perturbative analysis, the Bose factor due to the large number density of axion particles enhances the rate of the particle production [29].
Apart from the non-perturbative particle production that we discussed in this subsection, in principle, perturbative decay of the axion to the dark photons can also happen similarly to what happens for inflaton during perturbative reheating after non-perturbative preheating [28,29]. In appendix D, based on the analytic formalism that we extended in this section, we compute the decay rate of the axion to dark photons. Using the decay rate, we consider the possibility for the perturbative decay of the axion. As it is well known, for the favourable ranges of the axion mass and coupling, the perturbative decay may not happen in the time scales shorter than the age of the universe (see appendix D). In this regard, axion remnants can serve as part of the dark matter [6,42,59,61].
In the next subsection, we find semi-analytic solution of (2.11) for the mode function. The readers who are interested only in the results for the GWs may skip the next subsection and directly move to subsection 2.3.

Analytic computation of number density
Although the picture presented in the previous subsection gives a correct qualitative description, it is not very useful in obtaining the analytic solutions for the equations of motion for the mode functions (2.11). In this subsection, following the method of successive scattering introduced in Refs. [29,68], we find analytic solutions for the mode functions of the dark photons. The starting point is to consider the exact forms (2.12) and (2.17) which serve as natural basis for the adiabatic approximation. It is easy to see that in the domain of validity of the adiabatic approximation, the coefficients α(t), β(t), a(t) and b(t) are approximately constant in time and thus no particle production occurs. As a result, it is sufficient to seek regions where the adiabatic approximation breaks down and then compute the change in the particle number in each event. A caveat is that this procedure does not work in the narrow resonance regime where the adiabatic approximation almost always holds (see Ref. [29]). However, since in most part of the analysis we are interested in the broad resonance regime, this is enough for our purposes.
In the broad resonance regime, the adiabatic approximation holds for most of the time evolution of the mode function χ. As a result, we can approximately write the solutions for (2.11) as during the time in which ω 2 > 0. Note that the coefficients α and β are constant in adiabatic approximation. Comparing the above form with Eq. (2.12) we find that as long as the adiabatic approximation is valid the number density n = |β| 2 does not change. The integration in the exponential is only taken over times when ω 2 > 0. As we have seen before, taking κ < ϕ then ω 2 becomes negative for a time interval in each oscillation ofφ. If the adiabatic approximation holds during that interval, the solution would take the same form as Eq. (2.17) and we have

28)
. . where Ω 2 = −ω 2 > 0 and the integration is only taken over the period where ω 2 is negative. We call the regime, in which the approximate solution of the form Eq. (2.28) holds, tachyonic instability. As we already discussed, during the interval ω 2 < 0 the particle number is not welldefined but the mode function grows exponentially and results in typically large amount of new particles after ω 2 becomes positive once again.
Approximate solutions of the form (2.27) or (2.28) are valid for the regions ω 2 > 0 and ω 2 < 0 respectively and needed to be sewn together in regions where adiabatic approximation breaks down, i.e.
1 ω 2 dω dx 1 , or/and 1 ω 3 Let us write ω 2 in Eq. (2.21) as where A and q defined here are different from A m and q m defined under Eq. (2.22) for the standard Mathieu equation. The second term in the parenthesis is shown in Fig. 3 for λ = −1 which is only a sine function with decreasing amplitude q A = ϕ κ √ a as the universe expands. The j−th period ofφ oscillation is during the interval (x j−1 , x j ) with x j = 2πj. If we assume that q A > 1 then as shown by shaded area in Fig. 3, ω 2 becomes negative in j−th oscillation during the times between in which we have approximately set A and q to their values at x min j denoted respectively by A j and q j [69]. Here, x min j ≡ 2πj + (λ − 2)π/2 corresponds to the minimum of λ sin x in the j−th period. It is then evident that the adiabatic approximation breaks down around the points x ± j where ω 2 ≈ 0 which are specified by the red circles in Fig. 3. Thus we can use expressions (2.27) and (2.28) respectively for positive and negative ω 2 except around the turning points x ± j . Around the turning points, the approximate solution for the mode function can be found in terms of Airy functions Ai(x) and Bi(x) by Taylor expanding ω 2 up to linear order. We can use the asymptotic form of Airy functions to match the coefficients of Eqs. (2.27) and (2.28). The procedure is very similar to connecting adiabatic solutions for the wave function in classically allowed and forbidden regions in quantum mechanics. An equivalent approach would be to extend the definition of χ to the complex x−plane [68,70]. By repeating this procedure we can connect the coefficients of Eq. (2.27) during j−th oscillation (α j , β j ) to the coefficients of j + 1−th oscillation (α j+1 , β j+1 ) as [68] α and θ j ≡ θ 0 + j Θ j is the accumulated phase during the time intervals of positive ω 2 up to the point x − j in which θ 0 is some initial phase and for j ≥ 1 is the accumulated phase during j−th oscillation. If we set α 0 = 1 and β 0 = 0 initially, corresponding to the initial vacuum state, then in j−th oscillation the number density of particles is obtained to be [69] This can be used to compute the number density of dark photons as a function of time by computing the integrals (2.33) and (2.34). Fig. 4 shows the number density of dark photons for both polarizations obtained from numerical calculation of the mode function (blue curve). The red triangles show the number density obtained from successive multiplication of matrices in Eq. (2.32) (or Eq. (2.35)) which is in very good agreement with the numerical result. As is evident from Fig. 4, left-handed polarization corresponding to λ = −1 dominates over the righthanded one. The reason is that if we assume ϕ > 0, then the point x min j happens earlier for the left-handed polarization (λ = −1) compared to the right-handed one (λ = +1) by an amount ∆x = π in each oscillation (see Fig. 3). As a result, the magnitude of X j is relatively larger for the left-handed polarization during the period of exponential growth and the total number of particles is exponentially larger. Needless to say the situation is reversed if we choose ϕ to be negative 3 . The number density in the j−th oscillation is mostly controlled by the exponential factor in Eq. (2.35) so in the following we try to find a closed form expression for it. First, note that from Eq. (2.33) we can compute where E(φ; m) is the incomplete elliptic function of the second kind and we have used the approximation introduced in Ref. [68]. As mentioned above, A and q are approximately evaluated at x = x min j for which the scale factor is a j = (2πj + (λ − 2)π/2 + 1) 1/2 for both polarizations λ = ±1. Also for λ = −1, the scale factor is evaluated at ∆x = π earlier in each period. Indeed a j must also carry an index for polarization but to avoid involved notation we do not show that explicitly. After substitution of A j and q j from Eq. (2.30) in terms of a j , the summation in the exponential of Eq. (2.35) can be computed by using the definition of the Hurwitz-Zeta function ζ s (r) = ∞ j=0 1/(j + r) s [71]. After a couple of lines of algebra we obtain where we have defined r λ ≡ (2 + (2 + λ)π)/4π. In the final expression we have replaced the total number of oscillations, j, in terms of a j in the argument of the Hurwitz-Zeta function.
Since the process of production of dark photons is very short it is useful to find approximate expression for (2.37) for a few number of oscillations. In this case, we can obtain approximate expression for the Hurwitz-Zeta functions as for the range 1 a 10. As a result, we can approximately write Eq. (2.37) as between dark photon and the axion. This confirms our earlier claim that the particle production happens in the non-perturbative regime [29].
Having obtained the approximate number density from Eq. (2.40) we can estimate the amount of parity violation by comparing the exponent for two polarizations. Note that if we compare the values of Eq. (2.40) after a couple of oscillations, the dependency of a j on λ becomes unimportant so we can write Note that, as discussed in the last section, we must have κ < ϕ/ √ a for the tachyonic instability to occur, so the above expression is always negative in this regime confirming the subdominance of right-handed polarization compared to the left-handed. Evaluating this for κ * reveals that the amount of parity violation becomes stronger for larger ϕ as expected.
As we include more and more oscillations, the amplitude of the sine term in Eq. (2.30) decreases and the time interval in which ω 2 is negative shrinks until finally after the time x p computed in Eq. (2.24), ω 2 would be always positive and tachyonic instability disappears. However, as shown schematically in Fig. 3, even before x p is reached the time interval between x − j and x + j becomes so tiny that the adiabatic approximation breaks down even at the point x min j [68]. To check this, we compute 1 which reveals that the adiabatic approximation is valid as long as However, due to the expansion the inequality is soon violated which means that the adiabatic approximation is no longer true around the point x min j . Thus, although ω 2 is still negative around x min j , Eq. (2.28) can no longer be exploited. We denote the time when this happens by x et and the corresponding scale factor by a et since this specifies the end of the tachyonic instability regime. For κ * given in Eq. (2.26) using the condition of Eq. (2.42) we find that a et ≈ 0.6ϕ 2/3 . This is approximately the time when, according to the intuitive description of Ref. [59], the tachyonic band closes, although as mentioned before, ω 2 is still negative in each period hence we do not use this terminology here. 4 end of narrow res. Semi-tachyonic instability

Semi-tachyonic regime
After the end of tachyonic instability regime at x et , we can still have particle production due to the parametric resonance if parameters happen to be located in one of the instability bands. Since adiabatic approximation is not valid for the whole interval between x − j and x + j we cannot use Eq.  [29] for the broad parametric resonance regime. The difference is that here ω 2 is negative around the point of approximation, i.e. at x min j while in Ref. [29] the value of ω 2 is always positive. The details of the calculations are given in appendix C.
From the transfer matrix in Eq. (C.8) we can relate n j+1 to n j , which for n j 1, yields where θ tot ≡ ϑ+2θ j +arg β j −arg α j in which ϑ and θ j are defined in Eq. (C.9) and below Eq. (2.33). Note that in this regime we have ς < 1 (see Eq. (2.42)). In related works (Refs. [6,59]) this regime is usually called parametric resonance regime. However, parametric resonance is a generic name for the whole process of particle production due to an oscillating background including the tachyonic regime investigated before. Instead, we call this regime semi-tachyonic instability since particle production still occurs in the time interval where ω 2 < 0. Interestingly the amount of particle production in each oscillation in the semi-tachyonic regime is sensitive to the phase θ tot in Eq. (2.43). One can see that for the number of particles decreases as a result of interference of the solutions before and after scattering (see Ref. [29] for more details). This is the counterpart of the so-called stochastic resonance during preheating in an expanding universe mentioned in Ref. [29]. Although, as mentioned earlier, in our case ω 2 becomes negative around the point of particle production. Due to the complicated structure of θ tot it is not straightforward to find an analytic expression for the total amount of particle production in this regime. . From x et the semi-tachyonic regime starts until x en (the second vertical boundary between colored regions) when we leave the first instability band as discussed around Eq. (2.25). Note that while particle production in the tachyonic regime is very efficient and the number of particles grows exponentially within a short time interval, in the semi-tachyonic regime even more amount of particles are produced although within a much longer time. After x en there would be no instability regime and hence no efficient, nonperturbative particle production. After that we might have the conventional perturbative decay. This possibility is discussed in appendix D.

Backreaction and backscattering
Thus far we have neglected the effects of the produced dark photons on the dynamics of the axion field φ. As briefly discussed in the introduction and Fig. 1, the system includes a complicated interplay between dark photons, axion background and its fluctuations. The corresponding equations are solved numerically in the literature [6,59]. However, initially the energy density of dark photons and axion fluctuations are negligible and the dynamics can be studied linearly, corresponding to the production of dark photons (or similarly axion fluctuations) in a homogeneous background of oscillating axion field. This is the implicit assumption of the previous sections. However, this assumption breaks down after significant dark photons are produced by tachyonic instability and one cannot completely trust the analyses of previous sections afterwards. 5 The produced dark photons can affect the evolution of the axion field in at least two ways. The first is the effect of the produced dark photons on the evolution of the background axion fieldφ(t) that causes its amplitude to decrease faster than a −3/2 and as a result resonance is terminated earlier. This effect is usually called the backreaction. The second deals with the inhomogeneities that are induced by the dark photons through which axion acquires momenta even if the initial axion field was homogeneous. This effect is sometimes called the backscattering. In the following we try to estimate when these effects become important from the linear analysis point of view.
The equation of motion for the background axion field Eq. (2.8) can be written approximately asφ Here we have substituted the quantum expectation value on the right hand side in the Hartree approximation and also used the Weyl ordering in order to get a Hermitian result. By using the definition of the electric and magnetic fields for the dark photons in terms of the mode functions defined in Eq. (2.10), we obtain Hχ χ * + cc. (2.46) In the second line above we have used the solution of Eq. (2.12) with ω 2 > 0 to write the mode functions in terms of the number density of the dark photons where we also defined ψ ≡ 2 dt ω ph + arg β − arg α. The analysis of previous section is valid as long as the right hand side is negligible compared to the terms on the left hand side of Eq. (2.45). As a result, we can roughly equate them to find an estimate of the time when backreaction becomes important. In order to do this we approximate the integral in Eq. (2.46) by evaluating the integrand at k ≈ k * where k * = mκ * is the peak of the momentum of the spectrum of dark photons given by Eq. (2.26). Further, we neglect the subdominant polarization and only consider λ = −1. Thus, we approximately obtain where n * is the number density of dark photons at κ * and in addition we have neglected the effect of the oscillating factor in Eq. (2.46). We denote the time when backreaction effects are important by x br which corresponds to the time when αa br . By using the approximate expression Eq. (2.47) we obtain where in the second approximate equality we have set a br ∼ ϕ 2/3 and the fact that κ * ∼ ϕ 2/3 . 6 By using the above criteria we have numerically solved for the time when backreaction becomes significant. The result is shown in Fig. 6 by solid curves. As expected, for small values of ϕ (or correspondingly small couplings α a ) the backreaction effects are not important because particle production is not very efficient. In these cases, we expect that the expansion of the universe plays the leading role in terminating particle production. Further, as the initial energy of the axion background is higher we expect that the backreaction becomes important later. This is also This is shown by the dashed curves in Fig. 6 and seems to be fairly consistent with the numeric solution. Note that from the analysis of the previous section this should only be valid in the tachyonic regime. Finally, we should also take into account the backscattering effects in which the produced dark photons create axion particles, i.e. generate inhomogeneities in the axion field. By expanding the axion (quantum) field into its Fourier mode we havë As a result, one can write a formal solution by using the Green's function method where G φ is the Green's function of the left hand side of Eq. (2.50). Ignoring the effects of the axion fluctuations on the dynamics of dark photons, one can use the results of the previous section to compute the power spectrum of these fluctuations, φ k φ k = (2π) 3 The backscattering effects can be neglected as long as the energy in the fluctuations of axion is smaller than the energy in its homogeneous part. However, as dark photons are produced this assumption can break down. We denote the time when the backscattering effects become important by x bs . As a rough estimation of this time we use the criteria (2.52) The computation procedure for ∆ φ is very similar to the one that we will present in the following section for the spectrum of the GWs. 7 The final result is approximately the same order as the time when the backreaction becomes important which is given in Eq. (2.48).
It should be noted that the analytic results of the previous sections cannot be employed after the time when the backreaction/backscattering effects become significant and the system must be studied numerically [49,22,72,6,7]. As we have seen above, for some parameters one can neglect these effects and follow the linear analysis of the previous section. However, for the large values of the coupling it is the nonlinear effects that become important. Further, the production of dark photons probably continues for sometime after backreaction/backscattering effects become important. Since our analysis is blind to these effects, in what follows (i) we will assume that nonlinear effects (and not the expansion) terminate particle production and (ii) we use Eq. (2.48) to estimate the number density of dark photons. As we will see in the following, we obtain expressions for ∆N eff and Ω λ GW which are consistent with the numerical simulations.

Contribution of dark photons to N eff
Although dark photons are not visible, they contribute to the number of relativistic degrees of freedom. Their corrections to the effective number of relativistic degrees of freedom N eff before matter-radiation equality and around the recombination time, will be where ρ A is the energy density of dark photons produced non-perturbatively that we discussed in subsection 2.2 and ρ γ is the energy density of the standard radiation which is the dominant energy content of the universe in our scenario. We can find an estimate of the correction (2.53) in terms of the parameters of the model with our analytic results from previous subsections. The energy density of the produced dark photons can be written in terms of its number density as (2.54) Similar to the previous subsection, we approximate the integral by its value at the peak momentum k * and as before assume that the backreaction of dark photons terminate particle production. Thus we can approximately write where in the second semi-equality we have substituted the value of n * from Eq. (2.48). Substituting Eq. (2.55) in (2.53), the correction to the effective number of relativistic degrees of freedom at the time of backreaction x br is obtained to be GeV is the reduced Planck mass, and also we have set H br ∼ m ϕ −4/3 .
We have compared the approximate expression of Eq. (2.56) with the numerical results of Ref. [59] in Fig. 7. They have done the numerical simulation taking into account the backreaction effects but neglected any other effects including the fluctuations in the axion field. The parameters they have used for simulations is presented in the table of Fig. 7. This figure shows that Eq. (2.56) is in qualitative agreement with the numerical studies. The discrepancy is due to the many approximations we have made during the computation such as neglecting the subdominant dark photon polarization, neglecting dark photons with momenta far from k * and neglecting dark photons produced after the time of backreaction. Improvement of those approximations should result in even better agreement with the numerical results.
Before concluding this section we note that the current bound on the change in the effective number of relativistic degrees of freedom is ∆N eff < 0.3 [27]. Fig. 7 shows that the model under consideration can be consistent with observations.

Production of gravitational waves
In the previous section, we have seen that in the presence of Chern-Simons coupling, considerable amount of energy can be extracted from the background axion field by tachyonic instability, and dark photons then will be produced mostly with a specific energy (or equivalently momentum) denoted by k * = mκ * where κ * is given by Eq. (2.26). As the dark sector is minimally coupled to the gravity sector, the produced dark photons provide a nonlinear source for the linear tensor perturbations in the metric sector which are defined as δg ij = a 2 h ij subject to the traceless and transverse conditions h ii = 0 = ∂ i h ij (see [73] for a review of GW in cosmology).
One may also consider scalar perturbations in both dark and gravity sectors. Note that we have assumed a homogeneous configuration for the axion field. However, as we have discussed in subsection 2.3, inhomogeneities in the axion field will show up due to the backscattering of the produced dark photons even if we start with a completely homogeneous axion field at the beginning. Indeed these inhomogeneities in the axion field will not affect curvature perturbations during the RD era as the axion background energy density is negligible compared to the standard radiation energy density. In this respect, the inhomogeneities in the axion field can be thought of as isocurvature modes. Of more importance is however the nonlinear effects of these scalar modes on the linear evolution of the tensor modes which was the subject of some very recent studies [74,75,76,77]. Moreover, as was mentioned in subsection 2.3, we only focus on the case in which the backreaction terminates the process of particle production in a very short interval in the tachyonic resonance regime. Since from the lattice simulations we know that nonlinear backscattering effects are negligible in this regime, we will not take into account the scalar perturbations in the axion field. Indeed, the main features of the produced GWs will survive through our linear analysis. In this regards, in this section we only look for the effects of the dark photons as a source for GWs.
The tensor modes h ij can be expressed as a combination of right or left handed polarization as follows where e λ ij (k) are circular polarization tensors properties of which are presented in appendix A. Considering the dark sector Lagrangian density (1.1) minimally coupled to the Einstein-Hilbert Lagrangian density, it is straightforward to show that the equations of motion for h λ (k) is where Π λ ≡ e −λ ij T ij (k) and T ij (k) is the Fourier transform of the spatial part of the electromagnetic energy-momentum tensor given by The physically interesting object is the GW energy density per logarithmic wave number, dρ GW /d log k, related to the energy density ρ GW of GWs in the form where means suitable temporal and spatial averaging [73]. To simplify the derivations, in what follows we use conformal time η but in order to use the results of the previous section we rewrite all results at the end in terms of the cosmic time. In order to obtain the GW energy density, we look at the power spectrum ∆(k, λ, η) of the GWs defined by the two-point correlation function of the tensor polarizations as follows where we have used the translation and rotation symmetry of the FLRW universe to ensure the proportionality to delta functions in the absence of any net polarization [73,78]. Note that the parity violating interaction shows itself through the explicit polarization dependence of ∆. The solution for the source free tensor modes, i.e. Eq. (3.2) in the RD universe with Π λ = 0, is easily obtained to be where A λ and B λ are two constants. For this solution, it is easy to show that in the sub-horizon limit k aH, in which ∆ ∝ |A λ | 2 + |B λ | 2 where we have ignored oscillatory terms which vanish after taking the temporal average. For our case with the source, we can write the solution by means of the Green's theorem as where η i is the initial conformal time at which the source is turned on, i.e. corresponding to the time when the axion field starts to oscillate. The retarded Green's function of Eq. (3.2) in terms of the conformal time is given by G h (η, η ) = (a(η )/a(η)) sin(k(η − η ))/k. We assume that before η i the solution to the homogeneous equation vanishes 8 . Since the vector modes start to decay soon after the parametric amplification, we can safely assume that the source Π λ is nonzero only until the time η f . After that (η η f ), the solution can be expressed in the form of Eq. (3.6) and the coefficients A λ and B λ can be computed by matching the solutions. As a result, the amplitude of the two-point function of GWs at late time then takes the following form where a prime over the unequal-time average means that the delta functions are not included. To compute the unequal-time average, we need to express Π λ in terms of the mode functions which were obtained in the last section. By using Eq. (3.3), and from the definition of electric and magnetic fields for the dark photons, after some algebra we get where we have defined p ≡ k − p. The first part of the integrand above, P r,r λ , comes from the polarization vectors and tensors as P r,r λ = 1 4 (1 + rλ cos γ)(1 + r λ cos γ ) , (3.11) where cos γ ≡k.p and cos γ ≡k.p . The second part of the integrand above, M r,r , is constructed from the mode functions as M r,r = Ẋ r (p) − HX r (p)/2 Ẋ r (p ) − HX r (p )/2 + rr pp a 2 X r (p)X r (p ) , (3.12) where X r = χ r a k,r + χ * r a † −k,r and χ r (t) are the mode functions for the dark photons. The derivation of the mode functions was one of the main focuses of the previous section. Note that the explicit dependence on the GW polarization λ comes only through the term P r,r λ . By using the Wick's theorem and keeping only the connected part, we obtain (3.13) where we have defined υ p,r ≡χ p,r − Hχ p,r /2 (note that the dot is the derivative with respect to the cosmic time) and the prime over the mode functions means that the corresponding momentum and polarization are p and r respectively. The subscripts η and η of the parentheses denote that the expressions are evaluated at the corresponding conformal times.
Using the above relations and after some re-arrangements, we find where we have defined the dimensionless fractional GW energy density Ω λ GW and ρ c = 3M 2 Pl H 2 is the critical energy density of the universe. In the above expression for the fractional GW energy density we have defined cos kη sin kη υ p,r υ p ,r + rr pp a 2 χ p,r χ p ,r . (3.15) We are interested in obtaining an analytic expression for Eq. (3.14) at the present time. In this regard, we have to impose a couple of approximations in various steps. The first approximation we make is to replace the spectrum of dark photons by an exact delta function in momentum space located at k * (= mκ * ) and also to take into account only the left-handed polarization. This is a reasonable approximation according to the discussion of the last section that the spectrum of photons is sharply peaked around k * and that one of the two polarizations is dominant. This approximation helps us to get rid of the momentum integration and the sum over polarizations in Eq. (3.14). More explicitly let us write the integral measure in the momentum space as d 3 p = pp k dp dp dΦ , (3.16) where Φ is the azimuthal angle with trivial integration yielding 2π. Also after trivial algebra, we find 17) . .
k 2k * Figure 8: By momentum conservation we must have k = p + p where k is the momentum of GW and p and p are momenta of dark photons. The main contribution of dark photons come from momenta p = p = k * which are significantly produced in tachyonic regime. As a result, the momentum of GW can be at most 2k * corresponding to a folded triangle.
which can be used in Eq. (3.11). Following the above approximation, we can set p p k * and ignore the contribution from the dark photons with other momenta. In this respect, for a general momentum k of the GWs, we have In other words, the three vectors p, p and k, which must satisfy the momentum conservation, form an isosceles triangle (see Fig. 8). Note that based on this approximation we must have k 2k * and modes with higher momenta cannot be produced. Further, considering the fact that one of the polarizations is dominant (in our convention it is chosen to be the left-handed one) we set r = r = −1. As a result, Eq. (3.14) can be written approximately as where a star in the subscript of I c and I s indicates that we have set p = p = k * , the superscript "−−" means that we have set r = r = −1, and ∆k * is the width of the peak in the spectrum of the produced photons which can be approximated as ∆k * ∼ k * . The remaining part is to take the time integral appearing in Eq. (3.15). Based on the approximation p = p = k * and r = r = −1 that we use here, it takes the form where we have neglected a term with the Hubble expansion rate in υ r,p χ r,p . Plugging the general solution of the mode functions from Eqs.
where we have defined the rescaled conformal time y ≡ η/m and the dimensionless GW momentum κ = k/m. Further, κ * is given by Eq. (2.26) and by ω * we mean that we have set the momentum of the dark photon to κ * in Eq. (2.21). Note that the dependency on the momentum of GWs appeared only through the argument of sine/cosine function. From the definition of the number of particles we have |α| 2 Form of the integrand near different saddle points. |αβ| n * which, from the previous section, we know that it grows exponentially due to tachyonic and semi-tachyonic particle productions. As a result, we expect that the value of the integral is dominated by the late time behaviour of its integrand which is exponentially larger than its early values. That is the reason why in writing Eq. Note that the amplitude of I c or I s is mainly determined by the prefactor n * . The integrand is of oscillating nature and can be approximately computed by the method of stationary phase. In the regime where particle production is no longer active, the phases arg α and arg β are approximately constant. As a result, the term in the second line of Eq. The saddle points are determined as usual by the equation dΨ/dx which results in Note that from the momentum conservation we have κ < 2κ * and also κ * < ϕ, thus the right hand side of the above equation is always between zero and unity which means that the above equation has many solutions as the sine function oscillates. This is shown in Fig. 9. As a result, to compute the integral we need to sum over all the saddle points. Fortunately, most of these contributions cancel each other as justified below.
We label the solutions of Eq. (3.24) by x ± i where i enumerates the number of oscillations of the sine function while ± shows that we have two solutions located in either side of the maximum of the sine function (see Fig. 9). However most of these saddle points do not contribute to the final result. The reason is that for the early solutions the contribution from the point x − i cancels the contribution from the successive point x + i on the other side of the maximum since the sign of the term d 2 Ψ/dx 2 corresponding to x − i and x + i are opposite to each other so they contribute to the final result with different signs (see for example Ref. [79]). Here, we have neglected the difference in the value of the scale factor corresponding to to x − i and x + i (similar approximation was exploited around Eq. (2.31)). Thus, contributions from the two successive saddle points x − i and x + i cancel each other at early times and we do not need to consider them. In the course of time, the successive saddle points x − i and x + i becomes closer and closer since the amplitude decreases due to the expansion and at some point x − i and x + i become so close to each other (and to the maximum of the sine function). Besides, the width of the Gaussian factor by which we approximate the contributions from each saddle point is proportional to 1/ d 2 Ψ/dx 2 and grows with the scale factor. As a result, the two solutions x − i and x + i cannot be treated as two separate saddle points. Based on these considerations, we can only consider the solution to Eq. (3.24) for which we have sin x ≈ 1.
This is a unique solution denoted by x s and with good accuracy gives the dominant contribution to the integral with properties Then the value of the integral in Eq. (3.23) is computed to be where Ψ s is the value of the phase at the dominant saddle point. Plugging the above result into the integral of Eq. (3.22) and then using the outcome in Eq. (3.19), we obtain the following result for the GWs at the time of emission Ω λ GW,em where the subscript em denotes the value of the corresponding quantity at the time of emission. Since GWs are produced in RD era we can write a 4 em H 2 em a 4 os H 2 os m 2 /4 as we have set a os = 1 and m = 2H os . The amplitude of the GWs also depends on the number of produced photons n * . As we already asserted we use Eq. (2.48) as an estimate of n * at the time backreaction becomes important. Furthermore, the value of k * can be computed from Eq. (2.26) in terms of the parameters of the setup. Putting these all together, we find the following expression for the GW energy density at the time of emission (3.28) The spectrum of GW is depicted in Fig. 10. The solid blue and red curves are analytic results of Eq. show the analytic spectrum we obtained in Eq. (3.28) while the dashed curves show the corresponding numerical results found in Ref. [59]. The numeric data is extracted from Fig. 3 of Ref. [59] where the parameters are α a = 55, θ = 1.2 and f a = 10 17 GeV. This corresponds to φ os = θf a = 1.2 × 10 17 GeV and ϕ = α a θ = 66. In extracting the numerical data we have re-scaled the horizontal axes such that the peak of the analytic and numeric spectra happens at the same k/2k * . The numerical value needed for a perfect match is k num * = 1.8 mκ * where κ * is given in Eq. (2.26).
k 2k * is because we have included only photons with momenta ∼ k * . For comparison we have also extracted the numerical data of the spectrum from Ref. [59] for α a = 55, θ = 1.2 and f a = 10 17 GeV. For these parameters the analytic results of previous sections predict the value of k * = mκ * with κ * ≈ 12.4 obtained from Eq. (2.26). However, the peak of the numerical data is larger than this value by a factor less than 2 for this case. In fact the perfect match in Fig. 10 is obtained by re-scaling the numerical data with k num * = 1.8 mκ * . We trace this discrepancy back to the numerous approximations we have made during the computation. Fig. 10 shows that our analysis provides a good approximation for the amplitude of the GW, the position of the peak (up to an O(1) factor) as well as its spectrum specially in the vicinity of the peak. However, it should be noted that the analytic expression Eq. (3.28) cannot be trusted in the IR limit, i.e. for very long wavelengths. The main reason is that for very small momenta, the saddle point approximation fails as the resulting integral is divergent which can be seen from Eqs. (3.25) and (3.26). In fact, it is not difficult to see that in this limit the integral appearing in Eq. (3.14) is actually independent of k. As a result, one expects the momentum dependence of the IR tail be of the form Ω GW ∝ k 3 . This is a universal behaviour for any casual mechanism of GW production [80,81]. This behaviour applies for modes with wavelength k −1 T where T is the typical time or length scale of GWs production. In our setup, we have only focused on the growing modes for which the typical time scale for tachyonic growth is T ∝ m −1 as we explained in the paragraph after Eq. (2.26). As a result, the IR behaviour can be seen for modes k m (corresponding to (k/2k * ) ϕ −2/3 ). This is the region where our analytic expression deviates from the numerical results shown in Fig. 10. Note that in expression Eq. (3.28), the dependence on the coupling α a is very weak and only comes in through the value of k * . Instead, the main dependence is on the initial value of the axion field which is usually written in the form φ os = θf a . While Fig. 10 suggests that what we have here is fairly consistent with the numerical result, it is parametrically different from what has been suggested by Refs. [59,60]. It seems that both derivations have their own pitfalls and do not have enough numerical data to discriminate between them.
The energy density of GW at present time can be computed by taking into account the effects of decoupling of the relativistic particles during the course of expansion, and we have [73] Ω λ GW,0 = Ω λ GW,em g s,0 g s,em where g s and g ρ are the effective numbers of relativistic degrees of freedom associated to the entropy and the energy density respectively, the subscripts 0 and em show the corresponding values at the present time and at the time of emission respectively, and Ω γ,0 is the present value of the fractional dimensionless energy density of radiation. Considering the typical values of g s,0 = 3.91 [82], g s,em g ρ,em = 106.8, g ρ,0 = 2, Ω γ,0 = 5.38 × 10 −5 [83], the GW energy density today can be obtained from (3.28) and (3.29) as It should be noted that the analytic results for the GW spectrum in Eq. (3.28) or (3.30) can be used as templates in GW data analysis to extract signals. Further, one can easily relate the observed spectrum to the fundamental parameters of the theory by fitting the values of k * and φ os . Finally, Eq. (3.28) shows that the chirality of the dark photons has resulted in a polarized GW spectrum which is the consequence of the parity violating nature of the Chern-Simons interaction. Interestingly, astrophysical stochastic GWs are expected to be unpolarized. Hence it is in principle possible to discriminate the polarized GW signal we have obtained in this model from the unpolarized astrophysical sources by the GW detectors. In this regard, the parity violation can be measured following the method suggested in Refs. [84,85,86,87,88,89,90].

Summary and conclusions
For the range of masses m 10 −27 eV, axion-like fields are expected to oscillate in the radiation dominated era when the Hubble expansion rate drops below their mass scale. The natural interaction between the axion-like field and the gauge bosons is through the Chern-Simons coupling. If the axion-like field is coupled to the dark gauge bosons, the coupling constant can be large enough that sizeable amounts of energy can be transferred from the axion-like field to the dark gauge boson through the tachyonic resonance process. These types of scenarios are well studied in both the early universe setups like inflation and late time scenarios like axionic vector dark matter scenarios.
Some interesting results arise when we consider interaction of the axion-like field with the hidden sector fields such as the dark gauge boson that is studied in the present paper. An important question is whether we can detect these axion-like fields when they do not interact directly with the visible sector. Since any type of matter is universally coupled to gravity according to the equivalence principle, gravitational waves are the natural (if not the only) candidates to probe these axion-like fields. More precisely, in our setup, the amplified dark photons serve as a quadratic source for the linear equations of motion of the tensor perturbations leading to the production of gravitational waves. Due to the parity-violating nature of the Chern-Simons interaction, one helicity of tensor modes is more amplified than the other, leading to the chiral gravitational waves. This scenario was recently suggested in Refs. [59,60] and analyzed by means of numerical methods. In this paper, we have studied this scenario analytically.
Section 2 is devoted to the analytic study of dark photon production from axion-like field. Ignoring the nonlinear effects, such as backreaction and backscattering, we have studied the amplification of the dark photons due to the tachyonic and also semi-tachyonic resonances. Using the methods that had been adopted to study the non-perturbative preheating process after inflation, we have found an analytic expression for the number density of the produced particles in terms of the mass and coupling of the axion. We have also studied the effects of backreaction and backscattering and estimated the time when these effects become important and beyond which the linear analysis cannot be trusted anymore. Further, for the range of parameters where these effects become significant, we have obtained an analytic expression for the number density of dark photons. With this assumption we have obtained an analytic estimation for the correction to the effective number of the relativistic degrees of freedom induced by the produced dark photons. Comparison with numerical results of Ref. [59] shows the consistency of our result within an O(1) prefactor.
In section 3 we looked for the effects of amplified dark photons on the gravitational wave spectral density. Since the tachyonic resonance of the dark photon is the source of the tensor modes, it is not easy to solve the equations of motion for the gravitational waves analytically and that is the reason why this problem was studied numerically in the previous works. Using the saddle point approximations, we have found that there are many saddle points that may play roles in the integral of the dark photon contributions to the gravitational waves. We have shown that most of these saddle points neatly cancel one another so that we only need to take into account the effects of the most important saddle point. In this regard, we have found an analytical expression for the peak of the gravitational wave spectral energy density in terms of the parameters of the model. We have compared the analytic result for the gravitational wave spectrum with numerical simulations of Ref. [59] and showed that they are in a very good agreement.
The analytic results of our paper, besides giving insight to the process of particle production of the model, are useful for observational purposes. The analytic expression for the gravitational wave spectrum can be used as a template to extract signal from observations. By fitting to the data one can find the amplitude and the position of the peak which are directly related to the parameters of the model. Further, one can in principle distinguish between other gravitational wave production mechanisms by checking consistency of other features such as the chirality of gravitational waves or the change in the effective relativistic degrees of freedom.

Supplementary note
After completion of this work, the paper [91] appeared where the authors have performed precise lattice simulations of the model, i.e. they have solved for the coupled system of axion, dark photon and metric fluctuations. Their numerical analyses confirm our analytical results at the linear regime. On the other hand, new features are found at the nonlinear regime when the backreaction and backscattering effects become important. First, the axion remnant decays with slower rate in comparison with the linear analysis result. Second, the chirality of the gravitational waves can be washed out for some allowed regions of the parameter space.
where γ = cos −1 (k i k i /kk ) is the angle between the two vectors k and k . Similarly we can construct linear polarization tensors where superscripts p and c correspond to plus and cross polarizations. Furthermore, circular polarization tensor is defined to be with λ = ±1. A useful relation between polarization vector and tensor is the following such that all identities of Eqs. (A.3)-(A.5) can be translated straightforwardly to the case of tensors.

B Basic Floquet theory
In this appendix we review the classification of the solutions of a linear differential equation with periodic coefficients, namely the Floquet theory. For more details see for example Refs. [66,92]. If we have the linear system of differential equationṡ with A(t) a T -periodic n × n matrix and x a n × 1 vector, then a typical solution x(t) does not have to be periodic but instead is of the form where p(t) is a periodic function with period T and µ is one of the n characteristic exponents µ i which are intrinsic to the problem and satisfy More specifically, if we have n linearly independent solutions x (n) (t) of the problem, we construct the principal fundamental matrix X(t) such that X(0) = I. Then ρ i ≡ exp(µ i T ) are eigenvalues of B ≡ X(T ) and ρ i are sometimes called characteristic multipliers. Note that we have three categories • |ρ| < 1 which means Re(µ) < 0 then x(t) tends to zero at infinity.
The special interesting case for us is a second order linear differential equation of the form where a(t) is T -periodic. It can be put into the general form of Eq. (B.1) if we define Then we can find two linearly independent solutions with linearly independent initial conditions Form Eq. (B.3) and definition of ρ i we have If we define θ ≡ tr(B)/2 then we have In terms of characteristic exponents µ 1 and µ 2 we find out that µ 1 + µ 2 = 0 and cosh(µ 1 T ) = θ.
Then it is easy to see that the solution is stable for −1 < θ < 1 and unstable if |θ| > 1. The boundary between stable and unstable solution is given by the condition θ = 1 or θ = −1 which correspond to a T -periodic or 2T -periodic solutions respectively [92]. Thus we can find boundaries, in the parameter space of the problem, between stable and unstable solutions by making an ansatz which is T -periodic or 2T -periodic and demand that they be a solution of Eq. (B.5). This gives us four homogeneous sets of linear system of equations with a n s and b n s in Eqs. (B.11) and (B.12) as unknowns. A nonzero solution for a n s or b n s is achieved if the determinant of the coefficient matrix is zero. This gives us four different relations among parameters of the differential equation. The coefficient matrix is indeed infinite dimensional but we can approximately solve these equations by keeping only a N × N sub-matrix for large enough N . The interesting second order equation for us is the Mathieu equationẍ+(A m −2q m cos(2t))x = 0 in which T = π. In the limit that q m A m such that we can ignore the oscillating term, we can approximately construct the matrix B defined below Eq. (B.4) as follows As a result, we have θ = cos π √ A m . Then we must have A m = (2n) 2 for θ = 1 and A m = (2n − 1) 2 for θ = −1 with n = 1, 2, . . . . We can find periodic solutions as perturbative series in q m x(t) = x 0 + q m x 1 + . . . (B.14) for n = 1, 2, . . . . Following this analysis it can be shown that for the first instability band n = 1, we must have while for higher n the width of the instability band is of order O(q n m ) [93]. Also in the first band the characteristic exponent is approximately given by For larger q m there is no closed form solution and we must solve the system numerically to find the characteristic exponents. The boundaries between the stable and unstable solutions can be found by the ansatz in Eqs. (B.11) and (B.12). For Eq. (2.22) this procedure results in the stability/instability chart shown in Fig. 2.

C Semi-tachyonic regime
In this appendix we derive the transfer matrix in the semi-tachyonic regime. Around x min j (see shaded regions after x et in Fig. 3), the equation for the mode function Eq. (2.20) can be approximately written as where we have defined τ ≡ (q j /2) 1/4 (x − x min j ) and ς 2 is given by Eq. (2.42). The general solution is given by where ν ± ≡ −(1 ± iς 2 )/2 and D ν (τ ) is the parabolic cylinder function which is the solution of the standard differential equation y + (ν + 1/2 − τ 2 /4)y = 0 (see Ref. [94] for more details) with c 1 and c 2 are some constants. Far from the point τ = 0 (or equivalently x = x max j ), in the region where ω 2 > 0, adiabatic approximation is valid and Eq. (2.27) can be used. However, there is an overlapping region which is not too close to the point τ = 0 to break the adiabatic approximation and still not too far so that the solution of the form (C.2) is still valid. The asymptotic form of the parabolic cylinder functions for τ > 0 can be obtained to be and for τ < 0 In the overlapping region, we can write the exponent in Eq. (2.27) as Then in the limit |τ | ς and for τ < 0 we can write Eq. (2.27) as where θ j ≡ θ j + ς 2 4 1 + log 4 and θ j is the accumulated phase by the time τ = 0 defined below Eq. (2.33). Similarly for τ > 0 we write The solutions of Eqs. (C.6) and (C.7) are connected via Eq. (C.2) and we can obtain a transfer matrix to relate the coefficients (α j , β j ) to (α j+1 , β j+1 ) as in which we have used the fact that Γ 1 2 + i ς 2 2 2 = π cosh(πς 2 /2) . The transfer matrix in Eq. (C.8) is the counterpart of Eq. (2.32) obtained in the tachyonic regime. Note that the transfer matrix is the same as obtained in Ref. [29] if we replace κ 2 in their expression by −ς 2 everywhere except in the Logarithm of Eq. (C.9). The transfer matrix Eq. (C.8) is used to obtain a recursive relation for number of dark photons in the main text.

D Possible perturbative decay
The non-perturbative particle production that we studied in section 2 is all one needs to study production of the GWs by dark photons in our setup. However, as we mentioned at the end of subsection 2.1, the subsequent evolution of the remnant of the axion with energy density ρ rem φ after the end of the non-perturbative particle production is important to see whether axion remnant decays or survives. Indeed, depending on the values of the mass and coupling, the decay is possible. If axion remnant decays, it perturbatively produces again dark photons which contribute to the effective number of relativistic degrees of freedom. Otherwise, it survives and would contribute to the dark matter as it is well known in the context of the axionic vector dark matter models. In this appendix, we track analytically the evolution of the axion remnant after the non-perturbative particle production.

D.1 Decay rate of axion to photons
Let us first compute the axion decay rate Γφ →AA to the photons. Of course we know the answer from the standard perturbative analysis of quantum field theory. However, here we find it via the adiabatic approximation method that we presented in section 2 and through the Boltzmann equation following the method suggested in Ref. [95]. The readers who are not interested in this derivation can simply move to the next subsection.
Since we are in the perturbative regime, we can iteratively solve Eq. (2.14) for |β| 1 with the initial conditionsᾱ(x = x ep ) = 1 andβ(x = x ep ) = 0 which yields where the notationᾱ andβ show that they correspond to the perturbative particle production. We use this notation for the number density and energy density of the particles that produce through the perturbative decay as well. Moreover, the initial conditionβ(x = x br ) = 0 corresponds to the state without any particle while we know that some dark photons are already produced through the non-perturbative particle production process. By consideringβ(x = x br ) = 0 and the notationβ used in Eq. (D.1) we mean that we only look at particles that are produced through the perturbative process. In this regime, we approximately have e −2i dx ω(t ) ≈ e −2iκ dx /a(x ) and also from (2.20) we find 1 ω where we have also neglected time derivative of the scale factor as before.
Substituting the above results in (D.1) yields β(x) ≈ − λϕ 8κ where a s = a(x = x s ) is the scale factor at the saddle time x s . Using this result, we find the number density of dark photons during perturbative decay of axion as follows n k,λ = |β| 2 = πmϕ 2 64κ 3 H s , (D. 5) where H s = H(x = x s ) is the value of the Hubble expansion rate at the saddle time. The above expression gives the number of produced particles through the perturbative decay of the axion to the dark photons. Ignoring the vacuum energy density, the energy density of the dark photonsρ A that are produced perturbatively is given by where Γφ →AA is the decay rate and ρ rem φ is the energy density of the remnant of the axion. Taking the time derivative of (D.7) and then comparing the result with the Boltzmann equation (D.8), we find the following expression for the decay rate of the axion remnant which is in agreement with the standard result of quantum field theory for a decay through a trilinear interaction. We use this result to find the time scale of the decay of the axion remnant in the next subsection.

D.2 Fate of the axion remnant
As we already mentioned, some parts of the energy density of the initial homogeneous axion field are transferred to the dark photons during the non-perturbative process of tachyonic and also possibly semi-tachyonic and parametric resonances while the remaining part is the remnant of the axion which is dealt with in this appendix. During the nonperturbative production of the dark photons, the energy conservation equation between the axion field and dark photons can be written in the form where ρφ is the energy density of the initial homogeneous axion field and ρ A is the energy density of dark photon produced non-perturbatively through the resonance processes. The interaction between the axion field and the dark photons is given by the Chern-Simons term φFF . Non-perturbative particle production starts from a os and ends at min [a br , a en ]. Therefore, we should integrate the above equation from x os to min [x br , x en ]. As we already mentioned, for the favourable range of parameter space, backreaction terminates the process of non-perturbative particle production and therefore we assume that non-perturbative particle production happens in a very short time interval and min [a br , a en ] = a br . In that case, the above equation can be integrated approximately to find the axion remnant at the time of end of non-perturbative particle production as where the inequality appears when we take into account the effects of backreaction and backscattering. If the time duration [x os , x br ] is large so that a br aos ∼ O(10), the above approximation is no longer applicable and we need to solve Eq. (D.10) numerically 9 . However, as we will see, we do not need the explicit form of the remnant and we determined it just for the concreteness. After the tachyonic resonance terminates at the time x br , we are left with the axion remnant which approximately given by ρ rem φ (a br ) and also dark photons ρ A that are produced via non-perturbative resonance process.
The decay rate of the axion to the dark photons is given by (D.9). After non-perturbative particle production is terminated by the backreaction, the Hubble expansion rate decreases as H 2 ∼ a −4 in the RD era and, depending on the mass and coupling, it can approach the axion decay rate (D.9). When the Hubble expansion rate drops below the decay rate (D.9), we would have perturbative particle production for the remnant of the axion. This is similar to what happens for the inflaton during perturbative reheating after non-perturbative preheating. From the conservation of the energy after the time x br , we have The above result shows that the axion remnant will decay completely to the dark photons in time scale t decay ∼ Γ −1 φ→AA . To have axionic dark matter, we need to prevent the axion to decay, requiring t decay > H −1 or equivalently Γφ →AA < H .
(D.13) Therefore, from Eq. (D.9) we see that depending on the mass and coupling, the above condition may or may not be satisfied. In summary, the axion mass determines the time that axion starts to oscillate while the combination of the mass and coupling determines the decay time. If the mass and coupling are such that the axion remnant survives, it can contribute to the dark matter [6,42,43,96,97]. This is the necessary condition for the axion dark matter scenarios but it is not enough. For example, depending on the parameters, the remnant can be considered as a part of dark matter or possibly the whole dark matter. In the latter case, there is a low mass bound m ≥ 10 −22 eV [98,5,99]. Otherwise, it can decay into the dark photons which contribute to the number of relativistic degrees of freedom similar to what we considered for the dark photons produced non-perturbatively in subsection 2.4. In this case, axion decays completely in the short time scale Γ −1 φ→AA , as shown in (D.12), and we can approximately consider instantaneous decay. Therefore, the whole of axion remnant energy density is converted to the dark photons through instantaneous perturbative decay around the time x dec when H(x dec ) ∼ Γφ →AA . The shift in the effective number of relativistic degrees of freedom then will be ∆N eff | t dec ≈ 8 7 11 4 4 3 ρ rem φ (a dec ) ρ γ (a dec ) , (D. 14) and the weakest condition that has to be satisfied in order to respect the predictions of the standard big bang cosmology is ∆N eff | t dec < O(1).