The tensor bi-spectrum in a matter bounce

The tensor bi-spectrum in a matter bounce

Debika Chowdhury,    V. Sreenath111Current address: Department of Physics and Astronomy, Louisiana State University, Baton Rouge, LA 70803, U. S. A.. E-mail: sreenath@lsu.edu.    and L. Sriramkumar
today
Abstract

Matter bounces are bouncing scenarios wherein the universe contracts as in a matter dominated phase at early times. Such scenarios are known to lead to a scale invariant spectrum of tensor perturbations, just as de Sitter inflation does. In this work, we examine if the tensor bi-spectrum can discriminate between the inflationary and the bouncing scenarios. Using the Maldacena formalism, we analytically evaluate the tensor bi-spectrum in a matter bounce for an arbitrary triangular configuration of the wavevectors. We show that, over scales of cosmological interest, the non-Gaussianity parameter that characterizes the amplitude of the tensor bi-spectrum is quite small when compared to the corresponding values in de Sitter inflation. During inflation, the amplitude of the tensor perturbations freeze on super-Hubble scales, a behavior that results in the so-called consistency condition relating the tensor bi-spectrum and the power spectrum in the squeezed limit. In contrast, in the bouncing scenarios, the amplitude of the tensor perturbations grow strongly as one approaches the bounce, which suggests that the consistency condition will not be valid in such situations. We explicitly show that the consistency relation is indeed violated in the matter bounce. We discuss the implications of the results.

The tensor bi-spectrum in a matter bounce

• Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India

1 Introduction

Bouncing models correspond to situations wherein the universe initially goes through a period of contraction until the scale factor reaches a certain minimum value before transiting to the expanding phase. They offer an alternative to inflation to overcome the horizon problem, as they permit well motivated, Minkowski-like initial conditions to be imposed on the perturbations at early times during the contracting phase (see, for instance, Refs. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]; for reviews, see Refs. [19, 20, 21]). Interestingly, certain bouncing scenarios can lead to nearly scale invariant perturbation spectra (see, for example, Refs. [3, 11, 18]), as is required by observations [22, 23]. For instance, a bouncing model wherein the universe goes through a contracting phase as driven by matter—referred to as a matter bounce—is known to lead to an exactly scale invariant spectrum of tensor perturbations as in de Sitter inflation [1, 2, 3]. Clearly, it will be worthwhile to examine if non-Gaussianities can help us discriminate between such scenarios [24, 25, 26].

The most dominant of the non-Gaussian signatures are the non-vanishing three-point functions involving the scalars as well as the tensors [27, 28, 29, 30, 31, 32, 33, 34]. In order to drive a bounce, it is well known that one requires matter fields that violate the null energy condition. Therefore, analyzing the evolution of the scalar perturbations require suitable modelling of the matter fields [20, 21, 19]. In contrast, the tensor perturbations depend only on the scale factor and hence are simpler to study. For this reason, we shall focus on the tensor bi-spectrum in this work. Further, we shall assume a specific functional form for the scale factor and we shall not attempt to construct sources that can give rise to such a behavior.

An interesting aspect of the three-point functions is their property in the so-called squeezed limit wherein the wavelength of one of the three modes involved is much larger than the other two [27, 35, 36, 37, 38, 39, 40, 41, 42, 43]. In this limit, under fairly general conditions, it is known that the three-point functions can be expressed completely in terms of the two-point functions, a relation that is referred to as the consistency condition. We should mention that, while the scalar consistency relation has drawn most of the attention, it has been established that all the four three-point functions involving scalars and tensors satisfy similar relations under certain conditions [28, 29, 30, 44, 45]. It is interesting to examine if the three-point functions generated in the bouncing scenarios satisfy the consistency condition. In the context of inflation, it is well known that the consistency relations arise due to the fact that the amplitude of the long wavelength mode freezes on super-Hubble scales. In contrast, in a bouncing universe it can be readily shown that the amplitude of the long wavelength mode grows sharply as one approaches the bounce during the contracting phase. This behavior suggests that the consistency relation may not hold in bouncing models [24]. The primordial consistency conditions lead to corresponding imprints on the anisotropies in the cosmic microwave background (see, for instance, Refs. [46, 47, 48]; in particular, see Ref. [49] for the signatures of the tensor modes) and the large scale structure (see, for example, Refs. [50, 51, 52]). It is clear that the consistency condition, if it can be confirmed by the observations, can help us discriminate between models of the early universe.

The most comprehensive formalism to study the generation of non-Gaussianities in the early universe is the approach due to Maldacena [27]. In this work, we analytically evaluate the tensor bi-spectrum in a matter bounce using the Maldacena formalism. To arrive at the tensor bi-spectrum analytically, one requires not only the behavior of the tensor modes, one also needs to be able to evaluate a certain integral involving the scale factor and the tensor modes. We conveniently divide the evolution into three domains and use the analytic solutions available in these domains to carry out the integrals and obtain the tensor bi-spectrum.

This paper is organized as follows. In the following section, considering a specific form for the scale factor, we shall divide the period before the bounce into two domains, the first corresponding to early times during the contracting phase and the other close to the bounce. We shall describe the analytic solutions to the tensor modes during these two domains and also discuss the behavior of the modes after the bounce to arrive at the corresponding tensor power spectrum over wavenumbers much smaller than the wavenumber associated with the bounce. We shall also compare the analytical solutions for the tensor modes with the corresponding results obtained numerically. In Sec. 3, we shall quickly summarize the essential expressions describing the tensor bi-spectrum in the Maldacena formalism. We shall also introduce the tensor non-Gaussianity parameter , a dimensionless quantity that reflects the amplitude of the tensor bi-spectrum. In Sec. 4, we shall evaluate the tensor bi-spectrum using the analytic solutions to the modes and the behavior of the scale factor in the three domains. We shall calculate the bi-spectrum for an arbitrary triangular configuration of the wavevectors. In Sec. 5, we shall illustrate the results in the equilateral and the squeezed limits. We shall show that the non-Gaussianity parameter that characterizes the tensor bi-spectrum is very small for cosmological scales and that the consistency relation is violated in the squeezed limit. We shall conclude with a brief discussion in Sec. 6. In an Appendix, we shall briefly outline a proof of the consistency condition satisfied by the tensor bi-spectrum during inflation.

Note that we shall work with natural units wherein , and define the Planck mass to be .

2 The tensor modes and the power spectrum

We shall consider the background to be the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric that is described by the line-element

 ds2=a2(η)(−dη2+δijdxidxj), (2.1)

where denotes the scale factor and is the conformal time coordinate. We shall assume that the scale factor describing the bouncing scenario is given in terms of the conformal time coordinate by the relation

 a(η)=a0(1+η2/η20)=a0(1+k20η2), (2.2)

where is the minimum value of the scale factor at the bounce (i.e. when ) and denotes the time scale that determines the duration of the bounce. Note that, at very early times, viz. when , the scale factor behaves as in a matter dominated universe (i.e. as ) and, for this reason, such a bouncing model is often referred to as the matter bounce. We shall assume that the scale associated with the bounce, viz. , is of the order of the Planck scale . Therefore, the wavenumbers of cosmological interest are orders of magnitude smaller than the wavenumber .

Upon taking into account the tensor perturbations characterized by , the spatially flat FLRW metric can be expressed as [27]

 ds2=a2(η)(−dη2+[eγ(η,x)]ijdxidxj). (2.3)

Recall that the primordial perturbations are generated due to quantum fluctuations. On quantization, the tensor perturbation can be decomposed in terms of the corresponding Fourier modes as follows:

 ^γij(η,x) = ∫d3k(2π)3/2^γkij(η)eik⋅x (2.4) = ∑s∫d3k(2π)3/2(^askεsij(k)hk(η)eik⋅x+^as†kεs∗ij(k)h∗k(η)e−ik⋅x).

In this decomposition, the pair of operators represent the annihilation and creation operators corresponding to the tensor modes associated with the wavevector , and they satisfy the standard commutation relations. The quantity represents the polarization tensor of the gravitational waves with their helicity being denoted by the index . The transverse and traceless nature of the gravitational waves leads to the conditions . We shall work with a normalization such that  [27]. The tensor power spectrum, viz.  , is defined as follows:

 ⟨^γkij(ηe)^γk′mn(ηe)⟩ = (2π)22k3Πkij,mn4PT(k)δ(3)(k+k′), (2.5)

where the expectation values on the left hand sides are to be evaluated in the specified initial quantum state of the perturbations, and denotes a suitably late conformal time when the power spectrum is to be evaluated. The quantity is given by [28, 29, 30, 44, 45]

 Πkij,mn=∑sεsij(k)εs∗mn(k). (2.6)

The tensor spectral index is defined as

 nT=dlnPT(k)dlnk. (2.7)

The tensor modes satisfy the differential equation

 h′′k+2a′ah′k+k2hk=0, (2.8)

where the overprimes denote differentiation with respect to the conformal time . If we write , then the modes satisfy the equation

 u′′k+(k2−a′′a)uk=0. (2.9)

The quantity corresponding to the scale factor (2.2) is given by

 a′′a=2k201+k20η2, (2.10)

which has essentially a Lorentzian profile. Note that the quantity exhibits a maximum at the bounce, with the maximum value being of the order of . Also, it goes to zero as . For modes of cosmological interest, one finds that at suitably early times (i.e. as ). In this domain, the quantity oscillates and we can impose the standard initial conditions on these modes and study their evolution thereafter.

Let us divide the period before the bounce into two domains, a domain corresponding to early times and another closer to the bounce. Let the first domain be determined by the condition , where is a relatively large number, which we shall set to be, say, . The second domain evidently corresponds to . In the first domain, we can assume that the scale factor behaves as

 a(η)≃a0k20η2, (2.11)

so that . Since the condition corresponds to, say, , the initial conditions can be imposed when . The modes can be easily obtained in such a case and the positive frequency modes that correspond to the vacuum state at early times are given by [1, 2, 3]

 hk=√2MPl1√2k1a0k20η2(1−ikη)e−ikη. (2.12)

Let us now consider the behavior of the modes in the second domain, i.e. when . Since we are interested in scales much smaller than , we shall assume that , which corresponds to the condition . Therefore, in this domain, for scales of cosmological interest, the equation governing the tensor mode reduces to

 h′′k+2a′ah′k≃0. (2.13)

This equation can be immediately integrated to yield

 h′k(η)≃h′k(η∗)a2(η∗)a2(η), (2.14)

where is a suitably chosen time and the scale factor is given by the complete expression (2.2). On further integration, we obtain that

 hk(η)=hk(η∗)+h′k(η∗)a2(η∗)∫ηη∗d~ηa2(~η), (2.15)

where we have chosen the constant of integration to be . If we choose , we can make use of the solution (2.12) to determine and . Note that, in the domain of interest, the first term in the above expression is, evidently, a constant, while the second term grows rapidly as one approaches the bounce. Upon using the form (2.2) of the scale factor, we find that we can express the behavior of the mode in the second domain as

 hk=Ak+Bkf(k0η), (2.16)

where

 f(k0η)=k0η1+k20η2+tan−1(k0η), (2.17)

while the quantities and are given by

 Ak = √2MPl1√2k1a0α2(1+ik0αk)eiαk/k0+Bkf(α), (2.18) Bk = √2MPl1√2k12a0α2(1+α2)2(3ik0α2k+3α−ikk0)eiαk/k0. (2.19)

Let us now turn to the third domain, i.e. immediately after the bounce. In this case too, for modes such that , the solution to is given by Eq. (2.16). We should highlight the fact that, whereas the bounce (2.2) is a symmetric one, the solution (2.16) is asymmetric in . Moreover, one may have naively expected the amplitude of the long wavelength modes to freeze once the universe starts expanding. This is largely true, though not completely so, and the behavior can possibly be attributed to the specific form of the scale factor (2.2). Note that, during this domain, while the first term in decays, the second term actually grows, albeit extremely mildly. We shall assume that, after the bounce, the universe transits to the conventional radiation domination epoch at, say, , where we shall set . We should mention that this choice is somewhat arbitrary and we shall discuss the dependence of the tensor power spectrum and the bi-spectrum on in due course.

In order to understand the extent of accuracy of the approximations involved, it would be worthwhile to compare the above analytical results for the mode with the corresponding numerical results. Clearly, given the scale factor (2.2), it is a matter of integrating the differential equation (2.8), along with the standard Bunch-Davies initial conditions, to arrive at the behavior of . The conformal time coordinate does not prove to be an efficient time variable for numerical integration, particularly when a large range in the scale factor needs to be covered. In the context of inflation, it is often the e-fold , defined as , where is a suitable time at which the scale factor takes the value , that is utilized to integrate the equations of motion (see, for instance, Refs. [53, 54, 55, 56]). Due to the exponential factor involved, a small range in e-folds covers a large range in time and scale factor. However, since is a monotonically increasing function, while it is useful to describe expanding universes, e-folds are not helpful in characterizing bouncing scenarios. In order to characterize a bounce, it would be convenient to choose a variable that is negative during the contracting phase of the universe, zero at the bounce and positive during the expanding phase. We shall choose to perform the integration using a new variable , which we call the e-N-fold, in terms of which the scale factor is defined as  [57]. We shall assume that is zero at the bounce, with negative values representing the phase prior to the bounce and positive values after.

In terms of the e-N-fold, the differential equation (2.8) governing the evolution of the tensor modes can be expressed as

 d2hkdN2+(3N+1HdHdN−1N)dhkdN+(kNaH)2hk=0, (2.20)

where is the Hubble parameter. Given the scale factor (2.2), the corresponding Hubble parameter can be easily evaluated in terms of the conformal time . In order to express the Hubble parameter in terms of the e-N-fold, we shall require as a function of . Upon using the definition of the e-N-folds and the expression (2.2) for the scale factor, we obtain that

 η(N)=±k−10(eN2/2−1)1/2. (2.21)

Since the Hubble parameter is negative during the contracting phase and positive during the expanding regime, we have to choose the root of accordingly during each phase. From the expression for the Hubble parameter, we evaluate the coefficients of the differential equation (2.20) in terms of . With the coefficients at hand, we numerically integrate the differential equation using a fifth order Runge-Kutta algorithm. We should mention that we have also independently checked the numerical results using Mathematica. Recall that, the initial conditions are imposed in a domain during the contracting phase wherein . As is done in the context of inflation, we shall impose the initial conditions when corresponding to, say, the e-N-fold . It should be pointed out that the initial conditions on the different modes are imposed at different times. In terms of the e-N-folds, the standard Bunch-Davies initial conditions can be expressed as

 hk = 1a(Ni)√2k, (2.22a) dhkdN = −iNia2(Ni)H(Ni)√k2−Nia(Ni)√2k. (2.22b)

We impose these initial conditions well before the bounce and evolve the modes until a suitable time after the bounce. The tensor mode evaluated in such a manner has been plotted for a given wavenumber (such that ) in Fig. 1. The figure also contains a plot of the analytical results (2.12) and (2.16) for the same wavenumber. As is evident from the figure, prior to the bounce and immediately after, the analytical results match the exact numerical results exceedingly well.

The tensor power spectrum after the bounce can be calculated using the solutions we have obtained. Recall that the tensor power spectrum is defined as

 PT(k)=4k32π2|hk(η)|2, (2.23)

with the spectrum to be evaluated at a suitable time. If we evaluate the tensor power spectrum at , we find that it can be expressed as

 PT(k)=4k32π2|Ak+Bkf(β)|2. (2.24)

Note that, is a quantity that we have artificially introduced and the actual problem does not contain . For and a sufficiently large (as we had said, for or so), the above power spectrum reduces to a scale invariant form with a weak dependence on , if is reasonably larger than unity. If we further assume that is large enough (say, ), then the scale invariant amplitude is found to be: , as expected [1, 2, 3]. In Fig. 2, we have plotted the complete tensor power spectrum described by the expression (2.24) for a given set of parameters.

We should stress that the power spectrum is actually valid only for modes which satisfy the condition . It is evident from the figure that the power spectrum is strictly scale invariant over this domain. Moreover, we find that the spectrum indeed reduces to the above-mentioned scale invariant amplitude for small values of the wavenumbers. We have also evaluated the tensor power spectrum numerically using the method described above. We have computed the spectrum at a given time soon after the bounce (corresponding to ) for all the modes. We find that, for wavenumbers such that , the numerical analysis also leads to a scale invariant spectrum whose amplitude matches the above analytical result to about .

3 The tensor bi-spectrum and the corresponding non-Gaussianity parameter

As we have mentioned, the most comprehensive formalism to calculate the three-point functions generated in the early universe is the formalism due to Maldacena [27]. The primary aim of Maldacena’s approach is to obtain the cubic order action that governs the scalar and the tensor perturbations using the ADM formalism. Then, based on the action, one arrives at the corresponding three-point functions using the standard rules of perturbative quantum field theory.

The tensor bi-spectrum in Fourier space, viz.  , evaluated at the conformal time, say, , is defined as

 ⟨^γk1m1n1(ηe)^γk2m2n2(ηe)^γk3m3n3(ηe)⟩ ≡ (2π)3Bm1n1m2n2m3n3γγγ(k1,k2,k3)δ(3)(k1+k2+k3).

Note that the delta function on the right hand side implies that the wavevectors , and form the edges of a triangle. For convenience, hereafter, we shall set

 Bm1n1m2n2m3n3γγγ(k1,k2,k3)=(2π)−9/2Gm1n1m2n2m3n3γγγ(k1,k2,k3). (3.2)

The tensor bi-spectrum , calculated in the perturbative vacuum using the Maldacena formalism, can be written in terms of the modes as follows [27, 31, 32, 34, 44, 33]:

 Gm1n1m2n2m3n3γγγ(k1,k2,k3) = M2Pl[(Πk1m1n1,ijΠk2m2n2,imΠk3m3n3,lj (3.3) −12Πk1m1n1,ijΠk2m2n2,mlΠk3m3n3,ij)k1mk1l+five permutations] ×[hk1(ηe)hk2(ηe)hk3(ηe)Gγγγ(k1,k2,k3) +complex conjugate],

where is described by the integral

 Gγγγ(k1,k2,k3)=−i4∫ηeηidηa2h∗k1h∗k2h∗k3, (3.4)

with denoting the time when the initial conditions are imposed on the perturbations and representing the time when the bi-spectrum is to be evaluated. Also, we should mention that denote the components of the three wavevectors along the -spatial direction111Such a clarification seems necessary to avoid confusion between , and which denote the wavenumbers associated with the wavevectors , and , and the quantity which represents the component of the wavevector along the -spatial direction..

The dimensionless non-Gaussianity parameter that characterizes the amplitude of the tensor bi-spectrum is defined as [34]

 hNL(k1,k2,k3) = −(42π2)2[k31k32k33Gm1n1m2n2m3n3γγγ(k1,k2,k3)] (3.5) × [Πk1m1n1,m3n3Πk2m2n2,¯m¯nk33PT(k1)PT(k2)+five permutations]−1,

where the overbars on the indices imply that they need to be summed over all allowed values. Our aim in this work is to evaluate the magnitude and shape of the tensor bi-spectrum and the corresponding non-Gaussianity parameter and compare with, say, the results in de Sitter inflation. Therefore, for simplicity, we shall set the polarization tensor to unity. In such a case, the expression (3.3) for the tensor bi-spectrum above simplifies to

 Gγγγ(k1,k2,k3) = M2Pl[hk1(ηe)hk2(ηe)hk3(ηe)¯Gγγγ(k1,k2,k3) (3.6) +complex conjugate],

where the quantity is described by the integral

 ¯Gγγγ(k1,k2,k3) = −i4(k21+k22+k23)∫ηeηidηa2h∗k1h∗k2h∗k3. (3.7)

We shall choose to be an early time during the contracting phase when the initial conditions are imposed on the modes (i.e. when ), and to be a suitably late time, say, some time after the bounce, when the bi-spectrum is evaluated. If we ignore the factors involving the polarization tensor, the non-Gaussianity parameter reduces to

 hNL(k1,k2,k3) = −(42π2)2[k31k32k33Gγγγ(k1,k2,k3)] (3.8) ×[2k33PT(k1)PT(k2)+two permutations]−1.

4 Evaluating the tensor bi-spectrum

With the forms of the scale factor and the mode functions at hand, in order to arrive at the tensor bi-spectrum, it is now a matter of evaluating the integral (3.7) in the three domains.

Let us begin by considering the first domain. Upon using the behavior (2.11) of the scale factor and the mode (2.12) in the first domain, we find that the quantity can be expressed as

 ¯G1γγγ(k1,k2,k3) = −i(k21+k22+k23)4M3Pla0k20√k1k2k3[I2(kT,k0,α)+i(1k1+1k2+1k3)I3(kT,k0,α) (4.1) −(1k1k2+1k2k3+1k1k3)I4(kT,k0,α)−ik1k2k3I5(kT,k0,α)],

where and the quantities are described by the integrals

 In(kT,k0,α)=∫−α/k0−∞dηηneikTη. (4.2)

For , these integrals can be evaluated to yield

 In+1(kT,k0,α) = −1n(−k0α)ne−iαkT/k0+ikTnIn(kT,k0,α), (4.3)

while is given by

 I1(kT,k0,α)=iπ+Ei(−iαkT/k0), (4.4)

where is the exponential integral function [58].

Let us now turn to evaluating in the second domain. Upon using the behavior (2.2) of the scale factor and the mode (2.16), we find that the quantity can be expressed as

 ¯G2γγγ(k1,k2,k3) = −ia20(k21+k22+k23)4k0[A∗k1A∗k2A∗k3J0(α) (4.5) +(A∗k1A∗k2B∗k3+A∗k1B∗k2A∗k3+B∗k1A∗k2A∗k3)J1(α) +(A∗k1B∗k2B∗k3+B∗k1A∗k2B∗k3+B∗k1B∗k2A∗k3)J2(α) +B∗k1B∗k2B∗k3J3(α)],

where are described by the integrals

 Jn(α)=∫0−αdx(1+x2)2fn(x), (4.6)

with the function being given by Eq. (2.17). The integrals and can be readily evaluated to obtain that

 J0(α)=α+2α33+α55 (4.7)

and

 J1(α) = −12(α2+α42)−1160+2(1+α2)15+(1+α2)220−8α15tan−1α (4.8) −4α15(1+α2)tan−1α−α5(1+α2)2tan−1α+415ln(1+α2).

In contrast, the integrals and are more involved. The integral can be divided into three parts and written as

 J2(α)=J21(α)+J22(α)+J23(α), (4.9)

where the integrals and can be easily evaluated to be

 J21(α) = ∫0−αdxx2=α33, (4.10) J22(α) = 2∫0−αdxx(1+x2)tan−1x (4.11) = α2(1+α22)tan−1α−12(α−tan−1α)−α36.

The quantity is given by

 J23(α)=∫0−αdx(1+x2)2(tan−1x)2, (4.12)

and, upon setting , it reduces to

 J23(α)=∫0−tan−1αdyy2sec6y. (4.13)

The integral involved can be evaluated to be (see, for instance, Ref. [58])

 ∫dyy2sec6y = −y(cosy−2ysiny)10cos5y−4y(cosy−ysiny)15cos3y+(1130+8y215)tany (4.14) +tan3y30+1615∞∑n=1(−1)n22n(22n−1)y2n+1(2n+1)(2n)!B2n,

where are the Bernoulli numbers. Needless to add, this result can be used to arrive at . We should add that the infinite series in the above expression is convergent, and we find that it can be expressed as follows [59]:

 = y{ln[Γ(1+yπ)]+ln[Γ(1−yπ)] (4.15) −ln[Γ(1−2yπ)]−ln[Γ(1+2yπ)]} +π{−ζ′(−1,1+yπ)+ζ′(−1,1−yπ) +12ζ′(−1,1+2yπ)−12ζ′(−1,1−2yπ)},

where denotes the derivative of the Hurwitz zeta function with respect to the first argument and is the Gamma function.

Let us now evaluate the last of the integrals, viz. . It proves to be convenient to divide the integral into four parts as follows:

 J3(α)=J31(α)+J32(α)+J33(α)+J34(α). (4.16)

If we set , we find that the integrals , and can be easily evaluated to be

 J31(α) = ∫0−tan−1αdytan3y=−α22+12ln(1+α2), (4.17) J32(α) = 3∫0−tan−1αdyy2tanysec4y (4.18) = −34(1+α2)2(tan−1α)2+α2(1+α2)tan−1α−α24 +αtan−1α−12ln(1+α2), J33(α) = 3∫0−tan−1αdyytan2ysec2y=α22−α3tan−1α−12ln(1+α2). (4.19)

The integral is given by

 J34(α)=∫0−tan−1αdyy3sec6y, (4.20)

which can be evaluated using the result [58]

 ∫dyy3sec6y = −y2(3cosy−4ysiny)20cos5y−2y2(3cosy−2ysiny)15cos3y (4.21) +(y+8y315)tany+ln|cosy|+ysiny10cos3y−120cos2y +85∞∑n=1(−1