Exact solutions of primordial gravitational waves

The future detection projects of gravitational waves are expected to have the sensitivity of detecting the primordial gravitational waves, which may be useful to get new insight into the very early universe. It is essential to analyze the evolution equation of the gravitational waves to estimate the present field strength of the primordial gravitational waves. In this paper, we obtain analytic solutions of the gravitational wave equation in the presence of non-relativistic matters as well as the cosmological constant. Although it is difficult to obtain the solution directly, we find that the equation for the square of the amplitudes has a simple polynomial solution. This quantity, which is directly related to the energy density of the gravitational waves, turns out to be useful to construct analytic solutions for the amplitudes by using Weierstrass's elliptic functions.


Introduction
The direct detection of gravitational waves [1] opened an era of gravitational wave astronomy. The ground-based experiments, such as KAGRA [2], Advanced LIGO [3], and Advanced VIRGO [4] are expected to observe many gravitational wave signals from binary systems of black holes and neutron stars, which may help us investigate various merging processes. The future space-based experiments such as eLISA [5], DECIGO [6], and BBO [7] are expected to have the sensitivity to investigate the primordial gravitational waves [8], namely, the gravitational waves produced in the early stage of the universe. These primordial gravitational waves are excellent tools to investigate some cosmological phenomena such as 1st order phase transition [9] and the inflation [10].
As for the predictions of the gravitational waves, the production and the evolution are two essential factors of the analysis. However, in most cases where the transition of the era is concerned, it is difficult to get any analytic solutions of the evolution equations. For example, the transition from the radiation-dominant era to the matter-dominant era has been analyzed in ref. [10]. It was shown by numerical calculation that the transition function can be obtained in a rather simple form. Later, various numerical improvements and analytical approximations have been achieved [11,12,13,14]. However, no analytic solution has been found yet. The main reason for the difficulty is that the equation has too many singular points. If an ordinary second-order differential equation has three definite singular points or the confluent limit of them, it is well known that the solution is given by the hypergeometric (or confluent hypergeometric) function, where various analytic continuations are applicable. In the case of an equation with four definite singular points, which is called Heun's equation (see for example [15]), we cannot usually get any analytical expression, so that it is difficult to expect the matching beyond the region of the convergence in the power series solutions. There are several exceptions, however. For example, analytic solutions have been obtained for Lamé equation, which has four singular points. As we will see later, the generic evolution equation of the gravitational waves has six singular points, so that we do not expect to have any analytic solutions. In the case of radiation to matter transition, the equation is the confluent Heun's equation.
The technique for solving the equations appeared in 2004 [16]. The technique is basically for Lamé equations [17] and later extended to the case where the equation has more singularities [18]. The basic idea of these papers are the followings: we solve the third order equation for the square of the solutions of the original equation instead of solving the equation directly.
In this paper, we will show that for the transition from the matter-dominant era to the cosmological-constant-dominant era, the evolution of the square of the amplitude of gravitational waves can be expressed as a polynomial of the reciprocal of the scale factor a, which is used to construct exact solutions to the gravitational wave equation in terms of Weierstrass elliptic function ℘.
In the next section, we will review the evolution equation in the cosmological setting. In section 3, we will explain the formalism of solving the second order equation out of the corresponding third order equation. We will apply the formalism to the gravitational wave equation and get exact results in the presence of matters and the cosmological constant in section 4 and 5. The final section will be devoted to discussions.

The Equation for gravitational waves
We consider homogeneous and isotropic background, namely the Friedmann-Robertson-Walker metric of the following form: (2.1) Here, we use the conformal time η as a time coordinate and K represents spatial curvature. Then, time evolution of the scale factor a is determined by the Friedmann equation, where H 0 is the present (a = 1) Hubble constant, a ′ = da/dη, and Here, G is the Newton constant, ρ r,0 and ρ m,0 denote energy density of relativistic and nonrelativistic matter at present respectively, and Λ is the cosmological constant. Defining a function g(a) as g(a) = Ω r + Ω m a + Ω K a 2 + Ω Λ a 4 , The gravitational wave equation (tensor mode) is given by represents the initial condition of the polarization, and e A ij (k) are the spin-2 polarization tensor satisfying the normalization condition ij e A ij (e B ij ) * = 2δ AB . Then, the equation for χ(η, k) is given by To observe the generic structure of the equation, it is convenient to change the variable η to a. Using the equation (2.5), we have Since g(a) is a quartic polynomial, this equation has singular points at g(a) = 0, a = 0, and a = ∞, in total six singular points.
If the equation has three definite singular points, the solution can be written by hypergeometric functions and we can use various analytic continuation. If the equation has four definite singular points, the equation is known as Heun's differential equation where we cannot obtain any good analytic continuation because these are essentially power series solution although many attempts has been done. When we include the radiation, matter and the curvature, i.e. g = Ω r +Ω m a+Ω K a 2 , the above equation reduces to the Heun's equation and further to the confluent Heun's equation in the limit Ω K → 0. This is the reason why we do not have any analytic solution for the transition from radiation dominant to matter dominant era.
We are going to deal with matter + cosmological constant dominant universe, which covers most of the cosmological time. In this case, the equation (2.9) still has five singular points. Normally we do not expect the equation has any analytic solution.
3 The third order equation associated with the second order equation We will follow the argument of ref. [17]. Let us consider a generic second order differential equation of the following form: where p(x) and q(x) are known functions. We consider a product of solutions y(x) = χ 1 (x)χ 2 (x), where χ 1 (x), χ 2 (x) are any solutions of the equation (3.1). Then the derivative can be written as Differentiating the both sides and using the equation (3.1), we get Then it is straightforward to obtain the equation for y(x) as follows: y ′′′ + 3py ′′ + (p ′ + 4q + 2p 2 )y ′ + (2q ′ + 4pq)y = 0. (3.5) In the second order equation of the form (3.1), it is well known that we can construct a constant from the Wronskian of χ 1 (x) and χ 2 (x) as and C(χ 1 , χ 2 ) is nonzero if and only if the two solutions χ 1 and χ 2 are linearly independent. We expect that we can get a constant from two independent solution of the equation (3.5). For two solutions y 1 = χ 1 χ 2 , y 2 = χ 3 χ 4 , we can define the following constant L out of the Wronskian, It is easy to see that this constant (3.7) can be written in terms of y 1 and y 2 as This result can be used to find the C(χ 1 , χ 2 ) directly from the solutions of the third order equation (3.5). Actually, the square of the constant C(χ 1 , χ 2 ) can be obtained by equating y 2 = y 1 ≡ y = χ 1 χ 2 ; Given a solution of the third order equation (3.5), which is assumed to have a nonzero constant (3.9), 1 we can construct two independent solutions of (3.1) as follows. From (3.9), we can obtain the value of C(χ 1 , χ 2 ) so that Dividing this by y = χ 1 χ 2 , we get On the other hand, equation (3.2) gives Combining (3.11) and (3.12), we get Integration of the above leads to the following solutions, (3.14) Since, for wave-like equations, the constant C is usually pure imaginary as we will see later, we write C = iC so that Another choice of the solution is given by If we can find non-oscillatory solution y for the third order equation (3.5), the solutions (3.15) can be regarded as phase space representation of the solutions. In our application, the construction of the solution above has a significant advantage. In physics, the power of the wave can be written as the time average (x = η or x = a in the cosmological setting) of the square of the field strength.
In any form of the wave equation, the power is proportional to y. Therefore, the evolution of y is physically significant.
We are going to give simple examples. The first one is the equation for a harmonic oscillator In this case, we know the solutions very well, The generic solutions are arbitrary linear combinations of these two independent solutions. The associated third order equation (3.5) becomes whose independent solutions are as follows: If we take a solution y = A 2 (A : const.), then from (3.9) we can see, as mentioned before, the constant C becomes pure imaginary, The construction of the solution (3.15) tells us that Therefore, y = A 2 = |χ ± | 2 is the square of the amplitude χ ± and is a constant in this case.
As a next example, let us consider the Bessel's equation of order ν, Then the third order equation (3.5) can be written as It is easy to get three independent power series solutions as . The constants C 2 for these solutions are calculated as In order to find the solution with the form χ = |χ|e iθ where the amplitude |χ| is not oscillatory, we have to search for y whose constant C 2 is negative. In this case, we know that such solutions are the Hankel functions, , (3.29) are the Bessel's functions of the 1st and 2nd kind, respectively. Comparing the power of the series of solutions, we find that the square of the amplitude |H ν (x) is written in terms of y 0 and y ±ν as This result is listed in ref. [19]. The constant C 2 for this solution y = H ν H (2) ν turns out to be as desired. Especially, when ν is a half-integer (ν = 1/2 + N, N :integer), we have a polynomial solution, In this case, we have used analytic continuation since we know the Bessel functions very well. However, if we find the polynomial solution for the third order equation, we can construct the Hankel function from (3.15). In other words, (3.15) represents a new way of writing Hankel functions. Figure 1 shows N = 1 case. These examples tell us that non-oscillatory solutions of the third order equation can be regarded as the evolution of the square of the amplitudes. We will apply the construction of the solution to the gravitational wave equation in the next section. Figure 1: 3/2 (x) which is a simple polynomial representing the amplitude of the oscillations.

Solutions for matter-cosmological constant dominant universe
In the previous section, we have shown that the non-oscillatory solution of the third order equation can be regarded as the square of the amplitudes. Before applying the formalism to the evolution equation for gravitational waves, we first show that physical observable quantities are directly related to the solution of the third order equation.
The energy density of the relic gravitational wave can be obtained by linear perturbation of the Ricci tensor in flat space as where P (k) is the primordial power spectrum and the overdot represents the derivative with respect to t,χ ≡ dχ/dt. The spectrum of the gravitational waves is usually described by the fraction of the energy density per logarithic frequency interval, where ρ crit = 3H 2 /8πG is the critical energy density of the universe. Substituting equation (4.1) and expressing it in the conformal time, we get If we set y = |χ(η, k)| 2 , then by using equation (3.4) we can express Ω gw in terms of y, Although we have to take time average for an oscillatory solution y, such a procedure is not necessary if we can find a non-oscillatory solution, e.g. polynomial solution, for the third order equation (3.5) directly. Energy and momentum of gravitational waves are expressed by the quadratic of the amplitude, therefore the direct analysis of the quadratic is also related to the analysis of these physical quantities. Now, let us apply the construction of the solution to the equation (2.9). The corresponding third order equation (3.5) becomes 4 n=0 b n a n θ(θ + n + 2) θ + 1 + n 2 + 4k 2 a 2 H 2 0 (θ + 2) y = 0, (4.5) where In general we can get power series solutions. However, in the case where the radiation is negligible, namely b 0 = 0, equation (4.5) allows the following polynomial solution, We can calculate C 2 for this solution by using (3.9) as follows, Furthermore, if we set b 2 = 0, one can see that C 2 is negative for any k ∈ R. In that case, the solution (3.15) for the second order equation becomes , Note that the time (a) average of the power is proportional to y, which depends on the wavelength. χ + and √ y are shown in figure 2.

The Explicit solutions
We have constructed gravitational wave solutions in the epoch where radiation and curvature components are negligible 2 . In this section, we will show that the solutions (4.9) can be expressed by elliptic functions if we use the conformal time instead of the scale factor. We first rescale the scale factor a in such a way that a = 1 when the energy density of the matter and that of the cosmological constant component are equal. We call the Hubble constant at that time H eq . Then the Friedmann equation (2.5) reduces to Note that although y cannot be regarded as the square of the amplitude in the case of positive b2 because y could be negative, the expression (4.9) for the solution is still valid for b2 = 0. integration of (5.1) gives where we have set the integration constant so that a = 0 at η = 0. We can solve this relation for a by introducing the Weierstrass's elliptic function ℘(z) (see, for example, [20]), where ω m,n = 2mω 1 + 2nω 2 and ′ m,n sums over all integers m, n except for (m, n) = (0, 0). ω 1 and ω 2 are half periods, meaning ℘(z + ω m,n ) = ℘(z) for any z ∈ C and m, n ∈ Z. The inverse of ℘(z) is known as where g 2 and g 3 are constants determined by the half periods, The half periods are obtained by the formula where e j = −e 2πi(j−1)/3 (j = 1, 2, 3) are the roots of x 3 + 1 = 0. The fundamental periods are shown in figure 3. We can see that a = ∞ corresponds to a zero of ℘(z), where the conformal time η takes the following value,η which is exactly 2/3 of the real half periods Before considering the gravitational waves, let us briefly see the behavior of this scale factor in the matter-dominant and Λ-dominant eras. More detailed accounts are given in Appendix A. Since ℘(z) ∼ 1/z 2 around z ∼ 0, scale factor behaves as a(η) ∼η 2 in the matter-dominant era (η ≪ 1). In order to see the behavior in the Λ-dominant era (η ∼η f ), we use the following differential equation for ℘(z): (℘ ′ (z)) 2 = 4(℘ 3 (z) + 1). (5.11) Then we can expand ℘(z) aroundη =η f as which leads to the de-Sitter expansion .
It is also instructive to get the relation between the cosmic time t and the conformal time η. Since the Friedmann equation can be solved as we can get the following relation, Let us write the solutions (4.9) in terms of the conformal time and perform the integration in the phase factor explicitly. y is written as where ℘(c j ) (j = 1, 2, 3) are roots of the third order equation x 3 + 8k 2 x 2 + 4 = 0 and we introduced k = k/H eq . Then, the integral in the exponent of (4.9) is decomposed as In order to perform the integral, we introduce the following functions [20]: where product ′ m,n is taken over all integers m, n except for (m, n) = (0, 0). These functions are related to ℘ as There are some addition theorems for these functions, among which we use the following formula: Using this formula, we find Then, we obtain the following expression of the solution (4.9), and Θ 2,3 are given by the cyclic permutation of c j . For completeness, we note that ℘(c j ) and ℘ ′ (c j ) (j = 1, 2, 3) are explicitly given by We can further simplify the solution (5.22) by using another formula inside the bracket, Although this expression seems complicated, it turns out that this solution includes, as limiting cases, the solutions in matter-dominant and cosmological constant-dominant era. Since our convention for the scale factor a is such that a = 1 at matter-Λ equality, we have to introduce another scale factor A taking arbitrary value A = A eq at the equality which depends on the energy fraction Ω m and Ω Λ at A = 1 in taking such limits. These two scale factors are related as a = A/A eq . With this re-definition of the scale factor, conformal time η and comoving momentum k are also transformed into new ones as τ = A eq η + const., q = A eq k.
Appendix A is devoted to detailed calculation of the limits Ω m → 0 and Ω Λ → 0 showing These are the well-known solutions expressed by the Hankel functions up to normalization constants. Also, in the sub-horizon (k → ∞) and super-horizon (k → 0) regime, the solution (5.27) reduces to simple forms as shown in Appendix B. Equation (B.7) shows that leading order in the sub-horizon limit is This is a manifestation that gravitational waves in this regime can be regarded as a collection of massless particles, namely gravitons, because it follows that the energy density ρ gw scales as ρ gw ∝ a −4 in this regime [21]. To see the super-horizon behavior, we expand the solution (5.27) as χ ± (η) = f ± (η) +k 2 g ± (η) + O(k 4 ), wherek = k/H eq . Then, f ± and g ± are determined in Appendix B as Considering long wavelength mode that enters the Hubble horizon after the energy density of radiation becomes negligible, two independent solutions χ ± should be superposed to form the almost-constant mode χ C which behaves χ C ∼const. at the begining of matter-dominant era (η ∼ 0). From (A.19), we can see the imaginary part of χ ± is such a superposition. Therefore, we set χ C = Imχ + = (χ + − χ − )/2i. By using thek expansion of χ ± (5.28), one can read super-horizon behavior of χ C as Note that because comoving Hubble parameter H = a ′ /a has a minimum H min = √ 3H eq /2 5/6 at a = 1/2 1/3 , modes with low comoving wavenumberk < √ 3/2 5/6 ∼ 0.97 never enter the Hubble horizon. Therefore expansion (5.29), which is valid fork ≪ 1, is not applicable to modes that enter the horizon in matter-dominant era. For such modes, we can use the expression (4.9) to get the following expression:

Conclusions and Discussions
We have shown that gravitational wave equation can be solved exactly in the presence of matters and the cosmological constant. While the gravitational wave amplitude can be written in terms of elliptic functions, square of the amplitude, which is directly related to the energy density, can be written as a simple polynomial of 1/a (4.10). Because this solution can describe the propagation of gravitational waves of arbitrary wavelength in the epoch after energy density of radiation becomes ignorable, which includes the present, it could be used to precise calculations of stochastic backgrounds of gravitational waves, which are expected to be directly observed in the future, as well as long wavelength modes affecting anisotropy of the cosmic microwave background. In the former case, a crucial quantity is the intensity (4.4). Since our solutions are valid after entering matter-dominant era, once a matching condition at radiation to matter transition is given, we can complute Ω gw exactly. On the other hand, in the latter case, longer wavelength modes are of main interest. For such a purpose, expression (5.30) would be useful. Unfortunately, we cannot obtain the exact results for the equation in the presence of radiations and matters because the third order equation (3.5) do not allow any polynomial solutions and it is difficult to get non-oscillatory solutions. However, it was shown in ref. [10] that the transfer function of the amplitudes are obtained numerically as T (k/k eq ) = [1.0 + 1.34(k/k eq ) + 2.50(k/k eq ) 2 ] 1/2 , (6.1) where k eq is the scale that entered the horizon at matter-radiation equality. This result indicates that the analysis of the square of the amplitude may be convenient for obtaining the transfer functions analytically.

A Matter dominant & Λ dominant limit
Because our consideration treats both the non-relativistic matters and the cosmological constant, all the resuts obtained above should include matter-dominant case and Λ-dominant case. In order to see that, we have to change the scale factor a, which is normalized such that a = 1 at the matter-Λ equality, to another one A which takes general value A eq at the equality. The relation between these two scale factors is A = A eq a. Note that this relation induces the relation between the conformal time η and one associated with A, τ , as η = A eq τ + const.. Introducing the energy fraction of matter Ω m and cosmological constant Ω Λ at A = 1, the Friedmann equation for A is given as where H 0 = (dA/dτ ) A=1 . Comparing this equation with (5.1), we get A eq = (Ω m /Ω Λ ) 1/3 and H eq = H 0 2/(1 + A 3 eq ). The matter-dominant case and Λ-dominant case correspond to the limit Ω Λ → 0 and Ω m → 0, where A eq → ∞ and A eq → 0, respectively. By using the solution (5.7), A is written as where τ 0 is a constant. In the following subsections, we will give the behavior of the scale factor and the solution χ ± under the limit of Ω m → 0 and Ω Λ → 0.
A.1 Λ-dominant case : Ω m → 0 Because Λ-dominant era corresponds toη ∼η f , we expand ℘(η) around there. We set τ 0 so that η =η f at τ = 0. We can obtain the explicit form by using the differential equation ℘ ′ (z) = −2 ℘ 3 (z) + 1 and the definition ofη f , i.e. ℘(η f ) = 0. The result is By substituting this into (A.1), and taking the limit Ω m → 0, we obtain which represents exactly the de-Sitter expansion. Before taking the limit of wave function χ ± , we expand it aroundη =η f . By using equations (C.2) and (C.3), we get Summations are evaluated by using equations (C.4) -(C.7), which lead to dχ ± /dη|η f = 0 and Finally, we consider the limit Ω m → 0 of the gravitational wave (5.27). Square of absolute value of the wave function |χ ± | 2 = y is given by (5.16). We note that the comoving wavenumber k is the physical wavenumber at a = 1, which is related to the comoving wavenumber q for the scale factor A as k = q/A eq . Thus, what should be kept constant in this limit is not k but q. Then, y approaches 2y = A eq A 3 + 8 q A eq H eq 2 A eq A 2 + 4 → 4(1 + q 2 τ 2 ).
Next we pick up the phase factor of χ ± normalized by the value atη =η f , Because k → ∞ as Ω m → 0, we consider c n under the limit k → ∞. The limit of ℘(c j ) is easily obtained from the equation (5.25) as Since the elliptic function ℘(z) has second-order poles at z = 2nω 1 + 2mω 2 ( ∀ n, m ∈ Z), c 1 approaches one of these poles. However, at the same time, the phase factor (A.4) is invariant under c j → c j + 2nω 1 + 2mω 2 , therefore we can take c 1 such that c 1 → 0 with k → ∞. By the similar argument, we can take c 2,3 such that c 2,3 →η f . In order to get the expansion of c j around k = ∞, we use the expansion of ℘(η) aroundη = 0 andη =η f , the latter of which is given by (A.2). The former is known as These expansions of ℘(η) determines the expansion of c j as follows: , , We must convert thesek expansion into the Ω m expansion. To do so, we need the expansion of k around Ω m = 0, which is easily obtained by the definitionk = q/A eq H eq and the constraint Also, we need the expansion ofη provided that τ andη f are kept fixed. The result is Note thatk(η −η f ) = qτ /2 √ 2 is kept constant under this limit in full order. We first compute n = 1 factor in the product (A.4). As c 1 ±η → ±η f , we need the expansion of σ(z) around z = ±η f , which is derived from the relation (5.19) and the definition ofη f , namely Then, it is straightforward to get Combining the above results, we find (A.14) Similar calculation shows As a final result, we find that the exact solution (5.27) converges into which is exactly the known result for the case of de-Sitter expansion.
As in the previous case, we expand the wave function χ ± aroundη = 0. ℘(η) is expanded as where we used the relations (5.19). Therefore we obtain Let us calculate the limit Ω Λ → 0 of the exact solution (5.27). The procedure is almost the same as in the case of Ω m → 0 in the previous section. In this case, the value of the scale factor at the equality, A = A eq , diverges as A eq = O(Ω The treatment of the phase factor doesn't so differ from the case of Ω m → 0, therefore we show only the results: Thus, we obatin which, apart from the normalization, recovers the solution for matter-dominant era.

B Wavenumber dependence of the solution
In this appendix, we consider smallk and largek limits of the exact solution (5.27). First, let f ± = lim k→0 χ ± , then f ± should be a solution for the equation (C.1) withk = 0, This zero-mode equation can be solved as where A and B are integration constants. Integration in the first term is performed by using the fact ℘ ′′ (z) = 6℘ 2 (z), where we re-defined the constant A → 6A. We fix these constants by considering the behavior aroundη = 0. f ± is expanded by setting k = 0 in (A. 19) as whereas ℘ ′ (η) = −2/η 3 + O(η 3 ). From these equations, we conclude that In order to see the super-horizon behavior of χ ± , we expand as where g ± will be determined below. Expanding (5.27) in terms ofk and directly deriving g ± (η) are very hard task, therefore we use the differential equation (C.1). Inserting (B.3) into (C.1), the equation for g ± reads This equation is easily solved for dg ± /dη as Thus, g ± is given as Two integration constants implicitly included in this expression reflect the fact that the equation (B.4) can determine g ± only up to the zero modes A℘ ′ (η) + B. These constants can be determined by considering expansion aroundη = 0 again. The result is We can explicitly perform the first integral of the second term by using ℘ ′′ (η) = 6℘ 2 (η), followed by the changes of integration variables, In the second and third equalities, we have changed the integration variables as a = 1/℘(z) and x = 1/a respectively, and done partial integration to obtain the final expression.