Analytic integral solutions for induced gravitational waves
Jinn-Ouk Gong



























Korea Astronomy and Space Science Institute, Daejeon 34055, Korea
Department of Science Education, Ewha Womans University, Seoul 03760, Korea
Asia Pacific Center for Theoretical Physics, Pohang 37673, Korea














We present analytic integral solutions for the second-order induced gravitational waves (GWs). After presenting all the possible second-order source terms, we calculate explicitly the solutions for the GWs induced by the linear scalar and tensor perturbations during matter- and radiation-dominated epochs. 




























































§ INTRODUCTION










A series of detection of the gravitational waves (GWs) by the LIGO and Virgo collaborations <cit.> has opened the era of multi-messenger astronomy led by GWs. The observed GW signals are originated from merging black holes and/or neutron stars, but there should be other events energetic enough to generate observable GWs. Such cosmological origins include cosmic strings <cit.>, phase transition <cit.> and preheating <cit.>. But the non-linear nature of gravity tells us that there is a persistent source of GWs – they can be induced by other cosmological perturbations at non-linear order <cit.>: while the scalar-, vector- and tensor-type perturbations are decoupled at linear order, they couple to each other and thus can be generated at non-linear order <cit.>. Especially, the second-order tensor perturbations, or GWs, induced by linear scalar perturbations[
It should be noted that large scalar perturbations may well be induced by large tensor perturbations on small scales, leading to copious production of primordial black holes <cit.>. Such induced scalar perturbations subsequently can further induce tensor perturbations. This is an interesting possibility, but is beyond the scope of the present work.
] may well be sizable if on small scales the primordial curvature perturbation is enhanced during inflation <cit.> or the density perturbation grows during an early matter-dominated (MD) epoch <cit.>.








In most of the literature, however, the study of second-order induced GWs has been focused on MD in one particular gauge choice for the scalar perturbations. The reason could be twofold. First, the (geometric) scalar perturbation remains constant even on sub-horizon scales during MD, but it decays quickly if it enters the horizon during radiation-dominated (RD) epoch [see (<ref>) and (<ref>)]. Thus, even if the primordial curvature perturbation is enhanced on small scales, one would naively expect such an enhancement to disappear during an RD epoch so that there should be no sizable induced GWs from scalar perturbations [see, however, <cit.>]. Furthermore, there is an upper bound on the contribution of the primordial GWs from the observations on the cosmic microwave background (CMB): in terms of the so-called tensor-to-scalar ratio, r < 0.07 at the pivot scale k=0.05 Mpc^-1 <cit.>. This means on the CMB scales the amplitude of the primordial GWs should be about 1/10 or even smaller than that of the primordial curvature perturbation. Thus, naturally scalar perturbations should be the most dominant source for the induced GWs compared to the other types of cosmological perturbations. For these reasons, it is very sensible to consider the induced GWs by scalar perturbations during MD.








Nevertheless, this does not mean at all that we have a complete and satisfactory understanding of the second-order induced GWs. First, there is no a priori reason why scalar perturbations should be considered only in a particular gauge condition. Moreover, as the linear scalar perturbations depend on the choice of gauge, the second-order GWs induced by their quadratic combinations should be also dependent on the gauge conditions. This is obvious and indeed was noticed early <cit.> but was explicitly shown only recently in <cit.>. Second, while scalar perturbations are likely to be the most important source for the induced GWs, other types of perturbations need not be neglected from the beginning. Especially, the linear tensor perturbations should be persistent. Furthermore, it is possible that the contributions of the tensor perturbations can be enhanced <cit.>, thus, the induced second-order GWs from them may well be significant accordingly[
More exactly, if the tensor-induced GWs are to be dominant, the scalar perturbations should not increase. Otherwise, the scalar-induced GWs are very likely to be more prominent than the tensor-induced ones. This condition allows certain models viable [e.g. <cit.>], but generally excludes models that include a break during inflation [e.g. <cit.>].
]. This indeed happens in certain concrete models beyond the standard slow-roll inflation <cit.>. Third, as the universe has evolved through both RD and MD, for a complete description of the induced GWs we need a proper understanding of RD as it is occurring. This was recognized early as well <cit.>, but analytic approach has been taken only recently <cit.>.






In this article, we provide analytic integral solutions for the second-order induced GWs from both linear scalar and tensor perturbations. We also present the full second-order source terms with all three types of cosmological perturbations, so it should be straightforward to calculate the solutions from the sources with vector perturbations. This article is outlined as follows. In Section <ref>, we provide the full traceless evolution equation for the spatial metric tensor, including all the explicit second-order source terms. In Section <ref>, we solve the equations of motion for the linear cosmological perturbations that will be used in Section <ref> to compute the analytic integral solutions for the second-order induced GWs. We briefly summarize our results in Section <ref>. Some technical details are relegated to the appendix sections.




























§ SECOND-ORDER EQUATION









Our metric convention of a flat Friedmann universe including cosmological perturbations is


    ds^2 
    = 
    -a^2 (1+2α) dη^2 - 2a^2 B_i dη dx^i
    + a^2 [ (1+2φ) δ_ij + 2γ_,ij + 2C_(i,j) + 2h_ij]
    dx^i dx^j
      ,


where dη = dt/a is the conformal time and a(η) is the scale factor. The indices of the perturbation variables are raised and lowered by δ_ij. Further, the shear χ_i is written as


    χ_i 
    ≡
    a ( B_i + aγ̇_,i + aĊ_i )
      .








The non-linear equation necessary for the induced GWs can be obtained from the traceless evolution equation for the spatial metric. Up to second order, the full equation is given by (<ref>) in Appendix <ref>. Further, writing B_i = β_,i + B^(v)_i with B^(v)i_,i = 0, the shear χ_i, anisotropic stress Π_ij and peculiar velocity v_i of the perfect fluid, whose energy density and pressure are written as ρ and p respectively, can be decomposed in terms of the scalar gradient, transverse vector and transverse and traceless tensor components as


    χ_i 
       
    =
    a ( β + aγ̇)_,i
    + a ( B^(v)_i + aĊ_i )
    ≡χ_,i + χ_i^(v)  ,
    
    Π_ij    
    = 
    1/a^2( Π_,ij - δ_ij/3ΔΠ)
    + 1/aΠ^(v)_(i,j) + Π^(t)_ij + δ_ij/3Π^k_k
      ,
    
    
    v_i
        
    =
    -v_,i + v^(v)_i
      ,


where obviously the superscripts (v) and (t) denote respectively the transverse vector and transverse and traceless tensor. Note that the last term of Π_ij is added because δ^ijΠ_ij≡Π^i_i ≠ 0 at non-linear order but is given by


    Π^i_i
    =
    2h^ijΠ_ij + v^iv^jΠ_ij + ⋯









With these decompositions, the second-order traceless evolution equation (<ref>) can be written as


    ḧ_ij + 3Hḣ_ij - Δ/a^2h_ij - 8π GΠ^(t)_ij
    + 1/a^2( ∂_i∂_j - δ_ij/3Δ) 
    [ 1/ad/dt( aχ) - α - φ - 8π GΠ]
    + 1/a[ 1/a^2d/dt( aχ^(v)_(i,j)) 
    - 8π G Π^(v)_(i,j)]
    =
    s_ij  ,


where s_ij denotes all second-order terms, which serve as the second-order “source”. As there are scalar, vector and tensor perturbations, at second order 6 combinations are possible:





1. Scalar-scalar source s_ij^(ss) denotes the collection of the products of two scalar perturbations:


    s_ij^(ss)
    =
       1/a^3d/dt[
    a ( 2φχ_,ij + φ_,iχ_,j + φ_,jχ_,i)
    ]
    + 1/a^2(
    κχ_,ij - 4φφ_,ij  -   3φ_,iφ_,j)
    + 1/a^4χ^,k_,iχ_,jk
        
    + 1/a^2[ 2αχ̇_,ij - Hαχ_,ij  + α̇χ_,ij 
    - 2(α+φ)α_,ij - α_,iα_,j - 2φ_(,iα_,j)]
    + 8π G (ρ+p) v_,iv_,j - 16π G φΠ_,ij/a^2
        - δ_ij/3{1/a^3d/dt[ 
    a ( 2φΔχ + 2φ_,kχ^,k)
    ]
    + 1/a^2(
    κΔχ - 4φΔφ   - 3φ^,kφ_,k)
    + 1/a^4χ^,klχ_,kl
            
    + 1/a^2[ 
    2αΔχ̇- HαΔχ + α̇Δχ 
    - 2(α+φ)Δα - α^,kα_,k - 2α^,kφ_,k]
    + 8 π G (ρ+p) v^,kv_,k - 16π GφΔ/a^2Π}  .







2. Scalar-tensor source s_ij^(st) denotes the collection of the products of one of the scalar perturbations and tensor perturbations:


    s_ij^(st) =
       d/dt[ 
     ḣ_ijα
    + 2 ( φḣ_ij + φ̇h_ij + 1/a^2h_i^kχ_,jk)
    + χ^,k/a^2(  h_ik,j  + h_jk,i  - h_ij,k)
    ]
    
        + 3H [ 
      ḣ_ijα
    + 2 ( φḣ_ij + φ̇h_ij + 1/a^2h_i^kχ_,jk)
    + χ^k/a^2(  h_ik,j  + h_jk,i  - h_ij,k)
    ]
    
        + αd/dt(  ḣ_ij)
    - 1/a^2χ^,kḣ_ij,k  
    + κḣ_ij  
    + 1/a^2[
    - 2   h_i^k  α_,jk
    - (  h_ik,j  + h_jk,i  - h_ij,k) α^,k]
    
        - δ_ij/3[ d/dt(
       1/a^2 2h^klχ_,kl)
    + 3H (
       1/a^2 2h^klχ_,kl)
    + 1/a^2(
    - 2   h^klα_,kl) ]
    + 1/a^2χ_,i^,kḣ_jk  
    - 1/a^2χ_,j^,kḣ_ik
        
    + 1/a^2[
     2 ( -2φΔh_ij + h_j^kφ_,ik
      - h_ijΔφ)
    + φ^,k( h_ik,j + h_jk,i - 3h_ij,k)
    - 2/3 
    h^klφ_,klδ_ij - 2h^klφ_,klδ_ij]
    
        
    - 16π G [ φΠ_ij^(t) + 1/a^2( h_i^kΠ_,ij - 1/3h_ijΔΠ 
    - δ_ij/3 h^klΠ_,kl) ]
      .










3. Tensor-tensor source s_ij^(tt) denotes the collection of the products of two tensor perturbations:


    s_ij^(tt) =
       d/dt( 
     2 h_i^k ḣ_jk)
    + 3H (
     2 h_i^k ḣ_jk)
    - δ_ij/3[ d/dt(
    2 h^klḣ_kl)
    + 3H (
     2 h^klḣ_kl)
    ]
    
        + 1/a^2{
    2 h^kl( h_il,jk  + h_jl,ik 
     - h_ij,kl  - h_kl,ij)
    - 2   h_i^k    Δh_jk 
    - h^kl_,i h_kl,j 
    + 2   h_i^k,l(  h_jl,k  - h_jk,l)
    
           
    - δ_ij/3[
    -4   h^klΔh_kl 
    +   h^kl,m(   2h_km,l   - 3h_kl,m)
    ]
    }
    - 16π G ( h_i^kΠ_jk^(t) - δ_ij/3 h^klΠ_kl^(t))
      .









4. Scalar-vector source s_ij^(sv) denotes the collection of the products of one of the scalar perturbations and vector perturbations:


    s_ij^(sv) =
       1/a^3d/dt{ 
      a 
      [ χ^(v)_(i,j)α
    + 2χ^(v)_(i,j)φ
    +   φ_,jχ^(v)_i  + φ_,iχ^(v)_j   ]
    }
        + αd/dt(  1/a^2χ^(v)_(i,j))
    - 1/a^4χ^,kχ^(v)_(i,j)k 
    +   1/a^2χ^(v)_(i,j)κ
     + 1/a^4(
    2χ^(v)_k,(iχ_,j)^,k  + χ^,kχ^(v)_k,ij)
    
        - δ_ij/3[ 1/a^3d/dt(
    2aχ^(v)_kφ^,k)
     + 1/a^4(
    2χ^(v)_k,lχ^,kl + χ^,kΔχ^(v)_k
     ) ]
    + 1/a^4( χ^(v)_i,kχ_,j^,k + χ_,i^,kχ^(v)_(j,k)) 
    - 1/a^4( χ^(v)_k,jχ_,i^,k + χ_,j^,kχ^(v)_(i,k)) 
    
        
    - 16π G φ/aΠ_(i,j)^(v)
    - 16 π G (ρ+p) ( v_,(iv^(v)_j)  - δ_ij/3 v^,kv^(v)_k  )
      .








5. Vector-vector source s_ij^(vv) denotes the collection of the products of two vector perturbations:


    s_ij^(vv) =
        
    - 1/a^4χ^(v)kχ^(v)_(i,j)k 
    + 1/a^4(   χ^(v)kχ^(v)_k,ij + χ^(v)k_,iχ^(v)_k,j)
    - δ_ij/3[
     1/a^4( χ^(v)k,lχ^(v)_k,l + χ^(v)kΔχ^(v)_k )
    ]
    
        
    + 1/a^4χ^(v)_i^,kχ^(v)_(j,k) 
    - 1/a^4χ^(v)k_,jχ^(v)_(i,k) 
    + 8 π G (ρ+p) [ v^(v)_iv^(v)_j - δ_ij/3 v^(v)kv^(v)_k ]
      .








6. Vector-tensor source s_ij^(vt) denotes the collection of the products of one of the vector perturbations and tensor perturbations:


    s_ij^(vt) =
       d/dt{1/a^2[
     2h_i^kχ^(v)_(j,k) 
    +  χ^(v)k(  h_ik,j  + h_jk,i  - h_ij,k)
    ]
    }
    + 3H 1/a^2[
     2h_i^kχ^(v)_(j,k) 
    +  χ^(v)k(  h_ik,j  + h_jk,i  - h_ij,k)
    ]
    
        
    - 1/a^2χ^(v)kḣ_ij,k  
     - 1/a^2χ^(v)k_,kḣ_ij 
    - δ_ij/3{d/dt[
     2/a^2   h^klχ^(v)_(k,l)]
    + 3H [
     2/a^2   h^klχ^(v)_(k,l)]
     }
    + 1/a^2χ^(v)_i^,kḣ_jk  
    - 1/a^2χ^(v)k_,jḣ_ik
        
    - 16π G 1/a( h_i^kΠ_(j,k)^(v) - δ_ij/3 h^klΠ_(k,l)^(v))
      .







Having sorted out all possible second-order source terms, we can proceed to find the solution of the second-order induced GWs as follows. First, we solve the linear equations and obtain their solutions. Then these linear solutions can be used to obtain the explicit form of the sources. After the transverse-traceless projection of the source s_ij [see (<ref>)], we can solve the inhomogeneous equation for the tensor perturbations and obtain the analytic integral solutions.





















§ LINEAR SOLUTIONS












 §.§ Vector perturbations at linear order










The following linear equations for the vector-type perturbations are derived respectively from the momentum constraint, traceless evolution equation (<ref>) and momentum conservation equation <cit.>:


    Δ/2a^3χ_i^(v) + 8π G(ρ+p)v_i^(v)    =
    0
      ,
    
    1/a^2d/dt( aχ_i^(v)) - 8π GΠ_i^(v)    =
    0
      ,
    
    1/a^4(ρ+p)d/dt[ a^4(ρ+p)v_i^(v)] + Δ/2a^2Π_i^(v)/ρ+p    = 
    0
      .


A great simplification is made in the case of vanishing vector-type stress, Π_i^(v) = 0. Then all the linear vector perturbations always vanish: 


    χ_i^(v) = v_i^(v) = 0   .


So among the possible sources to the second-order GWs, scalar-vector, vector-vector and vector-tensor contributions are absent, and we have only scalar-scalar, scalar-tensor and tensor-tensor sources.











 §.§ Scalar perturbations at linear order









With the perturbation in the extrinsic curvature κ being written as


    κ = 3Hα - 3φ̇- Δ/a^2χ  ,


the complete set of the linear equations for scalar perturbations is


    4π Gδρ + Hκ + Δ/a^2φ    = 
    0 
      ,
    
    κ + Δ/a^2χ - 12π G (ρ+p) av 
        = 
    0 
      ,
    
    κ̇+ 2Hκ - 4π G (δρ + 3δp) + ( 3Ḣ + Δ/a^2)α    = 
    0 
      ,
    
    χ̇+ Hχ - φ - α - 8π GΠ    = 
    0 
      ,
    
    δ̇ρ̇ + 3H(δρ + δp) - (ρ+p) ( κ - 3Hα + Δ/av ) 
        = 
    0 
      ,
    
    1/a^4(ρ+p)d/dt[ a^4(ρ+p)v ] - 1/aα 
    - 1/a(ρ+p)( δp + 2/3Δ/a^2Π) 
        = 
    0
      .


While δp = c_s^2δρ + τδS with δS being the entropy perturbation, for a barotropic fluid in the absence of δS, simply c_s^2 = w. Then, assuming no anisotropic stress (Π = 0) for simplicity, these equations become even simpler and allow analytic solutions.











  §.§.§ Solutions for linear scalar perturbations during MD







The solutions of the scalar perturbations during MD are already given in <cit.> up to linear order, in <cit.> up to second order and in <cit.> up to third order respectively. The solutions can be written conveniently in terms of the curvature perturbation φ, which does not decay but remains constant during MD even on sub-horizon scales. Thanks to gauge transformations, the solutions in one gauge are enough to find those in other gauges. 





We can readily solve the linear equations for the scalar perturbations during MD (p=0) without anisotropic stress (Π = 0). Summarizing, in the comoving gauge for which v = γ = 0, we find


    φ_v =   , 
      α_v = 0   , 
      χ_v = 2/5H  , 
      κ_v = 2/5k^2/a^2H  , 
      δ_v = 2/5k^2/a^2H^2  .


Here, the subscript v means the solutions are written in the comoving gauge.
These solutions can be by gauge transformation used to obtain solutions in different gauge conditions, e.g. zero-shear gauge for which β = γ = 0 thus, as the name stands, χ = 0 (thus a subscript χ) and we find


    φ_χ = 3/5  , 
      α_χ = -3/5  , 
      κ_χ = -9/5H  , 
      
    v_χ = -2/5aH  , 
      δ_χ = 6/5( 1+k^2/3a^2H^2)   .


Note that we can read easily the well-known relation during MD between the initial amplitude of the curvature perturbation in the comoving gauge  and that in the zero-shear gauge, or the “gravitational potential” Φ = -α_χ = φ_χ, as


    Φ = 3/5  .
















  §.§.§ Solutions for linear scalar perturbations during RD







During RD, the linear equations of motion are most readily solvable in the zero-shear gauge. With w=1/3 and a ∝η during RD, the equation of motion for the curvature perturbation φ_χ is obtained from the trace evolution equation combined with the energy constraint and traceless evolution equation as [see e.g. <cit.>]


    d^2φ_χ/dη^2 + 4/ηdφ_χ/dη
    - Δ/3φ_χ = 0
      .


Then, with z ≡ kη/√(3) where 1/√(3) is the sound speed during RD, we can straightforwardly find the linear solutions for the scalar perturbations:


    φ_χ    =
    2j_1(z)/z  ,
    
    
    v_χ    =
    - 1/aH[ j_0(z) - 2 j_1(z)/z]   ,
    
    κ_χ    =
    -6H [ j_0(z) - 2 j_1(z)/z]   ,
    
    α_χ    =
    -φ_χ  ,
    
    δ_χ    =
    4[ - j_0(z) + 2j_1(z)/z + zj_1(z) ]
      ,


where j_n is the first-kind spherical Bessel function of order n. The solutions in other gauges can be obtained by appropriate gauge transformations, e.g. the curvature perturbation in the comoving gauge as φ_v = φ - aHv. In the comoving gauge, the linear solutions for the scalar perturbations are


    φ_v
        =
     j_0(z)
      ,
    
    χ_v
        =
    1/H[ j_0(z) - 2 j_1(z)/z]   ,
    
    κ_v
        =
    3H z^2 [ j_0(z) - 2 j_1(z)/z]   ,
    
    α_v
        =
    [ 2j_0(z) - 4 j_1(z)/z - zj_1(z) ]   ,
    
    δ_v
        =
    4[ - 2j_0(z) + 4j_1(z)/z + zj_1(z) ]
      .


Here, we have set the coefficients in such a way that the initial amplitude of the curvature perturbation in the comoving gauge is, as for the solution during MD, , i.e. lim_z→0φ_v(z) =.

Note that since lim_z→0 j_0(z) = 1 and lim_z→0 j_1(z)/z = 1/3, we can find the well-known relation during RD between the comoving curvature perturbation  and the gravitational potential Φ as Φ = 2/3.
















 §.§ Tensor perturbations at linear order









Decomposing the tensor perturbations in terms of the two polarization tensors e_ij^+ and e_ij^× in the Fourier space, 


    h_ij(t,x) = ∫d^3k/(2π)^3 e^ik·x h_ij(t,k)
    = ∫d^3k/(2π)^3 e^ik·x[ h_+(t,k)e_ij^+(k) + h_×(t,k)e_ij^×(k) ] .


Since the polarization tensors are orthogonal to each other, i.e.


    e^ij_+ e_ij^× = 0 
      and  
    e^ij_+e_ij^+ = e^ij_× e_ij^× = 1
      ,


we can invert this to find


    h_λ(t,k) = e_λ^ij(k) ∫ d^3x e^-ik·x h_ij(t,x)


for each polarization λ. Then the linear equation of motion for each polarization mode is identical as 


    ḧ+3Hḣ + k^2/a^2h = 0 
      ,


where we have omitted the polarization index λ. Introducing v ≡ ah and moving to the conformal time dη = dt/a, the equation becomes


    d^2v/dη^2 + ( k^2 - 1/ad^2a/dη^2)v = 0 
      .













  §.§.§ Solutions for linear tensor perturbations during MD








During MD, a ∝η^2 so that


    1/ad^2a/dη^2 = 2/η^2  .


Thus in terms of a new variable x ≡ kη, (<ref>) becomes


    d^2v/dx^2 + ( 1 - 2/x^2) v = 0   .


The general solution of this equation is


    v = c_1 xj_1(x) + c_2 xy_1(x)
      ,


where y_n is the second-kind spherical Bessel function of order n. Since lim_x→0 y_1(x)/x → -∞ and lim_x→0 j_1(x)/x = 1/3, we choose j_1(x)/x as the proper solution, with the value at x→0 being the primordial value for the tensor perturbation h^λ_0( k) for each polarization λ: 


    h_λ(η, k)
    = 3h_0^λ( k) j_1(kη)/kη  .












  §.§.§ Solutions for linear tensor perturbations during RD







The basic equation of motion for the tensor perturbations is essentially the same as that during MD. That is, with v ≡ ah, we have the same linear equation for v given by (<ref>). The only difference is that since a ∝η during RD, with x ≡ kη, (<ref>) is simply


    d^2v/dx^2 + v = 0   ,


and the general solution is 


    v = ah = c_1 xj_0(x) + c_2 xy_0(x)   .


Since lim_x→0 y_0(x) → -∞ and lim_x→0 j_0(x) = 1, we choose j_0(x) as the proper solution, with the value at x→0 being the primordial value for the tensor perturbation h^λ_0( k) for each polarization λ:


    h_λ(η, k) = h_0^λ( k) j_0(kη)   .






















§ SECOND-ORDER SOLUTIONS FOR INDUCED GWS













 §.§ Equation of motion for tensor perturbations with source






To extract only the tensor parts from (<ref>), we apply the transverse-traceless projection so that s_ij on the right-hand side becomes what only sources tensor perturbations, s_ij^(tensor) <cit.>:


    s_ij^(tensor)   ≡ 
    s_ij - 3/2( ∂_i∂_j - δ_ij/3Δ) Δ^-2∂_k∂_l s^kl
    - 2 Δ^-1∂_(i∂^ks_j)k + 2Δ^-2∂_i∂_j ∂_k∂_l s^kl
        = 
    s_ij - 2 Δ^-1∂_(i∂^k s_j)k 
    + 1/2Δ^-2( ∂_i∂_j + δ_ijΔ) ∂_k∂_l s^kl  .


Since the two traceless polarization tensors are orthogonal to each other, we can extract the individual equation of each polarization mode h_λ by multiplying the corresponding polarization tensor e_ij^λ. Moreover, since e^ij_λ k_i = e^ij_λ k_j = 0, we have in the Fourier space


    e^ij_λ(k) ∫ d^3x e^-ik·x s_ij^(tensor)(x)
        = 
    e^ij_λ(k) ∫ d^3x e^-ik·x∫d^3q/(2π)^3 e^iq·x s_ij^(tensor)(q)
    
        = 
    e^ij_λ(k) ∫ d^3q 
    [ s_ij - q_iq_k/q^2 s^k_j  + q_jq_k/q^2 s^k_i
    + 1/2( -q_iq_jq_kq_l/q^4 + q_kq_l/q^2δ_ij) s^kl]
    ∫d^3x/(2π)^3 e^-i(k-q)·x
        = 
    e^ij_λ(k)   s_ij(k) 
      .


Thus, for each polarization λ the equation of motion is


    ḧ_λ(t,k) + 3Hḣ_λ(t,k) + k^2/a^2h_λ(t,k)
    = e^ij_λ(k)   s_ij(k)   .







One more simplification is ahead. Since the source term s_ij in the above equation is multiplied by the traceless polarization tensor e^ij_λ, the terms proportional to δ_ij in s_ij identically vanish on the right-hand side of (<ref>). Thus, in the absence of the anisotropic stress, the source terms that survive (<ref>) are


    s_ij^(ss) 
    =    1/a^3d/dt[
    a ( 2φχ_,ij + φ_,iχ_,j + φ_,jχ_,i)
    ]
    + 1/a^2(
    κχ_,ij - 4φφ_,ij  -   3φ_,iφ_,j)
    + 1/a^4χ^,k_,iχ_,jk
        
    + 1/a^2[ 2αχ̇_,ij - Hαχ_,ij  + α̇χ_,ij 
    - 2(α+φ)α_,ij - α_,iα_,j - 2φ_(,iα_,j)]
    + 8π G (ρ+p) v_,iv_,j  ,
    
    
    s_ij^(tt) 
    =    d/dt( 
     2 h_i^k ḣ_jk)
    + 3H (
     2 h_i^k ḣ_jk)
    
        + 1/a^2[
    2  h^kl(  h_il,jk  + h_jl,ik 
     - h_ij,kl  - h_kl,ij)
    - 2   h_i^k    Δh_jk 
    - h^kl_,i h_kl,j 
    + 2   h_i^k,l(  h_jl,k  - h_jk,l) ]
      ,
    
    
    s_ij^(st) 
    =    d/dt[ 
     ḣ_ijα
    + 2 ( φḣ_ij + φ̇h_ij + 1/a^2h_i^kχ_,jk)
    + χ^,k/a^2(  h_ik,j  + h_jk,i  - h_ij,k)
    ]
    
        + 3H [ 
      ḣ_ijα
    + 2 ( φḣ_ij + φ̇h_ij + 1/a^2h_i^kχ_,jk)
    + χ^,k/a^2(  h_ik,j  + h_jk,i  - h_ij,k)
    ]
    
        + αd/dt(  ḣ_ij)
    - 1/a^2χ^,kḣ_ij,k  
    + κḣ_ij  
    + 1/a^2[
    - 2   h_i^k  α_,jk
    - (  h_ik,j  + h_jk,i  - h_ij,k) α^,k]
    + 1/a^2χ_,i^,kḣ_jk  
    - 1/a^2χ_,j^,kḣ_ik
        
    + 1/a^2[
     2 ( -2φΔh_ij + h_j^kφ_,ik
      - h_ijΔφ)
    + φ^,k( h_ik,j + h_jk,i - 3h_ij,k)
    ]
      .










To find the solution for this inhomogeneous equation, we use the Green's function solution [for earlier attempts to apply the Green's function solutions to cosmological perturbations, see <cit.>]. That is, let L be a linear second-order differential operator, and the equation we want to solve is of the form


    Ly(x) = r(x)   ,


with the two homogeneous solutions being y_1 and y_2. Then, the full solution is given by


    y(x) 
        =  
    (appropriate combination of y_1 and y_2 according to the boundary conditions)
    
         
    + ∫ dx̃ r(x̃) y_1(x̃)y_2(x) - y_2(x̃)y_1(x)/y_1(x̃)y_2'(x̃) - y_2(x̃)y_1'(x̃)_≡ G(x,x̃)  .


We have seen that the homogeneous solutions during MD are given by (<ref>), then the Green's function during MD becomes


    G_MD(η,η̃) 
    = 
    xx̃/k[ j_1(x̃)y_1(x) - j_1(x)y_1(x̃) ] 
      .


With the homogeneous solutions during RD being given by (<ref>), the Green's function during RD becomes


    G_RD(η,η̃)
    =
    xx̃/k[ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ]
      .


Comparing with (<ref>), the only difference is the order of the spherical Bessel functions inside the square brackets.






Thus, we expect the second-order induced GWs sourced by the product of two linear perturbations X and Y would be, during both MD and RD, of the general integral form:


    h_λ(η, k)
        =
    1/a∫_0^η dη̃ [ a^3(η̃)  e^ij_λ(k) s_ij(k) ]
    G(η,η̃)
    
        =
    ∫d^3q/(2π)^3[ e^ij_λ(k) (⋯)_ij]
    X_0( k- q) Y_0( q) 
    ∫_0^x dx̃ K(x̃, k, q)
      .


Here, the integral over an internal momentum q is because the source s_ij is the product of two perturbations X and Y, it is written as a convolution in the Fourier space. The terms inside the square brackets constitute the (dimensionless) projection of the polarization tensor e^ij_λ(k). X_0( k) and Y_0( k) denote respectively the initial amplitudes of X and Y, i.e. ( k) and/or h_0^λ( k). Finally, the integral over x̃ is the kernel which is a function of both momenta as well as time. Our main concern in finding this analytic integral solution is to compute this kernel. In the following, we proceed with the sources (<ref>), (<ref>) and (<ref>) to calculate the closed analytic form of the kernel. 


























































































































 §.§ Scalar-scalar induced GWs during MD








We first consider the scalar-scalar source (<ref>). The analytic integral solutions in various gauge conditions are given only very recently in <cit.> so we can check our results in this section. We consider only two gauges, comoving and zero-shear gauges in which the solutions of the linear scalar perturbations during MD are given respectively by (<ref>) and (<ref>). We first compute the Fourier component of s_ij^(ss) in the zero-shear gauge for which χ = 0 and α = -φ, so s_ij^(ss) is greatly simplified. After straightforward calculations we find in the zero-shear gauge the source s_ij^(ss) purely in terms of the initial perturbation  as


    e_λ^ij(k) s_ij^(ss)(k)
    = 
    1/a^26/5∫d^3q/(2π)^3[ e_λ^ij(k) q_iq_j ] 
    (k-q)(q)
      .


Note that other than the overall 1/a^2, there is no time dependence. 

Then for each polarization, the solution of the GWs induced by the scalar-scalar source in the zero-shear gauge is


    h_λ(η,k) 
        = 
    6/5[ 1 - 3j_1(kη)/kη] 
    ∫d^3q/(2π)^3[ e_λ^ij(k) q_iq_j/k^2] 
    (k-q)(q)
      .


Compared with the general form of the solution (<ref>), the kernel is a function of k and η in the specific combination kη, and thus can be pulled out of the internal momentum integral. 








Likewise, in the comoving gauge, we have v=0 and α=0 during MD, so that the scalar-scalar source (<ref>) becomes


    e^ij_λ(k) s_ij^(ss) 
    = 
    1/a^2∫d^3q/(2π)^3[ e^ij_λ(k) q_iq_j ]
    ( 1 - 2/25k^2/a^2H^2)
    (k-q)(q)
      ,


where the additional term comes from the spatial gradient of the shear. Then we can find trivially the solution as, using aH=2/η during MD,


    h_λ(η,k) 
        = 
    6/5[ 1 - 3j_1(kη)/kη  - (kη)^2/60]
    ∫d^3q/(2π)^3[ e_λ^ij(k) q_iq_j/k^2] 
    (k-q)(q)
      .


Again, the kernel is a function of only kη. The reason why we have such a simple kernel for the scalar-scalar source is because the scalar perturbations can be written in terms of the constant , so the time integral is greatly simplified. These scalar-scalar induced solutions in the zero-shear gauge (<ref>) and comoving gauge (<ref>) agree with <cit.>. In Figure <ref>, we show the kernels barring the factor 6/5. As kη becomes bigger, the zero-shear gauge kernel approaches 1 while that in the comoving gauge increases as (kη)^2. Thus on small scales the amplitude of the induced GWs in the comoving gauge is much bigger than that in the zero-shear gauge. This shows clearly the gauge dependence of the scalar-induced GWs.














































 §.§ Scalar-tensor induced GWs during MD








Next we consider the GWs induced by the scalar-tensor source during MD.

Considering first the zero-shear gauge in which χ = 0 and α = - φ, we find (omitting the subscript χ for the scalar perturbation φ)


    s_ij^(st)( k) 
        =
    ∫d^3q_1d^3q_2/(2π)^3δ^(3)( k- q_12)
    1/2[
    ( d/dt + 3H ) ( φ_1 ḣ_λ'2 
    + 2φ̇_1 h_λ'2)
    - φ_1 d/dt( ḣ_λ'2)
    + κ_1 ḣ_λ'2 
    + 2/a^2( 2 k^k q_2k + q_1^2 ) φ_1 h_λ'2] e_ij2^λ'
         
    + ( q_1 ↔ q_2) 
      ,


where the subscript 1 for a perturbation variable means it is a function of q_1, e.g. φ_1 = φ( q_1) and so on.

Using φ_χ = 3/5 and κ_χ = -9H/5 during MD with  = constant, the terms including the time derivative of the linear tensor perturbation vanish and we have the following simple expression for the source:


    e^ij_λ( k)
    s_ij^(st)( k) 
    =
    1/a^26/5∫d^3q/(2π)^3[ e^ij_λ( k) e_ij^λ'( q) ]
    ( k^2 + q^2 ) ( k- q) h_λ'(η, q) 
      .


It is very important to note that the linear order perturbation h_λ' does possess time dependence as given by (<ref>). Thus, unlike the scalar-scalar source case, the time integral for the Green's function solution includes another time-dependent function from h_λ', namely, j_1(kη)/(kη):


    h_λ(η, k)
        =
    18/5∫d^3q/(2π)^3[ e^ij_λ( k) e_ij^λ'( q) ]
    ( q^2 + k^2 ) ( k- q) h_0^λ'( q) 
    1/kqx[
    y_1(x) ∫_0^x dx̃x̃^2 j_1( q/kx̃) j_1(x̃)
    - j_1(x) ∫_0^x dx̃x̃^2 j_1( q/kx̃) y_1(x̃)
    ]
    
        = 
    18/5∫d^3q/(2π)^3[ e^ij_λ( k) e_ij^λ'( q) ]
    ( k- q) h_0^λ'( q) 
    [  j_1(kη)/kη - j_1(qη)/qη]
    q^2+k^2/q^2-k^2  .


The detail of the x̃-integrals is given in Appendix <ref>.








In the comoving gauge where v=0, using = constant and 


    d/dt( 1/a^2H)
    =
    -1/2a^2  ,


which follows from H = 2/(3t) during MD, straightforward calculations give


    e^ij_λ( k) s_ij^(st)( k)
    =
    2/a^2∫d^3q/(2π)^3[ e^ij_λ( k) e^λ'_ij( q) ]
    ( k- q)
    [
    1/5H( k^2 - q^2 ) ḣ_λ'( q) 
    + k^2 h_λ'( q) 
    ]
      .


Unlike the zero-shear gauge case, the time derivative of the linear tensor perturbation remains, which is given by


    ḣ_λ'( q)
    =
    q/ad/d(qη) h_λ'(η, q)
    =
    q/a 3 h_0^λ'( q) 
    [ -3j_1(qη)/(qη)^2 + j_0(qη)/qη]
      .


We can find the solution in a manner very similar to (<ref>) but with different s_ij^(st) as:


    h_λ(η, k)
        =
    6/k^2x∫d^3q/(2π)^3[ e^ij_λ( k) e_ij^λ'( q) ]
    ( k- q) h_0^λ'( q)
    ∫_0^x dx̃x̃^3
    [ 3q^2+7k^2/10j_1(qx̃/k)/qx̃/k
    + k^2-q^2/10 j_0( q/kx̃) ]
    
         ×[ j_1(x̃)y_1(x) - j_1(x)y_1(x̃) ]
    
        =
    6
    ∫d^3q/(2π)^3[ e^ij_λ( k) e_ij^λ'( q) ]
    ( k- q) h_0^λ'( q)
    [
    q^2+5k^2/5(q^2-k^2)j_1(kη)/kη 
    - 5q^2+7k^2/10(q^2-k^2)j_1(qη)/qη + 1/10 j_0(qη) 
    ]
      .


This is the solution of the induced h_λ( k,η) from the scalar-tensor source in the comoving gauge. In Figure <ref> we show the kernels in both gauges. For kη≳ 1 with sizeable q/k, the comoving gauge kernel exhibits more rapid oscillations. But for small q/k, both kernels are approximated by 1/3-j_1(kη)/(kη) + (q^2/k^2) so they behave similarly. Compared to the scalar-scalar induced GWs as we have seen previously, the difference in the zero-shear gauge and comoving gauge is not prominent. This is because the gradient of the shear, that gives rise to the huge difference on small scales for the scalar-scalar induced GWs, is highly suppressed by the exponentially decaying linear GWs on small scales.














































 §.§ Tensor-tensor induced GWs during MD










Next we consider the tensor-tensor source during MD. Multiplying the polarization tensor e^ij_λ( k) and separating the time-dependent part gives


    e^ij_λ( k) s_ij^(tt)( k)
        =
    1/a^2∫d^3q_1d^3q_2/(2π)^3δ^(3)( k- q_12)
    { 
    9 h_0^λ_1( q_1) h_0^λ_2( q_2)
    e^ij_λ( k) e_i1^kλ_1 e_jk2^λ_2_≡_λ( q_1, q_2)  q_1q_2 
    [ -3j_1(q_1η)/(q_1η)^2 + j_0(q_1η)/q_1η]
    [ -3j_1(q_2η)/(q_2η)^2 + j_0(q_2η)/q_2η]
    
       
    + 
    9 h_0^λ_1( q_1) h_0^λ_2( q_2) e^ij_λ( k)
    [
    e^kl_λ_11(
    -q_2jq_2k e^λ_2_il2 - q_2iq_2k e_jl2^λ_2 
    + q_2kq_2l e^λ_2_ij2)
    
       
    + 1/2 q_1iq_1j e_λ_1^kl1 e_kl2^λ_2 
    + ( q_1· q_2) e_i1^kλ_1 e_jk2^λ_2
    - q_1^lq_2k e_i1^kλ_1 e_jl2^λ_2] 
    _≡_λ( q_1, q_2) j_1(q_1η)/q_1ηj_1(q_2η)/q_2η}
         
    + ( q_1 ↔ q_2) 
      .


As can be seen, only the Bessel function terms contain time dependence. Then the solution for h_λ can be written as


    h_λ(η, k) 
        =
    ∫d^3q_1d^3q_2/(2π)^3δ^(3)( k- q_12)
    {[ _λ( q_1, q_2) + ( q_1 ↔ q_2) ]
    q̅_1q̅_2/x[
    9/q̅_1^2q̅_2^2∫_0^x dx̃x̃^-1 j_1(q̅_1x̃) j_1(q̅_1x̃)
    - 3/q̅_1^2q̅_2∫_0^x dx̃ j_1(q̅_1x̃) j_0(q̅_2x̃)
    
       
    - 3/q̅_1q̅_2^2∫_0^x dx̃ j_0(q̅_1x̃) j_1(q̅_2x̃)
    + 1/q̅_1q̅_2∫_0^x dx̃x̃ j_0(q̅_1x̃) j_0(q̅_2x̃)
    ]
    [ j_1(x̃)y_1(x) - j_1(x)y_1(x̃) ] 
    
       
    +
    [ _λ( q_1, q_2) + ( q_1 ↔ q_2) ]
    1/k^2x1/q̅_1q̅_2∫_0^x dx̃x̃ j_1(q̅_1x̃) j_1(q̅_1x̃)
    [ j_1(x̃)y_1(x) - j_1(x)y_1(x̃) ] 
    }  ,


where q̅_1 ≡ q_1/k and q̅_2 ≡ q_2/k respectively.








Now, using the recurrence relation (<ref>) for n=1, we find


    j_2(q̅_2x) 
    = 
    3/q̅_2x j_1(q̅_1x) - j_0(q̅_1x) 
    = 
    3/q̅_2[ j_1(q̅_2x)/x - q̅_2/3j_0(q̅_2x) ]
      ,


so that the terms multiplied by _λ( q_1, q_2) become very simple as


    [
    9/q̅_1^2q̅_2^2∫_0^x dx̃x̃^-1 j_1(q̅_1x̃) j_1(q̅_2x̃)
    - 3/q̅_1^2q̅_2∫_0^x dx̃ j_1(q̅_1x̃) j_0(q̅_2x̃)
    - 3/q̅_1q̅_2^2∫_0^x dx̃ j_0(q̅_1x̃) j_1(q̅_2x̃)
    + 1/q̅_1q̅_2∫_0^x dx̃x̃ j_0(q̅_1x̃) j_0(q̅_2x̃)
    ]
    f_1(x̃)
    
        =
    1/q̅_1q̅_2∫_0^x dx̃x̃ j_2(q̅_1x̃) j_2(q̅_2x̃) f_1(x̃)
      ,


where f_1(x̃) denotes both first- and second-kind of the spherical Bessel functions. Since lim_x→0y_1(x) = -1/x^2-1/2+⋯ and lim_x→0j_n(x) ∼ x^n, we always have converging results. Thus, we can write (<ref>) as


    h_λ(η, k) 
        =
    ∫d^3q_1d^3q_2/(2π)^3δ^(3)( k- q_12)
    {[ _λ( q_1, q_2) + ( q_1 ↔ q_2) ]
    1/x∫_0^x dx̃x̃ j_2(q̅_1x̃) j_2(q̅_2x̃)
    [ j_1(x̃)y_1(x) - j_1(x)y_1(x̃) ] 
    
       
    +
    1/k^2[ _λ( q_1, q_2) + ( q_1 ↔ q_2) ]
    1/q̅_1q̅_2x∫_0^x dx̃x̃ j_1(q̅_1x̃) j_1(q̅_2x̃)
    [ j_1(x̃)y_1(x) - j_1(x)y_1(x̃) ] 
    }  .


Both the x̃ integrals can be performed analytically, with the details given in Appendix <ref>. Performing the momentum integral using the delta function, finally (<ref>) becomes


    h_λ(η, k) 
        =
    ∫d^3q/(2π)^3 18 h_0^λ_1( q) h_0^λ_2( k- q)
    {
    e_λ^ij( k) e_i^kλ_1( q) e_jk^λ_2( k- q) 
    F_MD( k, q,η)
    
       
    +
    1/k^2(
    e_λ^ij( k)
    e^kl_λ_1( q) (
    q_jk_k e_il^λ_2( k- q) 
    + q_ik_k e_jl^λ_2( k- q)
    + k_kk_l e_ij^λ_2( k- q)
    )
    
       
    +
    e_λ^ij( k) 
    [
    1/2 q_i q_j e^kl_λ_1( q) e_kl^λ_2( k- q)
    +  q·( k- q) e_i^kλ_1( q) e_jk^λ_2( k- q)
    - q^lk_k e_i^kλ_1( q) e_jl^λ_2( k- q)
    ]
    ) 
    G_MD( k, q,η)
    }  ,


where the kernels F_MD and G_MD are given respectively by (<ref>) and (<ref>). In Figure <ref> we show F_MD and G_MD.









































































 §.§ Scalar-scalar induced GWs during RD








Until now we have considered the second-order solutions for the induced GWs during MD. Now we consider the solutions during RD. Especially, after the first version of the present work appeared in arXiv, there have been a number of papers published in the literature, e.g. <cit.>, claiming that the gauge dependence of the induced GWs during RD disappears. This is not the case, at least regarding the solutions in the two gauge conditions we examine here explicitly – the zero-shear gauge and comoving gauge: the solutions are clearly different in these two gauges as can be seen in (<ref>) and (<ref>) [see also <cit.>]. 







An important difference for the linear solutions for the scalar perturbations is that now the curvature perturbation φ does not stay constant, but decays once the mode enters the horizon. Thus, the scalar-scalar induced GWs during RD do not behave simply as we have seen during MD, but exhibit rapid oscillations. 

We first consider the zero-shear gauge, in which the scalar-scalar source reads rather simply as


    s_ij^(ss)(k) 
    = 
    ∫d^3q/(2π)^3[
    2/a^2 q_iq_j φ_χ(q) φ_χ(k-q) 
    + 8π G (ρ+p) q_iq_j v_χ(q) v_χ(k-q) 
    ]
      .


With the scalar solutions given by (<ref>) and defining


    qη/√(3)
    =
    q/√(3)k kη≡q̅_1x
      and  | k- q|η/√(3)
    =
    | k- q|/√(3)k kη≡q̅_2x 
      ,


we can find


    e_λ^ij( k) s_ij^(ss)(k) 
    =
    4/a^2∫d^3q/(2π)^3[ e_λ^ij( k)q_iq_j ]
    (q) (k-q)
    {
    6 j_1(q̅_1x)/q̅_1xj_1(q̅_2x)/q̅_2x
    - 2 j_0(q̅_1x) j_1(q̅_2x)/q̅_2x
    - 2 j_1(q̅_1x)/q̅_1x  j_0(q̅_2x)
    + j_0(q̅_1x)j_0(q̅_2x)
    }  .


Thus, the inhomogeneous solution during RD that contains the integral including G_RD(η,η̃) is


    h_λ(η, k)
        = 
    4
    ∫d^3q/(2π)^3[ e_λ^ij( k) q_iq_j/k^2]
    (q) (k-q)
    ∫_0^x dx̃x̃^2
    {
    6 j_1(q̅_1x̃)/q̅_1x̃j_1(q̅_2x̃)/q̅_2x̃
    - 2 j_0(q̅_1x̃) j_1(q̅_2x̃)/q̅_2x̃
    - 2 j_1(q̅_1x̃)/q̅_1x̃  j_0(q̅_2x̃)
    
       
    + j_0(q̅_1x̃)j_0(q̅_2x̃)
    }[ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ]
      .


The integral over x̃ corresponds to the kernel for the scalar-scalar induced GWs during RD and can be performed analytically, given 
in Appendix <ref>.





In the comoving gauge, with the same short-handed notations as (<ref>), straightforward calculations give


    e_λ^ij( k) s_ij^(ss)(k) 
        =
    1/a^2∫d^3q/(2π)^3[ e_λ^ij( k)q_iq_j ]
    (q) (k-q)
    {
    j_1(q̅_1x) j_1(q̅_2x) [ q̅_1q̅_2x^2 + 2 q̅_2/q̅_1
    - 2/q̅_1q̅_2( 1 - 3q̅_1^2 ) ]
    
       
    + [ j_1(q̅_1x)/q̅_1x j_0(q̅_2x) 
    + j_0(q̅_1x) j_1(q̅_2x)/q̅_2x]
    x^2 ( 1 - 3q̅_1^2 + q̅_2^2 )
    + j_0(q̅_1x) j_0(q̅_2x) [ 1 - x^2/2( 1 - 3q̅_1^2 + q̅_2^2 ) ]
    }  ,


so that the solution for h_λ is written as


    h_λ(η, k)
        = 
    ∫d^3q/(2π)^3[ e_λ^ij( k) q_iq_j/k^2]
    (q) (k-q)
    ∫_0^x dx̃x̃^2
    {
    j_1(q̅_1x̃) j_1(q̅_2x̃) [ - q̅_1q̅_2x̃^2 + 2 q̅_2/q̅_1
    - 2/q̅_1q̅_2( 1 - 3q̅_1^2 ) ]
    
       
    + [ j_1(q̅_1x̃)/q̅_1x̃ j_0(q̅_2x̃) 
    + j_0(q̅_1x̃) j_1(q̅_2x̃)/q̅_2x̃]
    x̃^2 ( 1 - 3q̅_1^2 + q̅_2^2 )
    
       
    + j_0(q̅_1x̃) j_0(q̅_2x̃) 
    [ 1 - x̃^2/2( 1 - 3q̅_1^2 + q̅_2^2 ) ]
    }[ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ]
      .


The analytic results for the x̃-integral terms are given in Appendix <ref>. Comparing (<ref>) and (<ref>), we see that they are clearly different and the gauge dependence of the scalar-induced GWs are persistent during RD.

In Figure <ref> we show the kernels in the both gauges.

















































 §.§ Scalar-tensor induced GWs during RD










Next we consider the scalar-tensor induced GWs during RD. Again, we work first in the zero-shear gauge for the scalar perturbations. From the linear solutions during RD, the time derivatives of φ_χ and h_λ are


    φ̇_χ( q)
        =
    -2H j_2 ( qη/√(3)) 
      ,
    
    φ̈_χ( q)
        =
    2H^2 [ 5j_2 ( qη/√(3)) 
    - qη/√(3) j_1 ( qη/√(3)) ]
      ,
    
    ḣ_λ( q)
       
    =
    -Hqη h_0^λ( q) j_1(qη)
      ,
    
    ḧ_λ( q)
       
    =
    q^2/a^2 h_0^λ( q) j_2(qη)
      .


Then with qη = q̅_1x and | k- q|η = q̅_2x, we can write


    s_ij^(st)( k)
        =
    k^2/a^2∫d^3q/(2π)^3{2/x^2[ 2 j_2 ( q̅_1x/√(3)) 
    - q̅_1x/√(3) j_1( q̅_1x/√(3)) ] j_0(q̅_2x)
    
       
    + 2 ( 1 + | k- q|^2/k^2) 
    j_1(q̅_1x/√(3))/q̅_1x/√(3)  j_0(q̅_2x)
    }( q) h^λ'_0( k- q)
    e_ij^λ'( k- q)
    
         
    +
    k^2/a^2∫d^3q/(2π)^3{2/x^2 j_0 ( q̅_1x/√(3)) [ 2 j_2(q̅_2x) - q̅_2x j_1(q̅_2x) ] 
    + 2 ( 1 + q^2/k^2) 
    j_0 ( q̅_1x/√(3)) j_1(q̅_2x)/q̅_2x}( k- q) h^λ'_0( q)
    e_ij^λ'( q)
      .


We note that the first integral can be made identical to the second one upon defining k -  q≡ p, and then renaming the dummy integration variable p as q.

Thus, the scalar-tensor source term in the zero-shear gauge is, after multiplying the polarization tensor e^ij_λ( k),


    e^ij_λ( k) s_ij^(st)( k)
        =
    4k^2/a^2∫d^3q/(2π)^3[ e^ij_λ( k)e_ij^λ'( q) ]
    ( k- q) h^λ'_0( q)
    
       ×{1/x^2 j_0 ( q̅_1x̃/√(3)) 
    [ 2 j_2(q̅_2x) - q̅_2x j_1(q̅_2x) ] 
    + ( 1 + q^2/k^2) j_0 ( q̅_1x̃/√(3)) 
    j_1(q̅_2x)/q̅_2x
    
    
    
    }
      .


Then the solution is


    h_λ(η, k)
        =
    ∫d^3q/(2π)^3[ e^ij_λ( k)e_ij^λ'( q) ]
    ( k- q) h^λ'_0( q)
    
         ×
    4 ∫_0^x dx̃x̃^2 
    {1/x̃^2 j_0 ( q̅_1x̃/√(3)) 
    [ 2 j_2(q̅_2x̃) - q̅_2x̃ j_1(q̅_2x̃) ] 
    + ( 1 + q̅_1^2 ) j_0 ( q̅_1x̃/√(3)) 
    j_1(q̅_2x̃)/q̅_2x̃
    
    
    
    }[ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ]
      .


The x̃-integral can be performed analytically and the 
individual integrations are given in Appendix <ref>.






In the comoving gauge, similarly upon changing the dummy integration variable we can find the source term as


    e^ij_λ( k) s_ij^(st)( k)
        =
    2k^2/a^2∫d^3q/(2π)^3[ e^ij_λ( k)e_ij^λ'( q) ]
    ( k- q) h^λ'_0( q)
    [ 2(1-q̅_2^2) j_0(q̅_1x/√(3)) j_0(q̅_2x) 
    + ( -1 + q̅_1^2/3 + q̅_2^2 ) q̅_1x/√(3)
    j_1(q̅_1x/√(3)) j_0(q̅_2x) 
    
       
    + 2 q̅_1^2/3q̅_2 x j_0(q̅_1x/√(3)) j_1(q̅_2x) 
    + 2 q̅_1/√(3)/q̅_2( 1 - q̅_1^2 ) 
    j_1(q̅_1x/√(3)) j_1(q̅_2x) 
    ]
      ,


so that the analytic integral solution for h_λ is 


    h_λ(η, k)
        =
    ∫d^3q/(2π)^3[ e^ij_λ( k)e_ij^λ'( q) ]
    ( k- q) h^λ'_0( q)
    
         ×
    2 ∫_0^x dx̃x̃^2 
    [ 2(1-q̅_2^2) j_0(q̅_1x/√(3)) j_0(q̅_2x) 
    + ( -1 + q̅_1^2/3 + q̅_2^2 ) q̅_1x/√(3)
    j_1(q̅_1x/√(3)) j_0(q̅_2x) 
    
       
    + 2 q̅_1^2/3q̅_2 x j_0(q̅_1x/√(3)) j_1(q̅_2x) 
    + 2 q̅_1/√(3)/q̅_2( 1 - q̅_1^2 ) 
    j_1(q̅_1x/√(3)) j_1(q̅_2x) 
    ]
    [ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ]
      .


The details of each integration are given in Appendix <ref>. In Figure <ref>, we show both kernels, (<ref>) and (<ref>).








































 §.§ Tensor-tensor induced GWs during RD









Next we consider the tensor-tensor induced GWs during RD. Working in a similar manner to MD, multiplying the polarization tensor e^ij_λ( k) gives


    e^ij_λ( k) s_ij^(tt)( k)
        
    =
    1/a^2∫d^3q_1d^3q_2/(2π)^3δ^(3)( k- q_12)
    { h_0^λ_1( q_1) h_0^λ_2( q_2)
    e^ij_λ( k) q_1 q_2 e_i1^kλ_1 e_jk2^λ_2_≡_λ( q_1, q_2)
    
    j_1(q_1η) j_1(q_2η)
    
       
    +
    h_0^λ_1( q_1) h_0^λ_2( q_2)
    e^ij_λ( k)
    [
    e^kl_λ_11(
    -q_2jq_2k e^λ_2_il2 - q_2iq_2k e_jl2^λ_2 
    + q_2kq_2l e^λ_2_ij2)
    
       
    + 1/2 q_1iq_1j e_λ_1^kl1 e_kl2^λ_2 
    + ( q_1· q_2) e_i1^kλ_1 e_jk2^λ_2
    - q_1^lq_2k e_i1^kλ_1 e_jl2^λ_2] _≡_λ( q_1, q_2)
    
    j_0(q_1η) j_0(q_2η)
    }
         
    + ( q_1 ↔ q_2)
      .


Thus, with x=kη, q̅_1 ≡ q_1/k and q̅_2 ≡ q_2/k,


    h_λ(η, k)
        =
    ∫d^3q_1d^3q_2/(2π)^3δ^(3)( k- q_12) 1/k^2{[ _λ( q_1, q_2) + ( q_1 ↔ q_2) ]
    ∫_0^x dx̃x̃^2
    j_1(q̅_1x̃) j_1(q̅_2x̃)
    [ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ]
    
       
    +
    [ _λ( q_1, q_2) + ( q_1 ↔ q_2) ]
    ∫_0^x dx̃x̃^2
    j_0(q̅_1x̃) j_0(q̅_2x̃)
    [ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ] 
    }  .


The x̃-integrals can be performed analytically and finally we can write


    h_λ( k,η) 
        =
    ∫d^3q/(2π)^3 h_0^λ_1( q) h_0^λ_2( k- q)
    {
    e_λ^ij( k) e_i^kλ_1( q) e_jk^λ_2( k- q) 
    F_RD( k, q,η)
    
       
    +
    1/k^2(
    e_λ^ij( k)
    e^kl_λ_1( q) (
    q_jk_k e_il^λ_2( k- q) 
    + q_ik_k e_jl^λ_2( k- q)
    + k_kk_l e_ij^λ_2( k- q)
    )
    
       
    +
    e_λ^ij( k) 
    [
    1/2 q_i q_j e^kl_λ_1( q) e_kl^λ_2( k- q)
    +  q·( k- q) e_i^kλ_1( q) e_jk^λ_2( k- q)
    - q^lk_k e_i^kλ_1( q) e_jl^λ_2( k- q)
    ]
    ) 
    G_RD( k, q,η)
    }  ,


where the kernels F_RD and G_RD are given respectively by (<ref>) and (<ref>). In Figure <ref> we show F_RD and G_RD.




























§ CONCLUSIONS











In this article, we have presented the equation of motion for the tensor perturbations up to second order in perturbations, including all possible quadratic combinations of different types of cosmological perturbations. These terms serve as sources to generate second-order GWs. Given that the universe is filled with a perfect fluid matter with vanishing anisotropic stress, only linear scalar and tensor perturbations contribute to the source terms. And we have found the analytic integral solutions for the second-order GWs during both MD and RD induced by the scalar-scalar, scalar-tensor and tensor-tensor sources. The transition between MD and RD can be considered by separating the time integral of (<ref>). That is, for the transition from the epoch A to the epoch B at η = η_⋆, we may write <cit.>


    h_λ(η, k)
    =
    1/a{∫_0^η_⋆ dη̃ a^3(η)/a^3(η_⋆)[ a^3(η̃)  e^ij_λ(k) s^A_ij(k) ]
    G_A→ B(η,η̃)
    +
    ∫_η_⋆^η dη̃ [ a^3(η̃)  e^ij_λ(k) s^A→ B_ij(k) ]
    G_B(η,η̃)
    }  ,


where the first term represents the change of the propagation of the GWs produced during A through B, and the second term denotes the modification of the source. These terms can be found by matching, for the first term, the solutions and, for the second term, the kernels. However given the complicated source terms, it is a formidable task to compute the effects of the transition analytically, so we do not proceed any further but are satisfied with the above schematic form, which can in principle be further manipulated. 






Since the primary tensor perturbations are persistent irrespective of the sources, it is interesting to discuss if the tensor-induced GWs could ever be observationally significant. To address this question quantitatively, it is necessary to study further the relation between the tensor-induced GWs and the primary GWs. For example, if the primary GWs exhibit a sharp peak around at a certain scale k, the corresponding tensor-induced GWs are likely to be more prominent than the primary ones over a broader range of the wavenumber other than the peak. In this regards, the tensor-induced GWs may serve as a useful tool to probe the shape of the primary GWs and the relevant underlying physics. It is also interesting to consider how to test the origin of the induced GWs. We may consider, for example, the scalar-scalar-tensor bispectrum, ⟨ h⟩. If the tensor mode is scalar-induced, i.e. h ∼, the bispectrum is directly related to the local-type scalar bispectrum. If, on the other hand, the tensor mode is induced by the primary GWs, ⟨ h⟩ is proportional to the product of the scalar and tensor power spectra. Higher-order correlation functions may thus be able to serve as a possible tool to test the origin of the induced tensor perturbations.






As the kernels that involve rapid oscillations are analytically specified, the solutions given in this article should be useful for analytic and/or numerical studies of the second-order induced GWs.











































 §.§ Acknowledgements






I thank Jai-chan Hwang, Donghui Jeong, Kazunori Kohri, Sachiko Kuroyanagi, Takahiro Terada and Jaiyul Yoo for helpful comments and discussions. I also thank the anonymous referee for her/his constructive feedback that improved the content of this article. 

This work is supported in part by the Basic Science Research Program (2016R1D1A1B03930408) and Mid-career Research Program (2019R1A2C2085023) through the National Research Foundation of Korea Research Grants. 

I also acknowledge the Korea-Japan Basic Scientific Cooperation Program supported by the National Research Foundation of Korea and the Japan Society for the Promotion of Science (2018K2A9A2A08000127 and 2020K2A9A2A08000097) and the Ewha Womans University Research Grant of 2020 (1-2020-1630-001-1). 

I am grateful to the Asia Pacific Center for Theoretical Physics for Focus Research Program “The origin and evolution of the Universe” where parts of this work were presented and discussed.





























































§ TRACELESS EVOLUTION EQUATION FOR THE SPATIAL METRIC









The traceless evolution equation for the spatial metric is written compactly using the Arnowitt-Deser-Misner formulation <cit.>. With the metric


    ds^2 = -N^2(dx^0)^2 + γ_ij (N^idx^0+dx^i) (N^jdx^0+dx^j)   ,


where N, N^i and γ_ij denote respectively the lapse function, shift vector and spatial metric, the dynamics of the space-time is described by the spatial metric γ_ij through the curvature variables of the spatial hypersurfaces. Along with the matter contents residing in the space-time, the geometric equations for the curvature variables constitute a complete set of the equations of motion.
The extrinsic curvature K_ij is introduced as


    K_ij≡1/2N( N_i|j + N_j|i - γ_ij,0)   ,


where a vertical bar denotes a covariant derivative with respect to γ_ij. The evolution equation for the traceless part K_ij = K_ij - γ_ijK/3 with K ≡ K^i_i is


    K^i_j,0/N - K^i_j|kN^k/N 
    + K^k_jN^i_|k/N - K^i_kN^k_|j/N 
    = 
    KK^i_j - 1/N( N^|i_|j 
    - δ^i_j/3N^|k_|k) + R^i_j
    - 8π G T^i_j   ,


where R_ij is the traceless part of the intrinsic curvature tensor R_ij constructed from γ_ij. This is fully non-linear, geometric equation. To incorporate cosmological perturbations, we expand it perturbatively up to desired accuracy <cit.>, or write the exact non-linear equation for cosmological perturbations <cit.>. Using (<ref>) and expanding up to second order in perturbations, this equation is written as 



    ḧ^i_j + 3Hḣ^i_j - Δ/a^2h^i_j 
    + 1/a^2{[ 
    a^2 d/dt( χ^(i_,j)/a^2) + 3H χ^(i_,j)] 
    - δ^i_j/3[ 
    a^2 d/dt( χ^k_,k/a^2) + 3H χ^k_,k] }
       
    - 1/a^2( ∂^i∂_j - δ^i_j/3Δ)
    (α + φ)
    - 8π G ( Π^i_j - δ^i_j/3Π^k_k )
    
    
    =    d/dt[ 
    ( ḣ^i_j + χ^(i_,j)/a^2) α
    + 2 ( φδ^ik + h^ik) ( ḣ_jk + χ_(k,j)/a^2) + 2h^ikφ̇δ_jk
    + χ^k/a^2( φ^,iδ_jk + φ_,jδ^i_k - φ_,kδ^i_j
    + h^i_k,j + h_jk^,i - h^i_j,k) 
    ]
    
       
    + 3H [
    ( ḣ^i_j + χ^(i_,j)/a^2) α
    + 2 ( φδ^ik + h^ik) ( ḣ_jk + χ_(k,j)/a^2) + 2h^ikφ̇δ_jk
    + χ^k/a^2( φ^,iδ_jk + φ_,jδ^i_k - φ_,kδ^i_j
    + h^i_k,j + h_jk^,i - h^i_j,k) 
    ]
    
       
    + αd/dt( ḣ^i_j + χ^(i_,j)/a^2)
    - χ^k/a^2( ḣ^i_j,k + 1/a^2χ^(i_,j)k)
    + κ( ḣ^i_j + 1/a^2χ^(i_,j))
    
       
    - δ^i_j/3{d/dt[ 
    χ^k_,k/a^2α 
    + 2 ( φδ^kl + h^kl) χ_(k,l)/a^2 
    + 2h^klḣ_kl - 1/a^2χ^kφ_,k]
    + 3H [ 
    χ^k_,k/a^2α 
    + 2 ( φδ^kl + h^kl) χ_(k,l)/a^2 
    + 2h^klḣ_kl - 1/a^2χ^kφ_,k]
    
       
    + αd/dt( χ^k_,k/a^2)
    - χ^k/a^21/a^2χ^l_,lk
    + κχ^k_,k/a^2}
    + 1/a^2χ^i_,k( ḣ^k_j + 1/a^2χ^(k_,j)
    -1/a^2χ^k_,j( ḣ^i_k + 1/a^2χ^(i_,k)
    
       
    + 1/a^2{ 
    - αα^,i_,j
    + 1/2( -α^2 + 1/a^2χ^lχ_l )^,i_,j
    - 2 ( φδ^ik + h^ik) α_,jk  
    - ( 
    φ_,jδ^i_k  + h^i_k,j + φ^,iδ_jk + h_jk^,i
    - φ_,kδ^i_j - h^i_j,k) α^,k}
        -
    δ^i_j/3{1/a^2[
    - αΔα
    + 1/2Δ( -α^2 + 1/a^2χ^kχ_k )
    - 2 ( φδ^kl + h^kl) α_,kl
    + φ^,kα_,k] }
       
    + 1/a^2[ 
    - 3φ^,iφ_,i - 4φφ^,i_,j  
    - 4φΔh^i_j - 2h^i_jΔφ + φ_,kh^ik_,j + φ^,kh_jk^,i - 3φ^,kh^i_j,k
    + 2φ^,ikh_jk
           
    - 2h^ikΔh_kj - 2h^klh_kl,j^,i
    - h^kl,ih_kl,j - 2h^ik,lh_jk,l  
    + 2h^ik,lh_jl,k - 2h^kl( h^i_l,jk + h_jl,k^,i - h^i_j,kl) 
    
           
    - δ^i_j/3(
    -3φ^,kφ_,k - 4φΔφ + 2φ^,klh_,kl
    - 4h^klΔh_kl - 3h^kl,mh_kl,m + 2h^kl,mh_km,l) 
    ]
    
       
    - 16π G [
    ( φδ^ik + h^ik) Π_jk
    - δ^i_j/3( φδ^kl + h^kl)Π_kl]
    + 8π G (ρ+p) ( v^iv_j - δ^i_j/3v^kv_k )
      .


















§ INTEGRALS OF BESSEL FUNCTIONS












Here, a, b and c denote arbitrary positive constants. 











 §.§ Integrals for scalar-scalar and scalar-tensor induced GWs during MD








For s_ij^(ss) during MD, we have for the zero-shear gauge the following integrals for the spherical Bessel functions:


    ∫_0^x dx̃x̃^3 j_1(x̃) 
        = 
    3x^2j_1(x) - x^3j_0(x)
      ,
    
    ∫_0^x dx̃x̃^3 y_1(x̃) 
        = 
    3 + 3x^2y_1(x) - x^3y_0(x) 
      ,


along with the following identities of the spherical Bessel functions:


    ( 1/xd/dx)^m [ x^n+1f_n(x) ] 
        = 
    x^n-m+1f_n-m(x)
      ,
    
    
    f_n-1+f_n+1    = 
    2n+1/xf_n
      .


Also, we have the following integrals for the comoving gauge:


    ∫_0^x dx̃x̃^5 j_1(x̃) 
        = 
    5 x^2 ( x^2 - 6 ) j_1(x) - x^3 ( x^2 - 10 ) j_0(x)
      ,
    
    ∫_0^x dx̃x̃^5 y_1(x̃) 
        = 
    - 30 + 5x^2 ( x^2 - 6 ) y_1(x) - x^3 ( x^2 - 10 ) y_0(x)
      .







For s_ij^(st) during MD, for the zero-shear gauge, 


    ∫_0^x dx̃x̃^2 j_1(ax̃) j_1(x̃)
        =
    1/a^2-1[ x^2 j_1(ax) j_0(x) - ax^2j_0(ax)  j_1(x) ]
      ,
    
    ∫_0^x dx̃x̃^2 j_1(ax̃) y_1(x̃)
        =
    1/a^2-1[ - a + x^2 j_1(ax) y_0(x) - ax^2 j_0(ax) y_1(x) ]
      .


These give a rather simple result:


    ∫_0^x dx̃x̃^2 j_1(ax̃) [ j_1(x̃) y_1(x) - j_1(x) y_1(x̃) ]
    
        = 
    1/a^2-1[ aj_1(x) - j_1(ax) ]
      .


For the comoving gauge, we have the following integrals:


    ∫_0^x dx̃x̃^3 j_0(ax̃) j_1(x̃)
        =
    x^2/(a^2-1)^2{ j_0(x) [ (a^2-1)xj_0(ax) - 2aj_1(ax) ]
    + j_1(x) [ -(a^2-3)j_0(ax) + (a^2-1)axj_1(ax) ] }  ,
    
    ∫_0^x dx̃x̃^3 j_0(ax̃) y_1(x̃)
        =
    1/(a^2-1)^2{ 3 - a^2 + x^2j_0(ax) [ (a^2-1)xy_0(x) - (a^2-3)y_1(x) ]
    
       
    + ax^2j_1(ax) [ -2y_0(x) + (a^2-1)xy_1(x) ] }  .


Again, we find a rather simple result from these integrals:


    ∫_0^x dx̃x̃^3 j_0(ax̃) [ j_1(x̃) y_1(x) - j_1(x) y_1(x̃) ]
    
        =
    1/(a^2-1)^2[ (a^2-3) j_1(x) + 2aj_1(ax) - (a^2-1)xj_0(ax) ]
      .
















 §.§ Integrals for tensor-tensor induced GWs during MD









They can be arranged as, for the first two integrals,


    ∫_0^x dx̃x̃ j_2(ax̃) j_2(bx̃) j_1(cx̃)
    
        =
    (a^2+b^2-c^2)(a^2-b^2-c^2)x^2/16ab^2c j_1(ax)j_0(bx)j_0(cx)
    - (a^2+b^2-c^2)(a^2-b^2+c^2)x^2/16a^2bc j_0(ax)j_1(bx)j_0(cx)
    
         
    - (a^2+b^2-c^2)^2x^2/16a^2b^2 j_0(ax)j_0(bx)j_1(cx)
    + x/2a j_1(ax)j_0(bx)j_1(cx) + x/2b j_0(ax)j_1(bx)j_1(cx)
    - cx/2ab j_1(ax)j_1(bx)j_0(cx)
    
         
    - (a^2+b^2-c^2)x^2+12/8ab j_1(ax)j_1(bx)j_1(cx)
    
         
    + (a-b-c)(a+b-c)(a-b+c)(a+b+c)(a^2+b^2-c^2)/64a^3b^3c^2
           ×{ Si[(a-b-c)x] -  Si[(a+b-c)x]
    -  Si[(a-b+c)x] +  Si[(a+b+c)x] }  ,
    
       ∫_0^x dx̃x̃ j_2(ax̃) j_2(bx̃) y_1(cx̃)
    
        =
    (a^2+b^2-c^2)(a^2-b^2-c^2)x^2/16ab^2c j_1(ax)j_0(bx)y_0(cx)
    - (a^2+b^2-c^2)(a^2-b^2+c^2)x^2/16a^2bc j_0(ax)j_1(bx)y_0(cx)
    
         
    - (a^2+b^2-c^2)^2x^2/16a^2b^2 j_0(ax)j_0(bx)y_1(cx)
    + x/2a j_1(ax)j_0(bx)y_1(cx) + x/2b j_0(ax)j_1(bx)y_1(cx)
    - cx/2ab j_1(ax)j_1(bx)y_0(cx)
    
         
    - (a^2+b^2-c^2)x^2+12/8ab j_1(ax)j_1(bx)y_1(cx)
    
         
    - (a-b-c)(a+b-c)(a-b+c)(a+b+c)(a^2+b^2-c^2)/64a^3b^3c^2
           ×{ Cin[(a-b-c)x] -  Cin[(a+b-c)x]
    +  Cin[(a-b+c)x] -  Cin[(a+b+c)x] }
         
    - 3a^4+3(b^2-c^2)^2-2a^2(b^2+3c^2)/48a^2b^2c^2  .


Note that with the definition


    Cin(x) ≡∫_0^x 1-cost/t dt = γ + logx -  Ci(x) 
      ,


where γ≈ 0.577216 is the Euler-Mascheroni constant, while both Ci(x) and logx are diverging as x→0, lim_x→0 Cin(x) = 0. Then we find  


    1/x∫_0^x dx̃x̃ j_2(ax̃) j_2(bx̃)
    [ j_1(x̃)y_1(x) - j_1(x)y_1(x̃) ]
    ≡ 
    F_MD(a,b,x)
    
        =
    (a^2+b^2-1)(-a^2+b^2+1)/16b^2j_1(ax)/ax j_0(bx)
    + (a^2+b^2-1)(a^2-b^2+1)/16a^2 j_0(ax) j_1(bx)/bx
    + 1/2j_1(ax)/axj_1(bx)/bx
         
    + 3a^4+3b^4+3 - 6a^2-6b^2 - 2a^2b^2/48a^2b^2j_1(x)/x
         
    + (a-b-1)(a+b-1)(a-b+1)(a+b+1)(a^2+b^2-1)/64a^3b^3
           ×{(  Si[(a-b-1)x] y_1(x)/x +  Cin[(a-b-1)x] j_1(x)/x) 
    - (  Si[(a+b-1)x] y_1(x)/x +  Cin[(a+b-1)x] j_1(x)/x) 
    
       
    - (  Si[(a-b+1)x] y_1(x)/x -  Cin[(a-b+1)x] j_1(x)/x) 
    + (  Si[(a+b+1)x] y_1(x)/x -  Cin[(a+b+1)x] j_1(x)/x) }  .








Likewise, for the next two integrals,


    ∫_0^x dx̃x̃ j_1(ax̃) j_1(bx̃) j_1(cx̃)
    
        =
    -x^2/4 j_1(ax)j_1(bx)j_1(cx)
    
         
    - (a^2+b^2-c^2)x^2/8ab j_0(ax)j_0(bx)j_1(cx)
    - (-a^2+b^2+c^2)x^2/8ab j_1(ax)j_0(bx)j_0(cx)
    - (a^2-b^2+c^2)x^2/8ab j_0(ax)j_1(bx)j_0(cx)
    
         
    + (a-b-c)(a+b-c)(a-b+c)(a+b+c)/32a^2b^2c^2{ Si[(a-b-c)x] -  Si[(a+b-c)x]
    -  Si[(a-b+c)x] +  Si[(a+b+c)x] }  ,
    
       ∫_0^x dx̃x̃ j_1(ax̃) j_1(bx̃) y_1(cx̃)
    
        =
    -x^2/4 j_1(ax)j_1(bx)y_1(cx)
    
         
    - (a^2+b^2-c^2)x^2/8ab j_0(ax)j_0(bx)y_1(cx)
    - (-a^2+b^2+c^2)x^2/8ab j_1(ax)j_0(bx)y_0(cx)
    - (a^2-b^2+c^2)x^2/8ab j_0(ax)j_1(bx)y_0(cx)
    
         
    - (a-b-c)(a+b-c)(a-b+c)(a+b+c)/32a^2b^2c^2{ Cin[(a-b-c)x] -  Cin[(a+b-c)x]
    +  Cin[(a-b+c)x] -  Cin[(a+b+c)x] }
         
    - a^2+b^2-c^2/8abc^2  .


Thus,


    1/abx∫_0^x dx̃x̃ j_1(ax̃) j_1(bx̃)
    [ j_1(x̃)y_1(x) - j_1(x)y_1(x̃) ]
    ≡
    G_MD(a,b,x)
    
        =
    (-a^2+b^2+1)/8b^2j_1(ax)/ax j_0(bx) 
    + (a^2-b^2+1)/8a^2 j_0(ax) j_1(bx)/bx
    + a^2+b^2-1/8a^2b^2j_1(x)/x
         
    + (a-b-1)(a+b-1)(a-b+1)(a+b+1)/32a^3b^3
           ×{(  Si[(a-b-1)x] y_1(x)/x +  Cin[(a-b-1)x] j_1(x)/x) 
    - (  Si[(a+b-1)x] y_1(x)/x +  Cin[(a+b-1)x] j_1(x)/x) 
    
       
    - (  Si[(a-b+1)x] y_1(x)/x -  Cin[(a-b+1)x] j_1(x)/x) 
    + (  Si[(a+b+1)x] y_1(x)/x -  Cin[(a+b+1)x] j_1(x)/x) }  .


Notice that comparing with (<ref>), F_MD(a,b,x) and G_MD(a,b,x) are related by


    a^2+b^2-1/2G_MD(a,b,x)
    =
    F_MD(a,b,x) - 1/2j_1(ax)/axj_1(bx)/bx
    + 1/6j_1(x)/x  .

















 §.§ Integrals for scalar-scalar induced GWs during RD









We find


    ∫_0^x dx̃x̃^2 j_0(ax̃)j_0(bx̃)j_0(cx̃)
        =
    -1/4abc{ Si[ (a-b-c)x] -  Si[ (a+b-c)x] -  Si[ (a-b+c)x] +  Si[ (a+b+c)x]
    }  ,
    
    ∫_0^x dx̃x̃^2 j_0(ax̃)j_0(bx̃)y_0(cx̃)
        =
    1/4abc{ Cin[(a-b-c)x] -  Cin[(a+b-c)x] +  Cin[(a-b+c)x] -  Cin[(a+b+c)x]
    }  ,



    ∫_0^x dx̃x̃
    j_1(ax̃)j_0(bx̃) j_0(cx̃)
    
       
    =
    -x/a j_0(ax) j_0(bx) j_0(cx)
    - 1/2x^2 j_1(ax) j_0(bx) j_0(cx)
    + b/2ax^2 j_0(ax) j_1(bx) j_0(cx)
    + c/2ax^2 j_0(ax) j_0(bx) j_1(cx)
    
         
    + 1/8a^2bc{
    - ( a^2-(b+c)^2 )  Si[(a-b-c)x] + ( a^2-(b-c)^2 )  Si[(a+b-c)x]
    
         
    + ( a^2-(b-c)^2 )  Si[(a-b+c)x] - ( a^2-(b+c)^2 )  Si[(a+b+c)x]
    }  ,
    
       ∫_0^x dx̃x̃
    j_1(ax̃)j_0(bx̃) y_0(cx̃)
    
       
    = 
    -x/a j_0(ax) j_0(bx) y_0(cx)
    - 1/2x^2 j_1(ax) j_0(bx) y_0(cx)
    + b/2ax^2 j_0(ax) j_1(bx) y_0(cx)
    + c/2ax^2 j_0(ax) j_0(bx) y_1(cx)
    
         
    + 1/8a^2bc{( a^2-(b+c)^2 )  Cin[(a-b-c)x] - ( a^2-(b-c)^2 )  Cin[(a+b-c)x]
    
         
    + ( a^2-(b-c)^2 )  Cin[(a-b+c)x] - ( a^2-(b+c)^2 )  Cin[(a+b+c)x]
    }
    - 1/2ac  ,



    ∫_0^x dx̃
    j_1(ax̃)j_1(bx̃) j_0(cx̃)
    
       
    = 
    - a^2+b^2/3abx j_0(ax) j_0(bx) j_0(cx)
    + a^2-b^2/8bx^2 j_1(ax) j_0(bx) j_0(cx)
    + -a^2+b^2/8ax^2 j_0(ax) j_1(bx) j_0(cx)
    
         
    + c^2/24bx^2 j_1(ax) j_0(bx) j_0(cx)
    + c^2/24ax^2 j_0(ax) j_1(bx) j_0(cx)
    + c(5a^2+5b^2-c^2)/24abx^2 j_0(ax) j_0(bx) j_1(cx)
    
         
    - 1/3x j_1(ax) j_1(bx) j_0(cx)
    + c/12x^2 j_1(ax) j_1(bx) j_1(cx)
    
         
    + 1/96a^2b^2c{
    (a-b-c)^2 (3a^2 + 6ab + 3b^2 - 2ac + 2bc - c^2)  Si[(a-b-c)x]
    
         
    - (a+b-c)^2 (3a^2 - 6ab + 3b^2 - 2ac - 2bc - c^2)  Si[(a+b-c)x]
    
         
    - (a-b+c)^2 (3a^2 + 6ab + 3b^2 + 2ac - 2bc - c^2)  Si[(a-b+c)x]
    
         
    + (a+b+c)^2 (3a^2 - 6ab + 3b^2 + 2ac + 2bc - c^2)  Si[(a+b+c)x]
    }  ,
    
       ∫_0^x dx̃
    j_1(ax̃)j_1(bx̃) y_0(cx̃)
    
       
    = 
    - a^2+b^2/3abx j_0(ax) j_0(bx) y_0(cx)
    + a^2-b^2/8bx^2 j_1(ax) j_0(bx) y_0(cx)
    + -a^2+b^2/8ax^2 j_0(ax) j_1(bx) y_0(cx)
    
         
    + c^2/24bx^2 j_1(ax) j_0(bx) y_0(cx)
    + c^2/24ax^2 j_0(ax) j_1(bx) y_0(cx)
    + c(5a^2+5b^2-c^2)/24abx^2 j_0(ax) j_0(bx) y_1(cx)
    
         
    - 1/3x j_1(ax) j_1(bx) y_0(cx)
    + c/12x^2 j_1(ax) j_1(bx) y_1(cx)
    
         
    + 1/96a^2b^2c{
    -(a-b-c)^2 (3a^2 + 6ab + 3b^2 - 2ac + 2bc - c^2)  Cin[(a-b-c)x]
    
         
    + (a+b-c)^2 (3a^2 - 6ab + 3b^2 - 2ac - 2bc - c^2)  Cin[(a+b-c)x]
    
         
    - (a-b+c)^2 (3a^2 + 6ab + 3b^2 + 2ac - 2bc - c^2)  Cin[(a-b+c)x]
    
         
    + (a+b+c)^2 (3a^2 - 6ab + 3b^2 + 2ac + 2bc - c^2)  Cin[(a+b+c)x]
    }
    - 3a^2+3b^2+c^2/24abc  .













































For the comoving gauge, we can find


    ∫_0^x dx̃x̃^3 j_1(ax̃) j_0(bx̃) j_0(cx̃)
        =
    x/4abc{ j_0[(a-b-c)x] - j_0[(a+b-c)x] - j_0[(a-b+c)x] + j_0[(a+b+c)x] }
         
    + 1/4a^2bc{ Si[(a-b-c)x] +  Si[(a+b-c)x] +  Si[(a-b+c)x] -  Si[(a+b+c)x] }  ,
    
    ∫_0^x dx̃x̃^3 j_1(ax̃) j_0(bx̃) y_0(cx̃)
        =
    x/4abc{ -y_0[(a-b-c)x] + y_0[(a+b-c)x] - y_0[(a-b+c)x] + y_0[(a+b+c)x] }
         
    + 1/4a^2bc{ Cin[(a-b-c)x] -  Cin[(a+b-c)x] +  Cin[(a-b+c)x] -  Cin[(a+b+c)x] }
          
    - a^2-b^2+c^2/ac(a-b-c)(a+b-c)(a-b+c)(a+b+c)  ,



    ∫_0^x dx̃x̃^4 j_0(ax̃) j_0(bx̃) j_0(cx̃)
        =
    x^2/4abc{ -j_1[(a-b-c)x] + j_1[(a+b-c)x] + j_1[(a-b+c)x] - j_1[(a+b+c)x] }  ,
    
    ∫_0^x dx̃x̃^4 j_0(ax̃) j_0(bx̃) y_0(cx̃)
        =
    x^2/4abc{ y_1[(a-b-c)x] - y_1[(a+b-c)x] + y_1[(a-b+c)x] - y_1[(a+b+c)x] }
         
    + 2 ( a^4 + b^4 - 3c^4 - 2a^2b^2 + 2a^2c^2 + 2b^2c^2 )/(a-b-c)^2(a+b-c)^2(a-b+c)^2(a+b+c)^2  ,



    ∫_0^x dx̃x̃^4 j_1(ax̃) j_1(bx̃) j_0(cx̃)
        =
    x^2/4abc{ - y_0[(a-b-c)x] - y_0[(a+b-c)x] + y_0[(a-b+c)x] + y_0[(a+b+c)x] }
         
    + x^2/4a^2b^2c{( a^2 - 3ab + b^2 - ac + bc ) [ j_1[(a-b-c)x] - y_0[(a-b-c)x] ]
    
       
    - ( a^2 + 3ab + b^2 - ac - bc ) [ j_1[(a+b-c)x] - y_0[(a+b-c)x] ]
    
       
    - ( a^2 - 3ab + b^2 + ac - bc ) [ j_1[(a-b+c)x] - y_0[(a-b+c)x] ]
    
       
    + ( a^2+ 3ab + b^2 + ac + bc ) [ j_1[(a+b+c)x] - y_0[(a+b+c)x] ]
    }
         
    + 1/4a^2b^2c{ -  Si[(a-b-c)x] +  Si[(a+b-c)x] +  Si[(a-b+c)x] -  Si[(a+b+c)x] }  ,
    
    ∫_0^x dx̃x̃^4 j_1(ax̃) j_1(bx̃) y_0(cx̃)
        =
    - x^2/4abc{ j_0[(a-b-c)x] + j_0[(a+b-c)x] + j_0[(a-b+c)x] + j_0[(a+b+c)x] }
         
    + x^2/4a^2b^2c{
    - ( a^2 - 3ab + b^2 - ac + bc ) [ y_1[(a-b-c)x] + j_0[(a-b-c)x] ]
    
       
    + ( a^2 + 3ab + b^2 - ac - bc ) [ y_1[(a+b-c)x] + j_0[(a+b-c)x] ]
    
       
    - ( a^2 - 3ab + b^2 + ac - bc ) [ y_1[(a-b+c)x] + j_0[(a-b+c)x] ]
    
       
    + ( a^2+ 3ab + b^2 + ac + bc ) [ y_1[(a+b+c)x] + j_0[(a+b+c)x] ]
    }
         
    + a^6 + (b^2-c^2)^3 - a^4(b^2+3c^2) + a^2(-b^4+14b^2c^2+3c^4)/abc(a-b-c)^2(a+b-c)^2(a-b+c)^2(a+b+c)^2  .




















 §.§ Integrals for scalar-tensor induced GW during RD









We find


    ∫_0^x dx̃x̃ j_1(ax̃) j_1(bx̃) j_0(cx̃)
    
       
    = 
    1/6ab j_0(ax) j_0(bx) j_0 (cx)
    - a/3bx^2 j_0(ax) j_0(bx) j_0 (cx)
    - b/3ax^2 j_0(ax) j_0(bx) j_0 (cx)
    + c^2/6abx^2 j_0(ax) j_0(bx) j_0 (cx)
    
         
    - 1/6bx j_1(ax) j_0(bx) j_0 (cx)
    - 1/6ax j_0(ax) j_1(bx) j_0 (cx)
    - 1/3x^2 j_1(ax) j_1(bx) j_0 (cx)
    - c/6abx j_0(ax) j_0(bx) j_1(cx)
    
         
    + c/6bx^2 j_1(ax) j_0(bx) j_1(cx)
    + c/6ax^2 j_0(ax) j_1(bx) j_1(cx)
    
         
    + 1/24a^2b^2c{
    (a-b-c) ( 2a^2 + 2ab + 2b^2 - ac + bc - c^2 )  Cin[(a-b-c)x]
    
         
    - (a+b-c) ( 2a^2 - 2ab + 2b^2 - ac - bc - c^2 )  Cin[(a+b-c)x]
    
         
    - (a-b+c) ( 2a^2 + 2ab + 2b^2 + ac - bc - c^2 )  Cin[(a-b+c)x]
    
         
    + (a+b+c) ( 2a^2 - 2ab + 2b^2 + ac + bc - c^2 )  Cin[(a+b+c)x]
    }
    - 1/6ab  ,



    ∫_0^x dx̃x̃ j_1(ax̃) j_1(bx̃) y_0(cx̃)
    
       
    = 
    1/6ab j_0(ax) j_0(bx) y_0(cx)
    - a/3bx^2 j_0(ax) j_0(bx) y_0(cx)
    - b/3ax^2 j_0(ax) j_0(bx) y_0(cx)
    + c^2/6abx^2 j_0(ax) j_0(bx) y_0(cx)
    
         
    - 1/6bx j_1(ax) j_0(bx) y_0(cx)
    - 1/6ax j_0(ax) j_1(bx) y_0(cx)
    - 1/3x^2 j_1(ax) j_0(bx) y_0(cx)
    - c/6abx j_0(ax) j_0(bx) y_1(cx)
    
         
    + c/6bx^2 j_1(ax) j_0(bx) y_1(cx)
    + c/6ax^2 j_0(ax) j_1(bx) y_1(cx)
    
         
    + 1/24a^2b^2c{
    (a-b-c) ( 2a^2 + 2ab + 2b^2 - ac + bc - c^2 )  Si[(a-b-c)x]
    
         
    - (a+b-c) ( 2a^2 - 2ab + 2b^2 - ac - bc - c^2 )  Si[(a+b-c)x]
    
         
    + (a-b+c) ( 2a^2 + 2ab + 2b^2 + ac - bc - c^2 )  Si[(a-b+c)x]
    
         
    - (a+b+c) ( 2a^2 - 2ab + 2b^2 + ac + bc - c^2 )  Si[(a+b+c)x]
    }  ,




    ∫_0^x dx̃x̃ j_1(ax̃) j_0(bx̃) j_0(cx̃)
    
       
    =
    - 1/ax j_0(ax) j_0(bx) j_0(cx)
    - 1/2x^2 j_1(ax) j_0(bx) j_0(cx)
    + b/2ax^2 j_0(ax) j_1(bx) j_0(cx)
    + c/2ax^2 j_0(ax) j_0(bx) j_1(cx)
    
         
    + 1/8a^2bc{
    - [ a^2 - (b+c)^2 ]  Si[(a-b-c)x]
    + [ a^2 - (b-c)^2 ]  Si[(a+b-c)x]
    
         
    + [ a^2 - (b-c)^2 ]  Si[(a-b+c)x]
    - [ a^2 - (b+c)^2 ]  Si[(a+b+c)x]
    }  ,



    ∫_0^x dx̃x̃ j_1(ax̃) j_0(bx̃) y_0(cx̃)
    
       
    =
    - 1/ax j_0(ax) j_0(bx) y_0(cx)
    - 1/2x^2 j_1(ax) j_0(bx) y_0(cx)
    + b/2ax^2 j_0(ax) j_1(bx) y_0(cx)
    + c/2ax^2 j_0(ax) j_0(bx) y_1(cx)
    
         
    + 1/8a^2bc{[ a^2 - (b+c)^2 ]  Cin[(a-b-c)x]
    - [ a^2 - (b-c)^2 ]  Cin[(a+b-c)x]
    
         
    + [ a^2 - (b-c)^2 ]  Cin[(a-b+c)x]
    - [ a^2 - (b+c)^2 ]  Cin[(a+b+c)x]
    }
    - 1/2ac  ,



    ∫_0^x dx̃ j_2(ax̃) j_0(bx̃) j_0(cx̃)
    
       
    = 
    b^2/2a^2x j_0(ax) j_0(bx) j_0(cx)
    + c^2/2a^2x j_0(ax) j_0(bx) j_0(cx)
    - 3/2a j_1(ax) j_0(bx) j_0(cx)
    - a/8x^2 j_1(ax) j_0(bx) j_0(cx) 
    
         
    + b^2/8ax^2 j_1(ax) j_0(bx) j_0(cx)
    + c^2/8ax^2 j_1(ax) j_0(bx) j_0(cx)
    + b/8x^2 j_0(ax) j_1(bx) j_0(cx)
    - b^3/8a^2x^2 j_0(ax) j_1(bx) j_0(cx)
    
         
    - 3bc^2/8a^2x^2 j_0(ax) j_1(bx) j_0(cx)
    + b/2ax j_1(ax) j_1(bx) j_0(cx)
    + c/8x^2 j_0(ax) j_0(bx) j_1(cx)
    - 3b^2c/8a^2x^2 j_0(ax) j_0(bx) j_1(cx)
    
         
    - c^3/8a^2x^2 j_0(ax) j_0(bx) j_1(cx)
    + c/2ax j_1(ax) j_0(bx) j_1(cx)
    - bc/4ax^2 j_1(ax) j_1(bx) j_1(cx)
    
         
    + 1/32a^3bc{
    - [ a^2-(b+c)^2 ]^2  Si[(a-b-c)x]
    + [ a^2-(b-c)^2 ]^2  Si[(a+b-c)x]
    
         
    + [ a^2-(b-c)^2 ]^2  Si[(a-b+c)x]
    - [ a^2-(b+c)^2 ]^2  Si[(a+b+c)x]
    }  ,



    ∫_0^x dx̃ j_2(ax̃) j_0(bx̃) y_0(cx̃)
    
       
    = 
    b^2/2a^2x j_0(ax) j_0(bx) y_0(cx)
    + c^2/2a^2x j_0(ax) j_0(bx) y_0(cx)
    - 3/2a j_1(ax) j_0(bx) y_0(cx)
    - a/8x^2 j_1(ax) j_0(bx) y_0(cx)
    
         
    + b^2/8ax^2 j_1(ax) j_0(bx) y_0(cx)
    + c^2/8ax^2 j_1(ax) j_0(bx) y_0(cx)
    + b/8x^2 j_0(ax) j_1(bx) y_0(cx)
    - b^3/8a^2x^2 j_0(ax) j_1(bx) y_0(cx)
    
         
    - 3bc^2/8a^2x^2 j_0(ax) j_1(bx) y_0(cx)
    + b/2ax j_1(ax) j_1(bx) y_0(cx)
    + c/8x^2 j_0(ax) j_0(bx) y_1(cx)
    - 3b^2c/8a^2x^2 j_0(ax) j_0(bx) y_1(cx)
    
         
    - c^3/8a^2x^2 j_0(ax) j_0(bx) y_1(cx)
    + c/2ax j_1(ax) j_0(bx) y_1(cx)
    - bc/4a^2x^2 j_1(ax) j_1(bx) y_1(cx)
    
         
    + 1/32a^3bc{[ a^2-(b+c)^2 ]^2  Cin[(a-b-c)x]
    - [ a^2-(b-c)^2 ]^2  Cin[(a+b-c)x]
    
         
    + [ a^2-(b-c)^2 ]^2  Cin[(a-b+c)x]
    - [ a^2-(b+c)^2 ]^2  Cin[(a+b+c)x]
    }
    - 5/24c + b^2/8a^2c + 3c/8a^2  .








































































































































 §.§ Integrals for tensor-tensor induced GWs during RD








We can perform the integrals analytically to find


    ∫_0^x dx̃x̃^2 j_1(ax̃)j_1(bx̃)j_0(cx̃)
    
        =
    c/2abx^2 j_0(ax) j_0(bx) j_1(cx) 
    - 1/2a x^2 j_0(ax) j_1(bx) j_0(cx)
    - 1/2b x^2 j_1(ax) j_0(bx) j_0(cx)
    
         
    - a^2+b^2-c^2/8a^2b^2c{ Si[ (a-b-c)x] -  Si[ (a+b-c)x] -  Si[ (a-b+c)x] +  Si[ (a+b+c)x]
    }  ,
    
       ∫_0^x dx̃x̃^2 j_1(ax̃)j_1(bx̃)y_0(cx̃)
    
        =
    c/2abx^2 j_0(ax) j_0(bx) y_1(cx)
    - 1/2ax^2 j_0(ax) j_1(bx) y_0(cx)
    - 1/2bx^2 j_1(ax) j_0(bx) y_0(cx)
    + 1/2abc
         
    + a^2+b^2-c^2/8a^2b^2c{ Cin[(a-b-c)x] -  Cin[(a+b-c)x] +  Cin[(a-b+c)x] -  Cin[(a+b+c)x]
    }
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
      .


Thus,


    ∫_0^x dx̃x̃^2
    j_1(ax̃) j_1(bx̃)
    [ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ]
    ≡1/2abF_RD(a,b,x)
    
        =
    1/2ab[ j_0(ax) j_0(bx) -  j_0(x)
    + a^2+b^2-1/4ab{( - Si[ (a-b-1)x] y_0(x) -  Cin[(a-b-1)x] j_0(x) ) 
    
       
    - ( - Si[ (a+b-1)x] y_0(x) -  Cin[(a+b-1)x] j_0(x) ) 
    
       
    - ( - Si[ (a-b+1)x] y_0(x) +  Cin[(a-b+1)x] j_0(x) ) 
    
       
    + ( - Si[ (a+b+1)x] y_0(x) +  Cin[(a+b+1)x] j_0(x) ) 
    }]
      ,
    
       ∫_0^x dx̃x̃^2
    j_0(ax̃) j_0(bx̃)
    [ j_0(x̃)y_0(x) - j_0(x)y_0(x̃) ] 
    ≡1/2G_RD(a,b,x)
    
        =
    1/4ab{( - Si[ (a-b-1)x] y_0(x) -  Cin[(a-b-1)x] j_0(x) ) 
    + (  Si[ (a+b-1)x] y_0(x) +  Cin[(a+b-1)x] j_0(x) ) 
    
       
    + (  Si[ (a-b+1)x] y_0(x) -  Cin[(a-b+1)x] j_0(x) ) 
    + ( -  Si[ (a+b+1)x] y_0(x) +  Cin[(a+b+1)x] j_0(x) ) 
    }  .


Note that


    a^2+b^2-1/2G_RD(a,b,x)
    =
    F_RD(a,b,x) - j_0(ax)j_0(bx) + j_0(x)
      .





























99





[Abbott et al.(2016a)]2016PhRvL.116f1102A Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, , 116, 061102





[Abbott et al.(2016b)]2016PhRvL.116x1103A Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, , 116, 241103





[Abbott et al.(2017a)]2017PhRvL.118v1101A Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, , 118, 221101





[Abbott et al.(2017b)]2017PhRvL.119n1101A Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, , 119, 141101





[Abbott et al.(2017c)]2017PhRvL.119p1101A Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, , 119, 161101

  
  
  
  
[Abbott et al.(2017d)]2017ApJ...851L..35A Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, , 851, L35





[Alabidi et al.(2012)]2012JCAP...09..017A Alabidi, L., Kohri, K., Sasaki, M., et al. 2012, , 2012, 017





[Alabidi et al.(2013)]2013JCAP...05..033A Alabidi, L., Kohri, K., Sasaki, M., et al. 2013, , 2013, 033





[Ananda et al.(2007)]2007PhRvD..75l3518A Ananda, K. N., Clarkson, C., & Wands, D. 2007, , 75, 123518





[Arnowitt et al.(2008)]2008GReGr..40.1997A Arnowitt, R., Deser, S., & Misner, C. W. 2008, General Relativity and Gravitation, 40, 1997





[Arroja et al.(2009)]2009PhRvD..80l3526A Arroja, F., Assadullahi, H., Koyama, K., et al. 2009, , 80, 123526





[Assadullahi, & Wands(2009)]2009PhRvD..79h3511A Assadullahi, H., & Wands, D. 2009, , 79, 083511





[Assadullahi, & Wands(2010)]2010PhRvD..81b3527A Assadullahi, H., & Wands, D. 2010, , 81, 023527





[Baumann et al.(2007)]2007PhRvD..76h4019B Baumann, D., Steinhardt, P., Takahashi, K., et al. 2007, , 76, 084019





[BICEP2 Collaboration et al.(2016)]2016PhRvL.116c1302B BICEP2 Collaboration, Keck Array Collaboration, Ade, P. A. R., et al. 2016, , 116, 031302





[Brandenberger et al.(1986)]1986NuPhB.277..605B Brandenberger, R. H., Albrecht, A., & Turok, N. 1986, Nuclear Physics B, 277, 605




[Brandenberger et al.(2007)]2007PhRvL..98w1302B Brandenberger, R. H., Nayeri, A., Patil, S. P., et al. 2007, , 98, 231302




[Cai et al.(2015)]2015NuPhB.900..517C Cai, Y.-F., Gong, J.-O., Pi, S., et al. 2015, Nuclear Physics B, 900, 517




[De Luca et al.(2020)]2020JCAP...03..014D De Luca, V., Franciolini, G., Kehagias, A., et al. 2020, , 2020, 014





[Espinosa et al.(2018)]2018JCAP...09..012E Espinosa, J. R., Racco, D., & Riotto, A. 2018, , 2018, 012





[Gong, & Stewart(2002)]2002PhLB..538..213G Gong, J.-O., & Stewart, E. D. 2002, Physics Letters B, 538, 213




[Gong(2014)]2014JCAP...07..022G Gong, J.-O. 2014, , 2014, 022





[Gong et al.(2017)]2017JCAP...10..027G Gong, J.-O., Hwang, J.-. chan ., Noh, H., et al. 2017, , 2017, 027





[Hogan(1986)]1986MNRAS.218..629H Hogan, C. J. 1986, , 218, 629





[Hwang(1994)]1994ApJ...427..533H Hwang, J.-C. 1994, , 427, 533





[Hwang, & Noh(2007)]2007PhRvD..76j3527H Hwang, J.-C., & Noh, H. 2007, , 76, 103527





[Hwang et al.(2012)]2012ApJ...752...50H Hwang, J.-. chan ., Noh, H., & Gong, J.-O. 2012, , 752, 50





[Hwang et al.(2017)]2017ApJ...842...46H Hwang, J.-. chan ., Jeong, D., & Noh, H. 2017, , 842, 46






[Inomata et al.(2019a)]2019JCAP...10..071I Inomata, K., Kohri, K., Nakama, T., et al. 2019, , 2019, 071





[Inomata et al.(2019b)]2019PhRvD.100d3532I Inomata, K., Kohri, K., Nakama, T., et al. 2019, , 100, 043532




[Inomata & Terada(2020)]2020PhRvD.101b3523I Inomata, K. & Terada, T. 2020, , 101, 023523




[Jain et al.(2009)]2009JCAP...01..009J Jain, R. K., Chingangbam, P., Gong, J.-O., et al. 2009, , 2009, 009




[Jain et al.(2010)]2010PhRvD..82b3509J Jain, R. K., Chingangbam, P., Sriramkumar, L., et al. 2010, , 82, 023509





[Jedamzik et al.(2010)]2010JCAP...04..021J Jedamzik, K., Lemoine, M., & Martin, J. 2010, , 2010, 021





[Khlebnikov, & Tkachev(1997)]1997PhRvD..56..653K Khlebnikov, S., & Tkachev, I. 1997, , 56, 653




[Kobayashi et al.(2010)]2010PhRvL.105w1302K Kobayashi, T., Yamaguchi, M., & Yokoyama, J. 2010, , 105, 231302





[Kohri, & Terada(2018)]2018PhRvD..97l3532K Kohri, K., & Terada, T. 2018, , 97, 123532





[Mollerach et al.(2004)]2004PhRvD..69f3002M Mollerach, S., Harari, D., & Matarrese, S. 2004, , 69, 063002





[Mukhanov(2005)]2005pfc..book.....M Mukhanov, V. 2005, Physical Foundations of Cosmology, by Viatcheslav Mukhanov, pp. 442. Cambridge University Press, November 2005. 




[Mukhanov, & Vikman(2006)]2006JCAP...02..004M Mukhanov, V., & Vikman, A. 2006, , 2006, 004




[Mylova et al.(2018)]2018JCAP...12..024M Mylova, M., Özsoy, O., Parameswaran, S., et al. 2018, , 2018, 024




[Nakama & Suyama(2015)]2015PhRvD..92l1304N Nakama, T. & Suyama, T. 2015, , 92, 121304




[Nakama & Suyama(2016)]2016PhRvD..94d3507N Nakama, T. & Suyama, T. 2016, , 94, 043507





[Noh, & Hwang(2004)]2004PhRvD..69j4011N Noh, H., & Hwang, J.-C. 2004, , 69, 104011





[Pi et al.(2019)]2019JCAP...06..049P Pi, S., Sasaki, M., & Zhang, Y.-. li . 2019, , 2019, 049





[Stewart, & Gong(2001)]2001PhLB..510....1S Stewart, E. D., & Gong, J.-O. 2001, Physics Letters B, 510, 1




[Tomikawa & Kobayashi(2020)]2020PhRvD.101h3529T Tomikawa, K. & Kobayashi, T. 2020, , 101, 083529





[Vachaspati, & Vilenkin(1985)]1985PhRvD..31.3052V Vachaspati, T., & Vilenkin, A. 1985, , 31, 3052





[Witten(1984)]1984PhRvD..30..272W Witten, E. 1984, , 30, 272




  
[Yoo, & Gong(2016)]2016JCAP...07..017Y Yoo, J., & Gong, J.-O. 2016, , 2016, 017




[Yuan et al.(2020)]2020PhRvD.101f3018Y Yuan, C., Chen, Z.-C., & Huang, Q.-G. 2020, , 101, 063018






























