General solutions to N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {N}$$\end{document}-field cosmology with exponential potentials

We construct the general analytical solution for the N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {N}$$\end{document}-field product-exponential potential in an expanding FLRW background. We demonstrate the relevance of this analytical solution in more general contexts for the derivation of estimates for the transitional time between an arbitrary initial state and the slow-roll solutions. In certain cases, these estimates can also be used to demonstrate the non-linear convergence towards the slow-roll solutions. In addition, we extend the solution to include stiff matter as well.


Introduction
Inflation is the leading paradigm for the origin of structure formation in an initially isotropic and homogeneous universe [1]. It was originally proposed as an attempt to evade the fine-tunings of the standard cosmological model by theorizing a period of accelerating expansion in the very early universe [2][3][4][5][6][7][8]. Despite the simplicity of this idea, the non-linear nature of scalar-field equations makes the quest for analytical solutions impossible; without further simplifications we have to rely entirely on numerical tools. Although efficient numerical codes have been developed during the last decade, for instance the publicly available codes [9][10][11][12][13], (approximate) analytical solutions are desirable for a number of reasons: if inflation solves the fine-tuning problems by utilizing an attractor mechanism, then this should be shown analytically rather than numerically in a case-by-case analysis. In addition, for multiple fields calculation of non-Gaussianities and other higher order correlators is computationally demanding and therefore simple analytical expressions may provide more insight on how to better search for viable models.
The aim of this work is the derivation of general solutions for an arbitrary number of fields. These solutions can be used to derive estimates for the transition time (or field displacement) from arbitrary initial states to inflation, or alternatively they can be used to demonstrate non-linear convergence towards the slow-roll solutions.
The paper is organized as follows: in Sect. 2 we derive the general solution for a product-exponential potential. In Sect. 3 we provide analytical estimates for the duration of the passage from an arbitrary initial state towards inflation. In Sect. 4 we continue the investigation of general solutions in the presence of a fluid. Finally, in Sect. 5 we summarize our main results.
Conventions: dots (˙) and primes ( ) indicate derivatives with respect to the cosmic time and number of e-folds respectively; Greek letters refer to spacetime indices while Latin ones to field-space indices; the spacetime metric has the mostly plus signature (−, +, +, +); we assume the Einstein's convention for repeated indices; N is the e-folding number and N denotes the number of fields.

Earlier works
Simple models of the early and late universe assume an action with minimal derivative couplings and canonical kinetic terms where the Latin indices label the number of fields, running from 1 to N , and κ = 8π G/c 2 (we will set κ = 1 from now on). Variation of the action with an FLRW ansatz for the metric yields the two linearly independent Einstein's equations for the unknown metric function a(t) (one evolution equation and one constraint) and the generalized Klein-Gordon equations for the scalar fields φ i (t). This system of equations can be be written in compact form with the definition of the Hubble parameter H = d(ln a)/ dt The latter 1 + N equations form a non-linear system for the unknown functions φ(t) and a(t) (or H (t)) which is in general unsolvable for an arbitrary potential function. In the cosmic time parameterization the ODE (5) becomes separable for a constant potential, which describes a de Sitter universe. In Ref. [15] the general solution for a two-field model consisting of one massless and one self-interacting field with an exponential potential, V (φ, ψ) = V 0 e pφ , was first constructed. Using the conservation of canonical momentum the problem was shown to be equivalent to a single exponential potential and stiff matter. Equation (4) becomeṡ and in combination with Eq. (3) provides the solution for a. It is worth noticing that the same solution holds for an arbitrary number of additional massless fields, whose influence can be absorbed in the definition of the constant c 0 . In the presence of a non-zero gradient for the second field this method can not be applied (at least without serious modification). 2 In the following we adopt a different method that allows the derivation of the general solution of a product-exponential potential, i.e.
for some constant numbers λ i and without assuming N − 1 of them to be zero. 3 2.2 Integrability of the product-exponential potential For N canonically normalized fields switching to the efolding number, defined from N ≡ ln a = H dt, the Klein-Gordon system of Eq. (5) can be written in first order form as where we defined p i ≡ (ln V ) ,i . The velocities in terms of the number of e-folds are related to the slow-roll parameter ε as In this formulation it becomes clear that the logarithmic derivatives of the potential (and not the gradients V ,i ) determine the evolution of ε. For the product-exponential potential of Eq. (7) p i 's are constant and Eq. (9) decouple from Eq. (8), so one can study the reduced system. Even though Eq. (9) remain coupled through the term 1 2 v k v k , the system is fully integrable as there exist integrals of motion which enable us to further reduce the dynamical variables of the problem. These can be found by dividing Eq. (9) for two fields e.g. i and j where we absorbed the initial conditions dependence in the definition of For N = 2 and p i = 0 this integral of motion is identical to Eq. (10) of Ref. [15]. Note that these integrals of motion make 0 (for j = 1) to the case of general gradients p j = 0 requires the Footnote 2 continued explicit construction of the rotation N × N matrix, a problem that is considerably more involved than just solving the system of equations without considering p j = 0, as we will show in the next section. 3 Note that observables evaluated on the late-time solutions (as one typical does in all inflationary models) do not satisfy the constraints from the latest Planck results, as they are identical to observables for singlefield exponentials. These models, therefore, may be more relevant as dark energy candidates.
sense if the initial conditions for the two fields are not chosen exactly on the critical points of (9). If for instance v k,0 = −p k then the unique solution for that field is v k (N ) = −p k and hence it is not dynamical. Discarding this equation one then solves the reduced system of N − 1 equations.

Asymptotic behaviour
Before deriving the solution we will study the asymptotic behaviour for N → ∞. The system of Eq. (9) with constant p i have the following set of critical points: where the ratio of kinetic to potential energy becomes constant. This critical point exists for p i p i < 6 and the critical value for ε is 1/2 p i p i . 2. Kinetic energy domination, i.e.
where the kinetic energy of the fields dominate their potential energy. The critical value for ε is 3, which defines an N -dimensional hypersphere and, hence, it is a collection of critical points.
The stability of these critical points trivially follows after linearization of the system around the critical points. For the scaling solution, translating the critical point to the origin by defining then the linearized equation of motion reads and hence we find negative eigenvalues when the condition p i p i < 6 is satisfied. We observe that whenever this critical point exists it is also stable. Regarding kinetic domination, the set of Eq. (9) is not very useful because more than one realizations can give ε = 3. This also implies that critical points are not isolated and one has to apply more sophisticated methods to infer stability. For this reason, we will use the evolution equation for the slow-roll parameter which is found by contracting Eq. (9) with Defining a new variable the linearized equation for x becomes and hence stable realizations of kinetic domination require We can view the second term as the dot product of the vector p i and the vector The condition (20) can be written as and since |r| = 1, the previous inequality implies that is necessary and sufficient condition for kinetic domination to be stable. The points of stable kinetic domination should satisfy the inequality (21). Considering N → −N , i.e. moving backwards in time, the stability of these critical points is reversed. An important note is in order: because eigenvalues for stable solutions are strictly negative numbers, the flow towards the critical point is monotonic; this means that the velocities for individual fields would never cross their critical value when they approach it from higher or lower values. This observation will be used later in Sect. 2.5.

Derivation of the solution
Having expressed all velocities in terms of e.g. the reference field v 1 then in principle we can solve the equation of motion for v 1 and obtain the relation v 1 (N ) and then using the integrals of motion we can find the expression for the rest fields. However, we find it more instructive to work with ε as it is a central variable in inflation and dark energy models and, as we will show later, it better highlights the behaviour of the solution. Therefore, in the following we will express all velocities in terms of the slow-roll parameter. Using the integrals of motion in the expression of the slow-roll parameter provides a quadratic equation for v 1 where we defined A 2 ≡ A i A i and p 2 ≡ p i p i . The two solutions are At this point we are unaware which root is the correct one and we need to compare this expression with another one; contracting Eq.
and, hence, the last term of Eq. (25) is equal to Thus, the sign of A i v i determines which root is the correct one and the relation of v 1 in terms of ε is where Using this expression all velocities are written in terms of the slow-roll parameter and the differential equation of ε becomes separable. Note, however, that we are able to find the analytical solution of the inverse function N (ε) When the initial velocities of every field, except for v 1 , take their critical value v j,0 = −p j , then these fields are nondynamical giving A i = 0 and one obtains the single-field solution.
The field displacement as a function of the slow-roll parameter can be obtained from Eq. (30) after a time redefi- where no sum in j is implied. This equation is separable and as previously the inverse function φ j (ε) can be integrated with solution With the expressions (31) and (33) we have fully determined the solution {a, φ j } in terms of ε.

Analysis of the solution
The solution (31) elegantly captures the past or future asymptotic behaviour of the Klein-Gordon equations. Here one has to distinguish between two cases, depending on the sign of s relative to A i p i .
1. When s ·(A i p i ) > 0, the second and fourth terms diverge to minus infinity for ε → 3, proving that kinetic energy domination is the past attractor. 2. When s · (A i p i ) < 0, there are two more subcases: (a) p 2 < 6. The first term diverges as ε → 1/2 p 2 and the scaling solution is the future attractor. (b) p 2 > 6. The scaling solution does not exist and kinetic domination is the future attractor.
The previous analysis indicates that if initial conditions belong to case 1 then s should change to −s in order for the system to reach one of the critical points described in case 2. Now we will briefly explain how one can use the full expression (31), while (33) can be treated similarly. We follow the following steps: fied initial conditions the solution will then be given as a function of ε 0 and s: N = N (s, ε 0 , ε). 2. We calculate the value of ε for which s changes. This would happen for a specific value of ε = ε * , found by setting Eq. (27) to zero: Although this value always exists (due to the Schwarz inequality (A i p i ) 2 ≤ A 2 p 2 ), it may not be relevant for the assumed set of initial conditions. This is because the system may relax at one of its critical points before crossing this value. 3. Next, we compare the sign of A i p i with s. If s = sgn (A i p i ) then we move to step 4 and otherwise to step 5. 4. For this set of initial conditions the full solution will consist of both branches, namely where ε cr corresponds to the logarithmic divergence of the solution. Of course, we know from the asymptotic analysis that the divergence corresponds to one of the two types of critical points, that is either kinetic domination or the scaling solution. 5. In this case we have to compare ε * with ε cr , that corresponds to the logarithmic divergence of the solution. If ε * < ε cr < ε 0 then ε would never reach ε * and the solution is In any other case the solution is In Fig. 1 we plot the analytical solution (31) for 50 fields, along with the numerical result, for a model with p 1 = 1, p 2 = 0.5, p j = 0.05 ( j = 3, . . . , 50) for two different sets of initial conditions. In the first one we consider v 1 = −1.5, v 2 = −1, v j = 0.01 yielding ε 0 = 1.685, ε * = 0.552, s = −1 and sgn(A i p i ) = 1 and, thus, we are in step 5. Comparing ε * with the asymptotic value 1/2 p i p i = 0.685 we observe that ε * < ε cr < ε 0 and the analytical solution is given by (37). This is depicted with the pale blue line in Fig. 1; moreover, we observe that the solution extrapolates from kinetic domination to the scaling solution located at ε = 0.685, but it never crosses ε * . For our second set of initial conditions we consider v 1 = 1.5, v 2 = 1, v j = 0.01 and so we find ε 0 = 1.685, ε * = 0.007, s = 1 and Finally, we provide some physical insight on the terms A i v i and A i p i . Their physical significance becomes more apparent if we examine the single-field case where A i = 1 and the two terms become simply v 1 and p 1 . When e.g 0 < p 1 < √ 6 then the asymptotic value for v 1 would be v 1 = − p 1 . If the initial velocity v 1 is positive then it first has to decrease to zero and then decrease further until it reaches the value − p 1 . Equation (28) and when s becomes zero we need to switch to the other branch of the solution. For more fields the distinction between the two case is more complicated because all velocities are mixed through the A i v i terms.

How many e-folds before inflation?
For the product-exponential potential the existence of an analytical solution allows us to calculate the number of e-folds the system spends before it approaches the asymptotic solution. For accelerating solutions specifically, p i should be chosen such that 1/2 p i p i < 1. If the system starts at a state with ε 0 > 1, then the solution (31) provides the necessary number of e-folds before inflation begins.
For an arbitrary potential, one can ask a similar question; assuming that there exists an "asymptotic" state (in practice the slow-roll solution) and that the system begins with initial conditions close to kinetic domination, then how many e-folds are necessary before inflation starts? The exact duration for the passage from a state with ε > 1 to a state with ε < 1 requires the full analytical solution, which seems an impossible task. To make some progress, we can construct an estimate based on the previously derived solution (31). This relies on the observation that for small field displacements any continuous function p i can be approximated by its constant part (the first term in the Taylor series) and this gives a first estimate on the total number of e-folds. However, since the value of p i changes we do not know a priori which value to choose. In the following we will show that in general the minimum and maximum values of p i give the minimum and maximum number of e-folds respectively. To prove this we will first derive some properties of the solutions for product-exponential potentials.
In Fig. 2 we plot the solution for three different exponential potentials (V p ∝ e p i φ i and V λ ∝ e λ i φ i ) satisfying 0 < p i < λ i and for positive (upper plot) or negative (lower plot) initial velocities that are close to kinetic domination. In the first case, v (i) p (i) , v (i) λ (i) > 0 holds for each i (with no sum implied), while in the second one v (i) p (i) , v (i) λ (i) < 0. In the upper plot, ε first decreases until the value ε * and then it increases again until it reaches its asymptotic value. During the decreasing period we observe that steeper gradients lead to faster decrease, whereas this behaviour is reversed once ε starts to increase. This makes sense physically because in the first case the Hubble friction cooperates with the gradient term (until the velocities change sign) and hence ε decreases faster for steeper gradients; similarly, when the velocities increase steeper gradients lead to larger terminal velocities and hence the system needs more time to reach those values. However, in both plots the system reaches its asymptotic state slower for steeper gradients.
In certain cases we can understand the previous behaviour analytically. For two different product-exponential potentials with 0 < | p i | < |λ i | we obtain 0 < |v i p i | < |v i λ i |. When each velocity v i has sign opposite to p i , λ i then and we obtain the following inequality Using equation (17) it is transformed into a differential inequality where we solved the differential equations with equal initial conditions ε(N 0 ) =ε(N 0 ). In the opposite case where velocities have all the same sign with p i , λ i then the previous inequality is reversed. In order for this inequality to be satisfied ε should adjust the rate of change according to the magnitude of the gradients. This is exactly what is depicted e.g. in the lower plot of Fig. 2; we observe that for the three sets of exponents ε(N ) follows the same inequalities and hence, schematically we have black > orange >blue. A different way to view the previous is that for a given set of p i we can find two different sets of exponents, e.g. λ i and κ i , whose solutions bound the solution for p i at all times. In the derivation of these relations the strict inequality between the different pairs of exponents is crucial; the condition on the norm of the gradient vector |κ| 2 < | p| 2 < |λ| 2 is not sufficient to provide a bound on ε at all times (see Fig. 2). The inequalities should instead be satisfied for each individual component. Moreover, we required that each velocity should have either the same sign as the corresponding gradient or opposite sign; if this is not satisfied this method may fail to provide definite inequalities for ε. This property of the solutions can provide some conclusions about the evolution of more general potentials. For field-dependent p i (φ j ), we can estimate the number of e-folds by considering the minimum and maximum values for p i (φ j ) in some Fig. 3 Evolution of 2 product-exponential potentials with ( p 1 , p 2 ) = (0.6, 0.6) and (λ 1 , λ 2 ) = (0.01, 0.9) (blue and green dashed accordingly) for negative supercritical velocities with the same ε 0 . Even though | p| 2 < |λ| 2 the two lines cross Fig. 4 Estimation of ε for a quartic potential using two exponential potentials with p min = 4/φ 0 and p max = 4/φ * for negative velocities, ε 0 = 2.99 and φ 0 = 10M pl field interval [φ j,min , φ j,max ]. Using these values we calculate N min and N max that correspond to the the solutions for product exponential potentials (e p i,min φ i and e p i,max φ i respectively) and then estimate the actual transitional time in the interval N ∈ (N min , N max ).
As an illustrative example we will study quartic inflation V ∼ φ 4 /4. We set the initial conditions to φ 0 = 10 and ε 0 = 2.99. The field flows to the origin and in the interval (0, φ 0 ) the gradient p = 4/φ is a monotonically decreasing function; therefore, the minimum value for p would be given by p(φ 0 ), whereas the maximum value corresponds to the value of φ which gives ε(φ * ) = 1. The latter value is unknown, as we do not know the full solution; in the following we will demonstrate how to estimate it. The value φ * would be the displacement one would find using the analytical solution (33) if the value p(φ * ) was used for the gradient term Note that the previous is a transcendental equation for φ * . Once this value is known, then p max = p(φ * ) can be used in (31) to find the maximum number of e-folds. The precision of this estimate is portrayed in Fig. 4 and we can predict a priori the field displacement between 2.9 and 3.2 in units of M pl , whereas the usual estimation resulting from considering kinetic domination with p = 0 gives φ = 2.35 M pl . This difference, although small, can potentially be important in the light of different conjectures which severely restrict the allowed superplankian field displacement, as well as the minimum value of p [46,47]. Conversely, at a given point φ 0 one can constrain ε 0 by the requirement of 50-60 e-folds of inflation and therefore estimate the "basin of attraction" of the initial conditions that provide sufficient inflation.

The gradient-flow approximation
For potentials with slowly varying gradients (in some field interval) the two estimated values for ε will be very close. This implies that the system behaves closely to a product exponential with field-dependent exponents, which, nevertheless, satisfy p i 1. Moreover, the scaling relations for the velocities v i ≈ −p i imply that ε varies slowly and the two terms should be small. In the latter equation we recognise the conditions for the validity of the gradient-flow approximation [48][49][50][51][52][53]. Therefore, models that are approximated well by gradient flow lie at the intersection between de Sitter and scaling solutions. The reason why these models asymptote to a state with almost vanishing accelerations can be explained by the fact that they are small deformations of de Sitter-scaling models, which have identically zero accelerations; since they are not exactly product-exponentials, during their evolution they move adiabatically between asymptotic solutions of exponentials. 4 In this way, one can deduce the global convergence towards the slow-roll solution which can not be done using linearised analysis. As a final remark, the scaling relations fix the fields' velocities, through v i ≈ −p i , but not their position. Therefore, under this type of approximation observables inherit an explicit dependence on (a subset of) initial conditions.

Generalizations
In addition to the scalar fields we will consider a noninteracting fluid with an equation of state P = wρ. The Klein-Gordon equations remain unchanged, while Eqs. (4) and (5) are modified and we obtain an additional equation, the conservation equation for the fluid, The slow-roll parameter is now given by Introducing a new variable z ≡ ρ/H 2 and using the system of evolution equations can be written as In this form one can immediately read off asymptotic solutions, generalizing the solutions for one field (see e.g. [55]): 1. scalar kinetic domination: z = 0, ε = 3. There exists an infinite number of field realizations describing kinetic domination, namely all points on the hypersphere satisfying v i v i = 6. 2. scalar potential domination: z = 0, v i = −p i with ε = 1/2 p i p i . 3. fluid domination: v i = 0, z = 3 with ε = 3/2(1 + w). 4. scaling solution: (at least one of) the fields and the fluid are non-zero with ε = 3/2(1 + w) (the expressions for v i and z are more complicated but can be found by solving quadratic equations).
In the case of a fluid with w = 1 5 equations simplify to 5 This fluid resembles the effective equation of state of a scalar field in the kinetic domination regime or in the absence of a potential (V = 0). However, the problem is not equivalent to a (N + 1)-field model, where the extra field is massless, due to the extra factor of 2 in Eq. (52) (see also [56]). and, hence, one obtains another integral of motion which can be used to express the fluid in terms of a reference field Plugging back into the slow-roll parameter (49) we obtain which has the same form as Eq. (24) after the substitution A 2 →Ã 2 = A 2 + B. Moreover, for w = 1 Eq. (17) remains unchanged and the general solution in terms of ε has exactly the same form (Eqs. (31) and (33)) modulo the substitution A 2 →Ã 2 .

Summary and discussion
Scalar fields in cosmology have been extensively studied over the last decades mainly due to their applications in both the early and the late universe. To date, only a limited number of exact solutions exists in the literature. In this paper we derived the general solution for an interacting product-exponential potential, that has not been considered in earlier works, in terms of the slow-roll parameter ε. In the latter form it better highlights properties of its critical points. We argued how the solution becomes relevant in the derivation of estimates for the transitional time to inflation for models that follow the potential gradient flow. For the nonlinear evolution towards the slow-roll solution, i.e. starting with arbitrarily large deformations, one needs the full solution. Since there is no general prescription on how to solve the Klein-Gordon equation, one can alternatively derive estimates for the time of this transition using the solution of exponential potentials. This method generalizes the existing estimates, as it does not assume kinetic domination or negligible potential gradients K V, V ,φ (that was considered in previous works, e.g. [57]). We showed that the actual transitional time will lie between two bounds that depend on magnitude of the potential gradient (ln V ) ,φ . Besides providing more accurate results the method can be applied in multifield models where linearized stability becomes obscure due to the coupling of fields' perturbations.
It is worth noticing that the general solutions presented in this work can also describe multi-field domain walls. The correspondence between cosmological and domain wall solutions can be found e.g. in the Refs. [58][59][60][61]. The beta function is the analogue of the slow-roll parameter ε and hence the RG flow of multi-field domain walls with product exponential potentials can be found using the solutions (31) and (33) and the dictionary between cosmology and domain walls.
Finally, we considered the case of an extra non-interacting fluid. When the fluid obeys an equation of state w = 1 (stiff matter) another integral of motion exists which makes the problem again fully integrable with almost the same solution as before (more precisely with a minor substitution). 6 It would be interesting to investigate how this solution can be generalized for problems with non-trivial field manifolds.