title: Analytic solutions of scalar perturbations induced by scalar perturbations
authors: Keisuke Inomata
text: 
  
  
  RESCEU-19/20
  
  
  Vibronic response of a spin-1/2 state from a carbon impurity in two-dimensional WS_2
      Bruno Schuler
      December 30, 2023
  ====================================================================================
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  
  
  The inflation theory predicts that the cosmological perturbations are generated from the quantum fluctuations of the fields, which means that valuable information about the early Universe is printed on the perturbations.
  To reveal the history of the Universe, many authors have studied the perturbations for several decades.
  Based on the behavior under the spatial coordinate transformation, the cosmological perturbations can be divided into three types: scalar, vector, and tensor perturbations.
  In particular, scalar perturbations have been playing important roles in the study of the Universe.
  For example, on large scales (≳ 1Mpc), the amplitude and tilt of power spectrum of scalar perturbations have been determined through the observations of the large scale structure and the anisotropies of cosmic microwave background (CMB) <cit.>.
  These observational results put constraints on parameters of inflation models <cit.>.
  Also, scalar perturbations on small scales (≲ 1Mpc) have been attracting a lot of interest because of their rich phenomenology. 
  If the amplitude of the small-scale perturbations are large enough, some unique compact objects, primordial black holes (PBHs) <cit.> and ultracompact minihalos (UCMHs) <cit.>, are produced.
  In particular, PBHs are one of the hottest topics in cosmology because they are good candidates of dark matter and the heavy BHs detected by LIGO-Virgo collaborations <cit.> (see also Refs. <cit.> for recent reviews).
  
  
  
  In the context of the small-scale scalar perturbations, many authors have recently been interested in the gravitational waves (GWs) induced by scalar perturbations because the induced GWs could be observed by the future observations depending on the amplitude of the small-scale perturbations <cit.>.
  Here, we draw attention to the fact that scalar perturbations also induce scalar perturbations on different scales through the non-linear interaction between scalar perturbations.[Scalar perturbations induced by tensor perturbations are discussed from the viewpoint of PBH production in Refs. <cit.>.]
  Since the non-linear interaction appears at the second or higher order in perturbations, the induced scalar perturbations are smaller than the source (or first-order) scalar perturbations. 
  However, the production of scalar perturbations on different scales can be important when we discuss the sharp power spectrum of the curvature perturbations, which are often considered in PBH scenarios <cit.>, because the induced scalar perturbations can be dominant on scales where the source scalar perturbations are very small.
  As concrete examples, we will see this important case in Fig. <ref> (Sec. <ref>) and Fig. <ref> (Sec. <ref>) with the delta-function power spectrum of the source perturbations and in Fig. <ref> (Appendix <ref>) with the power spectrum for PBH scenarios.
  
  
  
  Motivated by the above discussion, we shed light on the scalar perturbations induced by scalar perturbations in this paper.
  Phenomenologically, the induced scalar perturbations have been studied in the context of the CMB anisotropies <cit.> 
  and the large scale structures <cit.>.
  Also, there are works focusing on the formalisms of the second-order perturbations, including scalar, vector, and tensor perturbations <cit.>.
  The main aim of this paper is to derive the analytic solutions of the induced scalar perturbations in two eras, a radiation-dominated (RD) era and a matter-dominated (MD) era.
  Although Ref. <cit.> derived the analytic solutions with simplified evolution of perturbations, we take into account the contributions that were neglected in the reference for the first time. 
  On top of the reference, Refs. <cit.> derived the analytic solutions in a MD era with the Newtonian approximation, valid in the subhorizon limit <cit.>. 
  By contrast, in this paper, we derive the solutions that can describe the induced perturbations even on superhorizon scales.[Ref. <cit.> showed that if the C-gauge and the comoving condition are taken, the equations in the Newtonian approximation can be applied to the perturbations even on superhorizon scales. In this paper,  we derive the analytic solutions of the induced perturbations in the conformal Newtonian gauge to make a comprehensive analysis in both the eras. ]
  
  
  This paper is organized as follows.
  In Sec. <ref>, we summarize the basic equations for the scalar perturbations.
  In Sec. <ref>, we derive the analytic solutions of the second-order perturbations and discuss their evolution in a RD.
  In Sec. <ref>, we discuss the second-order perturbations in a MD era, similarly to Sec. <ref>.
  Finally, we devote Sec. <ref> to the conclusion.
  
  
  
  
  
  
  
  § BASIC EQUATIONS FOR SCALAR PERTURBATIONS
  
  
  
  
  In this section, we summarize the basic equations for first- and second-order scalar perturbations.
  Throughout this paper, we assume a flat FRW universe, follow the notation in Ref. <cit.>, and take the conformal Newtonian (longitudinal) gauge. 
  Then, we can express the metric perturbations up to the second order as 
  
  
      s^2    = g_μν x^μ x^ν
          = a^2 {- (1 + 2 Φ^(1) + Φ^(2))η^2  + V_i^η x^i + [ (1-2 Ψ^(1) - Ψ^(2)) δ_ij + 1/2 h^_ij]  x^i  x^j },
  
  where the superscript denotes the order in perturbations and we have neglected the first-order vector and tensor perturbations, V^_i and h^_ij, because the main focus of this paper is on scalar perturbations induced by the first-order scalar perturbations.
  In this paper, we assume a perfect fluid for simplicity. 
  Then, we can express the energy-momentum tensor as
  
  
      T^μ_ ν   =  (ρ + P) u^μ u_ν + P δ^μ_ ν,
  
  
  where ρ is the energy density, P is the pressure, and u_μ is the 4-velocity.
  We take the following notation for their perturbations up to the second order:
  
  
      ρ   = ρ^ + δρ^ + 1/2δρ^, 
      
      	P    = P^ + δ P^ + 1/2δ P^, 
      
      	u^i    = 1/a( δ v^,i + 1/2δ v^,i + 1/2δv^i_V^),
  
  
  where the superscript “” means the background value, δv^i_V^ is the second-order vector perturbation of the velocity perturbation, and the comma denotes the spatial derivative as A^,i≡δ^ij A_,j≡δ^ij ∂ A/∂ x^j.
  Note that the transformation between the superscript and the subscript in the 3-dimensional space is performed by δ_ij because we consider a flat universe.
  Besides, we can express u^0 and u_μ (≡ g_μν u^ν) from their definitions as
  
  
      u^0 
      	   ≡ x^0/√(-g_μν x^μ x^ν)
         = 1/a( 1 - Φ^ - 1/2Φ^ + 3/2(Φ^)^2 + 1/2δ v^,iδ v^_   ,i), 
      
      	u_0    ≡ g_0 μ  u^μ
         = -a ( 1 + Φ^ + 1/2Φ^ - 1/2(Φ^)^2 + 1/2δ v^,iδ v^_   ,i), 
      
      	u_i    ≡ g_iμ  u^μ
         = a ( δ v^_   ,i + 1/2δ v^_   ,i + 1/2δv_V i^  -2 Ψ^δ v^_   ,i + 1/2V_i^),
  
  
  where we have neglected the terms of the third order or higher in perturbations.
  
  
  
  In what follows, we consider the adiabatic perturbations for simplicity.
  Then, the perturbations in a perfect fluid satisfy the following equations: (see Appendix <ref> for derivation)
  
  
      Φ^” + 3 (1 + c_^2) ℋΦ^' + [2ℋ' + (3c_^2 + 1) ℋ^2 ] Φ^ - c_^2 Φ^,i_    ,i = 0, 
         δ^≡δρ^/ρ^ = - ( 2 Φ^ + 2/ℋΦ^' - 2/3ℋ^2Φ^  ,i_     ,i),
  
  
  where the prime denotes the derivative with respect to η, ℋ (≡ a'/a) is the conformal Hubble parameter, and c_^2 (≡δ P^/δρ^) is the square of the adiabatic sound speed. 
  Note that we assume c_ is constant for simplicity throughout this paper.
  Solving these equations, we obtain the following solutions of the Fourier modes of the perturbations in a RD and a MD era:
  
  
      Φ^ ( k, η)    =  -2/3ζ^( k) T_Φ,(x)     (RD era) 
      
       -3/5ζ^( k) T_Φ,(x)     (MD era)
      	 ,	 
      δ^ ( k, η)    =  -2/3ζ^( k) T_δ,(x)     (RD era) 
      
       -3/5ζ^( k) T_δ,(x)     (MD era)	 	 
      	 ,
  
  
  where x≡ k η with k ≡ | k| and the first-order curvature perturbation ζ^ and the transfer functions are defined as
  
  
      ζ^   ≡ -Ψ^ + δρ^/3(ρ^ + P^), 
      
      	T_Φ,(x)    ≡ 3√(3)j_1(x/√(3))/x
          = 
      	3 [sin(x/√(3)) - (x/√(3)) cos (x/√(3))]/(x/√(3))^3, 
      
      	T_Φ,(x)    ≡ 1, 
      
      	T_δ,(x)    ≡6x(-6 + x^2) cos(x/√(3)) - 12 √(3) (-3 + x^2) sin(x/√(3))/x^3,
      
      	T_δ,(x)    ≡ -2 - x^2/6.
  
  
  Note that j_1 is the spherical Bessel function of the first kind.
  
  
  For the second-order perturbations, we can express the equation of motion for Ψ^ in a perfect fluid as <cit.>[
  Note that the notation of Ref. <cit.> is different from that in this paper. They are related to each other as, Φ^_ref = Φ^_this - 2 (Φ^_this)^2 and Ψ^_ref = Ψ^_this + 2 (Ψ^_this)^2, where the subscripts, “ref” and “this”, denote the variables in the reference and this paper, respectively.] (see also Appendix <ref> for derivation)
  
  
      Ψ^” + 3 (1 + c_^2) ℋΨ^' + [2ℋ' + (3c_^2 + 1) ℋ^2 ] Ψ^ - c_^2 Ψ^,i_    ,i = S^,
  
  
  where we have used the relation δ P^/ δρ^ = c_s^2 (≡δ P^/ δρ^), which are valid in the case the perturbations are adiabatic and the sound speed is constant <cit.>.
  The source term S^ is defined as 
  
  
      S^≡   (3c_^2 -1/3) Φ^ ,iΦ^_   ,i + 8c_^2 Φ^Φ^ ,i_   ,i + (3c_^2 + 1) (Φ^' )^2 + [ (3c_^2 + 1)ℋ^2 + 2 ℋ' ]  N^j_ iB^i_ j^
          
      	 + ℋ  N^j_ i(B^i_ j^)' + 1/3  N^j_ i(B^i_ j^)^,k_ ,k + (1/3 - c_^2 ) 4/3(1+w) ℋ^2( ℋΦ^,i +  Φ^,i' )  (  ℋΦ^_,i + Φ^'_, i),
  
  
  where w (≡ P^/ρ^) is the equation-of-state parameter and we assume w is constant, similarly to c_.
  The B^i_ j^ and N^j_ i are defined as 
  
  
      B^i_ j^≡[4(5+3w)/3(1+w)Φ^ ,iΦ^_   ,j 
      + 8/3(1+w) ℋ(Φ^ ,iΦ^_   ,j)' 
       + 8/3(1+w) ℋ^2Φ^ ,i' Φ^_   ,j'],
  
  
  
  
      N^j_ i A^i_ j( x)    ≡3/2∇^-2( ∂^j ∂_i/∇^2 - 1/3δ^j_ i) A^i_ j ( x) 
         = ∫^3 k/(2π)^3( -3/2k^2) ( k^j k_i/k^2 - 1/3δ^j_ i) A^i_ j ( k),
  
  
  where A^k_ j is an arbitrary tensor.
  
  Besides, we can express the second-order energy density perturbation with Ψ^ and Φ^ as 
  
  
      δ^≡   δρ^/ρ^
      
      	=    -2 Ψ^ +2 N^j_ iB^i_ j^ -2/ℋΨ^' + 2/3ℋ^2Ψ^,i_    ,i 
      	+ 2/ℋ^2(Φ^' )^2 + 16/3ℋ^2Φ^Φ^ ,i_   ,i
          + 1/ℋ^2( 2 - 8/9(1+w)) Φ^ ,iΦ^_   ,i - 8/9(1+w) ℋ^3 (Φ^ ,iΦ^_   ,i )'  - 8/9(1+w) ℋ^4Φ^'_, iΦ^,i',
  
  
  where its derivation is given in Appendix <ref>.
  In the superhorizon limit (kη≪ 1), the perturbation becomes 
  
  
      δ^≃ -2 Ψ^ +2 N^j_ i(4(5+3w)/3(1+w)Φ^ ,iΦ^_   ,j),
  
  
  where we have used Φ^' = Ψ^' =0 in the superhorizon limit. 
  
  As a gauge-invariant variable, we introduce the second-order curvature perturbation in the uniform density gauge,[Strictly speaking, we need one more gauge condition on top of the uniform density one to define the second-order curvature perturbations uniquely. In this paper, we take the same gauge condition as that in Eq. (7.71) of Ref. <cit.> (see also Appendix <ref>).]  which are expressed in the Newtonian gauge as <cit.>[We have fixed some typos in Eq. (7.71) of Ref. <cit.>.]
  
  
      ζ^ 
      	=   
      	-Ψ^ + 1/3(1+w)δ^ + 1/9 ℋ (1+w)^2[(  δ^)^2]' - 1+3w/9 (1+w)^2( δ^)^2  
          
      	- 2/3 ℋ (1+w)δ^( Ψ^' + 2 ℋΨ^) - 1/18 ℋ^2 (1+w)^2[ δ^  ,kδ^_,k - ∇^-2(  δ^_,i δ^_,j)^,ij],
  
  
  where the expression in a general gauge is given in Appendix <ref>.
  In the superhorizon limit, the curvature perturbation in a perfect fluid with Φ^ = Ψ^ can be expressed as
  
  
      ζ^≃ - 5+3w/3(1+w)Ψ^ + 2/3(1+w) N^j_ i(4(5+3w)/3(1+w)Φ^ ,iΦ^_   ,j) - 1+3w/9 (1+w)^2( δ^)^2 - 4/3(1+w)δ^Φ^,
  
  
  where we have used Eq. (<ref>). 
  The initial condition of the second-order curvature perturbation depends on the primordial non-Gaussianity.
  Following Ref. <cit.>, we assume the local-type non-Gaussianity throughout this paper, which is parametrized as ζ^ = 2 a_NL (ζ^)^2 in the superhorizon limit.
  The value of a_NL depends on the scenario and a_NL=1 corresponds to the Gaussian perturbation <cit.>.
  With the parametrization of the second-order curvature perturbation, Ψ^ in the superhorizon limit can be expressed as 
  
  
      Ψ^≃ -6(1+w)/5+3w a_NL (ζ^)^2 + 8/3(1+w) N^j_ i(Φ^ ,iΦ^_   ,j) - 1+3w/3 (1+w) (5+3w)( δ^)^2 - 4/5+3wδ^Φ^.
  
  
  
  
  
  
  § INDUCED SCALAR PERTURBATIONS IN A RADIATION-DOMINATED ERA
  
  
  
  
  
  In this section, we discuss the induced scalar perturbations in a RD era using the basic equations in the previous section.
  
  
  
  
   §.§ Analytic solutions
  
  
  First, we derive the analytic solutions for Ψ^, δ^, and ζ^.
  Taking c_^2 = w = 1/3, we can express Eq. (<ref>) in Fourier space as
  
  
      Ψ^”( k,η) + 4/ηΨ^'( k,η) + k^2/3Ψ^ ( k,η) = S^_ ( k,η),
  
  
  where we have substituted ℋ = 1/η.
  The source term S^_ is given by
  
  
      S^_ ( k,η) =   ∫^3 k̃/(2π)^3{ -( 2/3k̃· ( k - k̃) + 4/3 (k̃^2 + | k - k̃|^2 ) ) Φ(k̃) Φ( k - k̃) + 2 Φ'(k̃) Φ'( k - k̃)  . 
                   
      	+ 3/2k^2( 1/k^2 (  k ·k̃ ) (  k · ( k - k̃ ) ) - 1/3k̃· ( k - k̃) ) 
                   ×[
      	1/η( 6 Φ(k̃) Φ( k - k̃) + 4 ηΦ'(k̃) Φ( k - k̃) + 2 η^2 Φ'(k̃) Φ'( k - k̃) )' . 
                 . . ( 2/3)
      	- k^2/3( 6 Φ(k̃) Φ( k - k̃) + 4 ηΦ'(k̃) Φ( k - k̃) + 2 η^2 Φ'(k̃) Φ'( k - k̃) ) ] },
  
  
  where we have omitted the superscript for the first-order perturbations. 
  Using Eqs. (<ref>) and (<ref>), we can rewrite the source term as 
  
  
      S^_ ( k,η)
      	=   ∫^3 k̃/(2π)^3{ -( 2/3k̃· ( k - k̃) + 4/3 (k̃^2 + | k - k̃|^2) ) T_Φ,(k̃η) T_Φ,(| k - k̃| η)  . 
                     
      	+ 2 T_Φ,'(k̃η) T_Φ,'(| k - k̃| η) - 3/2k^2( 1/k^2 (  k ·k̃ ) (  k · ( k - k̃ ) ) - 1/3k̃· ( k - k̃) ) 
               . ( 2/3)×[ 
      	( 2k^2 + 2/3(k̃^2 + | k - k̃|^2) ) T_Φ,(k̃η) T_Φ,(| k - k̃| η) . . 
                         
      	+ ( 2k^2/3 +  2| k - k̃|^2/3) η T_Φ,'(k̃η) T_Φ,(| k - k̃| η)  
                         
      	 + ( 2k^2/3 + 2k̃^2/3) η T_Φ,(k̃η) T_Φ,'(| k - k̃| η) 
                 .. ( k̃^2/3) 
      	 + ( 2 k^2η^2/3 +8 ) T_Φ,'(k̃η) T_Φ,'(| k - k̃| η) ] }( 2/3)^2 ζ(k̃) ζ( k - k̃) 
      
      	 =   ∫^3 k̃/(2π)^3 k^2 uv f_(u,v,x) ( 2/3)^2 ζ(k̃) ζ( k - k̃),
  
  
  where u ≡ | k - k̃|/k, v ≡k̃/k, and f_(u,v,x) is defined as 
  
  
      f_(u,v,x) ≡   [ - ( 3u^2 + 3 v^2 +1/3 u v) - ( 2(u^2+v^2) - 3(u^2 - v^2)^2 +1/8uv) ( 2 + 2/3( u^2 + v^2 ) ) ] 
      	T_Φ,(v x) T_Φ,(u x) 
           
      	- ( 2(u^2+v^2) - 3(u^2 - v^2)^2 +1/8uv) ( 2/3 + 2u^2/3) v x  T_Φ,(v x)/ (vx) T_Φ,( u x) 
           
      	- ( 2(u^2+v^2) - 3(u^2 - v^2)^2 +1/8uv) ( 2/3 + 2v^2/3) u x  T_Φ,( v x)  T_Φ,(ux)/ (ux)
           
      	+ [ 2/uv - (  2(u^2+v^2) - 3(u^2 - v^2)^2 +1/8uv) ( 2x^2/3 + 8 ) ] u v   T_Φ,(vx)/ (vx) T_Φ,(ux)/ (ux).
  
  
  
  We solve Eq. (<ref>) using the Green function method.
  Here, we define z ≡ a^2 Ψ^ and rewrite Eq. (<ref>) as
  
  
      z”( k,η) + ( k^2/3 - 2/η^2) z( k,η) = a^2 S^_( k,η).
  
  
  The Green function is defined as the solution of the following equation:
  
  
      G”_ + ( k^2/3 - 2/η^2) G_ = δ(η - η̅),
  
  
  where the prime denotes the derivative with respect to η, not η̅.
  Then, we obtain the Green function as 
  
  
      k G_ (k, η;η̅) = -Θ(η - η̅) x x̅/√(3)[ j_1(x/√(3)) y_1(x̅/√(3)) - j_1(x̅/√(3)) y_1(x/√(3)) ],
  
  
  where y_1 is the spherical Bessel function of the second kind.
  Using this Green function, we can solve Eq. (<ref>) as 
  
  
      Ψ^ ( k,η) =   Ψ^( k, 0) T_Φ,(x) + ∫^η_0 η̅(a(η̅)/a(η))^2 G_(k,η;η̅) S^_ ( k,η) 
      
      	=   Ψ^( k, 0) T_Φ,(x) + ∫^3 k̃/(2π)^3 uv I_Ψ,,(u,v,x) ( 2/3)^2 ζ(k̃) ζ( k - k̃),
  	
  
  where I_Ψ,, is defined as 
  
  
      I_Ψ,,(u,v,x) ≡∫^x_0 x̅(x̅/x)^2 G_(k,η;η̅) f_(u,v,x).
  
  
  We can rewrite I_Ψ,, as 
  
  
      I_Ψ,,(u,v,x) 
      	=    ∫^x_0 x̅ (x̅/x)^2 { -x x̅/√(3)[ j_1(x/√(3)) y_1(x̅/√(3)) - j_1(x̅/√(3)) y_1(x/√(3)) ] } f(u,v, x̅) 
      
      	=     𝒥(u,v,x) 3√(3)j_1(x/√(3))/x + 𝒴(u,v,x) 3√(3)y_1(x/√(3))/x,
  
  
  where 𝒥 and 𝒴 are defined as 
  
  
      𝒥(u,v,x) ≡   ∫^x_0 x̅ ℐ_j(u,v,x), 
      𝒴(u,v,x) ≡   ∫^x_0 x̅ ℐ_y(u,v,x), 
      ℐ_j(u,v,x̅) ≡   ( -x̅^3/9) y_1(x̅/√(3)) f_(u,v,x̅), 
      ℐ_y(u,v,x̅) ≡   x̅^3/9  j_1(x̅/√(3)) f_(u,v,x̅).
  
  
  Here, we expand the integrands ℐ_j and ℐ_y as 
  
  
      ℐ_a(u,v, x̅) =  ∑_m=0^7 ∑_n=1^8 sin(α_n x̅  + ϕ_n) M^a_nm/x̅^m,
  
  
  where the sub/superscript a denotes j or y and α_n, ϕ_n, and M^a_nm are independent of x̅.
  Specifically, α_n and ϕ_n are given as 
  
  
      α_n    = 1 - u - v/√(3)      n=1,5 
      1 + u - v/√(3)      n=2,6 
      1 - u + v/√(3)      n=3,7 
      1 + u + v/√(3)      n=4,8 
      		, 
      ϕ_n    = π/2       n ≤ 4 
      	
      		0       n > 4
      		.
  
  
  The explicit expressions of M^a_nm are given in Appendix <ref>.
  For convenience, we divide 𝒥 into two parts depending on the power of x̅ of its integrand as
  
  
      𝒥(u,v,x) =   𝒥_0(u,v,x) + 𝒥_h(u,v,x), 
      𝒥_0(u,v,x) ≡   ∫^x_0 x̅ ∑_n=1^8 sin(α_n x̅  + ϕ_n) M^j_n0, 
      𝒥_h(u,v,x) ≡   ∫^x_0 x̅ ∑_m=1^7 ∑_n=1^8 sin(α_n x̅  + ϕ_n) M^j_nm/x̅^m,
  
  
  and similarly for 𝒴. 
  Note that the pioneering work, Ref. <cit.>, focused only on the contributions from 𝒥_0 and 𝒴_0.
  Performing the integrals of 𝒥_0 and 𝒴_0, we obtain the following expressions:
  
  
      𝒥_0 (u,v,x) =   3(-1 + 3(u^2-v^2)^2 -2(u^2 + v^2))/16 u^2 v^2
           ×( cos( -1+u-v/√(3) x )/-1+u-v - cos( 1+u-v/√(3) x )/1+u-v - cos( -1+u+v/√(3) x )/-1+u+v + cos( 1+u+v/√(3) x )/1+u+v)  
         - 3(-1 + 3(u^2-v^2)^2 -2(u^2 + v^2))/2uv ((u+v)^2-1)((u-v)^2 -1) , 
      𝒴_0(u,v,x) =    -3(-1 + 3(u^2-v^2)^2 -2(u^2 + v^2))/16 u^2 v^2
             ×( sin( -1+u-v/√(3) x )/-1+u-v + sin( 1+u-v/√(3) x )/1+u-v - sin( -1+u+v/√(3) x )/-1+u+v - sin( 1+u+v/√(3) x )/1+u+v).
  
  
  Note that 𝒥_0 and 𝒴_0 do not converge even in the late-time limit (x ≫ 1).
  On the other hand, 𝒥_h and 𝒴_h converge to 𝒥_h,late and 𝒴_h,late in the late-time limit, which are given as 
  
  
      𝒥_h,late(u,v) =   1/8 u^3 v^3[-9(u^6 + u^4 -u^2) -9(v^6 + v^4 - v^2) + 9  . 
                 . + 6 u^2 v^2 (u^2 - v^2)^2 + 5 u^2 v^2(u^2 + v^2) + 8 u^2 v^2 ] 
           -3/32 u^4 v^4( 3(u^4-v^4)^2 - 6(u^4 +v^4) + 3 - 8u^2v^2(1 + u^2 + v^2) ) log( 1-(u-v)^2 /(u+v)^2 -1), 
      𝒴_h,late(u,v) =    3 π/32 u^4 v^4( 3u^4(u^4 - 2) + 3v^4(v^4 -2) + 3 - 8u^2v^2(1 + u^2 + v^2) -6 u^4 v^4 ).
  
  
  The full expressions of 𝒥_h and 𝒴_h are given in Appendix <ref>.
  From Eq. (<ref>), we can express the contribution from the initial second-order perturbations as 
  
  
      Ψ^ ( k, 0) T_Φ,(x) = ∫^3 k̃/(2π)^3 uv I_Ψ,,i(u,v,x) ( 2/3)^2 ζ(k̃) ζ( k - k̃),
  
  
  where I_Ψ,,i is defined as
  
  
      I_Ψ,,i(u,v,x)    ≡𝒥_i(u,v)T_Φ,(x), 
      𝒥_i(u,v)    ≡2(u^2+v^2) - 3(u^2 - v^2)^2 +5 -12 a_NL/4uv.
  
  
  Combining all the contributions, we finally obtain the following expression:
  
  
      Ψ^ ( k,η) =   ∫^3 k̃/(2π)^3 uv I_Ψ,(u,v,x) ( 2/3)^2 ζ(k̃) ζ( k - k̃),
  
  
  where I_Ψ,≡ I_Ψ,,i + I_Ψ,,.
  This is the analytic solution of Ψ^.
  Here, we mention the asymptotic behavior in the large-scale limit, u ≫ 1 and v ≫ 1[
  The word “large-scale limit” here indicates the scales much larger than those of the source perturbations. On the other hand, the word “late-time limit” indicates the case where the induced perturbations (not the source perturbation) are on deep subhorizon scales. ], in which 𝒥_i can be approximated as
  
  
      𝒥_i(u,v) ≃2(u^2+v^2) - 3(u^2 - v^2)^2/4uv∼𝒪(1)      (u ≫ 1, v ≫ 1).
  
  
  When the induced perturbations are on superhorizon scales (x ≪ 1), |𝒥_h| and |𝒥_0| are much smaller than |𝒥_i| and therefore the leading contribution comes only from 𝒥_i. 
  On the other hand, when the induced perturbations are on deep subhorizon scales (x ≫ 1), 𝒥_h asymptotes to 𝒥_h,late, given as 
  
  
      𝒥_h,late(u,v) ≃-9(u^6+v^6) + 5 u^2 v^2 (u^2 + v^2)/8u^3 v^3 + 3(u^2 -v^2)^2/4uv     (u ≫ 1, v ≫ 1).
  
  
  Given |u-v| < 1 by definition, we can easily see that the cancellation between 𝒥_i and 𝒥_h,late occurs as 𝒥_i + 𝒥_h,late∼𝒪(1/u) or 𝒪(1/v).[In the large-scale limit (u≫ 1 and v ≫ 1), 𝒥_0 ∼𝒪(1/u^2) or 𝒪(1/v^2) regardless of x and therefore 𝒥_0 is subdominant in the limit. ]
  In this sense, the contribution from 𝒥_h, which was neglected in the previous work <cit.>, is very important for the late-time (x ≫ 1) evolution of the perturbations especially on the scales much larger than the source perturbations (u ≫ 1 and v ≫ 1).
  
  Similarly to Ψ^, from Eq. (<ref>), we obtain the analytic solution of δ^ as 
  
  
      δ^ ( k,η) =   ∫^3 k̃/(2π)^3 uv I_δ,(u,v,x) ( 2/3)^2 ζ(k̃) ζ( k - k̃),
  
  
  where 
  
  
      I_δ,(u,v,x) ≡    -2 (1+ x^2/3) I_Ψ,(u,v,x) - 2 x  I_Ψ,(u,v,x)/ x + I_δ,,(u,v,x), 
      
      	I_δ,,(u,v,x) ≡   9 -4 x^2 + 6(u^2 + v^2)(3 - 2x^2) - 27 (u^2 - v^2)^2 /6 uv T_Φ,(ux) T_Φ,(vx) 
         + 3 + 2x^2 + 2(3 - x^2) (u^2 + v^2) -9 (u^2 - v^2)^2/6 uv
                 ×( ux  T_Φ,(ux)/(ux) T_Φ,(vx) + vx T_Φ,(ux)  T_Φ,(vx)/ (vx)) 
          + x^2 [15 + 2x^2 + 2(3 -x^2)(u^2 + v^2) -9 (u^2 - v^2)^2 ]/6 T_Φ,(ux)/ (ux) T_Φ,(vx)/ (vx).
  
  
  
  
  For the second-order curvature perturbation, Eq. (<ref>) leads to 
  
  
      ζ^ ( k,η) =   ∫^3 k̃/(2π)^3 uv I_ζ,(u,v,x) ( 2/3)^2 ζ(k̃) ζ( k - k̃),
  
  
  where 
  
  
      I_ζ,(u,v,x) ≡    - I_Ψ, (u,v,x) + 1/4 I_δ,(u,v,x)  + I_ζ,,(u,v,x), 
      
      	I_ζ,,(u,v,x) ≡    -1/uv{1/8 T_δ,(vx) T_δ,(ux) - 1/16(vx  T_δ,(v x)/ (v x) T_δ,(ux) + ux  T_δ,(u x)/ (u x)T_δ,(vx)   ) . 
                 
      	+ 1/4[ T_δ,(vx) ( ux  T_Φ,(ux)/ (ux) + 2 T_Φ,(ux) ) + T_δ,(ux) ( vx  T_Φ,(vx)/ (vx) + 2 T_Φ,(vx) ) ]  
                 . - x^2/32[ (u^2 - v^2)^2 - 2(u^2 + v^2) + 1/4]  T_δ,(vx) T_δ,(ux) }.
  
  
  Note that I_Ψ,, I_δ,, and I_ζ, are symmetric under the exchange u ↔ v.[
  Strictly speaking, we can define I_Ψ,, I_δ,, and I_ζ, so that they are not symmetric under the exchange u ↔ v. 
  However, after the integration over k̃, only the symmetric part contributes to the second-order perturbations.
  ]
  
  
  Figure <ref> shows the evolution of the transfer functions, I_Ψ,, I_δ,, and I_ζ,, with different a_NL, u, and v.
  The lines for u=v=1 show the contributions from the perturbations whose wavenumbers are of the same order as those of the induced scalar perturbations.
  On the other hand, the lines for u=v=10 show the contributions from the perturbations whose wavenumbers are larger than those of the induced scalar perturbations.
  The time x ∼ 1/u and 1/v corresponds to the horizon entry of the source (or first-order) perturbations and x ∼ 1 corresponds to that of the induced (or second-order) perturbations.
  For a_NL, we take 0 and 1 as fiducial examples, which correspond to the perturbations with and without the non-Gaussianities, respectively.
  
  From this figure, we can see that, before the source perturbations enter the horizon (x ≪ 1/u and 1/v), the induced perturbations remain constant except for ζ^ with a_NL = 0, which grows even before that.
  Note that this behavior of ζ^ for a_NL = 0 does not contradict Ref. <cit.> because the reference neglects the spatial derivative of the first-order perturbations.
  This behavior just shows the fact that, even if we take the initial condition ζ^=0 in the superhorizon limit, ζ^ can evolve through the terms suppressed by the gradient factors such as u k η and v k η (≪ 1).
  
  Ψ^ and δ^ start to decay around x ∼ 1/u and 1/v and oscillate around x ∼ 1.
  On the other hand, ζ^ starts to oscillate around x ∼ 1/u and 1/v.
  In the late-time (x ≫ 1), all the second-order perturbations in Fig. <ref> oscillate and their oscillation amplitudes are related to the scale factor as I_Ψ,∝ a^-2, I_δ,∝ a^0, and I_ζ,∝ a^2.
  For I_Ψ, and I_δ,, their scale factor dependences in the late-time limit are the same as those of the first-order perturbations, Ψ^∝ a^-2 and δ^∝ a^0 up to their oscillations.
  On the other hand, I_ζ, evolves differently from the first-order curvature perturbation in the late-time limit (c.f. ζ^∝ a^0 up to its oscillation).
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Auto-power spectra
  
  
  
  In this subsection, we calculate the auto-power spectra of the second-order induced perturbations to see the evolution of the second-order perturbations.
  Note that, to compare with observations, we also need to take into account the cross-power spectra between the first-order and the third-order perturbations, which is of the same order of magnitude as the auto-power spectra.
  Since the main focus of this paper is the analytic solution of the second-order perturbations, we leave the study of the cross-power spectra for future work.
  For this reason, we should keep in mind that the observable power spectra might be different from the auto-power spectra.
  The main purpose of the auto-power spectrum in this paper is to show the evolution of the second-order perturbations.
  
  
  The expectation values of the second-order perturbations can be written as 
  
  
      ℳ^ ( k_1,η) ℳ^( k_2,η)  = ∫^3 k̃_1/(2π)^3∫^3 k̃_2/(2π)^3   (u_1 v_1) (u_2 v_2)  I_ℳ,(u_1,v_1,x_1) I_ℳ,(u_2,v_2,x_2) 
         ×( 2/3)^4 ζ(k̃_1) ζ( k_1 - k̃_1) ζ(k̃_2) ζ( k_2 - k̃_2) ,
  
  
  where ℳ indicates Ψ, δ, or ζ. The arguments are defined as v_1 ≡k̃_1/k_1, u_1 ≡ | k_1 - k̃_1|/k_1, x_1 ≡ k_1η, and similarly for the arguments with the subscript 2.
  Here, we assume that the non-Gaussianity of the primordial curvature perturbation is negligibly small for simplicity. 
  Then, we obtain
  
  
      ζ (k̃_1)  ζ (  k_1- k̃_1) ζ (k̃_2) ζ (  k_2- k̃_2)   = ζ (k̃_1)  ζ (k̃_2)ζ (  k_1- k̃_1) ζ (  k_2 - k̃_2)
            + ζ (k̃_1)  ζ (  k_2 - k̃_2)ζ (  k_1- k̃_1 ) ζ (k̃_2)
         = (2π)^3δ(k̃_1 + k̃_2) (2π)^3δ( k_1 - k̃_1 +  k_2 - k̃_2) 
                 ×2π^2/k̃^3_1𝒫_ζ^(k̃_1) 2π^2/| k_1 - k̃_1|^3𝒫_ζ^(| k_1 - k̃_1|)
            + (2π)^3δ(k̃_1 +   k_2 - k̃_2) (2π)^3δ( k_1 - k̃_1 + k̃_2) 
                 ×2π^2/k̃^3_1𝒫_ζ^(k̃_1) 2π^2/| k_1 - k̃_1|^3𝒫_ζ^(| k_1 - k̃_1|),
  
  
  where we have restored the superscript of the curvature perturbation in the power spectrum, given as
  
  
      ζ ( k_1) ζ ( k_2) = (2π)^3 δ( k_1 +  k_2) 2 π^2/k^3_1𝒫_ζ^(k_1).
  
  
  Then, we can rewrite Eq. (<ref>) as 
  
  
      ℳ^ ( k_1,η) ℳ^( k_2,η)  = (2π)^3 δ( k_1 +  k_2) 2 π^2/k^3_1𝒫_ℳ^(k_1,η), 
         𝒫_ℳ^(k,η) ≡∫_0^∞ v  ∫_|v-1|^v+1 u   I_ℳ,^2(u,v,x) ( 2/3)^4 𝒫_ζ^( k v) 𝒫_ζ^(k u).
  
  
  
  
  As a simple example, we here consider the delta-function power spectrum of curvature perturbation
  
  
      𝒫_ζ^(k) = A_ζδ(log(k/k_*)).
  
  
  In this case, we can express Eq. (<ref>) as 
  
  
      𝒫_ℳ^(k,η) =  A_ζ^2 ( 2/3)^4 ( k_*/k)^2 I^2_ℳ,(k_*/k,k_*/k,x) Θ(2- k/k_*).
  
  
  Figure <ref> shows the evolution of Eq. (<ref>) for Ψ^, δ^, and ζ^.
  We assume the Gaussian distribution of the curvature perturbations taking a_NL = 1.
  The power spectra of the three perturbations remain constant when the source perturbations at k_* are on superhorizon scale (x_* (≡ k_* η) ≲ 1).
  Once the source perturbations enter the horizon, the power spectra start to evolve even on superhorizon scales at that time through the non-linear interaction between the perturbations with different wavenumbers.
  We can also interpret the evolution on superhorizon scale as the cancellation between 𝒥_i and 𝒥_h, mentioned below Eq. (<ref>).
  See also Appendix <ref> for the auto-power spectra in a somewhat more realistic situation, which is often considered in the PBH scenarios.
  
  
  
  
  
  
  
  
  
  
  § INDUCED SCALAR PERTURBATIONS IN A MATTER-DOMINATED ERA
  
  
  
  
  In this section, we discuss the induced scalar perturbations in a MD era (c_s^2 = w = 0).
  
  
  
  
   §.§ Analytic solutions
  
  
  
  
  First, we derive the analytic solutions.
  The equations in a MD era are somewhat more simple than those in a RD era because Φ^' = 0 during the era.
  From Eq. (<ref>), we have
  
  
      Ψ^”( k,η) + 6/ηΨ^'( k,η) = S^_m( k,η),
  
  
  where we have substituted ℋ = 2/η and the source term is given as 
  
  
      S^_ ( k,η) 
      	=   ∫^3 k̃/(2π)^3[ k̃· ( k - k̃) 
      	- 10/3k^2   (  k ·k̃ ) (  k · ( k - k̃ ) ) ] Φ(k̃) Φ( k - k̃) 
      
      	=   ∫^3 k̃/(2π)^3 k^2 uv f_(u,v) ( 3/5)^2 ζ(k̃) ζ( k - k̃).
  
  
  The source function f_ is defined as 
  
  
      f_(u,v) ≡    - 2 + 3 (u^2 +v^2) - 5 (u^2 -v^2)^2/6  u v.
  
  
  We define the Green function, which is the solution of the following equation:
  
  
      G”_ + 6/η G_' = δ(η - η̅).
  
  
  Concretely, we can write down the Green function as
  
  
      G_ = Θ(η - η̅) η̅^6/5[ 1/η̅^5 -1/η^5] 
      ⇒    k G_ =  Θ(η - η̅) x̅^6/5[ 1/x̅^5 -1/x^5].
  
  
  Using this Green function, we can solve Eq. (<ref>) as
  
  
      Ψ^ ( k,η) =   Ψ^( k, 0) + ∫^η_0 η̅  G_(k,η;η̅) S^_ ( k,η) 
      
      	=   Ψ^( k, 0) + ∫^3 k̃/(2π)^3 uv I_Ψ,,(u,v,x) ( 3/5)^2 ζ(k̃) ζ( k - k̃),
  
  
  where I_Ψ,, is defined as 
  
  
      I_Ψ,,(u,v,x) ≡   ∫^x_0 x̅   kG_ (k,η; η̅) f_(u,v) 
      
      	=   x^2/14 f_(u,v).
  
  
  Note that this solution of Ψ^ was already obtained in Ref. <cit.>.
  From Eq. (<ref>), we can obtain the initial condition for Ψ^ as
  
  
      Ψ^( k, 0) = ∫^3 k̃/(2π)^3 uv I_Ψ,,i(u,v) ( 3/5)^2 ζ(k̃) ζ( k - k̃),
  
  
  where 
  
  
      I_Ψ,m,i(u,v)    ≡2(u^2+v^2) - 3(u^2 - v^2)^2 +5 - 10 a_NL/3uv.
  
  
  Finally, we can rewrite Eq. (<ref>) as
  
  
      Ψ^ ( k,η) =   ∫^3 k̃/(2π)^3 uv I_Ψ,(u,v,x) ( 3/5)^2 ζ(k̃) ζ( k - k̃),
  
  
  where I_Ψ,≡ I_Ψ,,i + I_Ψ,,.
  
  For the energy density perturbation, Eq. (<ref>) leads to the analytic solution given as
  
  
      δ^ ( k,η) =   ∫^3 k̃/(2π)^3 uv I_δ,(u,v,x) ( 3/5)^2 ζ(k̃) ζ( k - k̃),
  
  
  where 
  
  
      I_δ,(u,v,x) ≡    -(2+ x^2/6) I_Ψ,(u,v,x) - x  I_Ψ,(u,v,x)/ x + I_δ,,(u,v,x), 
      
      	I_δ,,(u,v,x) ≡   1/uv[ - 2 x^2/3( u^2 + v^2 ) 
      	 - 5x^2/181-u^2-v^2/2 
      	 + 20/3( 2(u^2+v^2) - 3(u^2 - v^2)^2 +1/4)  ].
  
  
  In the late-time limit, we obtain I_δ,(u,v,x) ≃ -x^4/84 f_ (u,v), which is consistent with Refs. <cit.>.
  
  Similarly, from Eq.(<ref>), we obtain the analytic solution of the curvature perturbations as 
  
  
      ζ^ ( k,η) =   ∫^3 k̃/(2π)^3 uv I_ζ,(u,v,x) ( 3/5)^2 ζ(k̃) ζ( k - k̃),
  
  
  where 
  
  
      I_ζ,(u,v,x) ≡    - I_Φ, (u,v,x) + 1/3 I_δ,(u,v,x)  + I_ζ,,(u,v,x), 
      
      	I_ζ,,(u,v,x) ≡    -1/uv{1/9  T_δ_(vx) T_δ_(ux)  - 1/18(vx  T_δ_(v x)/ (v x) T_δ_(ux) + ux  T_δ_(u x)/ (u x)T_δ_(vx)   ) . 
               
      	+ 2/3( T_δ_(vx) + T_δ_(ux) ) 
               . - x^2/72[ (u^2 - v^2)^2 - 2(u^2 + v^2) + 1/4]  T_δ_(vx) T_δ_(ux) }.
  
  
  Note that we have substituted T_Ψ_(x)=1.
  After the source perturbations enter the horizon (ux ≫ 1 and vx ≫ 1), I_ζ, can be approximated as 
  
  
      I_ζ,(u,v,x) ≃x^6 uv/10368( (u^2 - v^2)^2 - 2(u^2 + v^2) + 1 ).
  
  
  
  Figure <ref> shows the evolution of I_Ψ,, I_δ,, and I_ζ,.
  We can see that I_Ψ, and I_δ, are constant on superhorizon scales (x ≪ 1) and start to grow when the induced perturbations enter the horizon (x ≃ 1). 
  In the late time (x ≫ 1), I_Ψ, and I_δ, are proportional to x^2 (∝ a) and x^4 (∝ a^2), respectively.
  We can also see that the values in the late time do not depend on u, v, and a_NL so much.
  As for I_ζ,, it is constant for the case of a_NL = 1 and on the other hand grows proportionally to x^2 for a_NL = 0 before the source perturbations enter the horizon (x ≪ 1/u and 1/v).
  These behaviors are similar to those in a RD era, discussed in the previous section.
  After the source perturbations enter the horizon (x ≳ 1/u and 1/v),  I_ζ, grows proportionally to x^6 (∝ a^3) regardless of the value of a_NL.
  From these observations, we conclude that the behaviors of the second-order perturbations, Ψ^, δ^, and ζ^, in the late time (x ≫ 1) are different from those of the first-order perturbations, Ψ^∝ a^0, δ^∝ a, and ζ^∝ a.
  
  
  
  
  
  
  
  
  
   §.§ Auto-power spectra
  
  
  
  Here, similarly to the case of a RD era, we calculate the auto-power spectra of the induced scalar perturbations.
  Note again that we just use the power spectra to show the evolution of the second-order perturbations.
  To compare with the observation, we must pay attention to the cross-power spectra between the first-order and the third-order perturbations.[Ref. <cit.> showed that the leading contributions from the auto- and the cross-power spectra are canceled out in the IR limit, where wavenumbers of the two perturbations contributing to the power spectrum are significantly smaller than those of the other two perturbations in the power spectrum.]
  The power spectra are given similarly to Eq. (<ref>) as 
  
  
      𝒫_ℳ^(k,η) ≡∫_0^∞ v  ∫_|v-1|^v+1 u   I_ℳ,^2(u,v,x) ( 3/5)^4 𝒫_ζ^( k v) 𝒫_ζ^(k u),
  
  
  where ℳ denotes Ψ, δ, and ζ. 
  Note that the normalization factor (3/5)^4 is different from the one in a RD era, given in Eq. (<ref>).
  In the case of the delta-function power spectrum given by Eq. (<ref>), Eq. (<ref>) becomes 
  
  
      𝒫_ℳ^(k,η) =  A_ζ^2 ( 3/5)^4 ( k_*/k)^2 I^2_ℳ,(k_*/k,k_*/k,x) Θ(2- k/k_*).
  
  
  Figure <ref> shows the evolution of Eq. (<ref>) for Ψ^, δ^, and ζ^ with a_NL = 1.
  From this figure, we can see that 𝒫_Ψ^ and 𝒫_δ^ are constant on superhorizon scales and start to grow once they enter the horizon (x ≳ 1).
  On the other hand, 𝒫_ζ^ grows even on superhorizon scales once the source perturbations enter the horizon, similarly to that in a RD era. 
  See also Appendix <ref> for the auto-power spectra with the top-hat power spectrum of the curvature perturbations.
  
  
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSION
  
  
  
  
  In this paper, we have discussed the second-order scalar perturbations induced by the first-order scalar perturbations.
  Specifically, we have considered the two eras, a RD and a MD era, and obtained the analytic solutions for the gravitational potential Ψ^ and the energy density perturbation δ^ in the conformal Newtonian gauge, and the curvature perturbation ζ^ in the uniform density gauge.
  We have also shown the evolution of the second-order perturbations from outside to inside horizon.
  Here, we briefly summarize the updates from the previous analyses.
  
  For the case of a RD era, we have taken into account the contributions of 𝒥_h and 𝒴_h, which is given in Eq. (<ref>) and were neglected in Ref. <cit.>.
  In particular, we have found that the leading contribution in the scales much larger than the source perturbations (u ≫ 1 and v ≫ 1), coming from 𝒥_i, is canceled by 𝒥_h in the late-time (subhorizon) limit.
  
  For the case of a MD era, while the analyses on δ^ in Refs. <cit.> were based on the Newtonian approximation, we have derived the analytic solutions of the induced perturbations without the approximation.
  Our solutions are valid on all scales, including superhorizon scales, and are consistent with the previous results in the subhorizon limit, where the approximation is valid.
  Apart from the references, Ref. <cit.> also derived the analytic solution of Ψ^ in a MD era and our solution is consistent with the reference. The updates from the reference for the MD era case is the derivation of the analytic solutions of δ^ and ζ^.
  
  
  Finally, we mention some issues for future work.
  First, future studies on the effect of the anisotropic stress of the fluid need to be performed.
  Although we have assumed a perfect fluid and neglected the anisotropic stress for simplicity, the fluid in a realistic situation can have the anisotropic stress originating from the free-streaming particles, such as photons and neutrinos.
  The anisotropic stress causes the diffusion damping of the scalar perturbations on subhorizon scales in a RD era <cit.> and makes a difference between Φ^ and Ψ^ in both the eras, which could modify the analytic solutions derived in this paper.
  Second, the effects of the transition between a RD era and a MD era need to be discussed.
  In reality, the Universe has experienced the transition between the two eras.
  Since the second-order scalar perturbations depend on the subhorizon evolution of the first-order perturbations, the perturbations that enter the horizon in the first era and experience the transition at subhorizon scales could make a non-negligible contribution to the second-order perturbations in the subsequent era.
  In particular, this transition effect is closely related to the Meszaros effect, which is already discussed in the context of the second-order perturbations in Ref. <cit.>.
  Third, future work on the third-order perturbations needs to be performed.
  As we mentioned in Sec. <ref> and <ref>, to compare with observation, we need to take into account the cross-power spectra between the first-order and the third-order perturbations, which are of the same order of magnitude as the auto-power spectrum of the second-order perturbations.
  In particular, although the third-order perturbations in a MD era are already discussed in the context of the large scale structure <cit.> and CMB anisotropies <cit.>, those in a RD era have never been discussed and need to be studied in future work.
  
  
  
  
  
  
  
  We thank Takahiro Terada for explaining how to perform the integrals in Ref. <cit.> and Wayne Hu for useful comments on the manuscript.
  We also acknowledge the anonymous referee for her/his helpful comments.
  We used the Mathematica package xPand <cit.> during our work.
  This work was supported by JSPS KAKENHI Grant Numbers, 15H02082 and 20H05248.
  
  
  
  
  
  
  
  § DERIVATION OF EQUATIONS FOR THE PERTURBATIONS
  
  
  
  
  In this appendix, we derive the equations given in Sec. <ref> from the Einstein equations.
  Note that we assume a perfect fluid and follow the notation of the perturbations in Sec. <ref> unless otherwise noted.
  
  
  
  
   §.§ First-order perturbations
  
  
  
  First, we discuss the perturbations at first order.
  The Einstein and the energy-momentum tensors at first order are given as 
  
  
      δG^0_ 0^ =   1/a^2( 6 ℋ^2 Φ^ + 6 ℋΨ^' 
        - 2 Ψ^,i_   ,i), 
      δG^0_ i^ =   1/a^2( -2 ℋΦ^_,i - 2 Ψ^'_, i), 
      δG^i_ j^ =   1/a^2[( 2 ℋ^2 Φ^ + 4 ℋ' Φ^ + 2 ℋΦ^' + 4 ℋΨ^' + 2 Ψ^” + Φ^ ,k_    ,k - Ψ^ ,k_    ,k) δ^i_ j. 
             . -  Φ^ ,i_    ,j + Ψ^ ,i_    ,j], 
      δT^0_ 0^ =    - δρ^, 
      δT^0_ i^ =    (ρ^ + P^) δ v^_   ,i, 
      δT^i_ j^ =   δ P^δ^i_ j.
  
  
  From the perturbed Einstein equations G^μ_ ν^ = T^μ_ ν^/M_^2, we can derive some equations.
  The (i,j) entry with i≠ j leads to Φ^ = Ψ^ and the (0,i) entry leads to
  
  
      δ v_i = 1/3 (1 + w) ℋ^2( -2 ℋΦ^_,i - 2 Φ^'_, i),
  
  
  where we have used the relation Φ^ = Ψ^.
  Using the (0,0) entry and the (i,j) entry with i=j, we can derive the equations for the gravitational potential and the energy density perturbation (Eqs. (<ref>) and (<ref>)), 
  
  
      Φ^” + 3 (1+c_^2) ℋΦ^' - c_^2 Φ^ ,i_   ,i + (2 ℋ' + (1 + 3 c_^2) ℋ^2) Φ^ = 0, 
         δ^≡δρ^/ρ^ = - ( 2 Φ^ + 2/ℋΦ^' - 2/3ℋ^2Φ^  ,i_     ,i).
  
  
  Solving these equations, we obtain Eqs. (<ref>)-(<ref>).
  
  
  
  
  
   §.§ Second-order perturbations
  
  
  
  
  Next, we discuss the perturbations at second order.
  The Einstein tensor at second order is given as 
  
  
      δG^0_ 0^ =    1/a^2[ 3 ℋ^2 Φ^ + 3 ℋΨ^' - Ψ^,i_    ,i -12 ℋ^2 (Φ^)^2 -12 ℋΦ^Ψ^'  +12 ℋΨ^Ψ^' . 
         .      
      	- 3 (Ψ^' )^2 - 8 Ψ^Ψ^ ,i_   ,i - 3 Ψ^ ,iΨ^_   ,i], 
      δG^i_ j^ =    1/a^2{( ℋ^2 Φ^ + 2 ℋ' Φ^ + ℋΦ^' + 2 ℋΨ^' + Ψ^” + 1/2Φ^ ,k_    ,k - 1/2Ψ^ ,k_    ,k) δ^i_ j. 
              .  - 1/2Φ^ ,i_    ,j + 1/2Ψ^ ,i_    ,j. 
              
      	+ [ -4 ℋ^2 ( Φ^)^2 - 8 ℋ' ( Φ^)^2 - 8 ℋΦ^Φ^' - 8 ℋΦ^Ψ^' - 2 Φ^' Ψ^'  . 
                  
      	+ 8 ℋΨ^Ψ^' + ( Ψ^' )^2 -4 Φ^Ψ^” + 4 Ψ^Ψ^” - 2 Φ^Φ^ ,k_   ,k
          . ( Φ^)^2 + 2 Ψ^Φ^ ,k_   ,k - 4 Ψ^Ψ^ ,k_   ,k - Φ^ ,kΦ^_   ,k - 2 Ψ^ ,kΨ^_   ,k] δ^i_ j
              
      	+ Φ^ ,iΦ^_   ,j - Ψ^ ,iΦ^_   ,j - Φ^ ,iΨ^_   ,j + 3 Ψ^ ,iΨ^_   ,j
         .   1/2
      	+ 2 Φ^Φ^ ,i_   ,j - 2 Ψ^Φ^ ,i_   ,j + 4 Ψ^Ψ^ ,i_   ,j} + ( terms with  V_i^ or  h^_ij).
  
  
  Substituting Φ^ = Ψ^ into these equations, we obtain
  
  
      δG^0_ 0^ 
      	=    1/a^2[ 3 ℋ^2 Φ^ + 3 ℋΨ^' - Ψ^,i_    ,i -12 ℋ^2 (Φ^)^2 
      	- 3 (Φ^' )^2 - 8 Φ^Φ^ ,i_   ,i - 3 Φ^ ,iΦ^_   ,i], 
      δG^i_ j^ 
      	=
      	   1/a^2{( ℋ^2 Φ^ + 2 ℋ' Φ^ + ℋΦ^' + 2 ℋΨ^' + Ψ^” + 1/3Φ^ ,k_    ,k - 1/3Ψ^ ,k_    ,k) δ^i_ j. 
              
      	- 1/2Λ^i   k_ j   lΦ^ ,l_    ,k +  1/2Λ^i   k_ j   lΨ^ ,l_    ,k
              
      	+ [ -4 ℋ^2 ( Φ^)^2 - 8 ℋ' ( Φ^)^2 - 8 ℋΦ^Φ^' - (Φ^')^2  . 
                 . - 8/3Φ^Φ^ ,k_   ,k - 7/3Φ^ ,kΦ^_   ,k] δ^i_ j
           . 1/2
      	+   2  Λ^i   k_ j   lΦ^ ,lΦ^_   ,k + 4  Λ^i   k_ j   lΦ^Φ^ ,l_   ,k} + ( vector and tensor modes),
  
  
  where Λ^i   k_ j   l is the projection operator onto the traceless longitudinal mode and explicitly expressed as[Note that ∂_i ∂_j/∇^2 is the projection operator onto longitudinal modes.]
  
  
      Λ^i   k_ j   l A^l_ k ( x) = 3/2( ∂^i ∂_j/∇^2- 1/3δ^i_ j) ( ∂^k ∂_l/∇^2 - 1/3δ^k_ l) A^l_ k ( x).
  
  
  With the notation of the perturbations in Eqs. (<ref>)-(<ref>), we obtain the second-order energy-momentum tensor, 
  
  
      δT^0_ 0^   = - ( 1/2δρ^ + (ρ^ + P^) δ v^ ,iδ v^_   ,i), 
      δT^i_ j^   =  1/2δ P^δ^i_ j + (ρ^ + P^) δ v^ ,iδ v^_   ,j.
  
  
  From the traceless longitudinal mode of the perturbed Einstein equation δG^μ_ ν^ = δ T^μ_ ν^/M_^2, we obtain the following relation:
  
  
      1/a^2 N^j_ i( -1/2Φ^,i_   ,j + 1/2Ψ^,i_   ,j +  2 Φ^ ,iΦ^_   ,j + 4 Φ^Φ^ ,i_   ,j) 
      	 = ρ^ + P^/M_^2 N^j_ i (δ v^ ,iδ v^_   ,j).
  
  
  Note again N^j_ i is defined in Eq. (<ref>).
  Substituting Eq. (<ref>) into this equation, we obtain
  
  
      N^j_ i[-1/2Φ^,i_   ,j + 1/2Ψ^,i_   ,j + 2(1+3w)/3(1+w)Φ^ ,iΦ^_   ,j - 4/3(1+w) ℋ(Φ^ ,iΦ^_   ,j' + Φ^ ,i' Φ^_   ,j) . 
                                         . - 4/3(1+w) ℋ^2Φ^ ,i' Φ^_   ,j' + 4 Φ^Φ^ ,i_   ,j] = 0 
      ⇒     Φ^ = Ψ^ + 4 (Φ^)^2 - N^j_ iB^i_ j^,
  
  
  where B^i_ j^ is defined as 
  
  
      B^i_ j^≡[4(5+3w)/3(1+w)Φ^ ,iΦ^_   ,j 
      + 8/3(1+w) ℋ(Φ^ ,iΦ^_   ,j)' 
       + 8/3(1+w) ℋ^2Φ^ ,i' Φ^_   ,j'].
  
  
  Here, we use the relation δ P^ = c_^2 δρ^, which is valid for the adiabatic perturbations with a constant sound speed. 
  Then, from Eqs. (<ref>) and (<ref>) and the diagonal part of Eqs. (<ref>) and (<ref>), we obtain (Eqs.(<ref>) and (<ref>))
  
  
      c_^2/a^2[ 3 ℋ^2 Φ^ + 3 ℋΨ^' - Ψ^,i_    ,i -12 ℋ^2 (Φ^)^2 
      	- 3 (Φ^' )^2 - 8 Φ^Φ^ ,i_   ,i - 3 Φ^ ,iΦ^_   ,i] 
          + 1/a^2{ℋ^2 Φ^ + 2 ℋ' Φ^ + ℋΦ^' + 2 ℋΨ^' + Ψ^” + 1/3Φ^ ,i_    ,i - 1/3Ψ^ ,i_    ,i. 
             .
      	- [ 4 ℋ^2 ( Φ^)^2 + 8 ℋ' ( Φ^)^2 + 8 ℋΦ^Φ^' + (Φ^')^2  + 8/3Φ^Φ^ ,i_   ,i + 7/3Φ^ ,iΦ^_   ,i] }
         
      	= ρ^ + P^/M_^2(1/3 - c_^2 ) δ v^ ,iδ v^_   ,i
      ⇒     Ψ^” + 3 (1 + c_^2) ℋΨ^' + [2ℋ' + (3c_^2 + 1) ℋ^2 ] Ψ^ - c_^2 Ψ^,i_    ,i = S^,
  
  	
  where S^ is the source term, defined as 
  
  
      S^≡   (3c_^2 -1/3) Φ^ ,iΦ^_   ,i + 8c_^2 Φ^Φ^ ,i_   ,i + (3c_^2 + 1) (Φ^' )^2 + [ (3c_^2 + 1)ℋ^2 + 2 ℋ' ]  N^j_ iB^i_ j^
          
      	 + ℋ  N^j_ i(B^i_ j^)' + 1/3  N^j_ i(B^i_ j^)^,k_ ,k + (1/3 - c_^2 ) 4/3(1+w) ℋ^2( ℋΦ^,i +  Φ^,i' )  (  ℋΦ^_,i + Φ^'_, i).
  
  
  
  
  
  From the (0,0) entry of the perturbed Einstein equation, δG^0_ 0^ = δ T^0_ 0^/M_^2, we can obtain the equation for the second-order energy density perturbation, 
  
  
      δρ^ =    -2M_^2/a^2[ 3 ℋ^2 Φ^ + 3 ℋΨ^' - Ψ^,i_    ,i -12 ℋ^2 (Φ^)^2 
      	- 3 (Φ^' )^2 - 8 Φ^Φ^ ,i_   ,i - 3 Φ^ ,iΦ^_   ,i] 
         -2(ρ^ + P^) δ v^ ,iδ v^_   ,i.
  
  	
  Substituting Φ^ = Ψ^ and Eq. (<ref>), we obtain (Eq. (<ref>))
  
  
      δ^≡   δρ^/ρ^
      
      	=    -2 Ψ^ +2 N^j_ i B^i_ j -2/ℋΨ^' + 2/3ℋ^2Ψ^,i_    ,i 
      	+ 2/ℋ^2(Φ^' )^2 + 16/3ℋ^2Φ^Φ^ ,i_   ,i
          + 1/ℋ^2( 2 - 8/9(1+w)) Φ^ ,iΦ^_   ,i - 8/9(1+w) ℋ^3 (Φ^ ,iΦ^_   ,i )'  - 8/9(1+w) ℋ^4Φ^'_, iΦ^,i'.
  
  
  
  Finally, we summarize the concrete expression of the second-order curvature perturbation in Ref. <cit.>.[We have fixed some typos in Ref. <cit.>.]
  To show the expression in a general gauge, we restore some metric perturbations that are killed in the conformal Newtonian gauge as
  
  
      s^2    = g_μν x^μ x^ν
          = a^2 {- (1 + 2 Φ^(1) + Φ^(2))η^2  + (2B_,i^ + B_,i^) η x^i . 
                 .+ [ (1-2 Ψ^(1) - Ψ^(2)) δ_ij + 2E^_,ij + E^_,ij]  x^i  x^j },
  
  
  where we have neglected vector and tensor perturbations. 
  One of the often-used gauge-invariant quantities is the curvature perturbation in the uniform density hypersurfaces (δρ = 0).
  Although the condition of uniform density is enough to define the gauge-invariant curvature perturbation at first order in perturbations, it is not enough to define the second-order curvature perturbation.
  In this paper, following Ref. <cit.>, we impose the additional condition E^ =0 to define the second-order curvature perturbation.
  Then, we can define the curvature perturbation as <cit.> 
  
  
      ζ^ = -Ψ^ + ℋα^ + 1/4𝒳^k_ k - 1/4∇^-2𝒳_ij^  ,ij,
  
  
  where the quantities in a general gauge are given as
  
  
      α^≡    - δρ^/ρ^' + δρ^δρ^'/(ρ^')^2 + δρ^_,k/ρ^'E^  ,k, 
      𝒳_ij≡      2 [ (2ℋ^2 + ℋ' - ℋρ^”/ρ^') ( δρ^/ρ^')^2 + ℋδρ^δρ^'/ (ρ^')^2 + ℋδρ^_,k/ρ^' E^  ,k] δ_ij
           + 4 [δρ^/ρ^' (Ψ^' + 2 ℋΨ^) + Ψ^_,k E^  ,k] δ_ij
          + 8 Ψ^ E^_,ij - 2  δρ^/ρ^'E^_,ij' -4 E^_,ik E^  ,k_     ,j - 2 (B^_,iδρ^_,j/ρ^' + B^_,jδρ^_,i/ρ^')  -2 δρ^_,i δρ^_,j/ (ρ^')^2
          - 2 E^_,ijk E^  ,k+  ( E^_,i' δρ^_,j/ρ^' + E^_,j' δρ^_,i/ρ^').
  
  
  Then, we obtain
  
  
      𝒳^k_ k 
      	=     6 [ (2ℋ^2 + ℋ' - ℋρ^”/ρ^') ( δρ^/ρ^')^2 + ℋδρ^δρ^'/ (ρ^')^2 + ℋδρ^_,k/ρ^' E^  ,k] 
          + 12 [ Ψ^_,k E^  ,k + δρ^/ρ^'( Ψ^' +2 ℋΨ^) ]    
         
      	- 2 E^  ,k∇^2 E^_ ,k  -2  δρ^/ρ^'∇^2 E^'  +  8 Ψ^∇^2 E^ -4 E^  ,kl E^_,kl
          - 2δρ^_ k/ρ^'( 2B^  ,k - E^  ,k' )  - 2 δρ^  ,kδρ^_,k/(ρ^')^2,
  	
  
  
  
      ∇^-2𝒳_ij^  , ij 
       =     2 [ (2ℋ^2 + ℋ' - ℋρ^”/ρ^')  ( δρ^/ρ^')^2 + ℋδρ^δρ^'/ (ρ^')^2 + ℋδρ^_,k/ρ^' E^  ,k]  
          + 4  [δρ^/ρ^' (Ψ^' + 2 ℋΨ^) + Ψ^_,k E^  ,k]  
         
       + ∇^-2[ 8  Ψ^ E^_,ij - 2 δρ^/ρ^'E^_,ij'  -4 E^_,ik E^  ,k_     ,j  - 4 B^_,iδρ^_,j/ρ^'  -2 δρ^_,i δρ^_,j/ (ρ^')^2. 
         .          
       - 2 E^_,ijk E^  ,k + 2 E^_,i' δρ^_,j/ρ^']^,ij.
  
  
  Note that if we take the gauge satisfying δρ^ = δρ^ = E^ = 0, we have ζ^ = Ψ^ by definition.
  Taking the Newtonian gauge (E = B =0), we finally obtain Eq. (<ref>).
  
  
  
  
  
  
  
  § COMPLETE EXPRESSIONS OF M^A_NM, 𝒥_H AND 𝒴_H
  
  
  
  
  In this appendix, we summarize the complete expressions of M^a_nm, given in Eq. (<ref>), and 𝒥_h and 𝒴_h, defined in Eq. (<ref>). 
  
  First, we show the expressions of M^a_nm(u,v) (1 ≤ n ≤ 8,   0 ≤ m ≤ 7):
  
  
      M^j_1 0(u,v) = 0, 
         M^j_1 1(u,v) = -3 (3 u^4-2 u^2 (3 v^2+1)+3 v^4-2 v^2-1) (u (v (u+v+1)-2)-2 v)/16 u^3 v^3, 
         M^j_1 2(u,v) = 0, 
         M^j_1 3(u,v) = 1/8 u^4 v^4 (9 (-3 u^7 (v-1)+u^6 (3 (7-3 v) v+3)+u^5 (v (3 v (v+5)+11)-11) . 
                      
      	 +u^4 (v (3 v (v (6 v-13)+1)-23)-2)+u^3 (v (v (v (3 (v-13)
         v-14)+2)-11)+11) 
                      
         +u^2 (v (v (v (3 v ((5-3 v) v+1)+2)-5)+17)-1) 
                      
         +u (v (v (v (v (v (11-3 (v-7) v)-23)-11)+17)-5)+5) 
                     .         
         +v (3
         v^2+1) (v (v (v^2+v-4)-1)+5))), 
         M^j_1 4(u,v) = 0, 
         M^j_1 5(u,v) = 1/8 u^4 v^4 (27 (21 u^6+54 u^5 (v-1)-u^4 (3 v (7 v+18)+23)-36 u^3 (v-1) (3 v^2+1) . 
                      
         +u^2 (v (v (3 (36-7 v) v-10)+36)+17)+18
         u (v-1) (3 v^4-2 v^2+1) 
                     .
         +v (v (v (v (3 v (7 v-18)-23)+36)+17)-18)+5)), 
         M^j_1 6(u,v) = 0, 
         M^j_1 7(u,v) = -729 (3 u^4-2 u^2 (3 v^2+1)+3 v^4-2 v^2+1)/4 u^4 v^4, 
         M^j_2 m(u,v) = -M^j_1 m(-u,v), 
         M^j_3 m(u,v) = -M^j_1 m(u,-v), 
         M^j_4 m(u,v) = M^j_1 m(-u,-v), 
         M^j_5 0(u,v) = √(3)(-3 u^4+u^2 (6 v^2+2)-3 v^4+2 v^2+1)/16 u^2 v^2, 
         M^j_5 1(u,v) = 0, 
         M^j_5 2(u,v) = 1/16 u^4 v^4 (3 √(3)(3 u^7 (v-2) v+3 u^6 (v ((v-6) v-2)+2) .  
                       
         +2 u^5 v (v (-3 (v-1) v-4)+11)+u^4 (2 v (v (v (5-3 (v-6)
         v)+3)+2)-4) 
                      
         +u^3 v (v (v (v (3 v (v+2)+10)-28)+3)-22) 
                      
         +u^2 (v (v (v (v (v (3 (v-6) v-8)+6)+3)-10)+2)-2)  
                     . 
         -2 u v (3
         v^2+1) (v (v (v^2+v-4)-1)+5)+6 v^6-4 v^4-2 v^2))), 
         M^j_5 3(u,v) = 0, 
         M^j_5 4(u,v) = -1/8 u^4 v^4 (9 √(3)(3 u^7+21 u^6 (v-1)+u^5 (3 v (5 v-18)-11) . 
                         
         +u^4 (v (3 (7-13 v) v-23)+23)+u^3 (v (v (3 (36-13 v)
         v+2)+36)+11)
                      
         +u^2 (v (v (v (3 v (5 v+7)+2)+10)+17)-17)
                      
         +u (v (v (v (v (3 v (7 v-18)-23)+36)+17)-18)+5)
                     .   
         +v (v (v (v (v (3 (v-7)
         v-11)+23)+11)-17)+5)-5)), 
         M^j_5 5(u,v) = 0, 
         M^j_5 6(u,v) = 243 √(3) (u+v-1) (3 u^4-2 u^2 (3 v^2+1)+3 v^4-2 v^2+1)/4 u^4 v^4, 
         M^j_5 7(u,v) = 0, 
         M^j_6 m(u,v) = -M^j_5 m(-u,v), 
         M^j_7 m(u,v) = -M^j_5 m(u,-v), 
         M^j_8 m(u,v) = M^j_5 m(-u,-v), 
         M^y_1 m(u,v) = -M^j_5 m(u,v), 
         M^y_2 m(u,v) = M^j_5 m(-u,v), 
         M^y_3 m(u,v) = M^j_5 m(u,-v), 
         M^y_4 m(u,v) = -M^j_5 m(-u,-v), 
         M^y_5 m(u,v) = M^j_1 m(u,v), 
         M^y_6 m(u,v) = -M^j_1 m(-u,v), 
         M^y_7 m(u,v) = -M^j_1 m(u,-v), 
         M^y_8 m(u,v) = M^j_1 m(-u,-v).
  
  
  
  Next, we show the expressions of 𝒥_h and 𝒴_h, defined in Eq. (<ref>).
  To perform the integral, we use the following relation:[This relation is used to calculate the GWs induced by scalar perturbations analytically <cit.>.]
  
  
      ∫^x_2_x_1x̅sin (αx̅ + ϕ)/x̅^m =   [ ∑^m-2_k=0(m-k-2)!/(m-1)!α^k sin( αx̅ + ϕ + (k+2)/2π) x̅^1+k-m]^x_2_x_1
          - α^m-1/(m-1)!∫^x_2_x_1x̅1/x̅sin( αx̅ + ϕ + (m+1)/2π).
  
  
  After tedious calculations, we finally obtain
  
  
      𝒥_h(u,v,x) =    -3√(3)/8 u^4 v^4 x^6[ x sin( x/√(3)) ( 2u x cos(  u x/√(3)) {1/2. . . 
          v x [ 
      	 3 (x^2 + 18) u^4 -2 (3 (x^2 + 18)v^2 + 2(x^2 +9)) u^2 . 
             .  -3(x^2 -6) -4 v^2(x^2 +9) + 3v^4(x^2 + 18) ] cos(  v x/√(3)) 
         
      	 + √(3)[ 6 x^2 v^6 - (7x^2 + 54) v^4 + 6(x^2 + 6) v^2 + 3(x^2 -6) + u^4( (6v^2 -3)x^2 -54) . 
         . 1/2.     
      	 + 2 u^2(-6 x^2 v^4 + (x^2 + 54)v^2 + 2(x^2 + 9)) ] sin(  v x/√(3)) }
         + sin(  u x/√(3)) { 2√(3) v x [ 6 x^2 u^6 - ((12 v^2 +7)x^2 + 54)u^4  .  
                             
      	  + 2 (3x^2 v^4 + (x^2 + 54)v^2 + 3 (x^2 + 6))u^2 
                             . + 3(x^2 -6) + 4v^2(x^2 +9)-3 v^4(x^2 + 18) ] cos(  v x/√(3)) 
          + [ ((6v^2 -3) x^4 - 36 x^2) u^6 + (-(12 v^4 +v^2 +3)x^4 + 6(6v^2 + 7)x^2 + 324) u^4  . 
          + (6x^4 v^6 -x^2 (x^2 -36) v^4 + 4(x^4 + 3x^2 - 162)v^2 + 3(x^4 - 12x^2 -72)) u^2 
          -3 (x^2(x^2 + 12) v^6 + (x^4 - 14 x^2 - 108)v^4 
         . + (-x^4 + 12 x^2 + 72) v^2 - x^4 + 6 x^2 -36) ] sin(  v x/√(3)) }
         + cos(  x/√(3)) ( -u x cos( u x/√(3)){[3x^4 u^6 + ((3-9v^2)x^4 + (54-36v^2)x^2+324) u^4 .  
          + 3(3x^2 (x^2+8) v^4 -2(x^4 +14 x^2 +  108)v^2 -x^4 - 16 x^2 - 72)u^2 - 3v^6 x^2(x^2 + 12) 
         . -3(x^4 + 2 x^2 - 36) + 3 v^4 (x^4 + 26x^2 + 108) - v^2(5x^4 + 60x^2 + 216) ] sin(  v x/√(3)) 
         -2√(3) v x [ 9(x^2 + 6) u^4 -2(9(x^2 + 6)v^2 + 4 x^2 + 18) u^2 - x^2 . 
         . + 9 v^4(x^2 + 6) -4v^2(2x^2+9) + 18  ] cos(  v x/√(3)) }
          - sin(  u x/√(3)) { v x  [ -3 x^2 (x^2 + 12)u^6 + 3((3v^2 + 1) x^4 + (24 v^2 + 26) x^2 + 108) u^4  . 
          - (9 x^2(x^2 + 4)v^4 + 6(x^4 + 14 x^2 + 108)v^2 + 5x^4 + 60x^2 + 216 ) u^2 
         . + 3(x^4 v^6 + (x^4 + 18 x^2 + 108) v^4 - (x^4 + 16 x^2 + 72) v^2 -x^4 - 2x^2 + 36) ] cos(  v x/√(3)) 
          + √(3)[ 3x^2 ((1-2v^2) x^2 + 12) u^6 + ((12 v^4 + v^2 -3) x^4 -6(6v^2 + 13)x^2 -324) u^4  . 
          + (-6x^4 v^6 + x^2 (x^2 -36)v^4 + (8x^4 + 60 x^2 + 648) v^2 + 5x^4 + 60 x^2 + 216) u^2 
         . +  3v^6 x^2 (x^2 + 12) + 3(x^4 + 2 x^2 -36) - 3v^4 (x^4 + 26 x^2 + 108) . 
         . + v^2 (5 x^4 + 60 x^2 + 216) ] sin(  v x/√(3)) }) ] 
         + 3/32 u^4 v^4( 3u^4(u^4 - 2) + 3v^4(v^4 -2) + 3 - 8u^2v^2(1 + u^2 + v^2) -6 u^4 v^4 ) 
             ×( (|-1+u-v|/√(3) x ) + (1+u-v/√(3) x) - (-1+u+v/√(3) x ) - (1+u+v/√(3) x )  ) 
         + 1/8 u^3 v^3[-9(u^6 + u^4 -u^2) -9(v^6 + v^4 - v^2) + 9 + 6 u^2 v^2 (u^2 - v^2)^2  . 
         .          
      	+ 5 u^2 v^2(u^2 + v^2) + 8 u^2 v^2 ] 
         -3/32 u^4 v^4( 3u^4(u^4 - 2) + 3v^4(v^4 -2) + 3 - 8u^2v^2(1 + u^2 + v^2) -6 u^4 v^4 ) log( 1-(u-v)^2 /(u+v)^2 -1),
  
  
  
  
      𝒴_h (u,v,x) =    3√(3)/8 u^4 v^4 x^6[ x cos( x/√(3)) ( 2u x cos(  u x/√(3)) {1/2. . . 
          v x [ 
      	 3 (x^2 + 18) u^4 -2 (3 (x^2 + 18)v^2 + 2(x^2 +9)) u^2 . 
         .      
      	 -3(x^2 -6) -4 v^2(x^2 +9) + 3v^4(x^2 + 18) ] cos(  v x/√(3)) 
         
      	 + √(3)[ 6 x^2 v^6 - (7x^2 + 54) v^4 + 6(x^2 + 6) v^2 + 3(x^2 -6) + u^4( (6v^2 -3)x^2 -54) . 
         . 1/2.     
      	 + 2 u^2(-6 x^2 v^4 + (x^2 + 54)v^2 + 2(x^2 + 9)) ] sin(  v x/√(3)) }
         + sin(  u x/√(3)) { 2√(3) v x [ 6 x^2 u^6 - ((12 v^2 +7)x^2 + 54)u^4  .  
                             
      	 + 2 (3x^2 v^4 + (x^2 + 54)v^2 + 3 (x^2 + 6))u^2 
                             . + 3(x^2 -6) + 4v^2(x^2 +9)-3 v^4(x^2 + 18) ] cos(  v x/√(3)) 
          + [ ((6v^2 -3) x^4 - 36 x^2) u^6 + (-(12 v^4 +v^2 +3)x^4 + 6(6v^2 + 7)x^2 + 324) u^4  . 
          + (6x^4 v^6 -x^2 (x^2 -36) v^4 + 4(x^4 + 3x^2 - 162)v^2 + 3(x^4 - 12x^2 -72)) u^2 
         . -3 (x^2(x^2 + 12) v^6 + (x^4 - 14 x^2 - 108)v^4 . 
         . + (-x^4 + 12 x^2 + 72) v^2 - x^4 + 6 x^2 -36) ] sin(  v x/√(3)) }
         + sin(  x/√(3)) ( u x cos( u x/√(3)){[3x^4 u^6 + ((3-9v^2)x^4 + (54-36v^2)x^2+324) u^4 .  
          + 3(3x^2 (x^2+8) v^4 -2(x^4 +14 x^2 +  108)v^2 -x^4 - 16 x^2 - 72)u^2 - 3v^6 x^2(x^2 + 12) 
         . -3(x^4 + 2 x^2 - 36) + 3 v^4 (x^4 + 26x^2 + 108) - v^2(5x^4 + 60x^2 + 216) ] sin(  v x/√(3)) 
         -2√(3) v x [ 9(x^2 + 6) u^4 -2(9(x^2 + 6)v^2 + 4 x^2 + 18) u^2 - x^2  . 
         . + 9 v^4(x^2 + 6) -4v^2(2x^2+9) + 18  ] cos(  v x/√(3)) }
          + sin(  u x/√(3)) { v x  [ -3 x^2 (x^2 + 12)u^6 + 3((3v^2 + 1) x^4 + (24 v^2 + 26) x^2 + 108) u^4  . 
          - (9 x^2(x^2 + 4)v^4 + 6(x^4 + 14 x^2 + 108)v^2 + 5x^4 + 60x^2 + 216 ) u^2 
         . + 3(x^4 v^6 + (x^4 + 18 x^2 + 108) v^4 - (x^4 + 16 x^2 + 72) v^2 -x^4 - 2x^2 + 36) ] cos(  v x/√(3)) 
          + √(3)[ 3x^2 ((1-2v^2) x^2 + 12) u^6 + ((12 v^4 + v^2 -3) x^4 -6(6v^2 + 13)x^2 -324) u^4  . 
          + (-6x^4 v^6 + x^2 (x^2 -36)v^4 + (8x^4 + 60 x^2 + 648) v^2 + 5x^4 + 60 x^2 + 216) u^2 
         . +  3v^6 x^2 (x^2 + 12) + 3(x^4 + 2 x^2 -36) - 3v^4 (x^4 + 26 x^2 + 108) . 
         . + v^2 (5 x^4 + 60 x^2 + 216) ] sin(  v x/√(3)) }) ] 
          -3/32 u^4 v^4( 3u^4(u^4 - 2) + 3v^4(v^4 -2) + 3 - 8u^2v^2(1 + u^2 + v^2) -6 u^4 v^4 )  
             ×( (-1+u-v/√(3) x ) - (1+u-v/√(3) x) - (-1+u+v/√(3) x ) + (1+u+v/√(3) x )  ),
  
  
  where Si and Ci denote the sine and cosine integrals, defined as 
  
  
      Si(x)    ≡∫^x_0  t sin t/t, 
      Ci(x)    ≡ -∫^∞_x   t cos t/t.
  
  
  In the late-time limit (x≫ 1), 𝒥_h and 𝒴_h become 
  
  
      𝒥_h(u,v,x (≫ 1)) ≃   1/8 u^3 v^3[-9(u^6 + u^4 -u^2) -9(v^6 + v^4 - v^2) + 9 . 
         .           + 6 u^2 v^2 (u^2 - v^2)^2 + 5 u^2 v^2(u^2 + v^2) + 8 u^2 v^2 ] 
         -3/32 u^4 v^4( 3u^4(u^4 - 2) + 3v^4(v^4 -2) + 3 - 8u^2v^2(1 + u^2 + v^2) -6 u^4 v^4 ) 
                 ×log( 1-(u-v)^2 /(u+v)^2 -1), 
      𝒴_h(u,v,x (≫ 1)) ≃   3 π/32 u^4 v^4( 3u^4(u^4 - 2) + 3v^4(v^4 -2) + 3 - 8u^2v^2(1 + u^2 + v^2) -6 u^4 v^4 ).
  
  
  These expressions correspond to 𝒥_h,late and 𝒴_h,late, defined in Eqs. (<ref>) and (<ref>).
  
  
  
  
  
  
  
  § POWER SPECTRA IN OTHER SITUATIONS
  
  
  
  
  
  In this appendix, we discuss the auto-power spectra of the second-order scalar perturbations in situations not considered in the main body of this paper.
  
  
  
  
  
   §.§ Radiation-dominated era
  
  
  
  
  First, we discuss the power spectra in a RD era.
  Throughout this paper, we assume a perfect fluid, in which there is no anisotropic stress.
  However, in realistic situations, the fluid deviates from a perfect fluid and the scalar perturbations can be diffused on deeply subhorizon scales in a RD era.
  This phenomenon is called the diffusion (or Silk) damping <cit.>.
  The CMB distortion, which can constrain the small-scale perturbations, is caused by the energy injection to CMB through the diffusion damping.
  Because of this, we need to be careful when we compare the second-order perturbations with the observational constraints.
  
  For the first-order perturbations, the CMB distortion is caused only by the perturbation on k ≲ 10^4Mpc^-1. 
  This is because the photon-number changing process such as the double Compton scattering, which prevents the CMB distortion, is effective at the moment when the smaller-scale (k ≳  10^4Mpc^-1) perturbations are erased through the diffusion damping (see e.g. Ref. <cit.> for the detail of the CMB distortion).
  
  For the second-order perturbations, the situation is a little bit different.
  The contributions of 𝒥_0 and 𝒴_0 continue to depend on the evolution of the source perturbations even on deep subhorizon scales.
  𝒥_0 and 𝒴_0 are modified (might be erased) on all scales at the moment when the source perturbation, whose wavenumber is k ≫ 10^4 Mpc^-1, is erased by the diffusion damping.
  Note that, even if the contributions of 𝒥_0 and 𝒴_0 on k ≲ 10^4Mpc^-1 are also erased at that moment, the distortion is not produced because the photon-number changing process is effective at that time.
  On the other hand, the other contributions from 𝒥_h, 𝒥_i, and 𝒴_h do not depend on the source perturbations when the induced perturbations are on deep subhorizon scales[Conversely, 𝒥_h and 𝒴_h do not converge until the induced perturbations enter the horizon (x ∼ 1), which means that if the source perturbations are diffused at η_d, the convergence values are different from 𝒥_h,late and 𝒴_h,late on the scales k ≲ 1/η_d. Here, we assume that the source perturbations are diffused well after entering the horizon and the convergence values are given by 𝒥_h,late and 𝒴_h,late on all scales for simplicity.] and therefore they are not affected by the diffusion damping of the source perturbations at least on the subhorizon scales.
  Similarly to the first-order perturbations, the contributions of 𝒥_h, 𝒥_i, and 𝒴_h are erased at the moment when their scales are smaller than the diffusion scale and, if their scales are k ≲ 10^4Mpc^-1, they produce the distortion, which is constrained by the CMB observation.
  Since the analysis of the diffusion damping effect on 𝒥_0 and 𝒴_0 is beyond the scope of this paper, we just focus on the other contributions from 𝒥_h, 𝒥_i, and 𝒴_h hereafter.
  Note that, even if we neglect 𝒥_0 and 𝒴_0, Ψ^(2) still depends on time because 𝒥 and 𝒴 are the coefficients in front of the first-order solutions, which are proportional to j_1(x/√(3))/x and y_1(x/√(3))/x (see Eq. (<ref>))[The solution of y_1(x/√(3))/x does not appear in the analytic solution of Φ^ (Eq. (<ref>)) because of its initial condition on superhorizon scales. ]. 
  
  
  To compare the induced perturbations with the primordial first-order perturbations, we define the following power spectra:
  
  
      𝒫_𝒥_const(k) ≡   ∫_0^∞ v  ∫_|v-1|^v+1 u  (𝒥_h,late(u,v) + 𝒥_i (u,v) )^2 ( 2/3)^4 𝒫_ζ^( k v) 𝒫_ζ^(k u), 
      	
      𝒫_𝒴_const(k) ≡   ∫_0^∞ v  ∫_|v-1|^v+1 u  𝒴_h,late^2(u,v) ( 2/3)^4 𝒫_ζ^( k v) 𝒫_ζ^(k u).
  
  
  Note that these power spectra do not depend on time because they only include the contribution independent of the source perturbations on deep subhorizon scale.
  At first order in perturbations, the coefficient 𝒥 corresponds to the amplitude of the curvature perturbations on superhorizon scale (see Eq. (<ref>)).
  For this reason, we compare the Eqs. (<ref>) and (<ref>) with 𝒫_ζ^ with the appropriate normalization, explained below. 
  At first order, the power spectra of Ψ^ and ζ^ are related as 𝒫_ζ^≃ (9/4)𝒫_Ψ^ on superhorizon scales in a RD era. 
  From this observation, we compare (9/16) 𝒫_𝒥_const and (9/16)𝒫_𝒴_const with 𝒫_ζ^, where we have additionally multiplied them by 1/4 because the normalization of Ψ^ is different from Ψ^ by 1/2 (see Eq. (<ref>)).
  For 𝒫_ζ^, in addition to the delta-function power spectrum, we also consider the top-hat power spectrum, given as 
  
  
      𝒫_ζ^(k) = A_ζΘ(k_* - k) Θ(k - k_*,l),
  
  
  where Θ(x) denotes the Heaviside step function.
  
  Figure <ref> shows the normalized spectra, (9/16) 𝒫_𝒥_const/A_ζ^2 and (9/16)𝒫_𝒴_const/A_ζ^2, with the delta-function and the top-hat power spectra. 
  From this figure, we can see that the peaks of the power spectra become higher for the broader 𝒫_ζ^ and the wavenumber dependences of 𝒫_𝒥_const and 𝒫_𝒴_const on large scales are given by ∼ k^3 for the top-hat power spectrum and ∼ k^2 for the delta-function one up to the logarithmic factor. 
  These wavenumber dependences are the same as those of the GWs induced by scalar perturbations <cit.>.
  
  
  
  
  
  
  
  
  Since the small-scale perturbations are often discussed in the context of PBHs, we here consider the power spectrum, which predicts a sizable amount of PBHs that can explain BHs detected by LIGO-Virgo collaborations.
  Specifically, we take A_ζ = 0.01 and k_* = 5 × 10^5Mpc^-1 as fiducial values <cit.>.
  Figure <ref> shows the spectra of the induced scalar perturbations with the fiducial values.
  Although the power spectra are consistent with the null detection of the CMB μ-distortion and the light element abundance determined by the big bang nucleosynthesis so far, they could possibly be investigated in the future by a PIXIE-like spectrometer <cit.>. 
  Note again that we have neglected the contributions from 𝒥_0 and 𝒴_0 and those from the cross-power spectrum between the first- and the third-order perturbations for simplicity.
  For these reasons, we should keep in mind that the observable spectra of the induced perturbations might be different from those in Fig. <ref>.
  
  
  
  
  
  
  
  
  
  
  
   §.§ Matter-dominated era
  
  
  
  
  Next, we discuss the power spectra in a MD era.
  Similarly to the previous subsection, we consider the top-hat power spectrum, given in Eq. (<ref>).
  Since the observations of CMB and the large scale structure have revealed that the Universe has experienced the MD era starting around z ≃ 3400 and the 𝒫_ζ^ is almost scale-invariant on the large scales (k ≲ 1 Mpc^-1), we especially focus on the case of a broad top-hat power spectrum with k_*,l≪ k_* as a fiducial example.
  Figure <ref> shows the evolution of the power spectra with k_*,l/k_* = 10^-3.
  For comparison, we also show the spectra with k_*,l/k_* = 10^-1 (black lines), from which we can see the spectra do not depend on k_*,l so much especially for x_* = 10^2 and 10^3.
  From this observation, we conclude that the main contribution comes from the perturbations on the smallest scales (∼ k_*).
  Similarly to the spectra in Fig. <ref>, 𝒫_Ψ^ and 𝒫_δ^ remain constant on superhorizon scales and start to evolve once they enter the horizon.
  Finally, 𝒫_Ψ^ and 𝒫_δ^ on subhorizon scales grow proportionally to η^4 ∝ a^2 and η^8 ∝ a^4, respectively.
  On the other hand, 𝒫_ζ^ finally grows proportionally to η^12∝ a^6 on all scales, including superhorizon scales, well after the perturbations around the smallest scales enter the horizon.
  In the late time (x_* ≫ 1), the wavenumber dependences of the spectra can be roughly summarized as follows:
  
  
      𝒫_Ψ^∝
      	k^-1     (x ≪ 1) 
      
      	k^3     (1 ≪ x ≲ x_*)
      	, 
         𝒫_δ^∝
      	k^-1     (x ≪ 1) 
      
      	k^7     (1 ≪ x ≲ x_*)
      	, 
         𝒫_ζ^∝
      	k^3   (x ≲ x_*).
  
  
  
  
  
  
  
  
  
  
  
  apsrev4-1
  
  
  