Peaky Production of Light Dark Photon Dark Matter

Prepared for submission to JCAP arXiv:2212.11516v1 [hep-ph] 22 Dec 2022 Peaky Production of Light Dark Photon Dark Matter Yuichiro Nakai,a Ryo Namba,b Ippei Obatac,d a Tsung-Dao Lee Institute and School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China b RIKEN Interdisciplinary Theoretical and Mathematical Sciences (iTHEMS), Wako, Saitama 351-0198, Japan c Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany d Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI),UTIAS, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba, 277-8583, Japan E-mail: [email protected], [email protected], [email protected] Abstract. We explore a mechanism to produce a light dark photon dark matter through a coupling between the dark photon field and a spectator scalar field which plays no role in the inflationary expansion of the Universe while rolling down its potential during the inflation. The motion of the spectator field efficiently produces dark photons with large wavelengths which become non-relativistic before the time of matter-radiation equality. The spectrum of the wavelengths is peaky so that the constraint from the isocurvature perturbation can be evaded. The correct relic abundance is then achieved over a wide range of the dark photon mass down to 10−13 eV. Our mechanism favors high-scale inflation models which can be tested in future observations. Furthermore, fluctuations of the dark photon field during inflation produce gravitational waves detectable at future space-based interferometers and/or pulsar timing array experiments. Contents 1 Introduction 1 2 A spectator model 2.1 Background dynamics 2.2 Evolution of the dark photon field 2 3 6 3 The energy density spectrum 9 4 Dark photon DM: constraints and results 4.1 Relic abundance 4.2 Constraints 4.3 Results 12 12 13 17 5 Generation of tensor modes 17 6 Conclusion 21 A Dark photon production with constant n 22 1 Introduction Ultra-light dark matter (DM) can in principle account for the issues around the structures of small scales (see [1] for a review). A light spin-1 vector boson is an intriguing candidate of DM in our Universe. Such a dark photon DM can have a wide range of mass, extending down to the fuzzy DM region [2, 3]. Accordingly, searches for the dark photon DM have been conducted from various directions (see refs. [4, 5] for reviews). The Lyman-α constraint puts a lower bound on the mass, mγ 0 & 10−20 eV [6]. The black hole superradiance constraint spreads over a range of the dark photon mass, 10−20 eV . mγ 0 . 10−11 eV [7–11]. For a light dark photon DM kinetically mixing with the photon, additional constraints such as galactic heating [12] and distortions of the cosmic microwave background (CMB) [13] give upper bounds on the size of the kinetic mixing. Moreover, a number of experiments are being programmed for further explorations. For instance, DM Radio [14, 15] will search for the dark photon DM with a range of mass, 10−12 eV . mγ 0 . 10−3 eV. Recently, it has been revealed that laser interferometers for the detection of gravitational waves (GWs) are also sensitive to the vector DM with mass 10−18 eV . mγ 0 . 10−11 eV [16–18]. Therefore the phenomenology of a new vector boson DM is an exciting arena for the DM searches. Depending on the dark photon mass mγ 0 , different mechanisms for the production of the dark photon DM have been proposed. Inflationary fluctuations can give the correct DM abundance for mγ 0 & µeV [19]. As in the case of axions, the dark photon DM can also be produced by a misalignment mechanism [20], while a highly tuned coupling to the curvature is required [13]. Refs. [21–23] have realized a light dark photon DM whose mass is below µeV by the oscillation of an axion-like field coupling to the dark photon or a scalar field charged under the dark U (1). A network of cosmic strings can also produce a light dark photon DM [24]. –1– A dark photon coupling to the inflaton provides another attractive possibility for the dark photon DM production. Ref. [25] has introduced a coupling ϕF 0 Fe0 where ϕ denotes the inflaton, F 0 is the dark photon field strength and Fe0 is its dual. It was discussed that a dark photon DM with mass mγ 0 & µeV is produced through tachyonic instability during inflation. A more generic coupling I(ϕ)2 F 0 F 0 with some function I has been explored in refs. [26, 27]. This type of coupling in a cosmological context was first introduced as a mechanism of inflationary magnetogenesis [28] and has been actively studied [29–46]. A crucial feature of the coupling is that, depending on the achieved time dependence of I, diverse shapes of the produced dark photon spectrum can be realized. Thanks to this nature, dark photons with large wavelengths can be efficiently produced by the inflaton motion. The spectrum of produced dark photons has a peak around the wavelength that exits the horizon at the earliest time. The main target of the present paper is the production of a light dark photon DM with a range of mass, 10−11 eV . mγ 0 . 10−6 eV, which is above the mass region constrained by the black hole superradiance and below the region in which inflationary fluctuations can give the correct relic abundance. This range of mass is covered by the future DM radio experiment. Instead of a direct coupling to the inflaton, we here consider a coupling between the dark photon field and a spectator scalar field which does not play a role in the inflationary expansion but whose motion along its potential during inflation leads to production of the dark photon DM. The motion of the spectator field can efficiently produce dark photons with large wavelengths which become non-relativistic before the time of matter-radiation equality. It has been known that the production of a light dark photon DM through a coupling to the inflaton would receive a stringent constraint from the isocurvature perturbation [47]. On the other hand, in the present scenario, the wavelength spectrum of produced dark photons is peaky so that the constraint from the isocurvature perturbation based on the CMB observations can be evaded. The correct relic abundance is then achieved over the target mass range and below, and we show that the mechanism favors high-scale inflation models to provide the lower dark photon mass. In addition, fluctuations of the dark photon field during inflation produce GWs which can be detected at future space-based interferometers such as DECIGO [48], BBO [49] or µAres [50] and also pulsar timing array experiments such as SKA [51]. The rest of the paper is organized as follows. In section 2, we discuss the production of a light dark photon DM through a coupling between the dark photon field and a spectator scalar field. In Section 3, we compute the power spectrum of the dark photon. Section 4 then explores the parameter space that the correct DM abundance is realized. In section 5, the generation of tensor modes is investigated. Section 6 is devoted to conclusions and discussions. 2 A spectator model We consider a model of a light dark photon DM in which a dark photon field A0µ is kinetically coupled to a spectator scalar field σ. The total action is given as follows:  2  Z MPl 1 I(σ)2 0 0µν 1 2 0 0µ 4 √ µ S = d x −g R + Linf − ∂µ σ∂ σ − V (σ) − Fµν F − mγ 0 Aµ A , (2.1) 2 2 4 2 where R is the Ricci scalar associated with the 4-D metric gµν , whose determinant is denoted 0 = ∂ A0 − ∂ A0 is the field strength of the dark photon, and L by g, Fµν µ ν ν µ inf is the Lagrangian –2– density of the inflaton. In what follows, we do not specify the inflaton dynamics except for the assumption that it can be well approximated by de Sitter for our purpose. In (2.1), V (σ) is the potential of the spectator field σ, and I(σ) is a function of σ which will be specified later. The last term denotes a Proca term of the dark photon with a constant mass mγ 0 . We assume that there is no direct coupling between the dark photon and the Standard Model (SM) sector and no mixing with the ordinary photon. A crucial point of our work is to turn on a nontrivial evolution of the dark photon production. This enables a non-monotonic spectrum of the produced dark photon, resulting in a spectral peak at a certain wavenumber kpeak . This additional freedom is essential to evade the stringent constraint on the CMB isocurvature modes. The evolution of the production is controlled by the function I(σ) in (2.1), whose behavior inherits the dynamics of σ. In the following subsections, we hence consider a slow but non-constant homogeneous motion of σ and then perform the computation of the corresponding production of the dark photon. 2.1 Background dynamics Let us first consider the background dynamics of the model. Throughout the present work, we assume that there is no vacuum expectation value of the vector field:1 hA0µ i = 0 . (2.2) Thanks to this assumption, and ensuring a posteriori that the energy density of the dark photon is a subdominant component, we can consider the flat Friedmann-Lemaı̂tre-RobertsonWalker (FLRW) metric, ds2 = −dt2 + a(t)2 dx2 = a(τ )2 (−dτ 2 + dx2 ) , (2.3) as the background spacetime. Here, a is the scale factor and τ denotes the conformal time, dτ = dt/a. Throughout the paper, we assume that inflation proceeds as de Sitter essentially, and the time dependence of the scale factor is well approximated by a ' exp (Ht) ' −1/(Hτ ), where H ≡ ∂t a/a is the Hubble rate during inflation. For the kinetic function I of σ, we adopt the following simple exponential form:   σ I(σ) = I0 exp , (2.4) Λ1 with some dimension-1 energy scale Λ1 . The normalization constant I0 is assumed to be I0 = 1 to realize I(σ → 0) = 1 after inflation ends when σ settles into its potential minimum and recovers the canonical kinetic term of the dark photon field. The crucial parameter that dictates production of the dark photon, especially of the transverse modes, is the rate of the time variation of I with respect to the number of e-foldings N = ln a, i.e., n≡ dI/dN dσ/dN = , I Λ1 dN = Hdt . (2.5) 1 This assumption in fact depends on the scales of interest. Even if the classical background of A0µ is absent, the dark photon field can be generated from the quantum vacuum, as in our current consideration. If the produced field has a wavelength larger than a certain scale with a nonzero variance, then each local patch smaller than that stochastically experiences a coherent field of dark photon wave, which, strictly speaking, breaks spatial isotropy [45]. However, we are here interested in generation of the dark photon on scales smaller than the CMB ones, and therefore the vanishing 1-point function is a relevant assumption, at least under the consideration of the isocurvature constraints which is the main observational obstacle for dark photon production of long wavelength modes. –3– The critical values of n for dark photon production are |n| = 2: for |n| > 2, the associated dark electric or magnetic field, or both, goes through exponential amplifications [52]. With the form of I fixed by (2.4), we take the freedom to choose the spectator potential V (σ) so that the dark photon production occurs with nontrivial time dependence. In this regard, we consider a toy model,   σ 4 2 V (σ) = µ tanh , (2.6) Λ2 where µ is a potential energy scale and Λ2 is another energy scale which characterizes the steepness of the potential. This kind of functional form is sometimes considered in the context of the α-attractor model [53], while in our current study we assume that another independent sector drives inflation. To characterize the energy density of the spectator field σ, playing no role in the inflationary expansion, we define the following ratio, R≡ µ4 2 H2  1 . 3MPl (2.7) The Hubble slow-roll parameter H ≡ −Ḣ/H 2 (dot represents derivative with respect to the physical time t) is given by H = φ + σ + γ 0  1 , (2.8) 2 H 2 ) (α = φ, σ, γ 0 ), where ρ and where each contribution is defined as α ≡ (ρα + pα ) /(2MPl α pα are the corresponding energy density and pressure, respectively, and reads φ̇2 σ̇ 2 ,  = σ 2 H2 2 H2 , 2MPl 2MPl      1 1 2 1 0 0 1 0 0 1 2 0 0 02 = 2 2 I F0i F0i + 2 Fij Fij + mγ 0 A 0 + 2 A i A i . 3a2 2a 2 3a MPl H φ = γ 0 (2.9) (2.10) Here, γ 0 denotes a contribution from the backreaction of the dark photon field and we assume that its magnitude is small in comparison with the other contributions. Note that, in (2.10), the expression assumes the physical time (not conformal time), and the bracket h•i denotes the vacuum average. For our analytical convenience, we take the Hubble parameter as a constant value H = const. to solve the dynamics of the spectator-dark photon system. The equation of motion for the spectator field is d2 σ dσ Vσ IIσ µν +3 + 2 =− hF 0 F 0 i , dN 2 dN H 2H 2 µν (2.11) where the subscript σ denotes derivative with respect to σ, and the time dependence of H is assumed to be negligible (|∂N H/H|  3). The right-hand side describes a backreaction effect on the motion of the spectator field. We consider the parameter space where this effect is negligible as we will discuss later. Then, assuming that the first term of the left-hand side in Eq. (2.11) is also negligible, we can use the slow-roll condition 3dσ/dN ' −Vσ /H 2 and obtain from (2.5) 2 MPl tanh(σ̃) Vσ n'− = −2cR . (2.12) 2 2 3H Λ1 Λ1 cosh2 (σ̃) –4– 1 1 V(σ )/μ4 (1) (2) (3) 0 0.0 0 0.5 1.0 1.5 2.0 2.5 σ Figure 1. A schematic dynamics of the spectator field on the tanh2 potential (2.6) in our scenario. (1) The initial value is located toward the top of the hill and the velocity monotonically increases. (2) At a certain time, the spectator field crosses an inflection point of the potential and the velocity gets maximized. (3) The velocity turns to decrease toward the minimum of the potential until inflation ends. Here, we have defined dimensionless quantities, c ≡ Λ1 /Λ2 and σ̃ ≡ σ/Λ2 . Regarding the slow-roll solution of σ, we find that σ̃ and n are described in terms of the Lambert W-function as i hp W (y) , (2.13) σ̃(N ) = arcsinh p W (y) M2 , (2.14) n(N ) = −2cR Pl 2 Λ1 [1 + W (y)]3/2 y ≡ −4R 2 MPl N + 2 ln [sinh(σ̃i )] + sinh2 (σ̃i ) , Λ22 (2.15) where σ̃i = σ̃(0) > 0 is an initial value of σ̃. Inversely, using the definition of the W-function W (zez ) = z, the number of e-foldings is written as     sinh(σ̃i ) Λ22 2 2 + sinh (σ̃i ) − sinh (σ̃) . (2.16) N= 2 ln 2 sinh(σ̃) 4RMPl Note that, since we are requiring the weak coupling of the dark photon to the hidden sector, i.e. I > 1, during inflation, we choose the branch σ > 0, provided Λ1,2 > 0. We are interested in a situation where the magnitude of |n| is initially small but becomes greater than 2 at some point during inflation. Namely, the dark photon production occurs at the intermediate stage of inflation, and its overproduction on large scales is avoidable. To realize such a situation, we consider the following scenario: the spectator field rolls down on the potential and its velocity experiences a peak value at an intermediate time of inflation (see Figure 1). At an initial stage, σ̃ is near the top of the potential and |n| < 2 is satisfied. Then, the velocity of the scalar field gradually increases and gets maximized when it goes through an inflection point σ̃ = σ̃IP on the potential, Vσσ (σ̃IP ) = 2 µ4 1 − 2 sinh2 (σ̃IP ) =0 Λ22 cosh4 (σ̃IP ) –5– ←→ σ̃IP ' 0.658 , (2.17) 1.2 2.2 1.0 |n|  σ 2.0 0.8 1.8 0.6 1.6 0.4 0 10 20 30 40 50 60 10 N 20 30 40 50 60 N Figure 2. The time evolution of σ̃ (left panel) and n (right panel). The solid and dashed lines represent the numerical and slow-roll approximate solutions, respectively. Here, we take R = 10−2 , σ̃i = 1.2525, Λ1 = 5 × 10−3 MPl , c = 7.525 × 10−3 . at a certain number of e-foldings N = NIP , assuming n(NIP ) < −2. After that, the velocity decreases and finally the scalar field reaches the bottom of the potential. We define the number of e-foldings N1 and N2 at which n(N1 ) = n(N2 ) = −2 and N1 < NIP < N2 are satisfied. Using Eq. (2.12), they are obtained as      1 tanh(σ̃p ) sinh(σ̃i ) 2 2 Np = + sinh (σ̃i ) − sinh (σ̃p ) (p = 1, 2) , (2.18) 2 ln 4c cosh2 (σ̃p ) sinh(σ̃p ) where σ̃p = σ̃(Np ). Therefore, the parameter c, characterizing the slope of the potential, partially determines a duration of e-folds for which the particle production occurs. To obtain Np = O(10), we may choose c = O(10−2 ) for σ̃i = O(1). Figure 2 shows the time evolution of σ̃ and that of the corresponding n. We set the inflation end at Nend = 60. The solid and dashed curves represent the numerical and slow-roll approximate solutions (2.13) and (2.14), respectively. One can find that the slow-roll solution is well-fitted to the exact one. At an initial stage of inflation, |n| is smaller than 2 and the vector field is not amplified. At a certain time, |n| becomes greater than 2 and the vector field starts to be amplified. In the next section, we will show that the dark photon production occurs at an intermediate stage of inflation. As a result, the sourced power spectrum of the dark photon exhibits a peaky structure on small scales, reflecting the peaky evolution of n. 2.2 Evolution of the dark photon field Let us analyze the evolution of the dark photon field kinetically coupled to the scalar field σ. Since the vector field possesses a nonzero mass, its spatial components are decomposed into transverse modes and a longitudinal mode: A0i = A0T i + ∂i χ , ∂i A0T i =0 . (2.19) Integrating out the non-dynamical component A00 by A00 = −I 2 ∂ 2 ∂τ χ , −I 2 ∂ 2 + a2 m2γ 0 –6– (2.20) the quadratic action of the transverse modes and that of the longitudinal mode are respectively obtained as Z h i 1 0T 0T 0T 2 2 0T 0T ST = dτ dx I 2 (∂τ A0T ∂ A − ∂ A ∂ A ) − a m A A , (2.21) 0 i j i j i τ i i γ i 2 ! # " Z 1 −I 2 ∂ 2 2 2 ∂τ χ − ∂i χ∂i χ , (2.22) SL = dτ dx a mγ 0 ∂τ χ 2 −I 2 ∂ 2 + a2 m2γ 0 with ∂ 2 ≡ ∂i ∂i . Modulations of the kinetic couplings in these actions are crucial for the evolution of the vector field. To see this, we decompose the transverse modes (2 degrees of freedom) and the longitudinal one (1 d.o.f.) in terms of their Fourier modes, Z  dk  X X Y Y 0T V̂ (τ )e ( k̂) + i V̂ (τ )e ( k̂) eik·x , (2.23) IAi (τ, x) = i i k k (2π)3 Z amγ 0 Ik dk X̂k (τ ) ik·x χ(τ, x) = , (2.24) e , zk ≡ q 3 (2π) zk I 2 k 2 + a2 m2γ 0 where k ≡ |k|, superscripts X and Y label the polarization states, and the orthonormal polarization vectors take the forms eX i (k̂) = (cos θ cos φ, cos θ sin φ, − sin θ) , eYi (k̂) = (− sin φ, cos φ, 0) , (2.25) with the unit wave vector k̂ = (sin θ cos φ, sin θ sin φ, cos θ). They obey the properties: X Y Y X Y σ∗ ki eσi (k̂) = 0 (σ = X, Y ), eX i (k̂)ei (k̂) = ei (k̂)ei (k̂) = 1, ei (k̂)ei (k̂) = 0, ei (k̂) = σ X X Y Y 2 σ ei (k̂), ei (−k̂) = ei (k̂), ei (−k̂) = −ei (k̂). The canonical variables V̂k and X̂k satisfy the reality conditions, X , V̂kX † = V̂−k Y , V̂kY † = V̂−k X̂k† = X̂−k . (2.26) Note that we put i in front of the V̂kY term in (2.23) so that V̂kY respects the above property of Hermitian conjugate. In terms of these variables, the actions (2.21) and (2.22) are rewritten as # " ! Z 2 m2 2I a 0 1 X ∂ dk γ σ† σ† ST = (2.27) dτ ∂τ V̂k ∂τ V̂kσ − k 2 − τ + V̂k V̂kσ , 2 (2π)3 I I2 σ=X,Y # " ! Z a2 m2γ 0 1 dk ∂τ2 zk † † 2 SL = dτ ∂τ X̂k ∂τ X̂k − k − + X̂k X̂k . (2.28) 2 (2π)3 zk I2 Then, decomposing V̂kσ , X̂k into creation/annihilation operators, σ† V̂kσ = Vkσ âσk + Vkσ∗ â−k , X̂k = Xk b̂k + Xk∗ b̂†−k , 2 0 0 [âσk , âσ−k†0 ] = (2π)3 δ σσ δ(k + k0 ) , [b̂k , b̂†−k0 ] = (2π)3 δ(k + k0 ) , (2.29) (2.30) Note that we define the parity change k̂ → −k̂ by the simultaneous operations θ → π − θ and φ → φ + π. –7– with the vacuum |0i defined by âσk |0i = b̂k |0i = 0, we obtain the equations of motion for the mode functions as ! 2 m2 2I a 0 ∂ γ τ ∂τ2 Vk + k 2 − + Vk = 0 , (2.31) I I2 ! 2 m2 2z a 0 ∂ γ k Xk = 0 , (2.32) ∂τ2 Xk + k 2 − τ + zk I2 where the index σ (= X, Y ) of the transverse mode functions has been omitted since the two modes obey the same equation. Historically, for analytical convenience, many of the studies on inflationary kinetic coupling models have assumed a special functional form of I to make n constant in time. Such a case, with a negligible mass term, is briefly summarized in Appendix A. In the case of small mass, the longitudinal mode obeys the equation of motion as in (A.11) and does not experience an exponential enhancement, in contrast to the transverse modes when |n| > 2 is realized. Thus the contribution to the dark photon production from the longitudinal mode is always subdominant in the parameter space of our interest and shall be disregarded in the following consideration. In order to find the solution of the transverse modes with a mild time dependence of n(τ ), we assume a negligible mass of the dark photon, and then Eq. (2.31) is rewritten as   n(τ )(n(τ ) + 1) + dn/dN 2 2 ∂τ V k + k − Vk = 0 , (2.33) τ2 with a constant H. Neglecting a small correction of the velocity term dn/dN , we obtain a differential equation of the same form as Eq. (A.1). Although in the case of time-dependent n(τ ) there would be no analytically closed form of the solution to eq. (2.33) in general, there exists a useful technique called uniform approximation [54–56], which we employ here. By introducing variables   Z 2/3  3 τ∗ p   dτ̃ −g(τ̃ ) , τ < τ∗ , − 2 ν 2 (τ ) τ 2 g(τ ) ≡ (2.34) −k , f (τ ) ≡  Z τ 2/3  3 p τ2    dτ̃ g(τ̃ ) , τ∗ < τ , 2 τ∗ with ν 2 ≡ (1 + 4n(n + 1))/4 and the turning point τ∗ < 0 defined by g(τ∗ ) = 0, the solution of Eq. (2.33) is well described by the following linear combination of Airy functions: VkUA  1/4  1/4 f f ≡ Ak Ai(f ) + Bk Bi(f ) , g g (2.35) where the superscript “UA” is to remind (the leading order of) the uniform approximation. To utilize VkUA as an approximate solution to (2.33), the coefficients Ak and Bk are chosen to realize the adiabatic initial condition in the sub-horizon regime and found to be r π iθ Ak = iBk , Bk = e , (2.36) 2 –8– where θ is an overall phase factor irrelevant in the following discussion. Then, using the asymptotic behavior of Airy functions, the solution on the super-horizon regime is approximately given by ! √   Z τ r 2 Bk 2 3/2 ν (τ̃ ) UA iθ −kτ Vk ' √ 1/4 exp 'e √ exp dτ̃ f − k2 (2.37) 3 τ̃ 2 πg 2νk τ∗ for −kτ  |ν(τ )|. One caution to note is that this leading-order solution of the uniform approximation captures the correct spectral behavior, but its amplitude can be slightly off. Even in the case of constant n, or equivalently constant ν, denoted by ν0 , the comparison to the known exact solution using the Hankel function (see e.g. [30]), Vkexact , reveals the difference VkUA eiθ ' √ e−ν0 k  2ν0 −kτ ν0 − 1 2 , Vkexact eiθ 2ν0 −1 Γ(ν0 ) '√ √ 1 , k π (−kτ )ν0 − 2 √ ν0 − 1 VkUA 2π ν0 2 ' 0.95 , ' ν0 Vkexact e Γ(ν0 ) (2.38) where the last numerical value is evaluated for the case n0 = −2, or ν0 = 3/2. This is due to the error of the truncation at the leading order. The uniform approximation can accommodate the calculations of higher-order terms in an iterative manner for an arbitrary ν(τ ) [57]. UA(2) For constant ν = ν0 , for example, the first sub-leading correction gives Vk = VkUA /(12ν0 ) UA(2) [55], for which the difference from the exact result becomes |(VkUA +Vk )/Vkexact | ' 0.9999 for ν0 = 3/2. Nonetheless, our leading-order solution (2.35), and its super-horizon limit (2.37), gather the correct time evolution and the spectral feature, as can be seen in Fig. 3. We thus use the leading order and admit the ∼ 5 % error in the amplitude in the following considerations. 3 The energy density spectrum Based on the nontrivial background dynamics and the resulting characteristic production of dark photons described in the previous sections, we compute their energy density and show that its power spectrum can be peaked at scales much smaller than those of CMB observations but with wavelengths large enough to realize small-mass cold dark matter. The energy density ργ 0 is the (0,0)-component of the energy-momentum tensor of the dark photon A0 = −2δS[A0 ]/δg µν , which is calculated as Tµν ργ 0 = −T A0 0 0      1 1 0 0 2 0 0 2 2 02 0 0 = 4 I F0i F0i + Fij Fij + a mγ 0 A 0 + A i A i , 2a 2 (3.1) where this expression is compatible with the conformal time (not the physical one). The vacuum averaged energy density is split into two parts:3 3 The energy density is quadratic in the field Aµ , and there in principle exist cross terms between the transverse and longitudinal modes. They do not contribute to the vacuum average for two reasons: the first is that the different polarization modes are decoupled from each other at the linear order, and thus their cross correlations vanish. The second is that the cross terms always appear as total (spatial) derivatives, which is a consequence of the background isotropy and homogeneity, together with a vanishing vector vev hAµ i = 0, and therefore vanish when vacuum average is taken. –9– n<-2 107 2 k 3 |Ek |/H2 2 k 3 |Ek |/H2 10 5 1000 10 0.100 0.001 10 -5 10 7 n<-2 1000 10 0.100 10 -14 10 -9 x 10 -4 10 10 -5 1000 10 0.100 0.001 0.001 10-19 n<-2 10 5 105 2 k 3 |Ek |/H2 10 7 10 -14 10 -10 10-6 x 0.01 100.00 10 -5 10-4 0.001 0.010 0.100 1 10 100 x Figure 3. The time evolution of the electric field (3.8) with a mode exiting the horizon at N = 10 (left panel), N = 25 (middle panel) and N = 50 (right panel). The horizontal axis x ≡ −kτ is a dimensionless time flowing from the right to the left. The blue solid lines denote the exact numerical solutions. Regarding the initial conditions, we took Bunch-Davies vacuum and set the initial time variable as xini = 102 . The red dashed lines are the solutions with the uniform approximation (3.10) which start to fit at around the horizon crossing N = N∗ . The blue dashed lines are the approximate Gaussian fitting functions (5.23) discussed in section 5. The black dotted lines denote a time when n(t) = −2. We have used the same parameter set as in Figure 2. hργ 0 i = hργ 0 ,T i + hργ 0 ,L i , " ! #    ∗ Z a2 m2γ 0 Vk dk Vk 1 2 2 ∗ ∂τ + k + I ∂τ Vk Vk hργ 0 ,T i = 4 a (2π)3 I I I2 Z ≡ d ln k Pγ 0 ,T (k) , # " !    ∗ Z 2 m2 a 0 X X 1 dk γ k k hργ 0 ,L i = 4 z 2 ∂τ Xk Xk∗ ∂τ + k2 + 2a (2π)3 k zk zk I2 Z ≡ d ln k Pγ 0 ,L (k) , (3.2) (3.3) (3.4) where we have defined the power spectrum of the energy density for each mode, i.e. Pγ 0 ,T /L . Note that the two transverse polarization modes are summed over. The gradient and mass terms are sub-dominant on the super-horizon scales in comparison with the kinetic terms, and we neglect their contributions to the energy density. Let us first evaluate the power spectrum for the longitudinal mode. In our focused parameter space, p  M always holds and therefore the solution of Xk is given by Eq. (A.13). In this case, the spectral shape is blue tilted and only a few UV modes contribute to the energy density of the longitudinal mode at the inflation end N = Nend . Its magnitude is evaluated as [27] H4 hργ 0 ,L i|Nend ∼ 2 , (3.5) 8π and found to be much smaller than the energy density of the transverse modes as we will see below. In considering the transverse modes on the occasion of negligible mass, it is convenient 2 to define the corresponding (dark) electric and magnetic fields as EiT = −I∂τ A0T i /a and T 0T 2 Bi = Iijk ∂j Ak /a , where ijk is the Levi-Civita symbol in the flat spacetime. Their mode – 10 – functions can be defined as, 4 I d Ek ≡ − 2 a dτ  Vk I  , Bk ≡ k Vk . a2 (3.8) To evaluate these quantities, we use the result of the uniform approximation for Vk obtained in (2.35), or the corresponding super-horizon expression (2.37). In the same way as in Appendix A, the dark electric field Ek is dominant over the magnetic counterpart Bk in the case where the dark photon production occurs in the branch n < 0, and we neglect Bk in comparison to Ek in the following discussions. In order to calculate Ek under the uniform approximation, we define an “averaged” value of quantity ν(N ) by RN   σ̃(N ) − σ̃(N∗ ) 1 N − N∗ N∗ dÑ ν(Ñ ) ν̄k (N ) ≡ , (3.9) =− + N − Nk c(N − Nk ) 2 N − Nk where N∗ ≡ Nk − ln(ν) is the number of e-foldings at which g(τ (N∗ )) = 0 and the definitions of I(σ) and n, respectively (2.4) and (2.5), have been used. Then the following electric mode function obtained by the uniform approximation (2.35) and (2.37),   I d VkUA UA Ek ≡ 2 a dτ I    √ 2 3 iθ 3H 2 ν 'e √ exp ν̄k − (N − Nk ) , (−kτ → 0) , (3.10) 2 2k 3 3 becomes a good approximation of the exact solution. Figure 3 shows the time evolution of the electric mode function for a few different momenta. Recalling the e-folding times N1 and N2 defined below (2.17) such that n(N1 ) = n(N2 ) = −2 with N1 < N2 , the left panel of Fig. 3 shows an evolution of the electric field exiting the horizon before N = N1 . For a while after leaving the horizon, its magnitude is suppressed until |n| becomes greater than its critical value 2. After |n| exceeds 2, it starts to grow on the super-horizon regime. Its growth persists until N = N2 and after that the amplitude starts to decrease since |n| becomes smaller than 2 again. The middle panel gives an evolution of the mode function exiting the horizon just around |n| = 2. The amplitude of the electric field with momentum modes of around this scale is mostly enhanced. The right panel shows an evolution of the mode function exiting the horizon when |n| is already greater than 2. While it starts to grow, its magnitude is not maximally enhanced because the fluctuation is still in the sub-horizon for a while after |n| crosses 2. We can see that the expression (3.10) is well fitted to the numerical solution. We now evaluate the time evolution of the power spectrum of the dark photon energy density with a certain model parameter set. The result is shown in Figure 4. The spectrum has a peak at an intermediate scale where the index |n| becomes greater than 2. This peak scale roughly corresponds to the number of e-foldings Nk = N1 , at which |n| becomes equal 4 With Ek and Bk , we can express EiT and BiT as Z    i dk ik·x h ∗ X† X Y ∗ Y† Y EiT = e Ek âX k + Ek â−k ei (k̂) + i Ek âk + Ek â−k ei (k̂) , 3 (2π) Z   i  dk ik·x h T X ∗ X† Y Bi = e Bk âYk + Bk∗ âY−k† eX i (k̂) + i Bk âk + Bk â−k ei (k̂) , 3 (2π) (3.6) (3.7) Y Y X where we have used ijk kj eX k (k̂) = kei (k̂) and ijk kj ek (k̂) = −kei (k̂). Be alert for the mixing of the polarizations {X, Y } between the operators and polarization vectors in the expression of BiT . – 11 – γ',T /H4 105 10 0.001 10 -7 0.01 100.00 106 10 10 1014 1018 10 22 -1 k[Mpc ] Figure 4. The power spectrum of the vector field Pγ 0 ,T at e-folds N = N1 (dotted), N∗ (dashed), N2 (dot-dashed), Nend (solid). The blue and red lines denote the numerical and approximate solutions, respectively. The rapid increase for small scales represents the Bunch-Davies vacuum which should be renormalized as a UV contribution and we cut off its growth in our calculation. to 2 for the first time during the evolution. On the other hand, the amplitudes on large scales are suppressed because of a finite interval with |n| < 2 on the super-horizon regime (see Figure 2). We will see in the next section that this feature makes it possible to avoid the overproduction of isocurvature perturbation. The present calculation assumes that the backreaction of the dark photon field on the background motion of the spectator field is negligible. As we will discuss in the next section, this assumption is justified in our focused parameter region. 4 Dark photon DM: constraints and results In this section, we estimate the relic abundance of the dark photon DM and discuss theoretical and observational constraints on the current scenario. Then, the viable parameter space of the dark photon DM is shown. 4.1 Relic abundance Let us first derive the relic abundance of the dark photon DM, Ωγ 0 = hργ 0 it=t0 , ρ(t0 ) (4.1) where the average density of the dark photon is evaluated at the present time, denoted by 2 H 2 . To subscript 0, and ρ(t0 ) is the critical density of our present universe: ρ(t0 ) = 3MPl 0 ensure the produced dark photon becomes non-relativistic some time after inflation ends at t = tend , we compare the dark photon mass mγ 0 and a physical, time-dependent momentum scale with a comoving wave number at which the spectrum is peaked, qpeak (t) ≡ kpeak /a(t). Using kpeak = a(tpeak )H(tpeak ), where the time tpeak is defined by this relation and is taken to be during inflation, and defining a time duration of e-folds ∆N ≡ log(a(tend )/a(tpeak )), the physical momentum scale qpeak (t) is written as qpeak (t) = H(tpeak )e−∆N – 12 – a(tend ) a(treh ) , a(treh ) a(t) (4.2) where treh denotes the time when the reheating completes. For simplicity, we assume an instantaneous reheating, a(tend ) = a(treh ), where the inflationary energy scale is related to the reheating temperature Treh as 2 ρ(tend ) = 3MPl H(tend )2 = π2 4 g∗ (treh )Treh . 30 (4.3) Here, g∗ is the number of relativistic degrees of freedom. The momentum scale of the dark photon at the reheating period is then evaluated as   30   e Hinf , (4.4) qpeak (treh ) ∼ 0.1 GeV 12 10 GeV e∆N where Hinf denotes the inflationary Hubble scale. In our scenario, we assume a high-scale, almost de-Sitter inflation and H(tpeak ) ' H(tend ) ' Hinf . Therefore, the dark photon with mass mγ 0  eV in our interest is relativistic right after inflation ends. Then, the current energy density of the dark photon evolves from the end of inflation as  hργ 0 it=t0 = hργ 0 it=tend areh a(tNR ) 4  a(tNR ) a0 3 , (4.5) where tNR is the time when the dark photon becomes non-relativistic, determined by qpeak (tNR ) = mγ 0 . (4.6) By using the entropy conservation law, the present abundance of the dark photon DM is obtained as  2 3/2 3 mγ 0 T03 hργ 0 it=tend π gS (t0 )g∗ (treh )3/2 Treh 1 Ωγ 0 = . (4.7) 3 2 2 4 90 gS (treh ) −kpeak τend MPl 3MPl H0 Hinf Therefore, the mass of the dark photon is evaluated as mγ 0 ' 8.2 × 10 −11  eV Ωγ 0 h2 0.14  107 4 hργ 0 it=tend /Hinf  e30 e∆N  1015 GeV Treh 3 , (4.8) where we have used H0 = 2.133 × 10−42 h GeV, T0 = 2.725 K = 2.348 × 10−13 GeV, MPl = 2.435 × 1018 GeV and gS (t0 ) = 3.91, and assumed g∗ (treh ) = gS (treh ) = 106.75. From this expression, we can observe that, in order for the dark photon that explains the total dark matter abundance to have a smaller mass, either its energy density ργ 0 during inflation , ∆N or the reheating temperature Treh must be larger. The choice of these parameters in (4.8) is rather an optimistic one, in favor of small mass, and thus it gives a rough lower bound for the dark photon mass mγ 0 in our model. Yet we remark that it is sensitive to the values of ∆N and Treh , and a mass smaller than ∼ 10−10 eV by a few orders of magnitude is still feasible. 4.2 Constraints Next, we discuss relevant theoretical and observational constraints on the present scenario. – 13 – (i) Isocurvature mode Let us consider the isocurvature mode of the dark photon in our scenario, as it could be one of the stringent constraints on the inflationary production of dark photon DM [47]. The isocurvature fluctuation of the dark photon is defined as the entropy perturbation caused by the non-adiabatic mode: δργ 0 ργ 0 − hργ 0 i S≡ = . (4.9) hργ 0 i hργ 0 i Then, in terms of the Fourier decomposition, Z ργ 0 (x) dk ≡ δ̂γ 0 (k) eik·x (4.10) hργ 0 i (2π)3 its power spectrum is given by Z i dkdk0 h 6 0 0 ik·x+k0 ·y 0 (k) δ̂γ 0 (k )i − (2π) δ(k) δ(k ) e hS(x)S(y)i = h δ̂ γ 6 (2π) Z dk sin(k|x − y|) = PS (k) , k k|x − y| (4.11) (4.12) where the power spectrum PS is defined through 2π 2 PS (k) . (4.13) k3 Note that the disconnected contribution is explicitly subtracted in the expression (4.13). Since the transverse electric mode is energy-dominant, δ̂γ 0 (k) is approximately given by Z   dp  X X 1 Y \ X \ Y Y Y X e ( k − p) . e ( k − p) + i Ê Ê e ( p̂) + i Ê e ( p̂) Ê δ̂γ 0 (k) ' p i p i k−p i k−p i 2hργ 0 ,T i (2π)3 (4.14) X† ∗ X X A hat on the electric field manifests that Êp ≡ Ep âp + Ep â−p is an operator. Then, noting Y \ Y X \ 5 eX i (p̂)ei (k − p) = ei (p̂)ei (k − p) = 0, the power spectrum is evaluated as Z 2 1 X 1 dp∗ \ PS (k) ' eσi (p̂)eσi (k − p) |E|p|∗ |2 |E|k−p|∗ |2 , (4.15) 2 2 3 4 π hργ 0 ,T i (2π) hδ̂γ 0 (k) δ̂γ 0 (k0 )i − (2π)6 δ(k) δ(k0 ) = (2π)3 δ(k + k0 ) σ=X,Y where we have introduced the dimensionless momenta: p∗ ≡ p/k, |k − p|∗ ≡ |k − p|/k. Note that the disconnected term in hδ̂γ 0 (k)δ̂γ 0 (k0 )i, which is proportional to δ(k)δ(k0 ), exactly cancels with the corresponding term in (4.13), and thus PS takes into account only the connected diagrams. The amount of the isocurvature perturbation is tightly constrained at CMB scales. The current Planck satellite observation puts a constraint on the magnitude of the isocurvature power spectrum evaluated at three different scales k = 0.002, 0.05, 0.1 Mpc−1 . As shown in Figure 4, the resultant spectral shape is blue-tilted at CMB scales in our model. To constrain the amplitude of the power spectrum, we make use of the limits evaluated at the lowest and highest scales, k = 0.002 Mpc−1 and k = 0.1 Mpc−1 , which are given by [58] PS |k=0.002Mpc−1 . 0.8 × 10−10 , PS |k=0.1Mpc−1 . 2.2 × 10−9 . (4.16) The spectrum in our model has a peaky feature, and thus imposing the above observational isocurvature constraints on the peak value should provide a conservative upper bound on the dark photon production. 5 X \ Y This is because the vectors eX i (p̂) and ei (k − p) are both on the plane defined by k and p, and ei (p̂) Y \ and ei (k − p) are both perpendicular to it, under our construction of the polarization vectors. – 14 – (ii) Backreaction The Friedmann equation, 1 1 2 3MPl H 2 = φ̇2 + U (φ) + σ̇ 2 + V (σ) + hργ 0 i , 2 2 (4.17) where U (φ) denotes the inflaton potential, includes the contribution from the energy density of the produced dark photon. Our calculations in this paper relies on the assumption that the inflationary quasi-de-Sitter background is driven only by φ, and the produced dark photon has negligible impact on the dynamics of σ, which is already subdominant to φ. To neglect the backreaction contribution, we need hργ 0 i  1 2 2 σ̇  3MPl H2 , 2 (4.18) during inflation. This constraint is not tighter than that from the backreaction to the motion of the spectator field in Eq. (2.11). Namely, it suffices that the following condition must be satisfied: 2hργ 0 ,T i dσ IIσ 0 0µν hF F i , (4.19) 3  − ' µν dN 2H 2 H 2 Λ1 where in the most right-hand side we have ignored the contribution of the longitudinal mode. To characterize this hierarchy, we define a ratio, Rb ≡ 2hργ 0 ,T i , 3HΛ1 σ̇ (4.20) and consider the parameter region realizing Rb  1. By using σ̇/(HΛ1 ) ' n in Eq. (2.12) and the condition of slow-roll parameters (2.8) with (2.9), σ < H ←→ H Λ21 2 < n2 , 2MPl (4.21) on CMB scales, we can translate the condition Rb  1 to the following constraint on the magnitude of the dark photon energy density: hργ 0 ,T i 3|n| Λ21 3|n|  < 2 2 ∼ 2 × 107 |n| , 4 2 H 2 H 8π nCMB Pζ,CMB (4.22) where the power spectrum of the curvature perturbation in the vacuum state is given by Pζ,CMB = H2 ' 2 × 10−9 . 2 8π 2 MPl H (4.23) Note that in (4.22) we distinguish n from nCMB , the latter denoting the value of n at the time when the CMB modes exit the horizon. In our mechanism, the value of n changes during inflation by an O(1) amount, and specifically we demand |nCMB | < 2 in order to avoid overproduction of dark photon spoiling the CMB predictions. Typical values of nCMB in our considerations are ∼ 1, which is taken for the final evaluation in (4.22). – 15 – (iii) Non-relativistic time vs. equality time After inflation ends, the dark photon behaves as a relativistic radiation component since the scale of the physical momentum at the spectral peak qpeak is much larger than the dark photon mass mγ 0 . For the dark photon to behave as a viable dark matter candidate, it must become non-relativistic by the time of matter-radiation equality t = teq . We define the cosmic temperature TNR at time t = tNR when the scale of the physical momentum equals to the dark photon mass, kpeak qpeak (tNR ) = = mγ 0 . (4.24) a(tNR ) Using the entropy-conservation law, gS a3 T 3 = const., the temperature TNR is obtained as6   gS (treh ) 1/3 Treh qpeak (treh ) gS (tNR )  1/4     mγ 0 ρ(tend ) 1/12 ρ(tend )1/4 gS (treh ) 1/3 30 , = π 2 g∗ (treh ) H(tend ) ρ(treh ) −kpeak τend gS (tNR ) TNR = mγ 0 (4.25) assuming that the equation of state is of matter domination after inflation and before reheating. On the other hand, the temperature Teq at the equality time is determined by solving the energy density equality between the radiation and matter components, Ωr Ωm = , a(teq )4 a(teq )3 (4.26) with the normalization of the scale factor now taken at present, a(t0 ) = 1. Then, the entropy2 H 2 ) give conservation law and Ωr = π 2 g∗ (t0 )T04 /(90MPl 0 Teq = 2 H2 90Ωm MPl 0 . π 2 g∗ (t0 )T03 (4.27) Imposing TNR > Teq , we obtain the following constraint: TNR (π 2 /90)1/2 g∗ (t0 ) mγ 0 T03 = e∆N Teq Ωm g∗ (treh )1/2 Treh MPl H02  gS (treh ) gS (tNR ) 1/3 >1, (4.28) where the instantaneous reheating is assumed. Henceforth we take gS (teq ) = gS (t0 ) = 3.909, g∗ (teq ) = g∗ (t0 ) = 3.363. Defining ∆N = ∆Neq at which TNR = Teq , Eq. (4.28) is rewritten as      m 0  Ωm h2 Treh γ ∆N > ∆Neq ' 10.5 + ln + ln − ln . (4.29) 0.143 1013 GeV 10−10 eV We have by now collected all the constraints that are to be imposed on the amount of the production of dark photon dark matter in our model. 6 Ref. [27] has missed the factor (gS (treh )/gS (tNR ))1/3 , while the effect is small. – 16 – 4.3 Results We numerically find a viable parameter region of the dark photon DM evaluated in (4.8) consistent with the constraints discussed above, which is shown in Figure 5. To make this plot, we have solved the (massless approximate) equation of motion for dark photon, (2.33), together with the equation of motion for σ, (2.11), with constant H, with several values of σ̃i and c, and used the resultant solution of ργ 0 ,T shown in Figure 4. The initial field range σ̃i characterizes a timing when |n| becomes greater than 2. For a small σ̃i , |n| crosses 2 at an early stage of inflation, which leads to a large ∆N and correspondingly a small mγ 0 for the dark photon DM. However, the effect of the isocurvature mode is severe for this case because the dark photon is amplified on large scales. For a large σ̃i , |n| crosses 2 at a late stage of inflation, which leads to a large mγ 0 . The effect of backreaction tends to reduce in this case because N2 becomes closer to Nend or even exceeds it, and therefore the time interval of particle production becomes shorter. On the other hand, the parameter c is related to the steepness of the potential slope. As c increases, the slope of the potential becomes steeper and the value of |n| gets larger. As a result, the effect of the backreaction and/or the contribution of the isocurvature mode become severer. However, a large c and a small σ̃i tend to derive the dynamics of |n| damping so early and therefore predict a small hργ 0 ,T i leading to a large mγ 0 . To find a viable region for the light dark photon DM, a high-scale inflation would be preferable because a high reheating temperature is required in Eq. (4.8). We can see from the figure that the present scenario predicts a mass window mγ 0 & 10−13 eV. Due to the exponential sensitivity to the variation of model parameters, the viable parameter space leading to a preferable range of mass mγ 0 for light dark photon is localized in a small region on the σ̃i - c plane. The figure also describes the constraint that the dark photon must be non-relativistic until the equality time. We find that in the region of our interest the condition (4.29) is satisfied. In our study, we have assumed that the Hubble parameter H is constant in the whole period of inflation and therefore neglected its dynamics near and after the end of inflation. However, when we take it into account, we might expect that the allowed region of the dark photon mass would get wider, and even a smaller mass would be available. This is because the amplitude of the energy density of the dark photon is enhanced due to the rapid motion of the spectator field. However, the backreaction effect may become dominant and the system can be completely non-linear. Such a near-end and post-inflationary amplification is expected to be severe if the potential of σ becomes flat after inflation, unlike our choice of V (σ) as in (2.6). The analysis of such a non-linear system is beyond the scope of the present work. 5 Generation of tensor modes In this section, we evaluate the power spectrum of tensor modes sourced by the dark photon field during inflation. This effect is an inevitable consequence of the dark photon production and is potentially led to observable signals, which we would like to evaluate in this section. The tensor perturbation is given by fluctuations of the spacial components of the metric, gij (t, x) = a(t)2 (δij + 12 hij (t, x)), which obey the following equation of motion at the leading order:   ∇2 4 2 ∂t + 3H∂t − 2 hij ' − 2 Πlm El Em , (5.1) a MPl ij – 17 – 0.0084 c (potential steepness) c (potential steepness) 0.0084 0.0082 10-18 0.0080 0.0078 10-15 0.0076 0.0074 10-12 30 0.0082 20 0.0080 10 0 0.0078 40 0.0076 35 30 25 0.0074 20 15 0.0072 0.0072 10-9 1.20 1.25 10-6 10-3 1.30 1.35 1 1.40 1.45 1.20 1.50 1.25 1.30 1.35 1.40 1.45 1.50  σ i (initial field range)  σ i (initial field range) Figure 5. Left panel : contour plots of the dark photon DM with mass mγ 0 = 1, 10−3 , 10−6 , 10−9 , 10−12 , 10−15 , 10−18 eV (black lines) in the parameter space of σ̃i and c. Right panel : contour plots of ∆Neq = 0, 10, 20, 30 (black solid lines) and ∆N = 15, 20, 25, 30, 35, 40 (black dashed lines). In both panels, we set R = 10−2 , Λ1 = 5 × 10−3 MPl and rv = 5 × 10−4 . The shaded regions are constrained by the backreaction, Rb > 0.1 (light red) and Rb > 1 (red), and the isocurvature perturbation with k = 0.002 Mpc−1 (dark green) and k = 0.1 Mpc−1 (purple). where Πlm ij is the transverse-traceless projector defined by 1 l m lm Πlm , ij ≡ Πi Πj − Πij Π 2 Πij ≡ δij − ∂i ∂j . ∇2 We decompose hij into the linear polarization tensors in Fourier space, Z X Z dk dk ik·x ĥ (k, t)e = es (k̂)ĥsk (t)eik·x , hij (t, x) = ij 3 ij (2π)3 (2π) s=+,× +/× where the transverse-traceless polarization tensors eij ucts of the polarization vectors (5.2) (5.3) (k̂) are given by the following prod- X/Y ei (k̂):  1  X X Y Y √ = ei (k̂)ej (k̂) − ei (k̂)ej (k̂) , 2  i  X e× ei (k̂)eYj (k̂) + eYi (k̂)eX , j (k̂) ij (k̂) = √ 2 e+ ij (k̂) (5.4) (5.5) so that they satisfy the property esij∗ (−k̂) = esij (k̂). Then, by using ĥsk = esij∗ (k̂) ĥij (k̂) and s s Πlm ij elm (k̂) = eij (k̂), we obtain  ∂x2  2 4a3 + 1 − 2 (aĥsk ) = − esij∗ (k̂) 2 2 x k MPl Z  dp  X X Y Y × Ê e ( p̂) + i Ê e ( p̂) p i p i (2π)3   X Y Y \ \ × Êk−p eX ( k − p) + i Ê e ( k − p) . j k−p j – 18 – (5.6) on the de Sitter assumption. The solution to this equation is a linear combination of two components, ĥsk = ĥsk,v + ĥsk,s . The homogeneous solution ĥsk,v is the usual Banch-Davies vacuum mode and its mode function is given by, at the leading order in slow roll,   e−ikτ i s √ ahk,v = , (5.7) 1− kτ MPl 2k and the corresponding tensor-to-scalar ratio reads rv ≡ 8Ph,v = 16 H . Pζ,CMB (5.8) The numerical factor 8 in the intermediate step comes from our definition hij = 2δgij /a2 and from the 2 polarization states of the tensor modes. On the other hand, the peculiar solution ĥsk,s is sourced by the second order of the dark photon field. We find the mode function by using the Green function’s method, Z ∞ 4 dy a3 (y) GR (x, y) aĥsk,s = − 2 2 esij∗ (k̂) k MPl −∞ Z   dp  X Y \ X \ Y X Y Y X × (y) e ( k − p) , (y) e ( k − p) + i Ê Ê (y) e ( p̂) + i Ê (y) e ( p̂) Ê j j p i p i k−p k−p (2π)3 (5.9) where GR (x, y) ≡ Θ(y − x) [(x − y) cos(x − y) − (1 + xy) sin(x − y)] /(xy) ' −Θ(y − x)(x3 − y 3 )/(3xy) is the retarded Green function, where the last approximate equality is in the limit of x ≡ −kτ  1 and y ≡ −kτ 0  1, with τ 0 an auxiliary time variable. This approximation is valid thanks to the fact that our mechanism of dark photon amplification takes place on the super-horizon scales. Regarding the products of the polarization tensors and vectors, we assume k is directed in the ẑ axis and use the following identities: cos θp̂ cos θk−p [ cos 2φp̂ √ , 2 cos 2φp̂ Y Y \ √ e+∗ , ij (k̂)ei (p̂)ej (k − p) = 2 cos θp̂ sin 2φp̂ X Y \ √ e+∗ , ij (k̂)ei (p̂)ej (k − p) = 2 cos θk−p [ sin 2φp̂ Y X \ √ e+∗ ( k̂)e ( p̂)e ( k − p) = , i j ij 2 cos θpb cos θk−p [ sin 2φp̂ X X \ √ e×∗ ( k̂)e ( p̂)e ( k − p) = i , i j ij 2 sin 2φp̂ Y Y \ e×∗ , ij (k̂)ei (p̂)ej (k − p) = √ 2i cos θp̂ cos 2φp̂ X Y \ √ e×∗ , ij (k̂)ei (p̂)ej (k − p) = i 2 cos θk−p [ cos 2φp̂ Y X \ √ e×∗ ( k̂)e ( p̂)e ( k − p) = i , i j ij 2 X X \ e+∗ ij (k̂)ei (p̂)ej (k − p) = − – 19 – (5.10) (5.11) (5.12) (5.13) (5.14) (5.15) (5.16) (5.17) \ with cos θp̂ = k̂ · p̂ and cos θk−p [ = k̂ · k − p. Then, the dimensionless power spectrum of tensor modes is defined as 0 0 0 hĥsk ĥsk0 i = hĥsk,v ĥsk0 ,v i + hĥsk,s ĥsk0 ,s i 0 ≡ (2π)3 δ ss δ(k + k0 )
2π 2 ss Ph,v (k) + Ph,s (k) , 3 k where the cross terms vanish for linear perturbations. One can find Z Z τ dp 8k 3 dτ 0 x3 − y 3 ++ ++ \ Ph,s (k)|τ = 2 4 4 F ( p̂, k − p) Ep E|k−p| 3 0 (2π) 3y 3 π H MPl τmin τ Z Z τ 8k 3 dp dτ 0 x3 − y 3 ×× ×× \ Ph,s (k)|τ = 2 4 4 F ( p̂, k − p) Ep E|k−p| 0 (2π)3 3y 3 π H MPl τmin τ (5.18) 2 , (5.19) , (5.20)     2 2 2 F ++ ≡ cos2 (2φp̂ ) cos2 θp̂ cos2 θk−p , [ + 1 + sin (2φp̂ ) cos θp̂ + cos θk−p [     2 2 2 F ×× ≡ sin2 (2φp̂ ) cos2 θp̂ cos2 θk−p , [ + 1 + cos (2φp̂ ) cos θp̂ + cos θk−p [ (5.21) 2 where (5.22) and we have approximated the lower bound of the time integral by |τmin | = min(1/p, 1/|k − p|), because the dark photon can grow only after the horizon crossing and we focus on the contribution from the super-horizon modes. For our convenience, we perform the numerical evaluation of the time integral and that of the momentum integral separately. To do this, we define the following Gaussian fitting function for the electric component: " # 1 dn H2 Apeak (k) exp − (ln(τ /τ2 ))2 , (5.23) Ek (x) ' Ek,fit (x) ≡ √ 2 dN N =N2 2k 3 where Apeak (k) is the maximum amplitude and τ2 is the conformal time at which N = N2 . Plugging Eq. (5.23) into Eqs. (5.19) and (5.20), at the super-horizon limit τ → 0, i.e. x  y  1, we find ++ ×× Ph,s (k)|τ →0 = Ph,s (k)|τ →0 Z  A2 (p)A2 (|k − p|)
 dp∗ d cos θp̂ 1 F 2H 4 peak peak 2 2 1 + cos θ ' 2 1 + cos θ , p̂ [ 4 k−p π 9MPl (2π)2 p∗ |k − p|3∗ (5.24) where p∗ ≡ p/k , |k − p|∗ ≡ |k − p|/k and F is the numerical factor obtained by the time integration, " # !−1/2 Z ∞ √ dτ 0 dn dn exp − (ln(τ 0 /τm ))2 = π . (5.25) F≡ 0 dN N =Nm dN N =Nm −∞ τ Here, Nm (τm ) is the time when the amplitude of the electric component is maximized. Note that the time integration can be extended to ±∞ (−∞ → 0 in the terms of τ 0 ) since the integral has its support almost around the peak of the electric mode function. The cross power +×/×+ spectrum Ph,s (k) vanishes due to the cancellation in integrating the periodic function sin(2φp̂ ) cos(2φp̂ ). – 20 – ΩGW h2 10-12 10-15 10-18 10-11 10 -9 10 -7 10 -5 0.001 0.100 10 f[Hz] Figure 6. The GW spectrum sourced by the dark photon field during inflation. The sensitivity curves of SKA [51] (gray), LISA [60] (pink), DECIGO [61] (orange), BBO [61] (green) and µAres [50] (purple) are also shown. We here take R = 10−2 , Λ1 = 5 × 10−3 MPl , rv = 5 × 10−4 , σ̃i = 1.2525, c = 7.525 × 10−3 , which correspond to mγ 0 ' 1.64 × 10−13 eV. Let us now evaluate the logarithmic energy density of GWs at present, ΩGW (k) ≡ 1 dρGW , ρc d ln k (5.26) 2 H 2 is the present critical energy density of the Universe. Using the entropy where ρc = 3MPl 0 conservation law, ΩGW is related to the power spectrum of primordial tensor modes as [59],7  g −1/3   ∗ ++ ×× ΩGW (k)h2 ' 6.85 × 10−7 Ph,v (k) + Ph,s (k) + Ph,s (k) . (5.27) 100 Figure 6 depicts the power spectrum of the GW energy density. On large scales, the contribution from the vacuum modes is dominant and the spectrum is scale-invariant. At intermediate scales, however, its magnitude gets amplified by the contribution of the sourced tensor modes and has a peak at a frequency around µHz. This spectral shape is determined by the background time evolution of the index n. Even if the vacuum tensor-to-scalar ratio is small, the sourced power spectrum is potentially testable with the next generation pulsar timing array measurement (SKA) and the projected space-based laser interferometers such as DECIGO, BBO and µAres. On the other hand, it would be challenging to test it by LISA. 6 Conclusion We have explored a mechanism to produce a light dark photon DM through a coupling between the dark photon field and a spectator scalar field that does not play a role in the inflationary expansion of the Universe but is rolling down its potential during the inflation. The 7 The factor 4 difference in the numerical factor compared to the expression in [59] comes from the fact that our definition of hij contains 1/2, as seen right above (5.1), which amounts to 1/4 in terms of Ph . – 21 – motion of the spectator field efficiently produces dark photons with large wavelengths which become non-relativistic before the time of matter-radiation equality. We have constructed the mechanism such that the spectrum of wavelengths is peaky so that the constraint from the isocurvature perturbation can be evaded. The correct relic abundance is then achieved over a wide range of the dark photon mass. Depending on the model parameter set, we have found that our result could provide the dark photon DM with mass down to mγ 0 ≈ 10−13 eV. Our mechanism favors high-scale inflation models which can be tested in future observations. Furthermore, fluctuations of the dark photon field during inflation produce GWs detectable at future space-based interferometers. It would be interesting to consider the evolution of the dark photon energy density after inflation. At the end of inflation, the speed of the spectator field increases and the provided dark photon may backreact to the evolution. This effect would more or less modify the allowed parameter region of the dark photon DM. Our result is sensitive to both the amount of the dark photon production and the reheating temperature, and therefore, if there were to be an efficient enhancement around or after the end of inflation, a mass of the dark photon DM that is smaller than the one in our current study may be within the reach. This, however, likely requires non-linear studies of the coupled system between the spectator scalar and the dark photon and is beyond our current scope. Regarding to the generation of the dark photon, we have assumed no homogeneous component of the dark photon because our scenario does not generate the dark photon on large scales. In this sense, the resultant power spectrum is statistically isotropic on large scales. However, it might be worthwhile to prove the statistical anisotropy of tensor modes on intermediate scales, which would serve as another venue of interesting signatures of the primordial Universe. We leave these issues to a future study. Acknowledgements We would like to thank Eiichiro Komatsu for fruitful discussions. YN is supported by Natural Science Foundation of China under grant No. 12150610465. RN is in part supported by RIKEN Incentive Research Project grant. IO acknowledges the support from JSPS Overseas Research Fellowship and JSPS KAKENHI Grant No. JP20H05859 and 19K14702. A Dark photon production with constant n In this appendix, we briefly summarize the calculation of dark photon production for constant n, denoting it by n0 . The existing calculations in the literature can be found in, e.g., [27, 30, 52]. In this case, I simply behaves as a function of the scale factor, I ∝ an0 , and Eq. (2.31) is reduced to " # 2 m2 a 0 n (n + 1) γ 0 0 ∂τ2 Vk + k 2 − + Vk = 0 . (A.1) τ2 I2 This equation is derived under the de-Sitter approximation, i.e. a ' −1/(Hτ ). If the third term in the square parentheses of the above equation is negligible in comparison with the p other terms, an instability could occur on the super-horizon regime −kτ ≤ n0 (n0 + 1) for a certain range of n0 . To realize this instability, we consider a negative branch of the index n0 < 0, where I becomes a decreasing function in time [27]. 8 The condition of the instability 8 This branch is also appropriate for avoiding a strong coupling problem of the dark photon interacting with other matter sectors during inflation [29]. As can be seen in the action, normalization of I is relevant for the – 22 – holds for a whole period of inflation if m2γ 0 I(tend )2 H 2  n0 (n0 + 1) (A.2) is satisfied at the end of inflation t = tend . Then, neglecting the mass term, we obtain a solution of Vk with the Bunch-Davies initial condition, r i −πkτ (1) Vk = √ H−n0 −1/2 (−kτ ) , (A.3) 2 2k given by the Hankel function of the first kind. We have chosen the arbitrary initial phase such that the mode function becomes real at the leading order outside the horizon −kτ  1. The mode functions of the dark electric and magnetic fields, defined in (3.8), are then given by   √ −i π 2 I d Vk (1) = Ek ≡ − 2 H (−kτ )5/2 H−n0 +1/2 (−kτ ) , (A.4) 3/2 a dτ I 2k √ k i π 2 (1) Bk ≡ 2 Vk = 3/2 H (−kτ )5/2 H−n0 −1/2 (−kτ ) . (A.5) a 2k Using the asymptotic form of the Hankel function,   Γ(ν) 2 ν (1) Hν (x) ' −i , (x → 0 , Re (ν) > 0) , π x the expressions (A.4) and (A.5) in the super-horizon limit (|kτ | → 0) are given by
  −Γ 21 − n0 H2 1 Ek ' n +1/2 √ 3/2 , n0 < , (−kτ )−n0 −2 2 2 0 πk
  Γ − 12 − n0 H2 1 Bk ' n +3/2 √ 3/2 , n0 < − . (−kτ )−n0 −3 2 2 0 πk (A.6) (A.7) (A.8) Noting −kτ = eNk −N  1, where Nk ≡ ln(k/H) is the number of e-foldings at which fluctuations with momentum mode k exits the horizon, their amplitudes are proportional to |Ek | ∝ e(−n0 −2)(N −Nk ) , |Bk | ∝ e(−n0 −3)(N −Nk ) , (N − Nk > 0) . (A.9) Therefore, on super-horizon scales the electric field is proportional to a−(n0 +2) and grows when n0 < −2 is satisfied. On the other hand, the magnetic field evolves as a−(n0 +3) and grows for n0 < −3. In the parameter domain n0 < 0 of our interest, it is hence clear that (i) sufficient production of the dark photon energy is achieved for n0 < −2, and (ii) the dark electric field dominates over the magnetic counterpart. When n is not constant but dynamically evolves in time, n 6= n0 , the solution (A.3) for a constant value of n0 does not necessarily hold. However, particle production still occurs when |n| becomes greater than a threshold value |n| = 2, provided the time variation of n magnitude of an effective coupling strength because it is inversely proportional to I. Since I is exponentially large at early times, the coupling strength is highly suppressed throughout the inflationary period. This can be observed already from the effective mass term, which is essentially an illustration of the coupling to matter, as mγ 0 /I is suppressed earlier during inflation for n0 < 0 and would be divergent if n0 > 0. – 23 – is sufficiently small. As we explore in the main text, this fact enables to validate a scenario of an evolving n to realize production of the dark photon on scales much smaller than the CMB scale. For completeness, we next consider the evolution of the longitudinal mode. It is characterized by the time variation of the function ∂τ2 zk /zk which is given by ∂τ2 zk 1 2p4 − (2n2 − 7n + 1 + τ dn/dτ )p2 M 2 + (n2 + n − τ dn/dτ )M 4 = 2 , zk τ (p2 + M 2 )2 (A.10) where we have defined p ≡ k/a and M ≡ mγ 0 /I. With a constant n ' n0 . and a negligible mass in Eq. (2.32), we obtain   2 (p  M ) , (A.11) ∂τ2 Xk + k 2 − 2 Xk ' 0 τ n0 (n0 + 1) ∂τ2 Xk − Xk ' 0 (p  M ) . (A.12) τ2 Then, for the Bunch-Davies initial condition, we find   i e−ikτ 1− (p  M ) , Xk ' √ kτ 2k C2 Xk ' C1 (−τ )n0 +1 + (p  M ) , (−τ )n0 (A.13) (A.14) where C1,2 are integration constants determined by connecting the two solutions at the conformal time τNR,k when p = M [27]: r (−τNR,k )n0 −1 k (−τNR,k )−n0 √ , C2 = C1 = . (A.15) 2 2n0 + 1 2k 3/2 In order to obtain C1 and C2 , we have assumed (i) −kτNR  1, and (ii) the continuity of the energy density associated with the longitudinal mode, not the continuity of Xk and its time derivative. Therefore, in either case of p  M and p  M , the contribution from the longitudinal mode is subdominant in comparison to the exponentially enhanced energy density of the transverse modes for n < −2, and we neglect the former in what follows. References [1] E. G. M. Ferreira, Astron. Astrophys. Rev. 29, 7 (2021), arXiv:2005.03254 [astro-ph.CO] . [2] W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astro-ph/0003365 . [3] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] . [4] R. Essig et al., in Community Summer Study 2013: Snowmass on the Mississippi (2013) arXiv:1311.0029 [hep-ph] . [5] D. Antypas et al., (2022), arXiv:2203.14915 [hep-ex] . [6] V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, Phys. Rev. Lett. 119, 031302 (2017), arXiv:1703.04683 [astro-ph.CO] . – 24 – [7] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] . [8] A. Arvanitaki and S. Dubovsky, Phys. Rev. D 83, 044026 (2011), arXiv:1004.3558 [hep-th] . [9] M. Baryakhtar, R. Lasenby, and M. Teo, Phys. Rev. D 96, 035019 (2017), arXiv:1704.05081 [hep-ph] . [10] V. Cardoso, P. Pani, and T.-T. Yu, Phys. Rev. D 95, 124056 (2017), arXiv:1704.06151 [gr-qc] . [11] V. Cardoso, O. J. C. Dias, G. S. Hartnett, M. Middleton, P. Pani, and J. E. Santos, JCAP 03, 043 (2018), arXiv:1801.01420 [gr-qc] . [12] S. Dubovsky and G. Hernández-Chifflet, JCAP 12, 054 (2015), arXiv:1509.00039 [hep-ph] . [13] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, JCAP 06, 013 (2012), arXiv:1201.5902 [hep-ph] . [14] S. Chaudhuri, P. W. Graham, K. Irwin, J. Mardon, S. Rajendran, and Y. Zhao, Phys. Rev. D 92, 075012 (2015), arXiv:1411.7382 [hep-ph] . [15] M. Silva-Feaver et al., IEEE Trans. Appl. Supercond. 27, 1400204 (2017), arXiv:1610.09344 [astro-ph.IM] . [16] Y. Michimura, T. Fujita, S. Morisaki, H. Nakatsuka, and I. Obata, Phys. Rev. D 102, 102001 (2020), arXiv:2008.02482 [hep-ph] . [17] S. Morisaki, T. Fujita, Y. Michimura, H. Nakatsuka, and I. Obata, Phys. Rev. D 103, L051702 (2021), arXiv:2011.03589 [hep-ph] . [18] H. Nakatsuka, S. Morisaki, T. Fujita, J. Kume, Y. Michimura, K. Nagano, and I. Obata, (2022), arXiv:2205.02960 [astro-ph.CO] . [19] P. W. Graham, J. Mardon, and S. Rajendran, Phys. Rev. D 93, 103520 (2016), arXiv:1504.02102 [hep-ph] . [20] A. E. Nelson and J. Scholtz, Phys. Rev. D 84, 103501 (2011), arXiv:1105.2812 [hep-ph] . [21] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, and F. Takahashi, Phys. Lett. B 801, 135136 (2020), arXiv:1810.07188 [hep-ph] . [22] J. A. Dror, K. Harigaya, and V. Narayan, Phys. Rev. D 99, 035036 (2019), arXiv:1810.07195 [hep-ph] . [23] R. T. Co, A. Pierce, Z. Zhang, and Y. Zhao, Phys. Rev. D 99, 075002 (2019), arXiv:1810.07196 [hep-ph] . [24] A. J. Long and L.-T. Wang, Phys. Rev. D 99, 063529 (2019), arXiv:1901.03312 [hep-ph] . [25] M. Bastero-Gil, J. Santiago, L. Ubaldi, and R. Vega-Morales, JCAP 04, 015 (2019), arXiv:1810.07208 [hep-ph] . [26] K. Nakayama, JCAP 10, 019 (2019), arXiv:1907.06243 [hep-ph] . [27] Y. Nakai, R. Namba, and Z. Wang, JHEP 12, 170 (2020), arXiv:2004.10743 [hep-ph] . [28] B. Ratra, Astrophys. J. 391, L1 (1992). [29] V. Demozzi, V. Mukhanov, and H. Rubinstein, JCAP 08, 025 (2009), arXiv:0907.1030 [astro-ph.CO] . [30] N. Barnaby, R. Namba, and M. Peloso, Phys. Rev. D 85, 123523 (2012), arXiv:1202.1469 [astro-ph.CO] . [31] T. Fujita and S. Mukohyama, JCAP 10, 034 (2012), arXiv:1205.5031 [astro-ph.CO] . [32] T. Fujita and S. Yokoyama, JCAP 09, 009 (2013), arXiv:1306.2992 [astro-ph.CO] . – 25 – [33] T. Fujita and S. Yokoyama, JCAP 03, 013 (2014), [Erratum: JCAP 05, E02 (2014)], arXiv:1402.0596 [astro-ph.CO] . [34] R. J. Z. Ferreira, R. K. Jain, and M. S. Sloth, JCAP 10, 004 (2013), arXiv:1305.7151 [astro-ph.CO] . [35] R. J. Z. Ferreira, R. K. Jain, and M. S. Sloth, JCAP 06, 053 (2014), arXiv:1403.5516 [astro-ph.CO] . [36] D. Green and T. Kobayashi, JCAP 03, 010 (2016), arXiv:1511.08793 [astro-ph.CO] . [37] H. Bazrafshan Moghaddam, E. McDonough, R. Namba, and R. H. Brandenberger, Class. Quant. Grav. 35, 105015 (2018), arXiv:1707.05820 [astro-ph.CO] . [38] T. Kobayashi, JCAP 05, 040 (2014), arXiv:1403.5168 [astro-ph.CO] . [39] T. Fujita and R. Namba, Phys. Rev. D 94, 043523 (2016), arXiv:1602.05673 [astro-ph.CO] . [40] S. Kanno, M. Kimura, J. Soda, and S. Yokoyama, JCAP 08, 034 (2008), arXiv:0806.2422 [hep-ph] . [41] M.-a. Watanabe, S. Kanno, and J. Soda, Phys. Rev. Lett. 102, 191302 (2009), arXiv:0902.2833 [hep-th] . [42] M.-a. Watanabe, S. Kanno, and J. Soda, Prog. Theor. Phys. 123, 1041 (2010), arXiv:1003.0056 [astro-ph.CO] . [43] A. E. Gumrukcuoglu, B. Himmetoglu, and M. Peloso, Phys. Rev. D 81, 063528 (2010), arXiv:1001.4088 [astro-ph.CO] . [44] R. Namba, Phys. Rev. D 86, 083518 (2012), arXiv:1207.5547 [astro-ph.CO] . [45] N. Bartolo, S. Matarrese, M. Peloso, and A. Ricciardone, Phys. Rev. D 87, 023504 (2013), arXiv:1210.3257 [astro-ph.CO] . [46] A. Naruko, E. Komatsu, and M. Yamaguchi, JCAP 04, 045 (2015), arXiv:1411.5489 [astro-ph.CO] . [47] K. Nakayama, JCAP 08, 033 (2020), arXiv:2004.10036 [hep-ph] . [48] N. Seto, S. Kawamura, and T. Nakamura, Phys. Rev. Lett. 87, 221103 (2001), arXiv:astro-ph/0108011 . [49] J. Crowder and N. J. Cornish, Phys. Rev. D 72, 083005 (2005), arXiv:gr-qc/0506015 . [50] A. Sesana et al., Exper. Astron. 51, 1333 (2021), arXiv:1908.11391 [astro-ph.IM] . [51] M. Kramer, D. C. Backer, J. M. Cordes, T. J. W. Lazio, B. W. Stappers, and S. Johnston, New Astron. Rev. 48, 993 (2004), arXiv:astro-ph/0409379 . [52] K. Bamba and J. Yokoyama, Phys. Rev. D 69, 043507 (2004), arXiv:astro-ph/0310824 . [53] R. Kallosh, A. Linde, and D. Roest, JHEP 11, 198 (2013), arXiv:1311.0472 [hep-th] . [54] S. Habib, K. Heitmann, G. Jungman, and C. Molina-Paris, Phys. Rev. Lett. 89, 281301 (2002), arXiv:astro-ph/0208443 . [55] S. Habib, A. Heinen, K. Heitmann, G. Jungman, and C. Molina-Paris, Phys. Rev. D 70, 083507 (2004), arXiv:astro-ph/0406134 . [56] L. Lorenz, J. Martin, and C. Ringeval, Phys. Rev. D 78, 083513 (2008), arXiv:0807.3037 [astro-ph] . [57] F. Olver, Asymptotics and Special Functions, Computer science and applied mathematics : a series of monographs and textbooks (Academic Press, 1974). [58] Y. Akrami et al. (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] . – 26 – [59] T. L. Smith, M. Kamionkowski, and A. Cooray, Phys. Rev. D 73, 023504 (2006), arXiv:astro-ph/0506422 . [60] T. L. Smith and R. Caldwell, Phys. Rev. D 100, 104055 (2019), arXiv:1908.00546 [astro-ph.CO] . [61] M. Braglia and S. Kuroyanagi, Phys. Rev. D 104, 123547 (2021), arXiv:2106.03786 [astro-ph.CO] . – 27 –