Parametric Gevrey asymptotics in two complex time variables through truncated Laplace transforms

The work is devoted to the study of a family of linear initial value problems of partial differential equations in the complex domain, dealing with two complex time variables. The use of a truncated Laplace-like transformation in the construction of the analytic solution allows to overcome a small divisor phenomenon arising from the geometry of the problem and represents an alternative approach to the one proposed in a recent work [9] by the last two authors. The result leans on the application of a fixed point argument and the classical Ramis-Sibuya theorem.


Introduction
This work is devoted to the study of a family of singularly perturbed partial differential equations in the complex domain of the form (1) Q(∂ z )u(t 1 , t 2 , z, ) = P (t k 1 +1 1 ∂ t 1 , t 2 , ∂ t 2 , ∂ z , z, ) + f (t 1 , t 2 , z, ), under initial data u(0, t 2 , z, ) ≡ u(t 1 , 0, z, ) ≡ 0, with Q(X) ∈ C[X] and P (T 1 , T 21 , T 22 , Z, z, ) being a polynomial in (T 1 , T 21 , T 22 , Z) with holomorphic coefficients w.r.t. (z, ) on H β ×D(0, 0 ). Here, H β and D(0, 0 ) stand for the horizontal strip {z ∈ C : |Im(z)| < β} and the disc at the origin and radius 0 , for some β > 0, 0 > 0, respectively. The forcing term f (t 1 , t 2 , z, ) is holomorphic on C × D(0, h ) × H β × E, for any open sector E centered at 0 and contained in D(0, 0 ), for some h > 0, and remains close to a polynomial in t 1 , analytic in t 2 on D(0, h ) and in z on H β , as becomes close to the origin in C. The variable acts as a small complex parameter. The concrete assumptions on the elements involved in the main problem (1) are to be described and analysed in detail throughout the work. The study of a problem of such form is motivated by the recent research [9] of the second and third authors. The main aim in the preceeding work was related to the description of the asymptotic behavior of the analytic solutions, with respect to the perturbation parameter, near the origin, of singularly perturbed equations (2) Q(∂ z )u(t 1 , t 2 , z, ) =P (t k 1 +1 1 ∂ t 1 , t k 2 +1 2 ∂ t 2 , ∂ z , z, ) +f (t 1 , t 2 , z, ), withP (T 1 , T 2 , Z, ) being a polynomial in (T 1 , T 2 , Z, z, ) with holomorphic coefficients w.r.t. (z, ) on H β × D(0, 0 ). Two main novelties are considered here with respect to it. On the one hand, the irregular singular operators related to the second time variable stay rigid in (2), as a polynomial function of the operator t k 2 +1 2 ∂ t 2 . In the present study, the irregular operators in this variable fit a more general scheme within the problem, under certain technical assumptions (see (5) and (6)). This, at first sight slight, variation on the form of the main problem varies its underlying geometry radically. On the other hand, the appearance of different types of solutions observed in [9], known as inner and outer solutions, which describe boundary layer expansions do not appear in the present situation, since we study local solutions in time t 1 , t 2 near the origin in the complex domain. It is worth mentioning that, despite the fact that the form of the main equation under study resembles that of [9], the nature of the singularities appearing in the problem require to appeal different approaches and apply novel techniques, to be briefly described below.
This work continues a line of research on the study of the asymptotic behavior of solutions of singularly perturbed PDEs in the complex domain, under the action of two time variables: dealing with a symmetric factorized (resp. asymmetric) leading term [8] (resp. [6]), the mentioned work [9], and the corresponding q−analog [10] in the framework of q−difference−differential equations.
The technique developed in the present work consists on searching for solutions of the main problem ( see (8) for its precise expression) in the form of a Fourier, truncated Laplace and Laplace transform of certain function, for every fixed value of the perturbation parameter : (3) u(t, z, ) := The integration path L 1, stands for the segment [0, h 1 ( )e √ −1θ 1 ], for some holomorphic function → h 1 ( ) on the domain of definition of the perturbation parameter, approaching to infinity when tends to 0, and some θ 1 ∈ R which does not depend on . The integration with respect to the path L 2 stands for a usual Laplace transform along certain half line [0, ∞)e √ −1d 2 , for some d 2 ∈ R, whereas the function ω belongs to certain Banach space which depends on the choice of the perturbation parameter. From this point, the main problem is replaced by an auxiliary convolution problem (see (17)) in the Borel plane with respect to the time variables (t 1 , t 2 ). The precise knowledge on the geometry of the auxiliary problem is crucial in order to understand the location and control of the singularities (see Section 3). As a matter of fact, the singularities of the auxiliary problem are always located outside (but remain close to) a product of discs in the domain of the transformed time variables, say (τ 1 , τ 2 ). The radius of the discs depend on tending to infinity regarding τ 1 and shrinking to the origin with respect to τ 2 . In addition to this, one can choose narrow finite (resp. infinite) sectors with vertex at the origin with respect to τ 1 (resp. τ 2 ), valid for all the values of the pertubation parameter, which avoid the singular points. In other words, the function ω can be extended to the product of such sectors w.r.t. (τ 1 , τ 2 ). Therefore, a small denominator problem regarding movable singularities to infinity and to zero at the same time (in each of the time variables) has to be analysed. This construction through a truncated Laplace transform is proposed in order that the solutions (3) remain close, as tends to 0, to a double usual Laplace transform in both variables t 1 , t 2 . For such a complete double Laplace representable solution, a direct analysis of the asymptotic behaviour w.r.t. is unfortunately not possible (as shown in our previous work [9]). However, such study turns achievable within the new approach regarding truncated Laplace transform solutions.
The use of truncated Laplace transform with respect to one of the time variables in the Borel plane is used successfully to control the growth of the solutions, via complex Banach spaces of functions not only subject to an exponential growth in the monomial variables, but also whose domain of definition depends on each value of (see Section 4). Given a finite family of finite sectors E = (E p ) 0≤p≤ι−1 which conform a good covering (see Definition 3), the first main result in the work, Theorem 1, states the existence of a solution of the main problem in the form (3) for every 0 ≤ p ≤ ι − 1, remaining holomorphic in a domain T 1 × T 2 × H β × E p , where T 1 , T 2 are finite sectors with vertex at the origin. Moreover, the exponential decrease of the difference of two solutions associated to consecutive sectors in E enables the application of the classical Ramis-Sibuya theorem (RS) in order to achieve the second main result of our study, namely the asymptotic relation of the analytic solutions and the formal solution of the main problem in powers of , with coefficients in some complex Banach space (see Theorem 3).
In recent years, several steps have been taken to contribute to the knowledge of the asymptotic behavior of analytic solutions of singularly perturbed partial differential equations in the complex domain. We first refer to the recent works [17,18], by H. Yamazawa and M. Yoshino, and M. Yoshino resp. in which the parametric Borel summability of semilinear systems of PDEs is studied, first in the case of fuchsian operators, and second combining both irregular and fuchsian operators. We refer to [1,14] as introductory texts on the classical theory of summability of formal solutions of differential equations in the complex domain.
The appearance of truncated Laplace transform is closely related to the classical theory of asymptotic approximation of analytic functions (examples of this situation is the classical proof of Ritt's Theorem for Gevrey asymptotics, see [1] Proposition 10, and also Lemma 1.3.2 in [14]). Truncated Laplace transform also appears as a recent object of study in the literature, related to differential operators [12,13], but also from the numerical point of view [11]. The choice of an integration path for Laplace transform which depends on each fixed value of the perturbation parameter has been inspired from [3,15].
Throughout the work, we use bold letters to indicate a vector of two variables: we write τ for the pair (τ 1 , τ 2 ), u for (u 1 , u 2 ), T for (T 1 , T 2 ), etc.
The paper is organized as follows.
In Section 2.1, we recall some properties on Fourier transform which allow to transform the main problem, stated in Section 2.2, in the form a convolution problem, described in Section 2.3. The geometry of the problem is an important matter in this work, which needs to be explained in detail. Section 3 is focused on this issue. The Banach spaces involved in the construction of the analytic solution of the auxiliary problem, and some of their main properties, are stated in Section 4. Such function is constructed in Section 5. The analytic solution of the main problem is obtained in Section 6 (Theorem 1), and the work concludes with the description of the parametric Gevrey asymptotic expansions of the analytic solution, obtained in Section 7 (Theorem 3).

Layout of the main and auxiliary problems
In this initial section, we describe in detail the main problem under study (8) (Section 2.2), and the conditions on the elements involved in it. The solution of this problem is reduced to a convolution auxiliary problem in the Borel plane (17) when inspecting solutions in the particular form of a triple Fourier, Laplace and truncated Laplace transform (see Section 2.3). We first give some words about inverse Fourier transform on certain Banach spaces which act on the transformation of the problem (Section 2.1).

Inverse Fourier transform on certain function spaces
The transformation of the main problem with respect to variable z requires recalling some basic facts about inverse Fourier transform when acting on certain Banach spaces of real functions of exponential decrease at infinity.
The next result will be needed in our reasoning. We refer to [4], Proposition 7, for its proof.
Proposition 1 Let β > 0 and µ > 1. The inverse Fourier transform satisfies the following properties acting on every f ∈ E (β,µ) : • The function F −1 (f ) is well defined in R and can be analytically extended to the set • Let g ∈ E (β,µ) and let ψ(m) = 1 (2π) 1/2 f * g(m) be the convolution product of f and g, for all m ∈ R. Then, ψ ∈ E (β,µ) and it holds that
We assume that . We assume that Remark: In Section 3 we assume further geometric conditions on these polynomials. In particular, observe that condition (18) We choose µ ∈ R with µ > max The main aim in this work is to study the following initial value problem: for the initial conditions u(t 1 , 0, z, ) ≡ u(0, t 2 , z, ) ≡ 0. Let us describe the form of the elements involved in the problem. Let 0 > 0 and β > 0. For all 1 ≤ 1 ≤ D 1 − 1 and 1 ≤ 2 ≤ D 2 − 1, the term c 1 2 (z, ) are holomorphic functions on H β × D(0, 0 ). We recall that H β stands for the horizontal strip (4). The function c 1 2 is defined by where m → C 1 2 (m, ) is continuous for m ∈ R and is subject to uniform exponentially flat upper bounds with respect to ∈ D(0, 0 ), i.e. there exists C c > 0 such that (9) sup Observe that m → C 1 2 (m, ) belongs to E (β,µ) with sup ∈D(0, 0 ) The forcing term f (t 1 , t 2 , z, ) is a holomorphic function in C × D(0, h ) × H β × E, for any given open sector E centered at 0, and contained in D(0, 0 ) \ {0}, for some positive number h .
The forcing term is constructed as follows. Let N 1 ≥ 0 and F n 1 ,n 2 (m, ) ∈ E (β,µ) under uniform bounds with respect to in the disc D(0, 0 ). More precisely, assume that which turns out to be a holomorphic function on C 2 with respect to the first two variables, with coefficients in E (β,µ) . We write where κh 1 ( ) is a holomorphic function on any open sector centered at 0 in the punctured disc D(0, 0 ) \ {0} (see (20)), θ 1 is a real number to be determined and γ(n, z) stands for the incomplete Gamma function which is an entire function w.r.t. z, when n is a fixed positive real number. Observe that the forcing term F depends in particular on the choice of θ 1 .
The following property related to the lower incomplete Gamma function will be crucial in the construction of the auxiliary equation of the problem. Namely, (11) We recall that the infinite Laplace transform satisfies for every positive natural numbers n, k. This property will be used with respect to the second time variable, whereas a truncated Laplace transform depending on each value of the perturbation parameter near the origin is applied on the first variable. Both, (11) and (12) give rise to adequate algebraic properties which allow to reduce the main equation in the form of an auxiliary problem.
Regarding (10), F is holomorphic w.r.t. T 1 on C , T 2 on the disc D(0, T 0 /2) and on H β w.r.t. z. Furthermore, according to (11), (12), we observe that as tends to 0, for (well chosen) fixed T 1 . Therefore, F is getting closer to a polynomial in T 1 as tends to 0. The function f defined by is holomorphic on C ×D(0, h )×H β ×E, for any given open sector E centered at 0 and contained in D(0, 0 ) \ {0}, with h > 0 such that 0 < h 0 < T 0 /2. From the remark above, we check in particular that f becomes close to a polynomial in t 1 as becomes closer to the origin.

Auxiliary problems
We search for solutions of (8) in the form of an inverse Fourier transform The classical properties of inverse Fourier transform, together with (5), lead to an auxiliary functional equation satisfied by the expression U (T 1 , T 2 , m, ), namely Let 0 < κ < 1 and θ 1 , d 2 ∈ R. Let → h 1 ( ) be a holomorphic function defined on the domain of definition of the perturbation parameter, to be detailed afterwards. For every fixed value of the perturbation parameter , we search for solutions of (13) written as the Laplace transform with respect to T 2 along direction d 2 and the truncated Laplace transform with respect to T 1 along direction θ 1 − λk 2δD 2 arg( ) applied to a second auxiliary function. More precisely, we search for solutions of (13) of the form for some constants A m , 1 ≤ ≤ m − 1.
The assumption (6) guarantees the existence of d 2 k 2 ∈ N such that The following result states a one-to-one correspondence between the solution of (13) and (17). Its proof, which is omitted, can be adapted with minor modifications from [6], Lemma 1. (14). Then, it holds that Lemma 1, (15) together with the shape of the solution in (14), the assumptions on the coefficients c 1 2 and the forcing term f , and Lemma 2 entail ω being a solution of the following auxiliary convolution equation in the Borel plane: So far, the solution is of symbolic nature. The geometry of the problem, detailed in the following section, together with Section 4 provide convergence and growth estimates of such solution.
3 On the geometry of the problem In this section, we preserve the objects and assumptions detailed in Section 2.2 on the elements involved in the construction of the main problem under study (8), giving rise to the auxiliary problem (17). This section is devoted to the study of the geometry of the problem, which is crucial in the asymptotic approximation of the solution.
We define for every m ∈ R the polynomial In the case that τ 1 = 0 one can factorize P m in the form We assume that the polynomials Q and R D 1 D 2 satisfy that where S Q,R D 1 D 2 stands for an unbounded sector Let λ be a real number which satisfies that Let E be a sector with vertex at the origin which is contained in the disc D(0, 0 ) and for every ∈ E we define and the quantities with r 11 := λk 2δD 2 , and r 22 := λk 1 δ D 1 . The next result summarizes the main properties of the geometric construction above, which will be used to state the asymptotic behavior of the solutions of the main problem.
• Provided that λ > 0 is small enough, for any couple of directions (θ 1 , d 2 ) which satisfy that there exist an unbounded sector S d 2 with bisecting direction d 2 and small opening, and a sector S d 1 , , with for all ∈ E.
• Let S d 1 , and S d 2 be as above. We put Then, there exists C P > 0 which does not depend on ∈ E such that for every m ∈ R, τ ∈ Ω 1 ( ) × Ω 2 ( ).
The second statement is a direct consequence of the fact that for all τ 1 ∈ C and m ∈ R one has The pair (θ 1 , d 2 ) can be chosen accordingly, provided that λ, η Q,R D 1 D 2 ,δ 1 > 0 are small enough.

Banach spaces of functions with exponential growth
In this section, we recall the definition and main properties of certain Banach spaces previously used by the authors in [4], and adapted to the several variable case in [6,8]. The dependence of the domains of definition involved in the norm with respect to the values of the perturbation parameter has previously been consider in [5]. Let E be a sector of finite radius in the complex plane. For every ∈ E we consider the following two domains: a finite sector Ω 1 ( ) with vertex at the origin, bisecting direction d 1 which depends on , and radius r 1 ( ); and the union of an infinite sector S d 2 with vertex at the origin, fixed bisecting direction d 2 and positive opening which do not depend on together with the disc D(0, r 2 ( )) for some r 2 ( ) > 0, say Ω 2 ( ), i.e. Ω 2 ( ) = S d 2 ∪ D(0, r 2 ( )).
In the following we write d = (d 1 , d 2 ).
The first result follows directly from the definition of the norm of the Banach space in Definition 2.
Some parts of the proof of the following result can be adapted from that of Proposition 2 in [4]. We decided to include it completely for the sake of completeness and a self-contained work.
Proof Let f ∈ F d (ν,β,µ,k, ) . It holds that Therefore, one has , m)ds 2 (ν,β,µ,k, ) The classical estimates for m 1 ≥ 0 and m 2 > 0, yield At this point, we provide upper bounds for C 2.3 ( ). Let g(y) = ay 1+cy b , for some a, c > 0 and b > 1. The function g attains its maximum at y 0 = (c(b − 1)) −1/b . We apply this result to the case a = |τ 2 | k 2 , b = σ 1 /σ 1 and c = |τ 2 | k 2 σ 2 to arrive at , for some C(σ 1 ,σ 1 ) > 0. We plug the bound (32) into (31) and make the change of variable h = | | k 2 h at the integral in C 2.3 ( ) to arrive at for some C 24 > 0 only depending onσ 1 , σ 1 . Let x 0 > 0. The previous sup. particularized for those τ 2 ∈ Ω 2 ( ) such that x := |τ 2 / | k 2 > x 0 reads as For x > x 0 , one can perform the change of variable h = xu in the integral of ∆(x), in order to get that A partial fraction decomposition allows to write F k 2 (x) = F 1,k 2 (x) + F 2,k 2 (x), where We observe that F 1,k 2 (x) ≤ F 1,k 2 4 + x 2 , F 2,k 2 (x) ≤ F 2,k 2 4 + x 2 , for some positive constants F 1,k 2 , F 2,k 2 , valid for all x ≥ x 0 . Under the second assumption in (29), we obtain that sup x>x 0 ∆(x) is upper bounded by a constant. We conclude that the expression in (34) is upper bounded by for some C 25 > 0. It only rests to provide upper bounds for C 2.3 ( ) regarding the set of τ 2 ∈ Ω 2 ( ) such that 0 ≤ x ≤ x 0 . We observe that We plug this last expression into (31) to arrive at for some C 26 , C 27 , C 28 , C 29 > 0. We conclude that the left-hand side in (37) is upper bounded by In view of (35) and (38), we derive that which concludes the result. 2 The proof of the following result can be reproduced under minor adjustments from that of Proposition 2, [7].

Solution of an auxiliary problem
In this section, we preserve the elements and assumptions made on the main problem under study described in Section 2.2. More precisely, we assume the conditions (5), (6), (7) are satisfied by the parameters and elements involved. The coefficients c 1 2 (z, ) and the forcing term f (t, z, ) are constructed accordingly. We also assume the geometry of the problem is set in accordance with the assumptions made in Section 3 (see condition (18)), and preserve the values of r 1 ( ), r 2 ( ) for each ∈ E and λ > 0 (see (19) and (21)).
We provide a solution of the auxiliary problem (17) by means of a fixed point method in the Banach spaces introduced in Section 4.
Let ∈ D(0, 0 ) \ {0} and let d 1 , d 2 be chosen as described in Section 3. We define the operator We consider the Banach space of Definition 2, when fixing the domains described in (22), in accordance with the geometric analysis of the problem, in Section 3.
Corollary 1 Under the assumptions made in Proposition 2, the function ω d k (τ , m, ) is a solution of the auxiliary equation (17). Moreover, for every ∈ E, it satisfies that for every τ ∈ Ω 1 ( ) × Ω 2 ( ) and m ∈ R. The constant C ω > 0 can be chosen uniformly for all ∈ E.

Analytic solutions of the main problem
The main aim in this section is to provide analytic solutions of (8) for each of the elements of a family of sectors with respect to the perturbation parameter in the form of a truncated Laplace, Laplace and Fourier transforms. We first fix the geometric elements in this construction.

Definition 3
Let ι be an integer number, ι ≥ 2. Let E := (E p ) 0≤p≤ι−1 , where E p stands for a finite open sector with vertex at the origin, radius smaller than 0 . We assume the intersection of three different elements in E is empty, and 0≤p≤ι−1 E p = U \ {0}, for some neighborhood of the origin U ⊆ C. For the sake of simplicity, we arrange the sectors in order that nonempty intersections of sectors in E correspond to consecutive indices in the ring of integers modulo ι.
Under this configuration, we say that E describes a good covering in C .

Definition 4
Let ι be an integer number, ι ≥ 2, and let E := (E p ) 0≤p≤ι−1 be a good covering in C . Let T j be an open sector with vertex at the origin in C and finite radius r T j > 0, for j = 1, 2. For all 0 ≤ p ≤ ι − 1 we consider two bounded sectors S d j,p of bisecting direction d j,p , and small opening. In the following statements, we identify the indices p = ι and p = 0.
We say that the set is admissible if there exists δ > 0 such that for j = 1, 2 one has for every 0 ≤ p ≤ ι − 1, ∈ E p , t j ∈ T j and ξ j ∈ R (which may depend on t j and ) such that Let ι ≥ 2 be an integer number. Let E = (E p ) 0≤p≤ι−1 be a good covering and consider an admissible set {T 1 , T 2 , E, (S d 1,p ) 0≤p≤ι−1 , (S d 2,p ) 0≤p≤ι−1 }, which is associated to the good covering E. We briefly discuss the feasibility of such a construction. Indeed, let 0 ≤ p ≤ ι − 1 be fixed. We can first choose the direction ν 1,p (related to a fixed direction θ 1 depending on p) such that (52) holds for j = 1. Then, select the direction ν 2,p = d 2 in order that (52) holds for j = 2 together with the condition stated in the second item of Lemma 3.
Let 0 ≤ p ≤ ι − 1. For each 0 ≤ p ≤ ι − 1, we consider the main problem under study (8) under the assumptions (5)-(7), and departing from the coefficients c 1 2 (z, ) and the forcing term f (t, z, ) defined in Section 2.2. In virtue of the geometry of the problem described in Section 3 and Corollary 1, in particular the assumption of condition (18) and the choice of λ and r j ( ) for j = 1, 2 in (19) and (21) resp., one has that for every ∈ E p there exist a vector of directions d p = (d p,1 ( ), d p,2 ), a bounded sector with vertex at the origin S d p,1 , and bisecting direction d p,1 , with S d p,1 , ⊆ D(0, r 1 ( )) and an infinite sector S d p,2 of bisecting direction d p,2 such that the problem (17) admits a solution, say ω dp k (τ , m, ). Let us write Ω p,1 ( ) := S d p,1 , and Ω 2 ( ) := D(0, r 2 ( )) ∪ S d p,2 .
Indeed, the construction of ω dp k (τ , m, ) and the definition of u p (t, z, ) in (54) allow to affirm that the function The properties of Fourier transform (see Section 2.1) and Laplace transform (see Lemma 2), together with the definition of the elements involved in the main equation guarantee that (55) represents a solution of the main problem (8).
From now on, we refer to consecutive solutions of (8) to solutions associated to consecutive sectors in the corresponding good covering, which have nonempty intersection.
The next property on the difference of two consecutive solutions will be crucial in order to provide the asymptotic behavior of the solution at 0 regarding the perturbation parameter.
Moreover, there exist K, M > 0 such that for every 0 ≤ p ≤ ι − 1, one has (56) sup Proof The first part of the proof is guaranteed from the construction of the function u p (t, z, ) for every 0 ≤ p ≤ ι − 1. Let 0 ≤ p ≤ ι − 1. For every ∈ E p ∩ E p+1 we distinguish different situations depending on the relative position of the directions d p,1 , d p+1,1 , and d p,2 , d p+1,2 .
It is worth mentioning that the following cases state three equivalence classes regarding each element in the good covering. A continuity argument yields that for all 0 ≤ p ≤ ι − 1, if there exists ∈ E p ∩ E p+1 such that one of the following mutually excluded cases holds for such , then the same case holds for every element in E p ∩ E p+1 .
Case 1: Assume that L d p,1 , ≡ L d p+1,1 , and L d p,2 differs from L d p+1,2 . This situation occurs in case that the first component of every singularity in the Borel plane does not fall between the directions d p,1 and d p+1,1 but at least the second component of one singular point in the Borel plane occurs within angles between d p,2 and d p+1,2 .
In the previous expression, we have extended in a natural manner the notation adopted for the integration paths in Case 1 and Case 2. Analogous bounds as those stated for the integral I 25 (resp. I 23 ) are also valid for I 33 (resp. I 34 ), in Case 2. For the expression I 35 (resp. I 36 ) one can consider the estimates used to study I 14 (resp. I 13 ), involved in Case 1. We conclude the existence of C ω k ,6 , C 24 > 0 such that |u p+1 (t, z, ) − u p (t, z, )| ≤ C ω k ,6 exp − C 24 | | α , for all ∈ E p ∩ E p+1 , t ∈ ((T 1 ∩ D(0, h )) × (T 2 ∩ D(0, h ))) and z ∈ H β , with α defined in (57). Figure 1 illustrates the deformation of the paths involved in the procedure. 2

Parametric Gevrey asymptotic expansions of the analytic solutions
In this section, we analyse the asymptotic behavior of the analytic solutions of the main problem (8) obtained in the previous section, regarding the perturbation parameter approaching the origin. The classical criterion for k−summability of formal power series with coefficients in a Banach space, known as Ramis-Sibuya Theorem (see [1], p.121, or Lemma XI-2-6 in [2]) will be used to describe the Gevrey asymptotic approximation of the solution. The assumptions made in Section 2.2 and construction of the elements related to the main problem under study (8) are maintained in this section.
We first give some words on this classical summability theory for the sake of completeness.

k−summable formal power series and Ramis-Sibuya Theorem
Let (E, · E ) be a complex Banach space.
Definition 5 Let k ≥ 1 be an integer number. A formal power seriesf ( ) = n≥0 f n n ∈ E[[ ]] is k−summable with respect to along direction d ∈ R if there exists a bounded holomorphic function f defined in a finite sector V d of bisecting direction d and opening larger than π/k, and with values in E, which admitsf as its Gevrey asymptotic expansion of order 1/k on V d , i.e. for every proper subsector V 1 of V d , there exist K, M > 0 such that for every integer N ≥ 1 and ∈ V 1 . Watson's lemma guarantees uniqueness of such function, known as the k−sum of the formal power series.
Theorem 2 (RS) Let ι ≥ 2 and let (E p ) 0≤p≤ι−1 be a good covering in C . For every 0 ≤ p ≤ ι − 1 we consider a holomorphic function G p : E p → E, and define the function Θ p ( ) := G p+1 ( ) − G p ( ) holomorphic in Z p := E p ∩ E p+1 . We assume the following statements hold: • G p is a bounded function for ∈ Z p , → 0 for all 0 ≤ p ≤ ι − 1.
• Θ p is an exponentially flat function of order k in Z p for all 0 ≤ p ≤ ι − 1, i.e. there exist K, M > 0 such that Θ p ( ) E ≤ K exp − M | | k , valid for all ∈ Z p , and each 0 ≤ p ≤ ι−1.
Then, each of the functions G p ( ), for 0 ≤ p ≤ ι − 1 admits a common formal power serieŝ as Gevrey asymptotic expansion of order 1/k on E p . In addition to this, if the opening of E p 0 is larger than π/k for some 0 ≤ p 0 ≤ ι − 1, then G p 0 ( ) is unique, being the k−sum ofĜ( ) on E p 0 .

Asymptotic behavior of the solutions of (8) in the perturbation parameter
We are in conditions to describe the asymptotic behavior of the analytic solutions of the main problem under study (8) with respect to the perturbation parameter, at the origin. For this purpose, we consider a good covering E = (E p ) 0≤p≤ι−1 , for some integer number ι ≥ 2. We also fix an admissible set {T 1 , T 2 , E, (S d 1,p ) 0≤p≤ι−1 , (S d 2,p ) 0≤p≤ι−1 }, which is associated