{'fraction_non_alphanumeric': 0.04280458658932839, 'fraction_numerical': 0.027148956738902764, 'mean_word_length': 3.7175069673169494, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lagrangian perturbations theory : Third-order solution for general dark energy models

We investigate the Lagrangian perturbation theory of a general dark energy models with a constant equation of state, \omega, and derive the fitting form of the fastest growing mode solutions up to the third order. These fitting forms are less than few % errors compared to the numerical calculation. We also correct the solutions of Einstein-de Sitter and the open universe with the proper initial conditions and the correct consideration for the growing mode solutions in order to compare with the general dark energy models. Even though these fitting forms are limited to the constant equation of state models, one can apply these to the time-varying \omega by interpolating between models with the constant \omega. These fitting forms can be used for the correct prediction for the two point correlation function and can give the proper predictions for the different dark energy models.


INTRODUCTION
The formation and the evolution of the large scale structures (LSS) in the Universe through the gravitational instability is one of the important probes of the cosmological parameters [1]. One might be able to constrain the nature of dark energy through future galaxy surveys aim at a percent precision on a broad range of scales. One can rely on the linear perturbation theory on large scales where the amplitude of the density fluctuations is small, whereas numerical simulations or phenomenological models, such as the halo model are necessitate to investigate the small scale LSS which incorporate non-linear effects [2]. The analytical investigations beyond the linear theory are suitable for the accurate theoretical predictions to the intermediate scales observations, such as the baryonic acoustic oscillations (BAO) [3].
There have been several analytic techniques for the quasi-linear perturbative expansion in order to solve the equations for an irrotational and pressureless fluid of cold dark matter (CDM) [4]. Two major methods are the standard (or Eulerian) perturbation theory (SPT) and the Lagrangian perturbation theory (LPT). In SPT, the density and the velocity field of the fluid are the perturbed quantities which need to be small in order to satisfy the perturbative expansion [5]. The perturbed quantity is the deviation of the particle trajectory field from the homogeneous background expansion in LPT [6,7]. Thus, LPT does not rely on the smallness of the density and velocity fields. However, most of these approaches are based on the single stream approximation and neglect the shell-crossing effects.
The solutions for higher orders of the LPT on Friedmann Lemaître Robertson Walker (FLRW) background have been investigated [6]. For irrotational flows, the general second order solution is found [8], the third order one with slaved initial conditions is obtained [9], and the fourth order one is also solved [10,11]. These solutions are obtained for the Einstein-de Sitter universe and/or the open universe. However the Universe is dominated by the dark energy at present epoch and thus one needs to find the higher order solutions based on the dark energy with CDM (ωCDM) universe. The third order solutions in Eulerian perturbation theory for the general dark energy model was investigated [12] Lagrangian resummation theory (LRT) has been developed which investigate both non-linear biasing and redshift space distortions [13][14][15][16]. From these method, one can predict BAO peaks in both real and redshift spaces [17][18][19][20][21]. Even though LRT improves the theoretical prediction for the quasi-linear observable, it is not so accurate compared to the realistic value because one adopts the LPT higher order solutions based on EdS universe instead of ωCDM.
In the present paper, we reexamine the fastest growing mode solutions for the LPT up to the third order for the EdS and open universe. Even though there have been the known solutions for these models, they are given by incorrect initial conditions and thus provides the improper coefficients. From solutions for the EdS model with the correct initial condition, one can obtain the exact numerical LPT solutions for flat ωCDM universe. We provide the approximate analytic solutions for ωCDM models up to the third order with less than 5 % errors. We are preparing the dark energy dependent the power spectrum and other statistical observables incorporating these solutions.
In section 2 we briefly review LPT. The differential equations order by order for the general ωCDM models are given. Then, in section 3 we investigate the solution for the each order for the EdS and the open universe with the correct initial conditions. We also obtain the analytic approximation solution of the displacement vector up to the third order for ωCDM models. We make conclusions in section 4.

LAGRANGIAN PERTURBATION THEORY
It is convenient to define the physical distance as r = a(t) x where a(t) is the expansion scale factor and x is the comoving Eulerian coordinates. The gravitational evolution of structure for a collisionless (cold) dust in subhorizon scales corresponding to those of the large scale structure at present is described by the Newtonian equations (mass conservation, Euler, and Poisson equation, respectively) where δ( x, t) is the density contrast, v( x, t) = a˙ x is the peculiar velocity, φ( x, t) is the peculiar gravitational potential,ρ m is the mean matter energy density, σ ij is the stress tensor [1,4,5]. Note that the mass conservation equation couples the zeroth moment (ρ) to the first one ( v) of the distribution function, the Euler equation couples the first moment ( v) to the second one (σ ij ), and so on. One can close the hierarchy by postulating an ansatz for the stress tensor σ ij , i.e. the equation of state of the cosmological fluid [5].
For this purpose, one can assume the irrotationality condition ∇ × v = 0 in addition to the above equations (2.1) -(2.3). In Eulerian perturbation theory, one expands the density contrast as where ǫ is just a book-keeping device with ǫ ≪ 1. From the second order, the behavior of δ is no longer local. The density perturbation at one spot depends on the initial perturbation at other places through the peculiar velocity.
Instead of working on Eulerian non-linear perturbation theory which studies the dynamics of density and velocity fields, it is possible to develop non-linear perturbation theory in the Lagrangian scheme by following the trajectories of particles or fluid elements [22][23][24]. This is possible because the choice of fields representing the cosmic quantities is not unique and Lagrangian approach is a change of the coordinate system itself. According to the Lagrangian point of view, the path of each fluid is followed during its evolution. Each fluid particle is labeled by its initial coordinate, q ≡ x(t 0 ) and the position of the q-particle at later time t may be described by x( q, t) = q + S( q, t) , (2.5) where S is the displacement vector. Thus, the motion of the fluid element may be completely described by the displacement S. A slight perturbation of the Lagrangian particle paths carries a large amount of non-linear information about the corresponding Eulerian evolved observable, since the Lagrangian picture is intrinsically non-linear in the density field. This overcomes the difficulty of the Eulerian approach by allowing the large density contrast δ. During the highly non-linear evolution, many particles coming from the different original positions will tend to arrive at the same Eulerian place. Since the initial mass density field is sufficiently uniform, one can assume that the Eulerian mass density ρ m ( x, t) at any given time t satisfies the continuity relation If one uses x(t 0 ) = q, then one obtains where ∇ = ∇ q and the Jacobian matrix J αβ is given by where I is the identity matrix and S is a 3 × 3 matrix whose elements are S αβ called the deformation tensor. In general, the deformation tensor is not symmetric, i.e. S αβ = S βα . S αβ is symmetric if and only if the displacement vector S is an irrotational field in the Lagrangian space [25]. Thus, one obtains the inverse Jacobian matrix as where S c αβ is an element of the cofactor matrix S c . If one takes the divergence on the Euler equation (2.2) using the above relations between the Eulerian and the Lagrangian space, then one obtains [25,26] ( In addition to the above Eq. (2.11), one can impose the irrotationality of the peculiar velocity [10] ǫ αβγṠγβ = ǫ αβγ S µβṠµγ + S αν ǫ νβγ S µβṠµγ −Ṡ γβ . (2.12) Now one expands the displacement vector S according to the Lagrangian perturbative prescription This explicit separation with respect to the spatial and temporal coordinates for each order is not an assumption but a property of the perturbative Lagrangian description for an Einstein-de Sitter universe [25]. The solutions at each order can be non-separable functions of t and q for the general model, but we will adopt the above ansatz for our models. Thus, if one adopts the above expression for S in Eq. (2.13), then one can expand the right hand side term of Eq. (2.11) as = Dµ 1 (S (1) ) + Eµ 1 (S (2) ) + D 2 µ 2 (S (1) , S (1) ) + F µ 1 (S (3) ) (2.14) are defined by Thus, one can obtain Lagrangian Poisson equation order by order (from the linear to the third orders) where dots represent the derivatives with respect to the cosmic time t and µ 2 (S (1) , S (2) ) = µ 2 (S (2) , S (1) ) is satisfied for any tensor [10]. Also, one can obtain the one more equation for the irrotationalitÿ

LPT SOLUTIONS
We now solve the dynamical equations for the temporal part of the displacement vector S order by order given by Eqs. (2.18) -(2.21). We investigate the solutions for three different models, Einstein-de Sitter (EdS), open universe, and the universe of the general dark energy with the constant equation of state (ωCDM). Even though the solutions for both EdS and open universe are well known, we reinvestigate those models on purpose. We can obtain the initial conditions of the fastest growing mode solutions for the general dark energy models from those of EdS because both models are close to each other at early epoch. We also check the so-called BJCP solutions for the open universe with the proper initial conditions to correct them [24]. We obtain the fitting forms of the fastest growing solutions for ωCDM. We show the errors of these fitting forms in the appendix.

A. First order time component
Even though the exact first order solutions for the considered models are well known, we reinvestigate these solutions by solving Eq. (2.18)D Now we show the solutions of the above equation for the different cosmological models.

i) case I : EdS
By replacing the cosmic time t with τ as dτ = α a 2 dt, the above Eq. (3.1) is rewritten by The solution of the above equation can be written as where c aD and c bD are the integral constants and g 1a is a growing mode solution and g 1b is the decaying one. One can set c aD = 1 from the initial condtion g 1a (a i ) = a i .
ii) case II : Open Universe One can repeat the above process for the open universe. Again by defining dτ = τi a 2 dt, one obtains where A ≡ Ωm0 1−Ωm0 . Thus, one gets and If one replaces the variable y with t, then the above equation (3.1) becomes The analytic solution of the above Eq (3.9) is well known as BJCP and given by 10) where L(y) = 1 2 ln y−1 y+1 and L(a) = 1 2 ln In the references, one separates g 1a and g 1b as a growing and a decaying mode solution, respectively. However, this is not true. In general case, we are not able to separate the growing mode from the decaying one. One can check this from the above solutions Eqs. (3.10) and (3.11). Both g 1a and g 1b are growing as a increases even though g 1b is negative. The growing mode solution of Eq. (3.9) is obtained by using two initial conditions.
where we use a i is the initial epoch. We adopt a i = 1 50 which gives the proper sub-horizon behavior for the growth factor. Thus, we can obtain the correct c aD and c bD from these initial conditions. In BJCP, they ignore this fact and just choose the g 1a as the growing mode solution with c aD = 1. We show the evolutions of D + and D BJCP and the difference between them in Fig. 1. We choose Ω m0 = 1 3 in this figure. The error is about 23 % at the present epoch. We also show the rescaled BJCP solution by normalized the coefficient of g 1a to satisfy (3.14) In this case, the error between the exact solution and the rescaled BJCP is less than 1 % for 1 50 ≤ a ≤ 1. Later, we will consider the fastest second order solution and in this case we need to use the r aD instead of c aD in order to get the correct evolution of the second order solution.
iii) case III : ωCDM In this case, one can rewrite the above Eq. (3.1) as where Y ≡ Aa 3ω . The solution of the above equation is found as [28,29] where F is the hypergeometric function. If one uses the transformation formulas for the hypergeometric function [30], then one obtains The above solution can be interpreted as the solution of Eq. (3.15) after replace the variable Y = −Z −1 .
With this replacement the above Eq. (3.15) becomes When ω = − 1 3 , this is the same as the open universe in the case II. With using the initial conditions of the growth factor given by Eqs. (3.12) and (3.13), one can obtain the integral constants. Thus, we can regard the open universe as the specific case of ωCDM. Even though we already obtain the exact analytic solution for the first order in ωCDM, we also obtain the fitting form of g 1a where We demonstrate the evolution of D + for the different cosmological models in the first row of Fig. 2. In the left panel of the first row, we show the evolutions of D + for the different values of ω when Ω m0 = 0.3. The dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. One obtains the larger D + for the smaller value of ω. This is due to the fact that there has been the more matter density in the past for the smaller values of ω to give the larger values of D + at present. At present, the value of D + varies from 0.83 to 0.6 for − 4 3 ≤ ω ≤ − 1 2 . In the right panel, we also show the evolutions , respectively. The errors are less than 1 % for all models up to z ≤ 6. In the second column, we check the errors with the different Ω m0 values for the ΛCDM model using the same notation as that of the first row. The errors are less than 1 % for Ω m0 = 0.30 and 0.35. The error can be about 2 % for Ω m0 = 0.25.
We also investigate the dependence of D + on the different cosmological parameters at the specific z. In the first row of Fig. 3, we show the values of D + as a function of ω and Ω m0 . In the left panel, we fix the redshift at z = 1 and check the dependence of D + on ω for the different values of Ω m0 . The dashed, solid, and the dotdashed lines correspond to Ω m0 = 0.25, 0.30, and 0.35, respectively. Again, the bigger the ω, the smaller the D + for each Ω m0 . Also we notice that the change rate is steeper for the smaller value of Ω m0 . In the right panel, we show the property of D + as a function of Ω m0 for the different ω models at z = 1. Again, the dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. One interesting point is that D + dependence on Ω m0 becomes weaker as ω decreases. For example, D + for ω = −1 model changes only 2 % when Ω m0 varies from 0.25 to 0.35. However, D + changes about 10 % for the same variation of Ω m0 when ω = − 1 2 . The errors of the fitting form at the specific z for the different cosmological models are also studied. In the second row of Fig. 3, we investigate the errors of the fitting form for the different values of ω and Ω m0 . In the first column, we check the dependence of errors on ω for the different values of Ω m0 . Except Ω m0 = 0.25, the errors are less than 1 % for −1.2 ≤ ω ≤ −0.6. In the second column, we show the errors as a function of Ω m0 for the different ω models at z = 1. The error can be as large as 2 % when Ω m0 = 0.25.

B. Second order time component
In general, the second order equation satisfies, where D is the solution of the first order equation. Unlike the first order solution, the second order solution constrains both the spatial and temporal dependence as shown in Eq. (2.19). In the appendix, we prove that the solution E can be obtained from the specific relation with the first order solution for EdS case. Again, we consider the solutions of this equation for the different cosmological models. We need these solutions to obtain the higher order solutions.

i) case I : EdS
Again, one can rewrite the above equation by using x After replacing D given by Eq. (3.3) into the above equation, one obtains where e a and e b are the integral constants for the homogeneous solution of the above Eq. (3.26). For the fastest growing solution, one ignores the term including g 1b . Thus, the fastest growing mode solution of the second order becomes (3.28) Now one can consider the initial condition of the above solution. At the early epoch, only the first order term contributes to the perturbation and thus one can put E+(a i ) = 0. From this initial condition, one can link the e a with the initial epoch as where we use c aD = 1. With the above relation, the second order solution is rewritten as Thus, one can obtain the initial condition of dE+ da as ii) case II : Open Universe In this case, the second order perturbation equation (3.25) becomes From the above solutions, one can find that The coefficients of the particular solution are determined from the integral constants of the first order solutions. Thus, one needs to determine the integral constants of the homogeneous solutions. If one just considers the contribution from the first part of the first order solution, then one can ignore the three terms including g 1b and g 2b in the above Eq. (3.33) and the solution becomes where c +aE = r 2 aD and r aD is given by Eq. (3.14). We replace the coefficient c aE into c +aE as we mentioned before. If we only consider the first part of the first order solution, then we need to rescaled the coefficient by using the initial condition. Thus, we can obtain e a by using the fact that Now we can compare the correct solution given by Eq. (3.38) with BJCP solution given as (3.40) As we already see in the first order solution g 1a , the correct solution should have the proper constant of integration to satisfy D + (a i ) = a i . In BJCP, they ignore this fact and this causes the improper behavior for the second order solution. We show the evolutions of E + with the correct rescaled coefficients and E BJCP and the differences between them in Fig. 4. We show the evolutions of both the correctly rescaled second order solution E res and the E BJCP in the left panel of Fig. 4. The solid and dashed lines correspond to E res and E BJCP , respectively. As we can see E BJCP is overestimated compared to the correct second order solution. We again adopt Ω m0 = 1 3 in this figure. We also show the difference between E res and E BJCP in the right panel of Fig. 4. The error between them is about 15 % at the present epoch. Even though we show the exact solution, we also provide the approximate solution for this model.
ii) case III : ωCDM In this case, one can rewrite the above Eq. (3.1) as where D is the first order solution. The homogeneous solution is same as the first order one with the different integral constants and the nonhomogeneous solution should be obtained from the first order solutions. There is no analytic solution for the above equation but we can obtain the approximate fitting form of the fastest growing mode solution. If we adopt the initial conditions Eqs. (3.30) and (3.31), we can obtain the numerical solution of the above Eq. (3.45). Similar to the approximate first order solution, we obtain the fitting form of E + for the fastest growing solution as where we abbreviate subscript + and From now on, we will drop the subscript + for the fastest growing solution.
From the above fitting form, we obtain several properties of the fastest growing solution E. First, the signature of E is opposite to that of D. Thus, it decreases as a function of time. Second, as ω decreases, so does E. One can understand this because the dark energy with the negative ω acts like the negative pressure. Third, E decreases as Ω m0 increases. These properties allow us to study the behaviors of E for the various cosmological parameters compared to those of EdS. We show the time evolution of E for the different cosmological models in the first row of We investigate the accuracy of this fitting form for the different cosmological models. In the second row of Fig. 5, we show the errors of the fitting form as a function of the redshift z for the different models. In the first column, we show the ∆E E for the different values of ω when we fix Ω m0 = 0.3. The errors are less than 2 % for all considered models when z ≤ 4. In the second column, we check the errors of the fitting form for the different Ω m0 values for the ΛCDM model. The errors are less than 3 % for all the considered Ω m0 values.
We also investigate the dependence of E on the cosmological parameters at the specific z. In the first row of Fig. 6, we show the values of E as a function of ω and Ω m0 at the specific redshift z. In the left panel, we fix the redshift z = 1 and check the dependence of E on ω for the different values of Ω m0 . The dashed, solid, and the dotdashed lines correspond to Ω m0 = 0.25, 0.30, and 0.35, respectively. For Ω m0 = 0.35, E varies from -0.078 to -0.098 when ω changes from -0.6 to -1.2. We also obtain −0.096 ≤ E ≤ −0.068 for −1.2 ≤ ω ≤ −0.6 when Ω m0 = 0.25. Thus, we can conclude that E dependence on ω becomes weaker as Ω m0 increases. In the right panel, we show the dependence of E on Ω m0 for the different ω models at z = 1. Again, the dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. For ω = − 1 2 , E varies from -0.057 to -0.069 when Ω m0 changes from 0.25 to 0.35. E changes from -0.096 to -0.098 when Ω m0 changes from 0.25 to 0.35 for ω = − 4 3 . This case E is almost constant for the different values of Ω m0 . Thus, we can conclude that E dependence on Ω m0 becomes weaker as ω decreases. The errors on the fitting form at z = 1 for the different cosmological models are shown in the second row of Fig. 6. We show the errors on the fitting form as a function of ω and Ω m0 . In the first column, we show the dependence of errors on ω for the different values of Ω m0 . The errors are about less than 3 % for −1.2 ≤ ω ≤ −0.7. In the second column, we show the errors as a function of Ω m0 for the different ω models at z = 1. The errors are less than 3 % for all Ω m0 except ω = − 1 2 .

C. Third order time component
In general, the temporal component of the third order equation splits into two systems, where D and E are the solutions of the first and the second order equations, respectively. Because we have the analytic forms of both D and E, we can obtain the approximate analytic forms of F a , F b , and F T for the general dark energy models. We investigate the solutions for the different cosmological models.

i) case I : EdS
First, we rewrite the equation for the F a by using x By using D, we obtain are the integral constants of the homogeneous solution of F a . If we consider the contribution of the fastest growing solution only, then the F a becomes We also adopt the previous consideration to obtain the initial conditions. From F a+ (a i ) = 0, one obtains where we use c aD = 1 again. With this f (a) a , one gets the initial condition for dFa+ da as We can repeat the same process for the F b . The differential equation for F b becomes The solution for F b is given by are the integral constants of the homogeneous solution of F b . Again the fastest growing solution is given by where we use e a = 3 7 a i . From this, initial condition for dF b+ da is given by We write the differential equation for the F T by using the same parameter x The solution for F T is given by (3.65) From this, the initial condition for dFF + da is given by ii) case II ; Open Universe The third order perturbation equations are rewritten as As we show in the previous subsection, this is the special case of general dark energy model with equation of state ω = − 1 3 . Thus, we will show the solutions of the above equations in the next case. iii) case III : ωCDM In this case, one can rewrite the Eq. (3.50) as We obtain the fitting form of the fastest growing solution F a as (3.75) From the above fitting form Eq. (3.72), we obtain several properties of the fastest growing solution F a . First, the signature of F a is opposite to that of D. Thus, it decreases as a function of time. Second, as ω decreases, so does F a . Third, F a decreases as Ω m0 increases. These properties are similar to those of E in the previous subsection. It is natural because E and F a have the source terms which is proportional to D 2 and D 3 , respectively. We investigate the behaviors of F a for the various cosmological parameters. We show the time evolution of F a for the different cosmological models in the first row of Fig. 7. In the left panel, we show F a (z) for the different values of ω when Ω m0 = 0.3. The dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. The present values of F a vary from −0.07 to −0.19 for ω = − 1 2 and − 4 3 , respectively. In the right panel, we show the evolutions of F a for the different values of Ω m0 when we consider the ΛCDM model. The dotted, solid, and dotdashed lines correspond Ω m0 = 0.25, 0.30, and 0.35, respectively. We obtain −0.18 ≤ F a ≤ −0.14 for 0.25 ≤ Ω m0 ≤ 0.35. In the second row of Fig. 7, we show the errors of the fitting form as a function of the redshift z for the different models. In the first column, we show the ∆Fa Fa for the different values of ω when we fix Ω m0 = 0.3. The errors are less than 4 % for all models when we consider z ≤ 3. In the second column, we check the errors of the fitting form for the different Ω m0 values when we consider the ΛCDM model. Again, the errors are less than 4 % for all the considered Ω m0 values up to z ≤ 3.
We also investigate the dependence of F a on the cosmological parameters at the specific z. In the first row of Fig. 8, we show the values of F a as a function of ω and Ω m0 at the specific redshift z. In the left panel, we fix the redshift z = 1 and check the dependence of F a on ω for the different values of Ω m0 . The dashed, solid, and the dotdashed lines correspond to Ω m0 = 0.25, 0.30, and 0.35, respectively. For Ω m0 = 0.35, F a varies from -0.023 to -0.040 when ω changes from -0.5 to -1.2. We also obtain −0.039 ≤ F a ≤ −0.018 for −1.2 ≤ ω ≤ −0.5 when Ω m0 = 0.25. Thus, we can conclude that F a dependence on ω becomes weaker as Ω m0 increases. In the right panel, we show the dependence of F a on Ω m0 for the different ω models at z = 1. Again, the dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. For ω = − 1 2 , F a varies from -0.018 to -0.023 when Ω m0 changes from 0.25 to 0.35. F a changes from -0.039 to -0.040 when Ω m0 changes from 0.25 to 0.35 for ω = − 4 3 . This case F a is almost constant for the different values of Ω m0 . Thus, we can conclude that F a dependence on Ω m0 becomes weaker as ω decreases. All of these properties are same as those of E. We also investigate the errors on the fitting form at z = 1 for the different cosmological models. In the second row of Fig. 8, we show the errors on the fitting form as a function of ω and Ω m0 . In the first column, we show the dependence of errors on ω for the different values of Ω m0 . The errors are about less than 2 % for −1.2 ≤ ω ≤ −0.7. In the second column, we show the errors as a function of Ω m0 for the different ω models at z = 1. The errors are less than 2 % for all Ω m0 except ω = − 1 2 . One can rewrite the Eq. (3.51) as Now we can repeat to the same process as before to obtain the fitting form of F b as where (3.80) From the above fitting form Eq. (3.77), we obtain several properties of the fastest growing solution F b . First, the signature of F b is same as that of D. Thus, it increases as a function of time. Second, F b increases as ω decreases. Third, as Ω m0 increases so does F b . These properties are similar to those of D in the previous subsection and opposite to those of E or F a . It is easy to understand because F b has the opposite sign of the source term compared to those of E or F a . We investigate the behaviors of F b for the various cosmological parameters. We show the time evolution of F b for the different cosmological models in the first row Fig. 9. In the left panel, we show F b (z) for the different values of ω when Ω m0 = 0.3. The dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. The present values of F b vary from 0.104 to 0.271 for ω = − 1 2 and − 4 3 , respectively. In the right panel, we show the evolutions of F b for the different values of Ω m0 when we consider the ΛCDM model. The dotted, solid, and dotdashed lines correspond Ω m0 = 0.25, 0.30, and 0.35, respectively. We obtain 0.198 ≤ F b ≤ 0.249 for 0.25 ≤ Ω m0 ≤ 0.35 at present. The fitting form of F b is given by Eqs. (3.77) -(3.80). We investigate the accuracy of this fitting form for the different cosmological models. In the second row of Fig.9, we show the errors of the fitting form as a function of the redshift z for the different models. In the first column, we show the ∆F b F b for the different values of ω when we fix Ω m0 = 0.3. The errors are less than 4 % for all models when we consider z ≤ 3. In the second column, we check the errors of the fitting form for the different Ω m0 values when we consider the ΛCDM model. Again, the errors are less than 4 % for all the considered Ω m0 values up to z ≤ 3.
We also investigate the dependence of F b on the cosmological parameters at the specific z. In the first row of Fig. 10, we show the values of F b as a function of both ω and Ω m0 at the specific redshift z. In the left panel, we fix the redshift z = 1 and check the dependence of F b on ω for the different values of Ω m0 . The dashed, solid, and the dotdashed lines correspond to Ω m0 = 0.25, 0.30, and 0.35, respectively. For Ω m0 = 0.35, F b varies from 0.032 to 0.055 when ω changes from -0.5 to − 4 3 . We also obtain 0.024 ≤ F b ≤ 0.054 for − 4 3 ≤ ω ≤ − 1 2 when Ω m0 = 0.25. Thus, we can conclude that F b dependence on Ω m0 becomes weaker as ω decreases. In the right panel, we show the dependence of F b on Ω m0 for the different ω models at z = 1. Again, the dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. For ω = − 1 2 , F b varies from 0.024 to 0.032 when Ω m0 changes from 0.25 to 0.35. F b changes from 0.054 to 0.055 when Ω m0 changes from 0.25 to 0.35 for ω = − 4 3 . F b is almost constant for the different values of Ω m0 when ω = − 4 3 . Thus, we can conclude that F b dependence on Ω m0 becomes weaker as ω decreases. The errors on the fitting form at z = 1 for the different cosmological models are shown in the second row of Fig. 10. In the first column, we show the dependence of errors on ω for the different values of Ω m0 . The errors are about less than 3 % for −1.2 ≤ ω ≤ −0.7. In the second column, the errors of the fitting form as a function of Ω m0 are shown for the different ω models at z = 1. The errors are less than 2 % for all Ω m0 including ω = − 1 2 . The Eq. (3.52) for the transverse mode in this model becomes The homogeneous solution for the above equation is given by Even though, one can obtain the homogeneous solution of the transverse mode, we need to find the fitting form of the fastest growing solution F T which can be given by which is proportional to D 3 . First, the signature of F T is opposite to that of D. Thus, it decreases as a function of time. Second, as ω decreases, so does F T . Third, F T decreases as Ω m0 increases. We investigate the behaviors of F T for the various cosmological parameters. We show the time evolution of F T for the different cosmological models in the first row of Fig. 11. In the left panel, we show F T (z) for the different values of ω when Ω m0 = 0.3. The dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. The present values of F T vary from −0.03 to −0.8 for ω = − 1 2 and − 4 3 , respectively. In the right panel, we show the evolutions of F T for the different values of Ω m0 when we consider the ΛCDM model. The dotted, solid, and dotdashed lines correspond Ω m0 = 0.25, 0.30, and 0.35, respectively. We obtain −0.18 ≤ F a ≤ −0.14 for 0.25 ≤ Ω m0 ≤ 0.35. We investigate the accuracy of this fitting form for the different cosmological models. In the second row of Fig.11, we show the errors of the fitting form as a function of the redshift z for the different models. In the first column, we show the ∆FT FT for the different values of ω when we fix Ω m0 = 0.3. The errors are less than about 4 % for all models when we consider z ≤ 5. In the second column, we check the errors of the fitting form for the different Ω m0 values when we consider the ΛCDM model. Again, the errors are less than 3 % for all the considered Ω m0 values up to z ≤ 5.
We also investigate the dependence of F T on the cosmological parameters at the specific z. In the first row of Fig. 12, we show the values of F T as a function of ω and Ω m0 at the specific redshift z. In the left panel, we fix the redshift z = 1 and check the dependence of F T on ω for the different values of Ω m0 . The dashed, solid, and the dotdashed lines correspond to Ω m0 = 0.25, 0.30, and 0.35, respectively. For Ω m0 = 0.35, F T varies from -0.01 to -0.017 when ω changes from -0.5 to -1.2. We also obtain −0.016 ≤ F T ≤ −0.008 for −1.2 ≤ ω ≤ −0.5 when Ω m0 = 0.25. In the right panel, we show the dependence of F T on Ω m0 for the different ω models at z = 1. Again, the dotted, dashed, solid, and dotdashed lines correspond to ω = − 1 2 , − 5 6 , −1, and − 4 3 , respectively. For ω = − 1 2 , F T a varies from -0.008 to -0.01 when Ω m0 changes from 0.25 to 0.35. F T changes from -0.017 to -0.016 when Ω m0 changes from 0.25 to 0.35 for ω = − 4 3 . This case F T is almost constant for the different values of Ω m0 . Thus, we can conclude that F a dependence on Ω m0 becomes weaker as ω decreases. All of these properties are same as those of F a . We also investigate the errors on the fitting form at z = 1 for the different cosmological models. In the second row of Fig. 12, we show the errors on the fitting form as a function of ω and Ω m0 . In the first column, we show the dependence of errors on ω for the different values of Ω m0 . The errors are about less than 3 % for −1.2 ≤ ω ≤ −0.7. In the second column, the errors are depicted as a function of Ω m0 for the different ω models at z = 1. The errors are less than 4 % for all Ω m0 except ω = − 1 2 .

DISCUSSION AND CONCLUSIONS
We reinvestigate solutions for the Lagrangian perturbation theory of an irrotational fluid up to the third order for the Einstein-de Sitter and open universe. With the correct initial conditions and the proper consideration for the fastest growing mode solutions, we correct the known solutions for these models. For the first time, we obtain the analytic approximate solutions for the general dark energy with the constant equation of state. These fitting forms have less than 5 % errors compared to the numerical solutions for all orders up to z ≤ 3.
So far, one have used the EdS solutions for the time component to incorporate the power spectrum or higher order moments even when one adopt the cosmology dominated by the dark energy at present epoch. This is an inaccurate approximation. Thus, with our analytic approximate solutions one can consider the correct dark energy dependent on observable [31]. Even though these solutions are obtained for the constant equation of state dark energy models, one can apply these solution to the time varying dark energy models by interpolating between models with the constant equation of states [29].
These analytic solutions provide us the tools to the systematic study for the dependence of the solutions for each order on both the matter energy density and the dark energy equation of state. We are also able to investigate the time dependence of models to scan for large parameter spaces with solutions.

ACKNOWLEDGEMENTS
We would like to thank Cornelius Rampf for fruitful discussion and comments on the manuscript. We specially thanks to the anonymous referee for the useful comments. We also thank KIAS Center for Advanced Computation for providing computing resources.

A. Einstein de Sitter Universe
Then, we can rewrite the Friedmann equation as If we use the fact that a i → ∞ as τ i → 0, then we obtain where we adopt a 0 = 1, then τ 0 = − 2 8πGρ m0 3α 2 .