Gravitational wave from warm inflation
Xi-Bin Li He Wang Jian-Yang Zhu





lxbbnu@mail.bnu.edu.cn
Department of Physics, Beijing Normal University, Beijing 100875, China


hewang@mail.bnu.edu.cn
Department of Physics, Beijing Normal University, Beijing 100875, China


Corresponding author
zhujy@bnu.edu.cn
Department of Physics, Beijing Normal University, Beijing 100875, China




A fundamental prediction of inflation is a nearly scale-invariant spectrum of gravitational wave. The features of such a signal provide extremely important information about the physics of the early universe. In this paper, we focus on several topics about warm inflation. First, we discuss the stability property about warm inflation based on nonequilibrium statistical mechanics, which gives more fundamental physical illustrations to thermal property of such model. Then, we calculate the power spectrum of gravitational waves generated during warm inflation, in which there are three components contributing to such spectrum: thermal term, quantum term and cross term combining the both. We also discuss some interesting properties about these terms and illustrate them in different panels. As a model different from cold inflation, warm inflation model has its individual properties in observational practice, so we finally give a discussion about the observational effect to distinguish it from cold inflation.


98.80.-k, 98.80.Bp, 98.80.Es

Gravitational wave from warm inflation
    Jian-Yang Zhu
    December 30, 2023
======================================







§ INTRODUCTION


It has taken for several years to search for gravitational waves on astrophysical experiments, such as LIGO <cit.>, VIRGO <cit.> and other observational experiments <cit.> to test the prediction of general relativity. Hulse and Taylor <cit.> first announced the observational evidence for the existence of gravitational waves. Then came the good news that gravitational waves have been detected by LIGO <cit.> in recent years generated by binary black hole. This discovery was a great achievement that opened a new window to better understand our universe both at early epoch and late epoch.

Warm inflation model was established as a candidate scenario to overcome some defects in cold inflation <cit.>. However, it was realized a few years after its original proposal that the idea of warm inflation was not easy to realize in concrete models and even simply not be possible in relevant work <cit.>. Some problems were mentioned to suspect such scenario. First, it is hard to couple the inflaton directly with light fields. Considering a Yukawa interaction gϕψ̅ψ, the slow roll condition typically requires an inflaton with large value, while the fermion obtain an extra mass with also a large value unless the coupling a quite small. As a result, the dissipative effect may be too small either which implies that it is hard to obtain a period of inflation with dissipation coefficient long enough to require the 60 e-ford to solve the horizon and flatness problem. Second, a direct coupling to light fields may lead to large thermal corrections to the inflaton mass m_ψ=gϕ, which could prevent slow roll for T>H. Shortly afterwards successful models of warm inflation have been established, in which the inflaton indirectly interacts with the light degrees of freedom though a heavy mediator fields instead of being coupled with a light field directly <cit.>. In such scenarios dissipation can sustain both the slow-roll dynamics of the inflaton field and the temperature of the radiation bath for a sufficiently long period. One can read a Lagrangian density of the generic form,

    ℒ[Φ,X,Y]=ℒ[Φ]      +ℒ[X]+ℒ[Y]
          +ℒ_int[Φ,X]+ℒ_int[X,Y]

where Φ is the inflaton field, X are any field or degrees of freedom coupled directly to the inflaton field, while Y can be any other fields not necessarily coupled to the inflaton, but are coupled to X. The ℒ_int[Φ,X] and ℒ_int[X,Y] give the interaction among these fields. The evolution of the inflaton field can be properly determined in the context of the in-in, or the Schwinger closed-time path functional formalism <cit.>. This equation displays both dissipation and non-Markovian stochastic noise terms and it is a generalized Langevin-like equation of motion <cit.>.

Compared with the predictions of cold inflation that primordial density fluctuations mostly from quantum fluctuation and thermal bath are only generated at the end of inflation <cit.>, warm inflation model suggests that our universe is hot during the whole inflation when inflaton fields couple with the thermal bath and the primary source of density fluctuations come from thermal fluctuations <cit.>. The equation of motion for warm inflation can be written as a stochastic Lengevin equation, in which there is a dissipation term to describe the inflaton fields coupling with thermal bath and there is also a fluctuation term described by a stochastic noise term <cit.>. The fundamental principles of warm inflation have been reviewed recently in <cit.>. Warm inflation with strong dissipation requires Υ>H, with which thermal fluctuations dominate the whole inflationary epoch i.e. T>H <cit.>. Generally, we can still have the weak dissipation condition with Υ≪ H <cit.>. Noncanonical warm inflation model is under studying recently <cit.>.

Although primordial gravitational waves generated during inflation have not been detected yet till now, the discovery of gravitational waves by LIGO has shed a bright light to this prediction. The measurement of the cosmic microwave background radiation and other observations have given a good constraint to the ratio of tensor to scaler with r<0.07 at 95 CL <cit.>. The study of primordial gravitational wave is a way to prove the inflationary programm. What's more, it also provides an effective method to distinguish among different inflationary model. Now more and more researchers have cast their eyes to primordial gravitational waves generated by various kind of sources, such as primordial density perturbation <cit.>, some kind of inflationary model <cit.>, and even during reheating epoch <cit.>; a great deal of predictions have been given in these works which may be observed in observation experiments in the near future. However, till now there has not been any calculation about the gravitational waves generated from warm inflation. This approach must be an effective way to differentiate cold inflation model from warm inflation model.

In this paper, we attempt to illustrate the existence of gravitational wave generated from warm inflation. Green's function method has been used to calculate the power spectrum of tensor perturbation with source as the form of transverse-traceless tensor. This spectrum can separate
into three terms: thermal component, quantum component and cross component, and each component has its individual and interesting properties. With these properties, we discuss the observational prediction from the view of tensor-to-scalar ratio.

This paper is organized as follows: in Sec. <ref>, we give a brief introduction to warm inflation and stochastic approach to  deal with warm inflation. In Sec. <ref>, based on nonequilibrium statistical mechanics, we recalculate the statistical properties of warm inflation model and give more results in detail. In Sec. <ref> three programs for describing the thermal, quantum and their cross term are calculated and relevant calculations and equations are derived. In Sec. <ref>, we discuss our result by numerical analysis. Finally, in Sec. <ref>, we conclude our work and give some further discussions about our results.



§ STOCHASTIC APPROACH FOR WARM INFLATION



 §.§ Stochastic approach


First, let's have a brief review of stochastic approach for cold inflation. The equation of motion for cold inflation is the standard one:
In warm inflation model, the equation of background field is often written as the Langevin equation

    [∂^2 /∂ t^2+3H∂/∂ t-1/a^2∇^2]Φ+∂ V(Φ)/∂Φ = 0,

where Φ is the inflaton field operator, a is scale factor in Friedmann-Robertson-Walker metric and H is the Hubble parameter defended as H=ȧ/a. The stochastic approach assumes that inflaton field separates into two parts, one is Φ_> which denotes the long wavelength part, and another is Φ_< which denotes the short wavelength part for quantum vacuum fluctuation, i.e. Φ→Φ_>+Φ_<. Usually, Φ_< is written in terms of a filter (window) function.

The stochastic inflationary approach suggests that quantum inflaton field Φ is composed in a short wavelength part Φ_<, which denotes the quantum vacuum fluctuations, and a long wavelength part Φ_> i.e. Φ(x,t)=Φ_<(x,t)+Φ_>(x,t). As previously study, quantum fluctuation dominates on short wavelength. With this condition, we assume a number smaller than unity such that the quantum

    Φ_<(x,t)≡ϕ_q(x,t)=      ∫d^3k/(2π)^3/2W(k,t)
          ×[ϕ_k(t)e^-ix·kâ_k+h.c.],

where ϕ_k(t) is the field in momentum space, â_k is the annihilation operator whose Hermitian conjugate operators is â^†_k. In  (<ref>), W(k,t) are the window function with sharp momentum cutoff

    W(k,t)=θ(k-ϵ aH).

where ϵ is a suitable number smaller than 1. Then fields ϕ_k(t) in (<ref>) satisfy

    [∂^2 /∂ t^2+3H∂/∂ t-k^2/a^2+⟨ V,_ϕϕ(Φ_>)⟩]ϕ_k(t) = 0.

Generally, ϕ_k(τ) is given by

    ϕ_k(τ)=H√(π)/2(-τ)^3/2H_μ^(1)(-kτ).

where μ=√(9/4-m^2/H^2) with m^2≃ constant is the average of V,_ϕϕ(Φ_>), τ is the conformal time in de Sitter space-time defined as τ=-1/aH, and H^(1) is the Hankel function of the first kind.

Then, with  (<ref>), the equation of motion for long wavelength part of the field Φ_> reads

    [∂^2 /∂ t^2+3H∂/∂ t-1/a^2∇^2]Φ_>+V,_ϕ(Φ_>) = ξ_q,

where ξ_q is given by

    ξ_q=-[∂^2 /∂ t^2+3H∂/∂ t-1/a^2∇^2+ V,_ϕϕ(Φ_>)]ϕ_q(t)

and its average vanishes ⟨ξ_q⟩=0. With  (<ref>) and  (<ref>), the two-point correlation function is

    ⟨ξ_q(x,t)ξ_q(x',t')⟩ = ∫d^3k/(2π)^3e^ik·(x-x')Re[f_k(t)f_k'^*(t')],

where

    f_k(t)=[∂^2 W/∂ t^2+3H∂ W/∂ t]ϕ_k(t)+2 ∂ W/∂ t∂ϕ_k(t)/∂ t.

Here we just give a brief introduction of stochastic inflationary approach, and more details have been studied in several work
<cit.>.




 §.§ Warm inflation

In warm inflation model, the equation of motion of background field is often written as the Langevin equation
 
    [∂^2 /∂ t^2+(3H+Υ)∂/∂ t-1/a^2∇^2]Φ+∂ V(Φ)/∂Φ = ξ_T,

where Υ is the dissipation coefficient and ξ_T is the thermal noise fluctuation. In this paper, we consider only in the case of de Sitter space-time, where a(t)=exp(Ht) and H is regarded as a constant. According to the fluctuation-dissipation theorem, dissipation coefficient Υ and fluctuation noise ξ_T have the relation
 
    ⟨ξ_T(x,t)ξ_T(x,t')⟩=2Υ T a^-3δ(t-t').

The Fourier transformation of  (<ref>) is

    ⟨ξ_T(k,t)ξ_T(k',t')⟩=2(2π^3)Υ T a^-3δ^3(k+k')δ(t-t').

Usually Υ is a function of both background homogeneous inflaton field Φ and temperature T <cit.> and we do not attend to discuss this question here.

The the inflaton field operator Φ(x,t) is often separated into the parts as follow

    Φ(x,t)=ϕ(t)+δφ(x,t),

where δφ(x,t) is the perturbed part of inflaton, and ϕ(t) is the background homogeneous inflaton field which is defended as

    ϕ(t)=1/Ω∫_Ωd^3x Φ(x,t).

Here, Ω is particle horizon size Ω=1/H. With this relation,  (<ref>) reads

    ∂^2 ϕ/∂ t^2+[3H+Υ]∂ϕ/∂ t+V,_ϕ(ϕ)    =    0,
    {∂^2/∂ t^2+[3H+Υ(ϕ)]∂/∂ t-1/a^2∇^2+
    Υ_ϕ(ϕ)ϕ̇
        +V_ϕϕ(ϕ)}δφ   =   ξ_T.


With the slow-roll condition, we write  (<ref>) as form of that in cold inflation,

    3H(1+Q)ϕ̇+V,_ϕ(ϕ)=0,

where Q is the ratio of dissipation coefficient Γ and  Hubble parameter H i.e. Q≡Γ/3H. It is also necessary to define some slow-roll parameter for warm inflation,

    ε   =   1/16 π G(V,_ϕ/V)^2 ≪ 1+Q , 
    η   =   1/8 π GV,_ϕϕ/V≪ 1+Q ,

and

    β   =   1/8 π GΥ,_ϕ V,_ϕ/Υ V≪ 1+Q .


Now, looking again  (<ref>), we consider only the fluctuation from thermal noise and neglect that from quantum noise. Obviously, although thermal effect play a significant role in warm inflation, quantum effect may also be non-negligible. Quantum noise dominates  still for short wavelength perturbation, so we can use stochastic inflationary approach to deal the quantum fluctuation in warm inflation. Thus  (<ref>) in momentum space by defining the variable z=k/aH (ranging from 0 to ∞) becomes


    δφ”(k,z)-1/z(3Q+2)δφ'(k,z)+[1+3(η-β Q/(1+q))/z^2] δφ(k,z)=1/z^2H^2[ξ_q(k,z)+ξ_T(k,z)],


where primes denote the derivatives with respect to the variable z and ξ_q(k,z) is the quantum noise term reads

    ξ_q(k,z)=-H^2z^2/k^2[      ∂ ^2/∂ z^2-2+3Q/z∂/∂ z+1
          +Υ,_ϕϕ̇/H^2z^2+V,_ϕϕ/H^2z^2]ϕ̂_q(k,z),

in which

    ϕ̂_q(k,z)=W(k,z)[ϕ_k(z)â_-k+h.c.].

Here ϕ̂_q(k,z) are still the quantum field modes that satisfy the equation

    [∂ ^2/∂ z^2-2+3Q/z∂/∂ z+1+V,_ϕϕ/H^2z^2]ϕ_k(z)=0,

with the solution in terms of z:

    ϕ_k(z) = H√(π)/2k^3/2z^3/2H_μ^(1)(z),

where μ=√(9/4-V,_ϕϕ/H^2)≈ 3/2-3η. Then, we can get the expression of correlation function of quantum noise

    ⟨ξ_q(k,z)ξ_q(k',z')⟩ = (2π)^3z^2z'^2H^4/k^2k'^2[2n(k)+1]
    δ^3(k+k')Re[f_k(z)f_k'^*(z')],

where z'=k/a(τ')H and f_k(z) is <cit.>

    f_k(z)=k^2[W”-       3Q+2/zW'-3β Q/(1+Q)z^2W]ϕ_k(z) 
           +2k^2W'ϕ'_k(z),

and n(k) is the distribution function for high-frequency quantum inflaton which satisfies Bose-Einstein distribution, i.e. n(k)≡⟨â^†_k'â_-k⟩=1/[exp((k-μ)/aT)-1] where μ is the chemical potential to eliminate the divergence of several integrals which will be discussed below.



§ GRAVITATIONAL WAVE GENERATED FROM WARM INFLATION


Fluctuations (including thermal and quantum) through the energy-momentum tensors T_μν generates tensor perturbation in Friedmann-Robertson-Walker metrics as

    g_μν=a^2(τ)[dτ^2+(δ_ij+h_ij)dx^idx^j],

where dots denote the derivative with respect to conformal time τ and H_ij is the tensor perturbation with transverse-traceless gauge. Tensor perturbations h_ij satisfy the equation of motion

    ḧ_ij+2ȧ/aḣ_ij-∇^2 h_ij = 2/M^2_pΠ^  kl_ijT_kl,

with M^-2_p≡ 8π G. Define h̃_ij≡ a h_ij, then (<ref>) with Fourier transformation becomes

    ḧ̃̈_ij+(k^2-2/τ^2)h̃_ij = 2a/M^2_pΠ^  kl_ij(k)T_kl(k,τ) .

Here, we have used the relation τ=-1/aH in de Sitter space-time and Π^  kl_ij are the transverse-traceless projectors which follow <cit.>:

    Π^  kl_ij(k)p_kp_l   =   (p_i -k_i(p·k)/k^2)
    (p_j -k_j(p·k)/k^2) 
          -1/2(δ_ij-k̂_ik̂_j)(p^2 -(p·k)^2/k^2),


    Π^  kl_ij(k)Π^  mn_ij(k)p_kp_lp_mp_n=
    1/2(p^2 -(p·k)^2/k^2)^2,

and

    Π^  kl_ij(k)k_l=0,

where k̂_i=k_i/k. The solution of (<ref>) is

    h̃_ij(k,τ)=2/M_p^2∫dτ'G_k(τ,τ')a(τ')Π^  kl_ij(k)T_kl(k,τ'),

where G_k(τ,τ') is the Green function. Two linear solutions of (<ref>) are 


    y_1(τ)=-kτ n^(1)_1(-kτ)=(1+i/kτ)e^-ikτ

and

    y_2(τ)=-kτ n^(2)_1(-kτ)=(1-i/kτ)e^ikτ,


where n^(1)_1(z) and n^(2)_1(z) are the spherical Bessel functions of the third kind <cit.>.
According to the method in Appendix <ref>, Green's functions as a function of z and z' are

    G_k(z,z')=1/kzz'[(1+        zz')sin(z'-z)+(z-z') 
          ×cos(z'-z)]θ(z'-z),


Energy-momentum tensor for warn inflaton field write <cit.>

    T_ab=∂_a Φ∂_b Φ-g_ab(1/2∂_μΦ∂^μΦ-V).

The part on the right hand containing g_ab is projected away by Π^  kl_ij, and ϕ(t) is the space average of inflaton field that contributes nothing to the perturbation. So δφ are the primary part in T_ab. Thus, the tensor spectrum reads

    ⟨ h_ij( k,τ )h_ij( k^',τ )⟩
       =   4/(2π )^3a(τ )^2M_p^4∫dz^'a(τ ^')G_k(τ ,τ ^')∫dz^''a(τ ^'')G_k(τ ,τ ^'')
          ×∫d^3pd^3p^'Π _ij^  kl( k)Π _ij^ mn( k^')p_k(k_l-p_l)p_m^'(k_n^'-p_n^')
    
          ×⟨Φ ( p,τ ^')Φ ( k- p,τ
    ^')Φ ( p^',τ ^'')Φ ( k
    ^'- p^',τ ^'')⟩ .

It's convenient to write the spectrum as a function of z:

    ⟨ h_ij(k,z)h_ij(k',z_k')⟩
       =   4/(2π)^3 k^4 M_p^4∫dz'G̃_k(τ,τ')∫dz”G̃_k(τ,τ”)
          ×∫d^3pd^3p'Π^  kl_ij(k)Π^  mn_ij(k')p_k p_l p'_m p'_n  
          ×⟨δφ(p,τ')δφ(k-p,τ')δφ(p',τ”)δφ(k'-p',τ”)⟩,

where the relation z=k/aH have been used and G̃ are defined as

    G̃_k(z,z') = 1/z'^2[(1-zz')sin(z-z')+(z-z') 
    cos(z'-z)]θ(z'-z).

Now, one should note that z'_p are defined as z'_p=p/a(τ')H. If we omit the subscript like z”, it represents a variable in terms of wave number k,i.e. z”=k/a(τ”)H.

Using Wick theory, the average of ⟨...⟩ becomes

    ⟨δφ(p,τ')δφ(k-p,τ')
    δφ(p',τ”)δφ(k'-p',τ”)⟩
          =⟨δφ(p,τ')δφ(k-p,τ')⟩⟨δφ(p',τ”)δφ(k'-p',τ”)⟩
          +⟨δφ(p,τ')δφ(p',τ”)⟩⟨δφ(k-p,τ')δφ(k'-p',τ”)⟩
          +⟨δφ(p,τ')δφ(k'-p',τ”)⟩⟨δφ(p',τ”)δφ(p',τ”)⟩.

The first term on the right hand can be ignored because it is a diagram containing the disconnected term proportional to δ(k)δ(k'). And it is not hard to find that the second term and the third term are equivalent to each other.

Now the first priority is to solve Eq. (<ref>). According to (<ref>) and Appendix <ref>, the solution of
(<ref>) is

    δφ(k,z)=∫dz'g_k(z,z')1/z'^2H^2[ξ_T(k,z')+ξ_q(k,z')],

where

    g_k(z,z')   =   z^ν z'^ν/z'^2ν(2/π z')[J_α(z)Y_α(z')
                -J_α(z')Y_α(z)]θ(z'-z),

with

    ν   =   3(1+Q)/2, 
    α   =   √(ν^2-3β Q/1+Q-3η).

Write perturbed field as

    δφ=δφ_T+δφ_q,

where δφ_T is the part including integral of thermal noise ξ_T while δφ_q for quantum noise. Then we can find that tensor spectrum contains three components: thermal term ⟨δφ_Tδφ_T⟩⟨δφ_Tδφ_T⟩, quantum term ⟨δφ_qδφ_q⟩⟨δφ_qδφ_q⟩ and cross term ⟨δφ_qδφ_q⟩⟨δφ_Tδφ_T⟩. Next we will calculate the spectrum in terms of three components above.



 §.§ Thermal term

Thermal term of tensor spectrum at the end of inflation (τ=0, z→ 0) reads

    ⟨ h_ij(k)h_ij(k')⟩_T 
       =   4/(2π)^3 k^4 M_p^4∫dz'G̃_k(τ,τ')∫dz”G̃_k(τ,τ”)
          ×∫d^3pd^3p'Π^  kl_ij(k)Π^  mn_ij(k')p_k p_l p'_m p'_n  
          ×⟨δφ_T(p,τ')δφ_T(k-p,τ')⟩⟨δφ_T(p',τ”)δφ_T(k'-p',τ”)⟩,

The fluctuation-dissipation relation of (<ref>) acts in terms of t. Now using t=H^-1ln(k/Hz) together with

    δ(f(x))=∑_{x_0}δ(x-x_0)/|f'(x_0)|,

where x_0 are zero point of f(x), we obtain

    ⟨ξ_T(k,t)ξ_T(k',t')⟩=2(2π^3)Υ T H^4/k^2 k' z^3z' δ^3(k+k')δ(z-z').


It is helpful to make such calculation below:

    ⟨δφ_T(p,z'_p)δφ_T(p',z”_p')⟩
       =   π^2/4H^2∫_z'_p^∞dz1∫_z”_p'^∞dz_2(z'_p z”_p')^ν (z_1z_2)^-1-ν
          ×g̃_p(z'_p,z_1)g̃_p'(z”_p',z_2)⟨ξ_T(p,z_1)ξ_T(p',z_2)⟩
       ≃   π^5Υ T δ^3(k+k')/p^3[z'^ν_p Y_α(z'_p)][z”^ν_p Y_α(z'_p')]  
          ×∫_z'_p^∞dz_1 z^2-2ν_1 J^2_α(z_1),

where we have used (<ref>) and z^ν J_α(z)∼0 (z<1). Considering the function z^2-2νJ_α^2(z), we find that this function almost equals to zero
except a narrow peak at z≳1, so it is convenient to treat it as delta function. Thus

    ∫_z'_p^∞dz_1 z^2-2ν_1 J^2_α(z_1)
          ≃θ(1-p̅z')∫_0^∞dz_1 z^2-2ν_1 J^2_α(z_1)
          =θ(1-p̅z')Γ(ν-1)Γ(α-ν+3/2)/2√(π)Γ(ν-1/2)Γ(α+ν-1/2),

where z_p'=p/a(τ')H=(p/k)(k/a(τ')H)=p̅z' and (<ref>) has been used. With (<ref>) and (<ref>), get α≈ν. Thus

    ⟨δφ_T(p,z'_p)δφ_T(p',z”_p')⟩
       =   π^5 Υ T/p^3δ^3(p+p')(2^νΓ(ν))^2  
       ×   Γ(ν-1)Γ(3/2)/2√(π)Γ(ν-1/2)Γ(2ν-1/2).


The spectrum damps out as (sin kτ-kτcos kτ)/k^3 in Green's functions at large value of k. On the other hand, k only appear in the sublimit of the integral and then they are absorbed in θ function, and k become not so important which act only as the form of k-p. Besides,  and the most important, (<ref>) indicates that if p<ak_F, δφ(p,t_0) has not thermalized during inflation. So the integral in tensor spectrum gets its main contributions at p≫ k. With these approximations,

    ⟨ h_ij(k)h_ij(k')⟩_T 
       =   π^3 Υ^2 T^2δ^3(k+k')/16k^4 M_p^4{[2^νΓ(ν)]^2Γ(ν-1)/4Γ(ν-1/2)Γ(2ν-1/2)}^2 
          ×∫_0^∞dz'sin z'-z'cos z'/z'^2∫_z'^∞dz'sin z”-z”cos z”/z”^2
          ×∫d^3pd^3p'(p^2-(p·k)^2/k^2)^2δ^3(p+p')/p^6θ(1-p̅z').

Using Eqs. (<ref>), (<ref>), and (<ref>) ∼ (<ref>), thermal term of tensor spectrum (<ref>) is finally simplified to


    ⟨ h_ij(k)h_ij(k')⟩_T=3 π^4 H^4/10k^3 M_p^4δ^3(k+k')
    (T/H)^2  
          ×[Q8^QΓ(3Q/2+3/2)^3/(3Q+1)Γ(3Q/2+1)Γ(3Q+5/2)]^2.


Here, we have assumed z_2>z_1  and this assumption has no affect on the final result.



 §.§ Quantum term


According to (<ref>), the correlation function of δφ_q(k,z) is

    ⟨δφ_q(p,z'_p)δφ_q(p',z”_p')⟩
          =∫_z'_p^∞dz_1∫_z”_p'^∞dz_2 g(z'_p,z_1)g(z”_p',z_2)
          ×1/H^4(z_1z_2)^2⟨ξ_q(p,z_1)ξ_q(p',z_2)⟩
          =(2π)^3/p^4[2n(k)+1]δ^3(p+p')F_p(z_p')F^*_p(z_p”),

where

    F_p(z_p)=∫_z_p^∞dzG(z_p,z)f_p(z).

Using (<ref>), (<ref>), (<ref>) and (<ref>), F_p(z_p) reads

    F_p(z_p)≃-iz_p^3/2-μH/√(2)p^3,

with 0<z<ϵ<1 <cit.>. Notice that z_p=p/aH=p̅z<ϵ and z>ϵ, so 0<p̅<1. Thus the
quantum term of two-point function can be written as

    ⟨ h_ij(k)h_ij(k')⟩_q
    ≃32π^3H^4/5k^3M_p^4I(H/T)δ(k+k'),

where

    I(H/T)=∫_0^1dp̅[∫_0^∞dzsin z-zcos z/z^2.   
          ×.((p̅-p̅_0)zH/2T)]^2.


Next, we will proof that I(H/T) is normalized no matter at high temperature or low temperature. Note that p=(p̅-p̅_0)H/T and consider the low temperature condition (T≪ H or p≫ 1):

    ∫_0^∞sin x-xcos x/x^2(px/2) 
       =   ∫_0^∞dxsin x-xcos x/x^2+∫_0^∞dx/e^px-1sin x-xcos x/x^2
       =   1+2∑_n=1^∞(-1)^n-12n/(2n+1)!∫_0^∞x^2n-1/e^px-1dx
       =   1-2∑_n=1^∞ζ(2n)/2n+1(ip)^-2n
       =   1-ln(π/ip/sin(π/ip))+ip∫_-1/ip^1/iplnΓ(1+z)dz 
       =   1-ln(π/p/sinh(π/p))+2lnΓ(1+z)|_z→ 0
       ≃    1.

In the calculations above, we have used relevant properties of Gamma function and Zeta function in Appendix <ref>. We obtain I is unit at low temperature. Next consider the high temperature condition (T≫ H) and define z=exp(p_0/aT), p̃=p̅H/T. High temperature means that p_0 is a negative number with sufficient large value i.e. z≪ 1, thus

    ∫_0^∞sin x-xcos x/x^2((p̅-p̅_0)xH/2T) 
       =   ∫_0^∞dxsin x-xcos x/x^2+∫_0^∞dxze^-p̃x/1-ze^-p̃xsin x-xcos x/x^2
       =   1+2∑_n,m=1^∞(-1)^n-12n/(2n+1)!∫_0^∞x^2n-1(ze^-p̃x)^mdx
       =   1+2∑_n,m=1^∞(-1)^n-1/2n+11/p̃^2nz^m/m^2n
       =   1+2∑_n=1^∞(-1)^n-1/2n+11/p̃^2ng_2m(z)
       ≃    1.


Thus we assume that I(T/H)=1. This assumption is reliable. In stochastic approach, we assume a field ϕ_q as the vacuum fluctuation for
short wavelength with a sharp momentum cutoff, which means  such a fluctuation as a noise always exists no matter at low temperature or high
temperature. It is just the condition for cold inflation. On the other hand, chemical potential is a physical variable that cannot be ignored
when a particles coupling with other fields (like thermal bath) especially for the condition at phase transition. Most importantly, it is just
the chemical potential that eliminates the singularity at p=0.



 §.§ Cross term


If we simulate the calculation above, it is not hard to get the expression of the cross term for the two-point correlation function:

    ⟨ h_ij(k)h_ij(k')⟩_cross≃3π^4H^4/5k^3M_p^4δ(k+k') 
          ×T/HQ8^QΓ(3Q/2+3/2)^3/(3Q+1)Γ(3Q/2+1)Γ(3Q+5/2).

This component has nothing new compared with thermal and quantum component.

In the calculations above, we have ignored the slow parameter as the index of variable z. With this approximation, we obtain the scale-invariant spectrum of each components.



§ RESULT AND NUMERICAL ANALYSIS


Power spectrum of tensor perturbation P_h(k) is defined as:

    ⟨ h_ij(k)h_ij(k')⟩=2π^2/k^3P_h(k)δ(k+k').

The spectrum of tensor modes with vacuum form reads <cit.>

    P_h,vac(k)=2/π^2H^2/M_p^2.

We illustrate the power spectrum for thermal component in Fig. <ref> and total spectrum of tensor perturbation for warm inflation in
Fig. <ref>.




As discussed above, the spectrum for quantum component is  a constant with order of H^4/M_p^4. This value is so small that there is almost no correction to the spectrum of primordial tensors, which agrees with other cold inflation models <cit.>. However, what we interest most is the thermal component in total spectrum. As is shown in (<ref>), there are two variables determining the thermal spectrum i.e. Υ/H and T/H. Considering the weak dissipation condition with Υ/H≪ 1 and T/H≪ 1 first, we can directly get P_h,T(k)→ 0. In other words, this condition is just the same with cold inflation. In contrast, the condition with strong dissipation (Υ>H and T>H) is quite different from that of weak dissipation. The amplitude of thermal spectrum increases with the increase of Υ/H and T/H. Taking two simple examples, P_h,T(k)=63.87 at Υ/H=10 and T/H=10, or especially, P_h,T(k)∼ 10^4 at
Υ/H=100 and T/H=100. It has been widely studied that there has a significant correction on the spectrum of primordial
fluctuations in the warm regime with strong dissipation for T>H and such a correction lowers the tensor-to-scalar ratio. Thus the curvature
power spectrum is modified into <cit.>

    Δ^2_R=(H/ϕ̇)^2(H/2π)^2[1+2n+2π QT/H].

The temperature at the end of inflation is <cit.>

    4.09× 10^13GeV≤ T≤ 2.216× 10^14GeV.

With this and using (<ref>), (<ref>) and (<ref>), we obtain the tensor-to-scalar ratio in warm inflation:

    r≃[10^-8f(Q)+8]|n_T|/1+2n+2π QT/H,

where f(Q) is function containing Q in (<ref>) and n_T is the tensor index. From (<ref>) we can find that large Q will upper the tensor-to-scalar ratio , but it does not act obviously until Q reaches the level of 100 (Q>100). In
<cit.>, the authors have tabled the parameters like Q, T/H, r and others in warm inflation with different potential
by the constraint by recent years' observational data. The result shows T/H∼1 during inflation which agrees with the analysis in
<cit.>. The nondetection of cosmic gravitational waves background strongly constrains that thermal component of tensor
perturbation P_h,T will not upper the ratio dramatically which also agrees with a various of works
<cit.>. In this way, the primordial tensor spectrum thus can be used to distinguish warm inflation from
cold inflation.




§ CONCLUSIONS AND DISCUSSION


In this paper we first derive the two point correlation function of tensor perturbation in warm inflation and prove that primordial scale-invariant power spectrum of warm inflation is achieved only if system evolutes near thermal equilibrium state, in which condition warm inflaton fields are stable. We also explain the physical meaning of both scalar index and slow-roll parameter β in (<ref>) that correspond to the nonequilibrium properties during inflationary epoch. Then we calculate the power spectrum for warm inflation by Green's function method and mainly discuss the thermal component in total spectrum. We consider it as a new method to distinguish warm inflation from cold inflation. At last, we illustrate our result by numerical analysis. Using existing cosmic observational data, we find that fluctuations from thermal noise does not raise the tensor-to-scalar ratio dramatically at order of 10^-4 although the temperature is high enough.

There are also many issues which deserve further discussion in this paper. The first one is the thermal properties of warm inflation model. For example, how initial condition determines the scalar index n_s and whether the stability of thermal inflaton fields still holds with arbitrary initial condition? Then, we find that the spectrum of quantum fluctuation is almost a constant because of momentum cutoff in stochastic approach. So whether quantum fluctuations still exist in high temperature in other models becomes an interesting and challenging work. Finally, but the most important, is whether we can observe such a fluctuation? Although the cosmic gravitational waves background has not been detected yet, with the
discovery of gravitational waves, more and more new method coming out <cit.> and a series of
observation installations working or being build <cit.>, it is believable that primordial
gravitational waves can be detected in the foreseeable future.



This work was supported by the National Natural Science Foundation of China (Grants No. 11575270,
No. 11175019, and No. 11235003). The authors would like to thank Zong-Kuan Guo and Zong-Hong Zhu for helpful comments and discussions.







§ THERMAL PROPERTIES OF WARM INFLATION


Correlation function is defined as

    𝒫_δφ(x-y,t_1,t_2)=⟨δφ(x,t_1)δφ(y,t_2)⟩,

whose Fourier transformation is

    𝒫_δφ(k,t_1,t_2)=∫d^3k'/(2π)^3⟨δφ(k,t_1)
    δφ(k',t_2)⟩,


Let's consider (<ref>) again. Using the slow-roll condition and strong dissipation condition Γ≫ H, (<ref>) approximately writes  <cit.>

    Υdδφ(k,t)/dt+[k^2_p+V”(ϕ)]δφ(k,t) ≈ξ_T(k,t),

where V”(ϕ)=d^2V(ϕ)/dϕ^2 and ϕ is defined in  (<ref>), k_p=k/a is the physical wave number and k is the conformal wave number. Strong dissipation means that we can ignore the change of parameter a, k_p and T within the time interval 1/H and fluctuation from quantum noise is negligible. Then the solution of  (<ref>) is

    δφ(k,t)   ≈   1/Υe^-(t-t_0)/τ(ϕ)∫^t_t_0e^t'-t_0ξ_t(k,t')dt'
          +δφ(k,t_0)e^-(t-t_0)/τ(ϕ)

where t_0 is any coordinate time during inflation and τ(ϕ)=Υ/[k^2_p+V”(ϕ)].

The correlation function of perturbed inflation field is


    ⟨δφ(k,t_1)δφ(k',t_2)⟩   =   δφ(k,t_0)δφ(k',t_0)
    e^-(t_1+t_2)/τ(ϕ)+2(2π)^3Tδ^3(k+k')/Υ a^3e^-(t_1+t_2)/τ(ϕ)
            ×∫^t_1_t_0∫^t_2_t_0e^(s_1+s_2)/τ(ϕ)δ(s_1-s_2)ds_1 ds_2

where we have used  (<ref>). The double integral in  (<ref>) contains a δ function, so we need to integrate first to the lager one in t_1 and t_2. Then

    ⟨δφ(k,t_1)δφ(k',t_2)⟩   =   δφ(k,t_0)δφ(k',t_0)
    e^-(t_1+t_2)/τ(ϕ)+2(2π)^3Tδ^3(k+k')/Υ a^3e^-(t_1+t_2)/τ(ϕ)
          ×∫^min(t_1,t_2)_t_0∫^max(t_1,t_2)_t_0e^(s_1+s_2)/τ(ϕ)δ(s_1-s_2)ds_1 ds_2
       =   δφ(k,t_0)δφ(k',t_0)e^-(t_1+t_2)/τ(ϕ)+2(2π)^3Tδ^3(k+k')/Υ a^3e^-(t_1+t_2)/τ(ϕ)∫^min(t_1,t_2)_t_0e^2s/τ(ϕ)ds
       =   δφ(k,t_0)δφ(k',t_0)e^-(t_1+t_2)/τ(ϕ)+
    (2π)^3Tτ(ϕ)δ(k+k')/Υ a^3{e^-[t_1+t_2-2min(t_1,t_2)]/τ(ϕ)-e^-(t_1+t_2)/τ(ϕ)}
       =   δφ(k,t_0)δφ(k',t_0)e^-(t_1+t_2)/τ(ϕ)+
    (2π)^3Tτ(ϕ)δ^3(k+k')/Υ a^3{e^-|t_1-t_2|/τ(ϕ)-e^-(t_1+t_2)/τ(ϕ)}.


Now, let's make more detailed calculation of the second term in final equation.

    (2π)^3Tτ(ϕ)/Υ a^3
    =(2π)^3TΥ/a^3Υ(k^2/a^2+V”)  
    =(2π)^3HT/a^3H^3(k^2/a^2H^2+V”/H^2)
    = (2π)^3 HT/k^3/z^3(z^2+3η)≃(2π)^3 HTz/k^3,

where we have used relation z=k/aH and Eq. (<ref>). Finally, the correlation function reads

    ⟨δφ(k,t_1)δφ(k',t_2)⟩=(δφ(k,t_0)δφ(k',t_0)
    -(2π)^3 HTz/k^3δ^3(k+k'))e^-(t_1+t_2)/τ(ϕ)+(2π)^3 HTz/k^3δ^3(k+k')e^-|t_1-t_2|/τ(ϕ).



Thermalization requires that initial condition must thermalize at physical scale k_p within 1/H. From  (<ref>), the first term in right hand must be negligible with in Hubble time, so

    k^2_p-V”/HΥ>1.

If V”>HΥ, the freeze-out number is the same with that of cold inflation, k_F=H (z=1). However, if V”<HΥ, the freeze-out number reads

    k_F=(HΥ)^1/2.


Obviously 𝒫(k,t_1,t_2) is dependent on the initial state δφ(k,t_0). Now make average on initial state
and assume t_1=t_2=t together with t_0=0, so correlation function can be written as

    𝒫_δφ(k',t_0)=(2π^2/k^3P_δφ(k',t_0)-HTz_*/k^3)
    e^-2t/τϕ
    +HTz_*/k^3,

with the definition of power spectrum

    P_δφ(k,t)=k^3/2π^2∫d^3k'/(2π)^3⟨δφ(k,t_0)
    δφ(k',t_0)⟩,

and z_* as the value at horizon crossing z_*=(Υ H)^1/2/H.
Things become quite interesting. If the power spectrum is scale-invariant i.e.  𝒫_δφ(k,t_0)=(Υ H)^1/2
T/2π^2 <cit.>, the spectrum P_δφ(k,t) is also scale-invariant and totally the same with
𝒫_δφ(k,t_0). In other words, P_δφ is stable in this condition which means system is on thermal
equilibrium during time interval t>t_0. This is quite similar with the condition of correlation function for particles with Brown motion
<cit.>. Now, consider that P_δφ(k,t_0) is not a scale-invariant spectrum and assume
P_δφ(k,t_0)=A(k/k_0)^n'-1 where n' is an arbitrary number. Then P_δφ(k,t)becomes also
dependent on k i.e. P_δφ(k,t)=A(k/k_0)^n-1. However, this term damps with the increase of time, which means n
tends to unit with the time evolution. This is an effect dominated by nonequilibrium mechanics. In this way, we can say that scalar index n_s
and slow-roll parameter β are parameters that illustrate the deviation from equilibrium state. Probe on Cosmic Microwave Background shows
that our universe is almost in nearly equilibrium <cit.> if considering our universe in early epoch as a model in thermal bath. The
relation (<ref>) indicates that P_δφ(k,t_0) (the initial condition of universe) becomes not so
important even though we cannot give an accurate description till now.




§ A BRIEF INTRODUCTION TO GREEN'S FUNCTION 


Let's consider the equation <cit.>


    d^2 y/d t^2+ω^2(t)y=f(t),

with boundary condition

    y(0)=y'(0)=0.


To solve this equation, we need to solve another relevant equation in terms of Green function G(t,t'):


    d^2 G(t,t')/d t^2+ω^2(t)G(t,t')=δ(t-t'),

with boundary condition

    G|_t=0=d G(t,t')/dt|_x=0=0.


Obviously, when t<t', G satisfy

    d^2 G(t,t')/d t^2+ω^2(t)G(t,t')=0,

with boundary condition (<ref>). Then assume the general solution of (<ref>) is

    G(t<t')=c_1(t')y_1(t)+c_2(t')y_2(t),

where y_1 and y_2 is the solutions of (<ref>). Considering the boundary condition, we obtain

    G(t<t')=0.

Next, assume

    G(t>t')=c_3(t')y_1(t)+c_4(t')y_2(t).

It's not hard to find that G is continuous at t=t' while G' not, and

    d G(t,t')/d t|_t<t'=0,  d G(t,t')/d t|^t+0_t-0=1.

Thus

    c_3(t')y_1(t')+c_4(t')y_2(t')=0 
    
    c_3(t')y'_1(t')+c_4(t')y'_2(t')=1

Then we obtain the solution of (<ref>)

    G(t,t')=y_1(t)y_2(t')-y_2(t)y_1(t')/y_1(t)y'_2(t)-y_2(t)y'_1(t)θ(t-t').

Finally, the solution of (<ref>) is

    y(t)=∫_0^∞G(t,t')f(t')dt'.




§ SPECIAL FUNCTIONS AND INTEGRALS

The equation

    d^2y/dz^2+1/zdy/dz+(1-ν^2/z^2)y=0

has two linear independent solutions: Bessel function J_ν(z) and Neumann function Y_ν(z), with properties
<cit.>

    J_ν(z) ≈1/Γ(ν+1)(z/2)^ν (ν>0, z→ 0^+),


    Y_ν(z) ≈ -Γ(ν)/π(2/z)^ν (ν>0, z→ 0^+),

and

    J_ν(z)Y'_ν(z)-J'_ν(z)Y_ν(z)=2/π z.

The solution of equation


    d^2u/dz^2+1-2α/zdu/dz+(β^2+α^2-ν^2/z^2)u = 0

is also in terms of Bessel function

    u=z^α Z_ν(β z),


where Z_ν is any kind of Bessel function. Integral with double Bessel function is Schafgeitlin integral formula:

    ∫_0^∞      J_μ(ax)J_ν(bx)/x^λdx=b^νΓ(ν+μ-λ+1/2)/2^λ a^ν-λ+1Γ(ν+1)Γ(μ-ν+λ+1/2)
          × F(ν+μ-λ+1/2,ν-μ-λ+1/2;ν+1;b^2/a^2),
              μ+ν+1>λ>-1,0<b<a  or
              μ+ν+1>λ>0, a=b .

Specially,

    ∫_0^∞dz_1 z^2-2ν_1 J^2_α(z_1)
           = (1/2)^2ν-3Γ(2ν-2)Γ(α-ν-3/2)/2[Γ(ν-1/2)]^2Γ(α+ν-1/2)
           = Γ(ν-1)Γ(α-ν+3/2)/2√(π)Γ(ν-1/2)Γ(α+ν-1/2),

where we have used (<ref>). Zeta function is defined as

    ζ(z)=1/Γ(z)∫_0^∞x^z-1/e^x-1dx,

where Γ(z) are Gamma functions with

    Γ(z+1)=      zΓ(z) (Re(z)>1), 
    Γ(2z)=2^2z-1      π^-1/2Γ(z)Γ(z+1/2).

Zeta function has the property <cit.>

    ∑_k=1^∞ζ(2k)z^2k/2k+1   =   1/2ln(π z/sinπ z) 
          -1/2z∫_-z^zlnΓ(1+z)dz,   |z|<1.

We have used the two integral formulas below frequently in this paper

    ∫_0^∞dxsin x-xcos x/x^2=1,
    ∫_0^∞dxsin x-xcos x/x^3=π/4

and

    ∫_0^∞sin x-xcos x/x^3∫_0^xsin y-ycos y/y^2dydx=π/12





References

PhysRevX.6.041015 B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. X 6, 041015 (2016).

VIRGO <http://www.ego-gw.it/virgodescription/pag_4.html>.

GEO600<http://www.geo600.org/>.

TAMA300howpublished = <http://tamago.mtk.nao.ac.jp/>.

1975ApJ...195L..51HR. A. Hulse and J. H. Taylor, Astrophys. J., 195, L51 (1975).

Abbott2016 B. P. Abbott et al. , Living Reviews in Relativity 19, 1 (2016).

PhysRevLett.116.061102B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett.
    116, 061102 (2016).

PhysRevLett.116.241103B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett.
    116, 241103 (2016).

doi:10.1142/S0217751X09044206M. Bastero-Gil and A. Berera, Int. J. Mod. Phys. A 24, 2207 (2009).

BARTRUM2014116 A. Berera and R. O. Ramos, Phys. Lett. B 732, 116 (2014).

PhysRevD.60.083509J. Yokoyama and A. D. Linde, Phys.Rev. D 60, 083509 (1999).

PhysRevLett.83.264A. Berera, M. Gleiser, and R. O. Ramos, Phys. Rev. Lett. 83, 264 (1999).

PhysRevD.84.103503M. Bastero-Gil, A. Berera, and J. G. Rosa, Phys. Rev. D 84, 103503 (2011).

PhysRevLett.117.151301M. Bastero-Gil, A. Berera, R. O. Ramos, and J. G.Rosa, Phys. Rev. Lett. 117, 151301 (2016).

1475-7516-2011-09-033M. Bastero-Gil, A. Berera, and R. O. Ramos, J. Cosmol. Astropart. Phys 09, 033 (2011).

978-0-521-87499-1J. Rammer, Quantum Field Theory of Non-equilibrium States, Cambridge University Press, Cambridge U.K. (2007).

PhysRevD.76.083520 A. Berera, I. G. Moss, and R. O. Ramos, Phys. Rev. D 76, 083520 (2007).

PhysRevD.91.083540S. Bartrum, A. Berera, and J. G. Rosa, Phys. Rev. D 91, 083540 (2015)

RevModPhys.78.537 B. A. Bassett, S. Tsujikawa, and D. Wands, Rev. Mod. Phys. 78, 537 (2006).

PhysRevD.69.083525L. M. H. Hall, I. G. Moss, and A. Berera, Phys. Rev. D 69, 083525 (2004).

PhysRevD.62.083517A. N. Taylor and A. Berera, Phys. Rev. D 62, 083517 (2000).

PhysRevD.64.063513 H. P. de Oliveira and S. E. Jorás, Phys. Rev. D 64, 063513 (2001).

PhysRevD.71.023513 A. Berera and R. O. Ramos, Phys. Rev. D 71, 023513 (2005).

PhysRevD.50.2441M. Gleiser and R. O. Ramos, Phys. Rev. D 50, 2441 (1994).

PhysRevD.58.123508 A. Berera, M. Gleiser, and R. O. Ramos, Phys. Rev. D 58, 123508 (1998).

PhysRevD.59.083503 W. Lee and L.-Z. Fang Phys. Rev. D 59, 083503 (1999).

PhysRevD.57.741 H. P. de Oliveira and R. O. Ramos, Phys. Rev. D 57, 741 (1998).

Zhang1 X. -M. Zhang and J. -Y. Zhu, Phys. Rev. D 90, 123519 (2014).

Zhang2 X. -M. Zhang and J. -Y. Zhu, Phys. Rev. D 91, 063510 (2015).

Planck Planck Collaboration, P. A. R. Ade et al., Astronomy and Astrophyscs 594, A13 (2016).

PhysRevD.76.084019 D. Baumann, P. Steinhardt, K. Takahashi, and K. Ichiki, Phys. Rev. D 76, 084019 (2007).

PhysRevD.75.123518 K. N. Ananda, C. Clarkson, and D. Wands, Phys. Rev. D 75, 123518 (2007).

PhysRevD.85.023534 J. L. Cook and L. Sorbo, Phys. Rev. D 85, 023534 (2012).

PhysRevD.88.103518 M. Biagetti, M. Fasiello, and A. Riotto, Phys. Rev. D 88, 103518 (2013).
 	
1475-7516-2016-02-023 S. Kuroyanagi, T. Hiramatsu, and J. Yokoyama, J. Cosmol. Astropart. Phys 02, 023 (2016).

2004JCAP...08..011L M. Liguori, S. Matarrese,M. A. Musso, A. Riotto, J. Cosmol. Astropart. Phys 08, 011 (2004).

PhysRevD.61.084008 S. Winitzki and A. Vilenkin, Phys. Rev. D 61, 084008 (2000).

1475-7516-2013-03-032 R. O. Ramos and L. da Silva, J. Cosmol. Astropart. Phys 03, 032 (2013).

gravitationalwaves A. Buonanno, arxiv:0709.4682.

Gravitational...Waves..1 M. Maggiore, Gravitational waves (Volume 1): Theory and experiments, Oxford University Press (2007).

1970hmfw.book.....A M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, 9th Edition, Dover Publications
    (New York, 2012).

2014JCAP...05..004BM. Bastero-Gil, A. Berera, I. G. Moss, and R. O. Ramos, J. Cosmol. Astropart. Phys 5, 004 (2014).

2003moco.book.....D S. Dodelson, Modern Cosmology, Elsevier Press (Singapore, 2008).

2014PhLB..732..116B S. Bartrum, M. Bastero-Gil, A. Berera, R. Cerezo, R. O. amos, and J. G. Rosa, Physics Letters B 732, 116 (2014)

2016JCAP...11..022GY. Gim and W. Kim, J. Cosmol. Astropart. Phys 11, 022 (2016),

2017PhRvD..95b3517B M. Biagetti, M. Fasiello, and A. Riotto, Phys. Rev. D 95, 023517 (2017).

PAGANO2016823 L. Pagano, L. Salvati, and A. Melchiorri, Physics Letters B 760, 823 (2016).

PhysRevLett.118.091102 X.-L. Fan, K. Liao, M. Biesiada, A. Piórkowska-Kurpas, and Z.-H. Zhu, Phys. Rev. Lett. 118,
        091102 (2017).
  	
1475-7516-2015-12-006 X. Ding, M. Biesiada, and Z.-H. Zhu, J. Cosmol. Astropart. Phys 12, 006 (2015).

2000NuPhB.585..666B A. Berera, Nuclear Physics B 585, 666 (2000).

2006stme.book.....S F. Schwabl, Statistical Mechanics (Advanced Texts in Physics), Springer (2006).

WMAP E. Komatsu et al.The Astrophysical Journal Supplement Series, 192:18 (47pp), 2011.

MathemicalMethed G. B. Arfken, H. J. Weber and F. E. Harris, Mathematical Methods for Physicists (Seventh Edition) ,
        Academic Press, (2012).

SpecialFunctions Z.X. Wang, and D.R Guo, Special Functions (World Scientific Press, 1989).

SRIVASTAVA1988129 H. Srivastava, Journal of Mathematical Analysis and Applications 134, 129 (1988).




