Asymptotic Solutions of a Generalized Starobinski Model: Kinetic Dominance, Slow Roll and Separatrices

: We consider a generalized Starobinski inﬂationary model. We present a method for computing solutions as generalized asymptotic expansions, both in the kinetic dominance stage (psi series solutions) and in the slow roll stage (asymptotic expansions of the separatrix solutions). These asymptotic expansions are derived in the framework of the Hamilton-Jacobi formalism where the Hubble parameter is written as a function of the inﬂaton ﬁeld. They are applied to determine the values of the inﬂaton ﬁeld when the inﬂation period starts and ends as well as to estimate the corresponding amount of inﬂation. As a consequence, they can be used to select the appropriate initial conditions for determining a solution with a previously ﬁxed amount of inﬂation.


Introduction
The inflationary cosmology is an important ingredient for the description of the early universe [1][2][3]. Many inflationary models have been considered which are formulated in terms of an homogeneous spatially flat Friedman-Lemaître-Robertson-Walker (FLRW) spacetime with scale factor a(t), and a time dependent single real field φ(t), known as the inflaton field [4][5][6]. These models are described by the nonlinear ordinary second order differential equationφ where V = V(φ) is a given potential function and H =ȧ/a is the Hubble parameter which depends on the inflaton field according to the equation Here m Pl = √h c/8πG is the Planck mass and dots indicate derivatives with respect to the cosmic time t.
One of the most interesting inflationary models is the Starobinski model [1,7], which obtains from a modification of the Einsten-Hilbert action by an extra curvature quadratic term. It is associated with the potential function with Λ being constant. Several modifications and generalizations of this model have been proposed, such as the α-attractor models [8][9][10], leading to power law inflation [11], or models with exponential power law inflation in modified gravity [12]. Some of these extensions (see for example [10] or [12]) can be written as particular cases of the potential where Λ j , j = 1, 2, 3, λ and µ are positive coefficients with λ > µ. In particular, by setting in µ = λ 2 = 1 3 , Λ 1 = Λ 3 = Λ 2 2 = Λ 4 , the potential (4) reduces to (3). Figure 1 shows the graphic of the potential (4) for several values of the parameters.  (3)) and µ taking the values 1 3 (black line, also corresponding to the Starobinski model), 0.65 (red line) and 0.15 (blue line). The axes origin is (φ min , V min ) with φ min being the inflaton value at which the potential takes its minimum and V min := V(φ min ).
In this work we use the Hamilton-Jacobi formalism to determine asymptotic solutions of the model (4) for two different situations: the kinetic dominance (KD) and the slow-roll (SR) stages.
The KD stage [13][14][15] is the period when the kinetic energy of the inflaton field dominates over its potential energy.φ 2 V(φ).
It is a non-inflationary or pre-inflationary stage that is followed by a short fast-roll inflation phase [16] and afterwards by the traditional SR inflation stage [5,[17][18][19][20][21][22][23][24]. Recently, Handley et al [14,15,[25][26][27] have shown the relevance of the KD period (5) for setting initial conditions. Indeed, as it was stablished in [14,16], for models with a smooth and positive potential V the Hubble parameter is a positive monotonically decreasing function of t. As a consequence, the functions φ(t) and H(t) corresponding to arbitrary finite initial data are not singular forward in t. Nevertheless they may develop singularities backwards in t. Thus, it can be proved [14] that for generic initial conditions the evolution backwards in time starts by a KD regime followed by an inflationary regime. The solutions of Equation (1) manifest generically branch point singularities of logarithmic type related to the condition (5), and their presence is associated with the so-called psi-series [28] asymptotic solutions. Handley et al. developed in [15] a general method to compute psi-series expansions for the solutions of Equation (1) and the generalization of Equation (2) for FLRW spacetimes with curvature. They referred to their series as logolinear series. As it was discussed in Section 2-C of [16], the mathematical singularity t = t * does not mean a real physical phenomenon. It is an extrapolation of the classical treatment of the inflaton models which is outside the validity region of the classical treatment. At this point quantum cosmology has to be used.
In a previous work [29] we have presented an alternative method for computing expansions of the solutions of (1)- (2) in the KD stage for wide classes of potentials. The starting point of the method is the Hamilton-Jacobi formalism, we look for solutions of the Hamilton-Jacobi equations as psi-series in e ϕ (ϕ is the rescaled inflaton field ϕ ≡ 3 2 φ m Pl ) with polynomial coefficients. The potential (4) is not a particular case of the models studied in [29], nevertheless we extend here our previous method by considering polynomial coefficients depending on the two variables e −2λϕ and e −2µϕ . We refer to these series as expolinear series. At this point we need to impose the exponents λ, µ and the quotient λ µ to be irrational numbers, however a simple limit operation shows that the results also apply to rational exponents.
On the other hand, the SR stage is the period in which the potential energy dominates over the kinetic energyφ It is known [19,30] that many inflationary models exhibit separatrix solutions, these solutions behave as "attractor solutions" in the sense that they determine the asymptotic behaviour of all the solutions of the model as the cosmic time t is large enough. In particular they provide us with asymptotic expansions in the Hamilton-Jacobi formalism for the SR stage. In [31], we have studied the existence of separatrix solutions and have discussed their asymptotic expansions for wide classes of inflationary models. The results of [31] imply that the potential model (4) has a separatrix solution defined for φ > φ min for all values of the parameters λ and µ, however, the existence of a separatrix solution defined for φ < φ min requires that the condition λ < 1 is fulfilled. The paper is organized as follows: In Section 2, we briefly introduce the Hamilton-Jacobi formalism for the model (4). Section 3 describes the method for determining psiseries solutions for the Hamilton-Jacobi equations. We defer to Appendix A the discussion of how to derive from our psi-series in the inflaton field ϕ, the logolinear series for ϕ(t) involving integer and irrational powers of (t − t * ). We devote Section 4 to obtain asymptotic expansions of the separatrix solutions. The question of how to use the series solutions obtained in Sections 3 and 4 to provide suitable approximations to the Hamilton-Jacobi equations in the KD and SR stages is discussed in Section 5.1.
These approximations are applied in Section 5.2 to provide analytical approximations of several relevant quantities. Note that we need to take the parameter Λ 3 in (4) as so that it is satisfied that V min = 0 and consequently, the model exhibits solutions which leave the inflation period and overlap a reheating phase [5,6,21,[32][33][34]. Thus, we use the psi-series to determine the value of the inflaton field at the initial moment of the inflation period, and the separatrices to obtain the value of the inflaton field at the end of the inflation period. Using both approximations we provide a formula for the amount of inflation. The paper ends with a summary of conclusions in Section 6.

Hamilton Jacobi Formulation
Let us begin with the Hamilton Jacobi formulation of the inflationary models (1) and (2). We use the rescaled variables and rewrite equations (1) and (2) as and with c j = 3 m 2 Pl Λ j , j = 1, 2, 3. Furthermore, due to (7) the parameter c 3 is defined as In this work we use Plank units (G = c =h = 1). Thus, for the Starobinski model (3) we have that (see for example [15]) Λ 4 = 10 −10 m 2 Pl , so the coefficients c j , j = 1, 2, 3 are of order 10 −10 .

The Hamilton-Jacobi Equations
In the Hamilton-Jacobi formulation of inflationary models the reduced Hubble parameter is determined as a function h(ϕ) of the inflaton field, consequently, we consider two subsets in the phase space (ϕ,φ) whereφ has a constant sign.
The diffeomorphims and enable us to describe the dynamics of (9)-(10) on the regions R ± of the (ϕ, h) plane (see Figure 2 below).
Thus, the parts of each solution ϕ(t) of (9) in D ± satisfy the equations Here primes denote derivatives with respect to ϕ and the reduced Hubble function h is assumed to be the positive root Equations (17) and (18) are referred to as the Hamilton-Jacobi formalism for inflationary models [6,19,35,36].
The sets R ± plays the role of the phase space of the formalism. Thus each solution h = h(ϕ) of (17) determines a corresponding implicit solution ϕ(t) of (18) given by

The KD Period. Asymptotic Solutions
As a consequence of (16), it follows that the reduced Hubble parameter h is a positive monotonically decreasing function of t. This property implies that for smooth and positive potential functions v, the solutions ϕ(t) of (9) with arbitrary finite initial data do not have singularities forward in the cosmic time t. Nevertheless, the function h(t) increases without bound backwards in time, so that we may expect that h(t) and ϕ(t) may develop singularities. Thus, if the KD condition (5) holds then we may neglect v and v in the inflaton equations and from (9) we havë and we obtain two families of approximate solutions where t * and ϕ p are arbitrary constants. Consequently, the asymptotic form of the reduced Hubble parameter (19) is These approximate solutions of the inflaton equations are the dominant terms of the psi series expansions that we will consider below. The corresponding approximate solutions of (17) take the form where b is an arbitrary strictly positive parameter (b = e ±ϕ p ). We will denote by h (+) (ϕ) the solutions of (17) which emerge from the KD region with asymptotic behaviour Note that h (+) (ϕ) is indeed solution of the equation We also use the notation ϕ (+) (t) for the associated solution of (18) which blows up at a finite time t * given by and denote h (+) (t) = h (+) (ϕ (+) (t)). The asymptotic forms of ϕ (+) (t) and h (+) (t) near the blow up time are given by and Analogously, we denote by h (−) (ϕ) the solutions of (17) which emerge from the KD region with asymptotic behaviour We point out that only if λ < 1 the Equation (17) with potential (11) may admit solutions with asymptotic behaviour (30). Otherwise the potential function v(ϕ) is a dominant term with respect to h(ϕ) 2 and h (ϕ) 2 in (17) as ϕ → −∞. We assume that λ < 1, then it is clear that h (−) (ϕ) is a solution of the equation Now ϕ (−) (t) stands for the associated solution of (18) which blows up at a finite time given by and h (−) (t) = h (−) (ϕ (−) (t)). The asymptotic behaviours of ϕ (−) (t) and h (−) (t) near the blow up time are given by and

Slow-Roll Stage and Separatrix Solutions
In our previous work [31] we proved that for a potential v(ϕ) such that the Equation (17) has a unique solution h (+) Furthermore, using the symmetry (17) and (18), we have that provided the Equation (17) has a unique solution h The trajectories in the (ϕ,φ) phase space associated to these solutions are boundaries of the regions filled by the solutions with asymptotic behaviour given by either (28) or (33) (see Figure 3 below). Moreover, h s (ϕ))) is the boundary in the phase spaces R + (resp. R − ) between solutions of (17) defined for all ϕ > ϕ min (resp. ϕ < ϕ min ) and solutions of (17) which leave R + for a certain ϕ * > ϕ min (resp. which leave R − for a certain ϕ * < ϕ min ). These important solutions are referred to as separatrix solutions and it turns out that, for wide sets of initial conditions, the solutions of the model reduce asymptotically to these special solutions [19,30]. Thus, they provide us with accurate approximations to the solutions in the SR stage.
For the inflationary model with the generalized Starobinski potential (11) we have that Consequently, a separatrix solution arises with asymptotic behaviour Furthermore, provided that λ < 1 there is another separatrix solution such that

The Amount of Inflation
The inflation period of the universe evolution is characterized by an accelerated universe expansionä > 0. From the identity it follows that this period is determined by the constrainṫ Then, the inflation regions (43) in R ± are characterized by It is well known that a successful solution to the horizon and flatness cosmological problems requires that the amount of inflation should be close to N ∼ 60 [6,[37][38][39]. Thus, given a solution h = h(ϕ) of (17) it is essential to determine the values ϕ in , ϕ end for which inflation starts and ends. Now, due to (44) both values satisfy As we will show in Section 5.2, we obtain faithful approximations of ϕ in if we use KD approximations for h(ϕ) in (46) supplied by appropriately truncated psi-series expansions. Moreover, ϕ end is accurately approximated in terms of the SR approximations determined through separatrix expansions. Both approximations will be applied to determine the amount of inflation N as a function of the free parameter b.

Expolinear Series for the Reduced Hubble Parameter h(ϕ) in the KD Period
We consider the KD period and look for solutions of (17) with asymptotic behaviour (25). Thus we substitute the psi-series expansion in (17), where γ (+) n (ϕ) are polynomials on c 1 e −2 λ ϕ and c 2 e − 2 µ ϕ . By inserting (47) into (17) we obtain where we have introduced the variables Henceforth, we will assume that λ, µ and λ µ are irrational numbers.
Nevertheless a simple limit operation shows that the results below also apply to rational exponents, or irrational exponents λ and µ such that λ/µ is rational. From (50) we have that the powers of x, y and e −ϕ are linearly independent functions; consequently we may identify the coefficients of e −2nϕ for n = 0, 1, . . . on both sides of Equation (48). Thus, for n = 0 we obtain that γ verifies the linear first order PDE which can be solved by the characteristic method [40] and, after the change of variables it reduces to the linear first order ordinary differential equation where P 1 is an arbitrary function. In terms of the variables (x, y), the general solution can be written asγ Thus, Equation (51) has a unique polynomial solution given by Furthermore, identifying the coefficients of e −2nϕ , n ∈ N in both sides of Equation (48) leads to the recursion relation The Equations (54) are nonhomogeneous linear first order PDE, whose nonhomogeneous terms depend on the coefficients γ (+) j with j = 1, . . . , n. We now use induction to prove that the coefficients γ By using the induction hypothesis, we can write the right-hand side of Equation (54) as where the coefficients q l,j,k for l = 1, . . . , n and j, k such that j + k ≤ l. The change of variables (52) transforms (54) and (57) into with general solution Here P n+1 denotes an arbitrary function on ξ. Thus, in terms of the variables (x, y) the unique polynomial solution of Equation (54) is given by For example, we have that the first Equation (54) are and the second coefficient in (47) is given by Recalling Equations (47), (49) and (55) we get that the reduced Hubble parameter is On the other hand, if λ < 1 we can use the symmetry (17) and (18) to find that the psi-series corresponding to solutions with asymptotic behaviour (30) are given by where γ (−) n,j,k = γ (+) n,j,k λ→−λ,µ→−µ n ∈ N, j = 0, . . . , n, k = 0, . . . , n − j.

Asymptotic Series for Separatrix Solutions
In this section we determine asymptotic expansions for the separatrix solutions of Equation (17). According to (36) where the dots stand for lower order terms as ϕ → ∞. We next prove that the asymptotic expansion for the separatrix solution is of the form where the coefficients S n,l can be recursively determined. Equivalently, in terms of the variables (x, y) (49) we have Substitution of (66) into (17) leads us to the equation Since the powers of x and y are linearly independent, we can identify the coefficient of each monomial x n y l in both sides of Equation (67). For example, by setting (n, l) = (0, 0), (1, 0) and (0, 1), it follows that which is in agreement with (64). On the other hand, if (n, l) / ∈ {(0, 0), (1, 0), (0, 1)}, equating the coefficient of x n y l in both sides of (67) provides us with As a consequence, the Equations (68) and (69) supply a recursion formula to determine all the coefficients in (65). Analogously, if we assume that λ < 1, taking into account the asymptotic behaviour (41) of the separatrix solution h where the dots stand for lower order terms. Then, we look for an expansion of the form where the coefficients r n (y) are polynomials in y. Substitution of (71) into (17) leads us to the equation As the powers of x and y are linearly independent, we can equate the coefficient of each power of x in both sides of (72). Thus, from the coefficient of x 0 we have that The only polynomial solution of (73) is given by The coefficient of x 1−n with n > 1 implies the equation We now use induction in n to prove that equations (75) determine recursively r n (y) as a polynomial of degree at most n. First, from (74) we have that r 1 (y) is a polynomial of degree 1 in y. Next, we have that the right hand sides of equations (75) depend only on r j j = . . . , n − 1, moreover using the induction hypothesis, these terms are polynomials in y of degree less or equal than n. Consequently r n (y) is determined by (75) as a polynomial of degree at most n. For example, for n = 2 (75) takes the form 2 λ µ r 2 − (3 λ 2 + 1) r 2 = −2 µ 2 y 2 (r 1 ) 2 + 2 λ µ y r 1 r 1 whose polynomial solution is In general, it can be written that r n (y) = n ∑ j=0 r n,j y j , where the coefficients r n,j depend on λ, µ and c 3 . In this way, the expansion of the separatrix solution h Note that as µ < λ the exponents in (79) satisfy (2n − 1)λ − 2jµ > −λ so the terms in the sum are subdominant terms with respect to e −λϕ as ϕ → −∞.

Aplications
In this section we apply the expansions obtained in Sections 3 and 4 to get suitable approximations to the solutions of the Hamilton-Jacobi Equation (17) with v(ϕ) given by (11) in both the KD and the SR stages, as well as to the determination of the duration of the inflation period and the amount of inflation.

Approximate Solutions
We define the m-order approximation h s (ϕ)) as the truncated series (61) (resp. (65)) in which we keep all the terms in e −cϕ with c ≤ m − 1 and remove the terms with c > m − 1. Thus, as the expansion (61) for h (+) (ϕ) contains terms in e (−2jλ−2kµ−2n+1)ϕ with j and k being greater than or equal than zero, only the terms with n ≤ [ m 2 ] determine the m-order approximation. Here [.] stands for the integer part of a real number. More precisely, taking into account that 0 ≤ j, k ≤ n and j + k ≤ n it follows that where I (+) m is the finite set of indeces Analogously, the expansion (65) for h  Figure 4 shows the 6-order approximations for h (+) (ϕ) and h (+) s (ϕ) together with the numerical solution. We observe that h (+) approx,6 (ϕ) represents a very good approximation in the KD region while h (+) s,approx,6 (ϕ) gives a very good approximation in the SR region. To get suitable approximations of h (−) (ϕ) and h (−) s (ϕ) for λ < 1, we can proceed in a similar way. Thus we define the m-order approximation h s,approx,m (ϕ) to h s (ϕ)) as the truncated series (62) (resp. (79) in which we keep the terms in e cϕ for c ≤ m − 1 and remove the terms with c > m − 1. The expansion (62) for h (−) (ϕ) contains terms in e (−2jλ−2kµ+2n−1)ϕ for n = 1, . . . , j, k, ≤ n, thus by recalling that 0 < µ < λ < 1, it is clear that and consequently, only terms with n ≤ In order to compare the approximations h In terms of h = h(ϕ), the Equation (31) reads where The Equation (89) allows us to avoid the numerical instabilities due to the exponential growth of the solutions of (31) as ϕ → − ∞. Figure 5 shows the graphics of h (−) approx,6 (ϕ, b)/ v(ϕ) and h (−) s,approx,6 (ϕ)/ v(ϕ) together with the numerical solution h of (89). We point out that h (−) approx,6 (ϕ, b)/ v(ϕ) provides a very good approximation in the KD region while h (−) s,approx,6 (ϕ)/ v(ϕ) provides a very good approximation in the SR region.  It is worth noticing that in terms of the modified Hubble parameter the inflation region is given by and as we will show next, it takes place for ϕ along an interval about (−100, 0). Due to the exponential growth of v(ϕ) as ϕ → −∞, it is not possible to appreciate the inflation region (44) in the (ϕ, h) plane.

Aplications to the Inflation Period
Let us now look for appropriate approximations to the values ϕ in (b) and ϕ end of the inflation period as well as to the amount of inflation N(b). For a solution h (+) (ϕ) of (17) satisfying (25), its m-order approximations h We denote the solutions of (91) and (92) by ϕ in,m (b) and ϕ end,m , respectively. Then, to calculate N(b) we need a value ϕ * m (b) such that the KD approximation holds for the interval ϕ ∈ (ϕ * (b), ϕ in,m (b)), while the SR approximation holds for the remaining part of the inflation period. To this end, we choose the value of ϕ for which the function s,approx,m (ϕ) 2 , reaches its minimum. Thus, we can estimate the amount of inflation as a function of b as The approximation (93) can be used to select an appropriate value of the parameter b and, consequently, to determine the initial conditions, corresponding to a solution of Equation (17) with a previously fixed value of N. Figures 6 and 7 illustrate the application of this method. The values for the parameters of the potential are the same as those chosen in Section 5.1 and we are using again the six-order approximations. In particular, for N(b) = 60, we obtain b 60 ≈ 2.28186 × 10 7 while for N(b) = 70, we obtain b 70 ≈ 2.80532 × 10 7 . In Figure 7, we plot the numerical solution of (17) with initial conditions ϕ 0 = 12, h 0 = h (+) approx,6 (ϕ 0 , b 60 ) (left) and ϕ 0 = 12, approx,6 (ϕ 0 , b 70 ) (right). Numerical computation of the amount of inflation for the solution leads to N num ≈ 59.044 and N num ≈ 68.995 respectively. In the case λ < 1 we can proceed in the same way for solutions h (−) (ϕ) of (17) with asymptotic behaviour (30). For these solutions ϕ in,m (b) and ϕ end.m (b) are taken as the solutions of the equations The value ϕ * m (b) is taken as the minimum of the function s,approx,m (ϕ) 2 , and the amount of inflation is estimated by The gray dots correspond (from left to right) to ϕ in,6 (b 60 ), ϕ * 6 (b 60 ) and ϕ end,6 (left figure) and ϕ in,6 (b 70 ), ϕ * 6 (b 70 ) and ϕ end,6 (right figure).

Conclusions
We have analysed the asymptotic properties of the solutions of a generalization of the Starobinski inflationary model determined by the potential (4). We have considered both the kinetic dominance and the slow roll periods within the Hamilton-Jacobi formalism, which provides a natural framework to determine generalized asymptotic expansions of the Hubble parameter as a function of the inflaton field.
We have found that the Hamilton-Jacobi Equation (17) in the KD period admits asymptotic expansion solutions as ϕ → ∞ which can be written as series in the variable e −2ϕ with polynomial coefficients in the variables e −2λϕ and e −2µϕ . Moreover, if λ < 1 we have proved that Equation (17) in the kinetic dominance period admits asymptotic expansions solutions as ϕ → − ∞, which can be written as series in e 2ϕ with have also polynomial coefficients in e −2λϕ and e −2µϕ . These asymptotic expansions depend on a free parameter b, and the polynomial coefficients can, in both cases, be recursively obtained (see Equation (58)), and indeed this recursion formula is easily programmed.
On the other hand, the solutions of (17) in the SR period behaves as separatrix solutions [19,30], when these separatrix solutions exist. Using the results of a previous work [34] we know that the inflation model (4) admits separatrix solution defined for ϕ > ϕ min for all λ > 0 (we have denoted this separatrix solution by h (+) s (ϕ)), and admits also a separatrix solution defined for ϕ < ϕ min provided that 0 < λ < 1 (we have used here the notation h s (ϕ) can be written in powers of e 2λϕ with polynomial coefficients in e −2µϕ . Again the coefficients in these expansions can be recursively determined and the recursive method is easily programmed.
All these expansions prove to be useful to obtain approximations to the solutions of (17) in the KD and SR periods. Thus given m ∈ N only a finite number of terms in the expansions as ϕ → ∞ (resp. ϕ → −∞) are proportional to e −cϕ (resp. e −cϕ ) with c ≤ m − 1.
In particular, we have considered the model for the value of the parameter Λ 3 in (4) corresponding to V(φ min ) = 0. The corresponding expansions have been applied to determine the values of the inflaton field when the period of inflation starts ϕ in and ends ϕ end , as well as to compute the amount of inflation as a function of the free parameter b in the KD period. Thus, this function can be used to select the appropriate value of the free parameter, and consequently the appropriate initial conditions, for a solution with a previously fixed amount of inflation.
Finally, we mention that many of the results in this work could be easily formulated for models with an exponentially potential well with 0 < λ, µ < 1 and Λ 1 , Λ 2 , Λ 3 > 0 (see for instance the models consider in [41]).

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
Then, for n ∈ N and j, k such that j + k ≤ n we have that where we are introducing the vectorial functions Substitution of (A4) and (A6) into (A1) leads to n,j,k z j w k C r (Φ (+) r,n,j,k )(t − t * ) 2r+2n−1 . (A7) Since we assume that λ, µ and λ µ are irrational numbers, the powers of z, the powers of w and the powers of (t − t * ) 2n are linearly independent functions, so that the coefficients of (t − t * ) 2n for n = 0, 1, . . . on both sides of Equation (A7) must be equal. Consequently, ϕ n satisfies the first order linear PDE s,j,k z j w k C n−s (Φ (+) n−s,s,j,k ).

(A8)
As the right-hand side of (A8) depends only on ϕ j , j = 1, . . . , n − 1, this equation shows that the coefficients ϕ n can be recursively determined. Furthermore, using induction and applying the method of characteristics we can prove that the coefficients ϕ (+) n (z, w) can be determined as polynomials in the variables z and w of degree at most n. To prove this, let us start with the first Equation (A8) which in terms of the variables reduces to the ODE Its general solution is given bŷ with R 1 being an arbitrary function. Taking into account (56) and coming back to the variables (z, w) we have that the general solution of (A9) iŝ so that the only polynomial solution takes the form If we now use the induction hypothesis in (A8) we find that its right-hand side can be written as a polynomial in the variables (z, w) of degree n, i.e. ∑ j+k≤n,j,k≥0 s (+) n,j,k z j w k .