Moulay A Barkatou
email: moulay.barkatou@unilim.fr

Suzy S Maddah
email: suzy.maddah@inria.fr

Formal Reduction of Singularly-Perturbed Linear Differential Systems

Keywords: Singularly-perturbed linear differential systems, turning points, cyclic vectors, rank reduction, singularities, exponential part, formal solutions, computer algebra

In this article, we discuss formal invariants of singularly-perturbed linear differential systems in neighborhood of turning points and give algorithms which allow their computation. The algorithms proposed are implemented in the computer algebra system Maple.

ε 2 ∂ 2 f 2 -(x 3 -ε) f = 0
over the (x, ε)-region D : |x| ≤ α 0 , 0 < |ε| ≤ ε 0 for some positive constants α 0 and ε 0 . Setting

F = f ε∂ f
, this equation can be rewritten as the following first-order system:

ε ∂F = A(x, ε) F = 0 1 x 3 -ε 0 F. (1)
One can observe that the Jordan form of A 0 (x) := A(x, 0) is not stable in any neighborhood of x = 0. We refer by turning points to such points where the Jordan form of A 0 (x) changes, i.e., either the multiplicity of the eigenvalues or the degrees of the elementary divisors change in a neighborhood of such points [40, p. 57]. For example, x = 0 might be a turning point for [START_REF] Abbas | Formal Reduction of a class of Pfaffian Systems in Two Variables[END_REF].

Our goal is to compute a formal solution in a neighborhood of x = 0 as ε → 0. • In the region β 0 |ε| 1/3 ≤ |x| ≤ α 0 where β 0 is a positive constant, the transformation

F = T 1 G
where T 1 = 1 0 0 x 3/2 transforms system (1) into (x -3 ε)x 3/2 ∂G = { 0 1 1 0

+ (x -3 ε) 0 0 -1 -3 2 x 1/2 }G,
which can be rewritten, by a self-explanatory change of notation x -3 ε = ξ, as

ξx 3/2 ∂G = { 0 1 1 0 + ξ 0 0 -1 -3 2 x 1/2 }G. A transformation G = T 2 U where T 2 = 1 2 1 2 -1 2 1 2
+ O(ξ) would then result in the blockdiagonalized system

ξx 3/2 ∂U = { -1 0 0 1 + O(ξ)}U.
Consequently, the system splits into two first-order linear differential scalar equations. We can then write down the following fundamental matrix of formal solutions for [START_REF] Abbas | Formal Reduction of a class of Pfaffian Systems in Two Variables[END_REF] in this region:

F outer = T 1 T 2 exp( -2 5 x 5/2 ε -1 + O(x -1/2 ) 0 0 2 5 x 5/2 ε -1 + O(x -1/2 )
).

• In the region |x| ≤ β 0 |ε| 1/3 , we shall perform a so-called stretching transformation. Namely we set τ = xε -1/3 and ∂ τ = d dτ . Then, for all τ such as |τ| < ∞, except possibly for neighborhoods of the roots of τ 3 -1 = 02 , the transformation F = L 1 G where

L 1 = 1 0 0 ε 1/2
reduces system (1) to

ε 1/6 ∂ τ G = 0 1 τ 3 -1 0 G. (2)
And the transformation G = L 2 U where

L 2 =        1 2 + 1 8 τ 3 + O(τ 6 ) 1 2 + 1 8 τ 3 + O(τ 6 ) -i 2 + i 8 ω 3 + O(τ 6 ) i 2 -i 8 τ 3 + O(τ 6 )        ,
reduces system (2) to

ε 1/6 ∂ τ U = { -i + i 2 τ 3 + O(τ 6 ) 0 0 i -i 2 τ 3 + O(τ 6 )
+ O(ε 1/6 )}U.

A fundamental matrix of formal solution of system ( 1) is then given by:

F inner = L 1 L 2 exp(        -iτ+(i/8)τ 4 +O(τ 7 ) ε 1/6 0 0 iτ-(i/8)τ 4 +O(τ 7 ) ε 1/6        ).
We call F outer (resp. F inner ) an outer (resp. inner) solution as it is obtained in the outer (resp. inner) region around x = 0. The corresponding differential systems are sometimes referred to as outer and inner differential systems as well. However, unlike the above particular example, intermediate regions might be encountered.

In this article, we study the construction of formal solutions (the formal reduction), of systems of the following form:

ε h ∂F = A(x, ε)F, (3)
where the entries of the matrix A(x, ε) are formal power series in ε, whose coefficients are formal power series in x; and h is an integer which we call the ε-rank of the system. In [? , p. 74], Iwano summarized as follows the problems needed to be resolved to obtain the complete knowledge about the asymptotic behavior of the solutions of a singularly-perturbed linear differential system: (P1) Divide a domain [D] in (x, ε)-space into a finite number of sub-domains so that the solution behaves quite differently as ε tends to zero in each of these sub-domains;

(P2) Find out a complete set of asymptotic expressions of independent solutions in each of these sub-domains;

(P3) Determine the so-called connection formula; i.e. a relation connecting two different complete sets of the asymptotic expressions obtained in [START_REF] Abbas | On the Reduction of Singularly-Perturbed Linear Differential Systems[END_REF].

In order to resolve (P1) and (P2), Iwano proceeded by associating a convex polygon to firstorder systems in analogy with the scalar case, after imposing a precise triangular structure on A(x, ε). (P3) remains generally unresolved and it is out of the scope of this article as well (see, e.g. [START_REF] Fruchard | Composite Asymptotic Expansions[END_REF][START_REF] Hulek | Systémes d' équations différentielles linéaires singuliérement perturbées et développements asymptotiques combinés[END_REF][START_REF] Wasow | Linear Turning Point Theory[END_REF] and references therein). The hypotheses on A(x, ε) were eventually relaxed for special types of systems in a series of papers by Iwano, Sibuya, Wasow, Nakano, Nishimoto, and others (see references in [START_REF] Wasow | Topics in the Theory of Linear Ordinary Differential Equations Having Singularities with respect to a Parameter[END_REF]). A famous and prevalent method among scientists is the Wkb method (see, e.g. [START_REF] Langer | Turning points in linear asymptotic theory[END_REF]), which, roughly speaking, reduces the system at hand into one whose asymptotic behavior is known in the literature. However, as pointed out in [START_REF] Wasow | Simplification of Turning Point Problems For Systems of Linear Differential Equations[END_REF], in view of the great variety and complexity of systems given by (3), it is rarely expected that a system can be reduced into an already investigated simpler form. The methods proposed in the literature for the symbolic resolution of such systems are not yet fully algorithmic, exclude turning points (see [START_REF] Chen | Solutions Formelles de Systemes d'Equations Differentielles Lineaires Ordinaires Homogenes[END_REF] and references therein), treat systems of dimension two only, rely on Arnold-Wasow form [START_REF] Arnold | On Matrices Depending on Parameters[END_REF], and/or impose further restrictions on the structure of the input matrix. Moreover, there exists no computer algebra package dedicated to the symbolic resolution of neither system (3) nor the scalar n th -order scalar equation. The widely-used softwares Maple and Mathematica content themselves to the computation of outer formal solutions for scalar equations, and so does [START_REF] Macutan | Formal Solutions of Scalar Singularly-Perturbed Linear Differential Equations[END_REF].

In this article, we attempt to resolve (P2) algorithmically, that is give an algorithm to construct a fundamental matrix of formal solutions of an input system. We also attempt to give an insight into (P1) within some examples.

As illustrated by the introductory example, the formal reduction of a singularly-perturbed system given by (3) leads inevitably to the consideration of more general systems. These considerations will be discussed in Section 1. In Section 2, we give necessary preliminaries including results on n th -order scalar equations, the base cases, and the well-known Splitting lemma which splits the system into subsystems of lower dimensions, whenever possible. We discuss our contribution in the remaining sections, as we deviate from the classical reduction to give algorithms which act on the system directly without resorting to an equivalent scalar equation or to the Arnold form: First we refine in Sections 3.1 and 3.2 the algorithms which we proposed in [START_REF] Abbas | On the Reduction of Singularly-Perturbed Linear Differential Systems[END_REF] for the resolution of turning points and reduction of the ε-rank h to its minimal integer value. Then in Section 4, we associate an ε-Newton polygon to a given system and show its role in retrieving formal invariants. And finally, we make use in Section 5 of these notions introduced to propose an overall formal reduction algorithm.

As we have already indicated, the literature on this problem is vast. In the hope of keeping this article concise but self-contained, we restrict our presentation to constructive results which contribute directly to our process of formal reduction. The omitted proofs are collected in the appendix along with illustrating examples.

In this paper, we content ourselves to the process of formal reduction. Any reference to the asymptotic interpretation of formal solutions will be dropped in the sequel. One may consult in this direction the method of matched asymptotic expansions (see, e.g. [START_REF] Skinner | Singular Perturbation Theory[END_REF][START_REF] Wasow | Linear Turning Point Theory[END_REF] and references therein) or composite asymptotic expansions (see, e.g. [START_REF] Fruchard | Composite Asymptotic Expansions[END_REF][START_REF] Hulek | Systémes d' équations différentielles linéaires singuliérement perturbées et développements asymptotiques combinés[END_REF] and references therein).

In the sequel, we adopt C as the base field for the simplicity of the presentation. However, any other commutative field F of characteristic zero (Q ⊆ F ⊆ C) can be considered instead. The algorithms presented herein require naturally algebraic extensions of the base field. However, such extensions can be restricted as described in [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF]Section 7.2].

Since we are leading a local investigation, we assume that the input system has at most one singular point, otherwise the region of study can be shrunk. Moreover, we place this singular point at the origin, otherwise a translation in the independent variable can be performed (x → 1/x for a singularity at ∞).

Notations

• C is the field of complex numbers, Q is the field of rational numbers, and

Q -= {q ∈ Q | q ≤ 0}.
• We denote the dimension of the systems under consideration by n. We use r for the algebraic rank of specified matrices (leading matrix coefficients) and h for the ε-rank.

• We use upper case letters for algebraic structures, matrices, and vectors; the lower case letters i, j, k, l for indices; and some lower case letters (greek and latin) locally in proofs and sections, such as u, v, σ, ρ, , etc...

• We use x, t, τ as independent variables and ∂ to signify a derivation. Since x dominates this article, we drop it in the derivation. We set ∂ = d dx and ∂ τ = d dτ .

• We use F for the unknown n-dimensional column vector (or n × n-matrix); we use G, U, W, and Z in this context as well. For scalar equations, we use f and g.

• I d×n (resp. I n ) stands for the identity matrix of dimensions d × n (resp. n × n); and O d×n (resp. O n ) stands for the zero matrix of dimensions d × n (resp. n × n). If confusion is unlikely to arise, we simply use 0 to denote O d×n or O n .

• We say that M ∈ M n (R) whenever the matrix M is a square matrix of size n whose entries lie in a ring R.

• GL n (R) is the general linear group of degree n over R (the set of n × n invertible matrices together with the operation of matrix multiplication).

• C[[x]

] is the ring of formal power series in x whose coefficients lie in C; and C((x)) is its fraction field. The usual valuation of an element g(x) of C((x)), i.e. the order of g in x, is denoted by val x (g), with val x (0) = +∞. The ring C[[x]] (resp. the field C((x))) is endowed with a differential ring (resp. field) structure by considering the derivation ∂.

• C[[ε]

] is the ring of formal power series in ε whose coefficients lie in C; and

C((ε)) is its fraction field. C[[x]][[ε]] (resp. C[[x]]((ε))) is the ring of formal (resp. meromorphic) power series in ε whose coefficients lie in C[[x]]. Similarly, C((x))[[ε]] (resp. C((x))((ε)))
is the ring (resp. field) of formal (resp. meromorphic) power series in ε whose coefficients lie in C((x)).

• Any element f ∈ C((x))((ε)) can be written as f = ∞ k∈Z f k (x)ε k , where f k (x) ∈ C((x)) for all k ∈ Z. The derivation ∂ extends naturally to C((x))((ε)) by the formula ∂ f = ∞ k∈Z (∂ f k (x))ε k . Moreover, the map val ε : C((x))((ε)) → Z ∪ ∞ defines a valuation over C((x))((ε)), satisfying the following properties for all f (x, ε), g(x, ε) in C((x))((ε)): -val ε ( f ) = ∞ if and only if f = 0; -val ε ( f g) = val ε ( f ) + val ε (g); -val ε ( f + g) ≥ min (val ε ( f ), val ε (g)).
• We give the blocks of a matrix M with upper indices, e.g.

M = M 11 M 12 M 21 M 22 ,
and the size of the different blocks is dropped unless it is not clear from the context.

Ring of coefficients

In this first part, we investigate a suitable ring of coefficients to treat first-order singularlyperturbed linear differential systems and n th -order scalar singularly-perturbed differential equations, in a neighborhood of a turning point (see also [START_REF] Wasow | Topics in the Theory of Linear Ordinary Differential Equations Having Singularities with respect to a Parameter[END_REF]Chapter 2] and references therein). For (U 0 , V 0 ) ∈ R 2 , we put

P(U 0 , V 0 ) = {(U, V) ∈ R 2 | U ≥ U 0 and V ≥ V 0 }. Let f = ∞ k∈Z f k (x)ε k ∈ C((x))((ε))
and let P f be the union of P(k, val x ( f k )), k ∈ Z. Then the Newton polygon of f , denoted by N f , is the boundary of the convex hull in R 2 of the set P f . We now consider, for σ, p ∈ Q with σ ≤ 0, the half-plane

H σ,p = {(U, V) ∈ R 2 | V ≥ σU + p}.

One can verify that

K = { f ∈ C((x))((ε)) | N f ⊂ H σ,p for some σ ∈ Q -, p ∈ Q }
is a differential field endowed with the derivation ∂ = d/dx and it is the field of fractions of the ring R = K ∩ C((x)) [[ε]]. Moreover, the elements of K can be characterized geometrically as follows:

f = k∈Z f k (x)ε k ∈ K ⇐⇒ ∃ σ ∈ Q -, p ∈ Q | ∀ k, val x (f k ) ≥ σk + p.
In fact, given a non-zero element f of K, we set ν f := val ε ( f ). Among the half-planes which contain N f , we consider the half-planes which are bounded below by a straight line passing through the point (ν f , val x ( f ν f )), i.e. the half-planes H σ,p f with p f := val x ( f ν f ) -σν f . Then, among the H σ,p f 's, we pick the half-plane of maximal slope, which we denote by σ f . Evidently, H σ f ,p f is determined by a straight line which passes through the point (ν f , val x ( f ν f )) and point(s) (k, val x ( f k )) for at least one k > ν. We denote by (S f ) this straight line whose equation is given by V = σ f U + p f , and we say that it is associated to f . For f = 0, we set σ f = 0 and p f = 0. Example 1. Given a nonzero element f of C((x))((ε)). If we can construct a straight line of finite slope which passes through (ν f , val x ( f ν f )), and stays below the points (k, val x ( f k )) for all k > ν, then f ∈ K. We illustrate the following examples in Figure 1:

• Let a = x 3 ε + x -8 ε 2 = x 3 ε(1 + x -11 ε) then σ a = -11. • Let g = ∞ k=0 x -3k ε k then σ g = -3. • Let f = x 2 + ∞ k=1 x -3k ε k then σ f = -5. • Let s = x -2 + ∞ k=1 x -3k ε k then σ s = -3.
The following two lemmas give further insight into K:

Lemma 1. If f (x, ε) ∈ K is non-zero, then there exists e ∈ Z -such that f (x, x e ε) ∈ C[[x]]((ε)) . Proof. Let f = ∞ k=ν f f k (x)ε k . We set µ k = val x ( f k ) and f k = x µ k fk .
Then for any e ∈ Z -with e ≤ σ we have:

f (x, ε) = ε ν f ( f ν f + f ν f +1 ε + f ν f +2 ε 2 + f ν f +3 ε 3 + . . . ) = ε ν f x µ ν f ( fν f + x µ ν f +1 -µ ν f fν f +1 ε + x µ ν f +2 -µ ν f fν f +2 ε 2 + x µ ν f +3 -µ ν f fν f +3 ε 3 + . . . ) = ε ν f x µ ν f ( fν f + x µ ν f +1 -µ ν f -e fν f +1 (x e ε) + x µ ν f +2 -µ ν f -2e fν f +2 (x e ε) 2 + x µ ν f +3 -µ ν f -3e fν f +3 (x e ε) 3 + . . . ). But, µ k ≥ σ f (k -ν f ) + µ ν f for all k ≥ ν f . Hence, µ k ≥ e(k -ν f ) + µ ν f for all k ≥ ν f .
We can thus choose e to be the largest integer which is less than or equal to σ f . Lemma 2. The field of constants of the differential field

(K, ∂) is C((ε)). Proof. Let f ∈ K and suppose that ∂ f = 0. Then, ∞ k=ν f (∂ f k (x))ε k = 0 which yields ∂ f k (x) = 0 for all k ≥ ν f . Thus, for every k ≥ ν f , f k (x) = c k for some c k ∈ C, which yields f = ∞ k=ν f c k ε k .
Henceforth, algorithmically, we can use the following handy representation of any element f ∈ K:

f (x, ε) = ∞ k=ν f f k (x)ε k = ∞ k=ν f x -kσ f f k (x)(x σ f ε) k = x p f ∞ k=ν f x -kσ f -p f f k (x)ξ k f , with ξ f = x σ f ε.
(4) Evidently, val x (x -σ f k-p f f k (x)) ≥ 0 for all k ≥ ν f . Moreover, under the notation ξ f , we have:

∂( k=ν f f k ξ k f ) = k=ν f (∂( f k ) + σ f k x f k )ξ k f .
1.1. Input systems and input equations Motivated by the introductory example, and with the help of the above notations, we will investigate the following system rather than (3):

∂F = A(x, ε)F where A ∈ M n (K).
(

) 5
The valuation val ε (and the usual valuation val x of C((x))) extends naturally to matrices A = (a i, j ) ∈ M n (K) by val ε (A) = min i, j val ε (a i, j ). We can then define the Newton polygon of A ∈ M n (K), ν A ∈ Z, p A ∈ Q, and σ A ∈ Q -, in an analogous manner to those of f ∈ K. Hence, we can express A(x, ε) as follows:

A = ∞ k=ν A A k ε k A , A ν A 0, p A = val x (A ν A ) -σ A ν A , and val x (A k ) ≥ σ A k + p A ∀ k ≥ ν A . (6)
Hence by setting ξ A = x σ A ε, we can write in analogy with (4):

A(x, ε) = x p A ∞ k=ν A x -kσ A -p A A k (x) ξ k A = ξ ν A A x p A ∞ k=0 x -kσ A -val x (A ν A ) A ν A +k (x) ξ k A where x -val x (A ν A ) A ν A | x=0 0.
Thus, we can rewrite system (5) as follows:

ξ -ν A A x -p A ∂F = ( ∞ k=0 x -kσ A -val x (A ν A ) A ν A +k (x) ξ k A ) F, where x -kσ A -val x (A ν A ) A ν A +k ∈ M n (C[[x]]).
Similarly, if we consider a n th -order differential equation whose coefficients are elements of K:

a n (x, ε)∂ n f + a n-1 (x, ε) ∂ n-1 f + . . . + a 1 (x, ε) ∂ f + a 0 (x, ε) f = 0, ( 7
)
then setting F = [ f, ∂ f, ∂ 2 f, . . . , ∂ n-1 f
] T , we can express (7) equivalently as a system ∂F = A(x, ε)F where A(x, ε) ∈ M n (K) is a companion matrix, and we define ν a := ν A , p a := p A , σ a := σ A , and ξ a := ξ A .

1.1.1. Notations Thus, in this paper, we treat system (5) in the following form:

[A σ A ] ξ h A x p A ∂F = A(x, ξ A )F, A(x, ξ A ) = ∞ k=0 A k (x)ξ k A ∈ M n (R), A 0 (x = 0) 0, ( 8
) where ξ A = x σ A ε, σ A ∈ Q -, p A ∈ Q, and h ∈ Z. Under the definition of ξ A , A(x, ξ A ) = k=0 A k (x)ξ k A with A k (x) ∈ M n (C[[x]]
) for all k ≥ 0. We refer to A 0 (x) as the leading coefficient matrix and to A 0,0 := A 0 (x = 0) as the leading constant matrix.

We also consider n th -order scalar linear differential equations of the form:

[a σ a ] ∂ n f + a n-1 (x, ξ a )∂ n-1 f + . . . + a 1 (x, ξ a )∂ f + a 0 (x, ξ a ) f = 0, (9)
where

a i (x, ξ a ) = ∞ k∈Z a i,k (x) ε k ∈ K for all i ∈ {0, . . . , n}, with a i,k (x) ∈ C[[x]
] for all k ≥ 0. We remark that a system [A σ A ] given by ( 8) can be also expressed as a scalar equation of the form [START_REF] Barkatou | A Rational Version of Moser's Algorithm[END_REF]. The theoretical possibility stems from the work of [START_REF] Churchill | Cyclic vectors[END_REF] which discusses cyclic vectors. Moreover, an algorithm to compute a companion block diagonal form of [A σ A ] (and hence equivalent scalar form) can be obtained by generalizing the work in [START_REF] Barkatou | An Algorithm for Computing a Companion Block Diagonal Form for a system of Linear Differential Equations[END_REF] developed for unperturbed systems (see Appendix A).

In the sequel, for the clarity of the presentation, the index A (resp. a) will be dropped from ν A , p A , σ A , and ξ A (resp. ν a , p a , σ a , and ξ a ) whenever confusion is unlikely to arise.

Equivalent systems

Consider system [A σ A ] given by ( 8). Let T ∈ GL n (K) then the transformation F = TG (also called gauge transformation) yields a system

[ Ãσ à ] ξ à h x p∂G = Ã(x, ξ Ã)G, where à ∈ M n (R). ( 10
)
for some σ Ã ∈ Q -, p ∈ Q, and h ∈ Z. We say that systems [A σ A ] and T [A σ A ] := [ Ãσ Ã ] are equivalent. Examples of such transformations and their applications are the transformations T 1 , T 2 , L 1 , L 2 applied within the introductory example.

In the sequel, we seek at each step a transformation which yields for an input system [A σ A ], an equivalent system [ Ãσ Ã ], so that either h < h or [ Ãσ Ã ] can be decoupled into systems of lower dimensions. Using a recursive process which employs this approach, we aim to break down any input system [A σ A ] into system(s) which can be treated with known methods (see basic cases of Section 2.1).

Evidently, σ Ã might differ from σ A and hence ξ A will be updated to ξ Ã after the application of T (see the introductory example), which brings us to the next subsection.

Restraining index and (P1)

As we are lead naturally to the study of system [A σ A ] with σ A ∈ Q -within the study of system (8), we will need to understand the inevitable growth of the order of the poles in x within the reduction. For instance, in the introductory example, we start with an input system for which σ A = 0 and arrive at an decoupled system with σ U = -3.

As described in [40, Chapter 2] and as we have observed, these poles grow at worst linearly with k. From here initially stems our main motivation for our choice of the ring of coefficients and the representation in terms of ξ A : the orders of the poles in x which might be introduced in an input system within the process of reduction, are stored in p and σ, which allows investigating them at any step within the reduction process. We can thus talk about a restraining index ρ = -1/σ. This is the first step towards determining the inner and intermediate regions. In fact, the complicated behavior anticipated in the neighborhood of turning points can be investigated with the help of a sequence of positive rational numbers:

[ρ] 0 = ρ 0 < ρ 1 < ρ 2 < • • • < ρ m . (11)

Preliminaries

A detailed discussion of the scalar case (n = 1) and regularly-perturbed systems (h ≤ 0), is given in [18, Chapter 3, Section 1, p. 52 -56] for σ A = p A = 0. In this section, we give a brief discussion which discards this restriction on σ A and p A .

The base cases 2.1.1. The case h ≤ 0

If h ≤ 0 then the solution of system [A σ A ] can be sought, up to any order µ, upon presenting the solution as a power series in ξ A = x σ A ε, i.e.

F = ∞ k=0 F k (x) ξ -h+k A .
The latter can then be inserted into [A σ A ] ξ h x p ∂F = AF and the like powers of ξ A equated. This reduces the symbolic resolution to solving successively a set of inhomogeneous linear singular differential systems in x solely:

x p ∂F 0 = (A 0 (x) + h σ A x p-1 I) F 0 x p ∂F 1 = (A 0 (x) + (h -1) σ A x p-1 I) F 1 + A 1 (x) F 0 . . . x p ∂F µ = (A 0 (x) + (h -µ) σ A x p-1 I) F µ + A 1 (x) F µ-1 + . . . + A µ (x) F 0 .
For any k ≥ 0, the solution of the homogeneous system singular in x can be sought via the Maple package Isolde [START_REF] Barkatou | Integration of Systems of Ordinary Linear Differential Equations[END_REF] or the Mathemagix package Lindalg [START_REF] Maddah | Lindalg : Mathemagix Package for Symbolic Resolution of Linear Systems of Differential Equations[END_REF] (which are both based on [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF]). Afterwards, the solution of the inhomogeneous system can be obtained by the method of variation of constants (see, e.g. [START_REF] Balser | Formal Power Series and Linear Systems of Meromorphic Ordinary Differential Equations[END_REF]Theorem 3,p. 11]). We remark that a transformation F k = x -kσ A G reduces the inhomogeneous system in the dependent vector F k to:

x p ∂G = (A 0 (x) + h σ A x p-1 I n ) G(x) ,
which coincides with that of F 0 (x). Thus we have:

F k (x) = x -kσ A F 0 (x) F -1 0 (x) x kσ A [A 1 (x) F k-1 (x) + . . . + A k (x) F 0 (x)] dx.

The scalar case

The n th -order singularly-perturbed scalar differential equations are treated in [START_REF] Macutan | Formal Solutions of Scalar Singularly-Perturbed Linear Differential Equations[END_REF] upon introducing an analog to the Newton polygon and polynomials (see e.g. [6] and references therein for the unperturbed counterparts of these equations): the ε-polygon and ε-polynomials. In this subsection, we adapt this treatment to the more general equation [a σ a ] given by ( 9):

[a σ a ] ∂ n f + a n-1 (x, ξ a )∂ n-1 f + . . . + a 1 (x, ξ a )∂ f + a 0 (x, ξ a ) f = 0.
We first give an analog to [27, Definition 2.1] and then generalize [27, Proposition 4.1], whose proof is easily adaptable from [START_REF] Macutan | Formal Solutions of Scalar Singularly-Perturbed Linear Differential Equations[END_REF] (see Appendix B). Definition 1. Consider the scalar equation [a σ a ] given by [START_REF] Barkatou | A Rational Version of Moser's Algorithm[END_REF]. For (u 0 , v 0 ) ∈ R 2 , we put

P(U 0 , V 0 ) = {(U, V) ∈ R 2 | U ≤ U 0 and V ≥ V 0 }.
Set ν i = val ε (a i ) for i ∈ {0, . . . , n} and let P a be the union of P(i, ν i ). Then the ε-polygon of [a σ a ], denoted by N ε (a), is the intersection of R + ×R with the convex hull in R 2 of the set P a . We denote the slopes of the edges of N ε (a) by {e 1 , . . . , e }. These slopes are non-negative rational numbers. And, for every j ∈ {1, . . . , }, consider the algebraic equation given by

(E j ) k=0 a i k , val ε (a i k ) (x) X (i k -i 0 ) = 0,
where 0 ≤ i 0 < i 1 < • • • < i = n denote the integers i for which (i, ν i ) lie on the edge of slope e j of the ε-polygon, and a i,ν i (x) = ξ -ν i a a i (x, ξ a )| ξ a =0 . We say that (E j ) is the ε-polynomial associated to the slope e j .

Proposition 1. Consider a nonzero [a σ a ] given by (9) and its ε-polygon N ε (a) of slopes {e 1 , . . . , e }. Let f (x, ξ a ) = exp( q(x, ξ a )dx) with q(x, ξ a ) 0

∈ s∈N * C((x))((ξ 1/s a )). If f (x, ξ a ) satisfies [a σ a ] formally then q(x, ξ a ) = 1 ξ e j a (X(x) + O(ξ a ))
for some j ∈ {1, . . . , }, and X(x) is one of the non zero roots of the ε-polynomial (E j ) associated to e j .

The full expansion of q(x, ξ a ) can be obtained with successive substitutions of the form f = exp( X(x) Example 2. Consider the following linear differential equation with σ a = 0:

[a σ a ] ∂ 3 f - x ξ 2 a ∂ 2 f - 1 ξ 5 a f = 0. ( 12
)
We wish to study this equation according to its ε-polygon (see Figure 3). We have two slopes:

• The slope e 1 = 3 2 for which (E 1 ) x X 2 + 1 = 0 whose nonzero solutions are X = ± i √ x .

• The slope e 2 = 2 for which (E 2 ) X 3 -xX 2 = 0 whose nonzero solution is X = x.

Thus, the leading terms (of the exponential part) of solutions of (12) are given by exp ( x

ξ 2 a dx) = exp ( x ε 2 dx) and exp (± i √ xξ 3/2 a dx) = exp (± i √ xε 3/2 dx).
One can also refer to [START_REF] Macutan | Formal Solutions of Scalar Singularly-Perturbed Linear Differential Equations[END_REF]Section 5] for examples on the harmonic oscillator with damping and Orr-Sommerfeld equation.

Based on the above, we give the following definition:

Definition 2. Consider the scalar equation [a σ a ] given by [START_REF] Barkatou | A Rational Version of Moser's Algorithm[END_REF]. We call the largest slope of the ε-polygon N ε (a), the ε-formal exponential order (ε-exponential order, in short) :

ω ε (a) = max {e j , 1 ≤ j ≤ }.
Clearly, under the notations of this subsection, we have:

ω ε (a) = n-1 max i=0 (0, - val ε (a i ) n -i ).

Consequences of the equivalence between a scalar equation and a system

Due to their equivalence, the formal solutions of a first-order system [A σ A ] given by (8) can be computed from an equivalent scalar equation [a σ a ] given by (9), using the ε-polygon of [a σ a ].

Moreover, [START_REF] Barkatou | On k-simple forms of first-order linear differential systems and their computation[END_REF] can be computed as well using the Iwano-Sibuya's polygon (which collects valuations with respect to both x and ε [START_REF] Iwano | Reduction of the Order of a Linear Ordinary Differential Equation Containing a Small Parameter[END_REF]). However such a treatment is unsatisfactory (see e.g. [START_REF] Macutan | Formal Solutions of Scalar Singularly-Perturbed Linear Differential Equations[END_REF]Conclusion]) since it overlooks the information that we can derive from the system directly and demands an indirect treatment. Nevertheless, it plays a key role in the theoretical basis of the algorithm which we will develop in later parts. Most importantly, the equivalence between [A σ A ] and [a σ a ], and the invariance of the ε-exponential order under gauge transformations, allows the definition of the ε-exponential order of system [A σ A ] as follows: Definition 3. Consider a system [A σ A ] given by (8) and an equivalent scalar equation [a σ a ] given by [START_REF] Barkatou | A Rational Version of Moser's Algorithm[END_REF]. The ε-polynomials and ε-formal

exponential order ω ε (A) of [A σ A ] are those of [a σ a ].
It follows that, given a system [A σ A ], we have:

ω ε (A) = max 0≤i≤n (0, -val ε (a i ) n -i ). ( 13
)
where [a σ a ] is a scalar equation equivalent to

[A σ A ].
In this paper, we give a direct treatment of the system which will lead eventually to the computation of ω ε (A), the ε-polynomials, and consequently a basis of the space of formal solutions, without resorting to the equivalence of an input system to a scalar equation. Our treatment relies on the properties of the eigenvalues of the leading coefficient of the input system. We distinguish between three cases (distinct, unique, or zero eigenvalues), the first of which is classical and is recalled in the next subsection.

ε-Block Diagonalization

Consider a system [A σ A ] given by [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF]. A classical tool in the perturbation theory is the so-called splitting which separates off the existing distinct coalescence patterns. Whenever the leading constant matrix A 0,0 := A 0 (x = 0) admits at least two distinct eigenvalues, the system can be decoupled into subsystems of lower dimensions:

Theorem 1. Consider system [A σ A ] given by (8) [A σ A ] ξ h A x p ∂F = A(x, ξ A )F, A(x, ξ A ) ∈ M n (R), A 0,0 0. If A 0,0 = A 11 0,0 O O A 22 0,0
such that A 11 0,0 and A 22 0,0 have no eigenvalues in common, then there exists a unique transformation T

(x, ξ T ) = ∞ k=0 T k (x)ξ k T ∈ GL n (R), given by T (x, ξ T ) = I + ∞ k=0 O T 12 k (x) T 21 k (x) O ξ T k ,
such that the transformation F = TG gives

[ Ãσ à ] ξ h à x p ∂G = Ã(x, ξ Ã) G = Ã11 (x, ξ Ã) O O Ã22 (x, ξ Ã) G where Ã0,0 = A 0,0 . Moreover σ T , σ à ≥ σ where σ =        σ A + p-1 h if p < 1 σ A otherwise .
Well-known proofs proceed in two steps: Block-diagonalizing A 0 (x), followed by successive block-diagonalization of coefficients (see, e.g., [21, Theorem XII-4-1, p. 381] and [START_REF] Chen | Solutions Formelles de Systemes d'Equations Differentielles Lineaires Ordinaires Homogenes[END_REF]Chapter 3, for the particular case of σ A = p A = 0, and [START_REF] Wasow | Topics in the Theory of Linear Ordinary Differential Equations Having Singularities with respect to a Parameter[END_REF]Theorem 8.1, for the general case). We give a constructive proof in Appendix B. We remark that the specified form of A 0,0 in Theorem 1 is non-restrictive. It suffices that A 0,0 has at least two distinct eigenvalues. Such a block-diagonal form can be then attained by Jordan constant transformation. In case the base field is not algebraically closed, the weaker transformation given in [15, Lemma A.1] can be applied instead, since it does not require any algebraic field extensions.

After splitting the input system, we proceed in our reduction for each of the decoupled subsystems in parallel. We can thus assume, without any loss of generality, that the leading constant matrix of the input system [A σ A ] has a unique eigenvalue γ ∈ C. Thus, upon applying the socalled eigenvalue shifting, i.e.

F = EG = G exp( ξ -h A x -p γdx), (14)
one can verify that the resulting system E[A σ A ] has a nilpotent leading constant matrix. Thus, it remains to discuss this case of nilpotency. This is the goal of the two following sections.

Two-fold rank reduction

Without loss of generality, we can now assume that system [A σ A ] given by ( 8) is such that A 0,0 is nilpotent. At this stage of reduction, our approach diverges from the classical indirect ones which require Arnold-Wasow form, a cyclic vector, or reduction to a companion form. Two cases arise:

• A 0,0 is nilpotent but A 0 (x) is not: In this case, the eigenvalues of A 0 (x) might coalesce in some neighborhood of x = 0, which would cause the so-called turning point. In Subsection 3.1, we propose an algorithm to compute a transformation which reduces the input system to a system for which both A 0,0 and A 0 (x) are nilpotent.

• Both A 0,0 and A 0 (x) are nilpotent: In this case, we propose in Subsection 3.2 an algorithm to compute a transformation which can reduce the ε-rank h of system [A σ A ] to its minimal nonzero integer value. This minimality is an intrinsic property of the system and its solutions. Its computation paves the way for the retrieval of the ε-exponential order of the system, as we show later in Section 4.

Resolution of turning points

We illustrate the source of turning points in the following example:

Example 3. [38, p. 223] Let a(x), b(x) be holomorphic for x in a region Ω ∈ C and consider M(x) =           x 1 a(x) 0 x b(x) 0 0 0           . The transition matrix T (x) =           α(x) β(x) γ(x)(-a(x)x -1 + b(x)x -2 ) 0 α(x) -γ(x)b(x)x -1 0 0 γ(x)          
, with arbitrary scalar functions α(x), β(x), γ(x) such that α(x) 0 and γ(x) 0 for all x ∈ Ω, yields for x 0:

J(x) = T -1 M(x)T (x) =           x 1 0 0 x 0 0 0 0           . If α(x), β(x), γ(x) are holomorphic, so is T , provided that zero is not a point of Ω. If 0 ∈ Ω, then T (x)
has a pole at x = 0 unless the following conditions are satisfied: b(0) = 0 and ∂b(0) = a(0).

Therefore, M(x) is in general not holomorphically similar to J(x) in regions that contain x = 0. Moreover, if the second condition is not satisfied, two possibilities arise:

• if b(0) = 0 then there exists a constant invertible matrix S such that S -1 M(0)S =           0 1 0 0 0 0 0 0 0           .
Then M(x) is pointwise similar to J(x), even in regions containing x = 0.

• if b(0) 0 then the Jordan form of M(0) is J(0) =           0 1 0 0 0 1 0 0 0           so that J(x) is not holomor- phic at x = 0.
Hence, it might occur that the Jordan matrix is not holomorphic in some region or the Jordan matrix is itself holomorphic but, nevertheless, not holomorphically (although perhaps point-wise) similar to M(x). Below is an example which exhibits a system whose leading coefficient matrix has such complications: Example 4. [40, p. 57] Consider the following system with σ A = 0, i.e. ξ A = ε:

[A σ A ] ξ 2 A ∂F =           0 1 0 0 0 1 ξ A 0 x           F where A 0 (x) =           0 1 0 0 0 1 0 0 x           .
It follows from the form of A 0 (x) that the origin is possibly a turning point. We remark that this system is equivalent to [START_REF] Barkatou | Formal Solutions of Linear Differential Systems with Essential Singularities in their Coefficients[END_REF], as it is obtained from the latter by setting

F = [ f, ξ 2 A ∂ f, ξ 4 A ∂ 2 f ] T .
We recall the following proposition which we first gave in [START_REF] Abbas | On the Reduction of Singularly-Perturbed Linear Differential Systems[END_REF], and we refine its proof:

Proposition 2. [2, Section 3] Consider system (8) given by [A σ A ] ξ h A x p ∂F = A(x, ξ A ) F = ∞ k=0 A k (x)ξ k A F with ξ A = x σ A ε.
Suppose that the leading constant matrix A 0,0 is nilpotent and A 0 (x) has at least one nonzero eigenvalue. Then there exists an invertible polynomial transformation T in a root of x such that the transformation F = TG and a re-adjustment of the independent variable x results in a system whose leading constant matrix has at least one nonzero constant eigenvalue.

Example 5. Consider the matricial form of Weber's equation (σ A = 0, p A = 0) given by:

[A σ A ] ξ A dF dx = 0 1 x 2 0 F, with ξ A = x σ A ε. Let T (x) = 1 0 0 x then F = TG yields: [B σ A ] ξ dG dx = 0 x x -ξ A x -1 G.
Upon factorizing x, we have:

ξ A x -1 dG dx = 0 1 1 -ξ A x -2 G.
We compute σ Ã = -2 and thus set ξ Ã = x -2 ε. The above system can be rewritten as:

[ Ãσ Ã ] ξ Ã x dG dx = 0 1 1 -ξ Ã G.
Clearly, Ã0,0 has 2 constant eigenvalues and so the splitting lemma can be applied to decouple this system into two scalar equations.

Proof. (Proposition 2) The eigenvalues of A 0 (x) admit a formal expansion in the fractional powers of x in the neighborhood of x = 0 (see, e.g., [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]). We are interested only in their first nonzero terms. Let µ(x) = ∞ j=0 µ j x j/s be a nonzero eigenvalue of A 0 (x) with s ∈ N * , j ∧ s = 1, and whose leading exponent, i.e. smallest j/s for which µ j 0, is minimal among the other nonzero eigenvalues. Without loss of generality, we can assume that s = 1, otherwise we set x = t s . By [START_REF] Barkatou | A Rational Version of Moser's Algorithm[END_REF][START_REF] Jeannerod | A Reduction Algorithm for Matrices Depending on a Parameter[END_REF], there exists T ∈ GL n (C(x)) such that for B 0 (x) = T -1 A 0 (x)T = B 0,0 + B 1,0 x + . . . :

θ B 0 (x) (λ) = x rank(B 0,0 ) det(λI + B 0,0 x + B 1,0 )| x=0
does not vanish identically in λ. Let ν > 0 denote the valuation of B 0 (x) in x. By [24, Proposition 1], there are ndeg(θ) eigenvalues of B 0 (x) whose leading exponents lie in [ν, ν + 1[, and deg(θ) eigenvalues for which the leading exponent is equal to or greater than ν + 1. Then, applying the transformation

F = T (x)G to the system [A σ A ] yields [B σ A ] ξ h A x p ∂G = B(x, ξ A ) G = ∞ k=0 B k (x)ξ k A G where B 0 (x) := x ν (B 0,0 + B 1,0 x+. . . ) and B 0,0 has n-deg(θ) eigenvalues whose leading exponents lie in [0, 1[. Let Ã(x, ξ A ) = x -ν B(x, ξ A ) then the order of poles introduced in x in Ã(x, ξ A ) is at worst (span x (T )+ν)
, where span x (T ) is the difference between the valuation and the degree of the polynomial transformation T (x) in x. Hence, σ Ã ≥ σ Aspan x (T ) -ν and [B σ A ] can be rewritten as

[ Ãσ à ] ξ h à x p ∂G = Ã(x, ξ Ã) G = ∞ k=0 Ãk (x)ξ k à G where ξ à = x σ à ε, p ∈ Z, Ãk (x) ∈ M n ([[x]]), Ã0 (x) = x -ν B 0 (x)
, and Ã0,0 = B 0,0 . Then µ(x) is an eigenvalue of Ã0 (x) with a minimal leading exponent and hence it is among those whose leading exponents lie in [0, 1[. By our assumption s = 1, and hence the leading exponent of µ(x) is zero and µ 0 0. Since µ 0 is a nonzero eigenvalue of Ã0,0 , it follows that the Ã0,0 is non-nilpotent.

We remark that the eigenvalues of A 0 (x) are the roots of the algebraic scalar equation f (x, λ) = det(A 0 (x)-λI n ) = 0 and can be computed by Newton-Puiseux algorithm. The sought polynomial transformation T can be computed via Isolde, miniIsolde, or Lindalg.

Remark 1. Proposition 2 leads to the following observation about the detection of a turning point: Consider an input system [A σ A ] given by (3) with σ A = 0, p A = 0, and restraining index

ρ A = -1/σ A . If [A σ A ]
has a turning point at x = 0 then, by the end of the process of formal reduction of this system (i.e. whenever all decoupled subsystems are either of non-positive εrank or are scalar equations), at least one of its decoupled subsystems has a finite restraining index.

Example 6. Consider the following system whose σ A = 0 and ξ A = ε:

[A σ A ] ξ 2 A ∂F =           0 1 0 0 0 1 ξ A x 0           F where A 0 (x) =           0 1 0 0 0 1 0 x 0           . We first compute s = 2. We set x = t 2 and compute T =           0 0 1 0 t 0 0 0 t 2          
. Or equivalently, we consider

T =           0 0 1 0 x 1/2 0 0 0 x          
. Then F = TG yields:

ξ 2 A ∂G = x 1/2 {           0 1 0 1 0 0 0 1 0           +           0 0 0 0 0 0 1 0 0           x -3/2 ξ A +           0 0 0 0 -1 0 0 0 -2           x -1 ξ 2 A } G,
which we rewrite as:

[ Ãσ Ã ] ξ 2 Ã x 5/2 ∂G = {           0 1 0 1 0 0 0 1 0           +           0 0 0 0 0 0 1 0 0           ξ Ã +           0 0 0 0 -x 2 0 0 0 -2x 2           ξ 2 Ã } G,
where ξ Ã = x -3/2 ε. Now that the leading term is a constant matrix with three distinct eigenvalues, we can proceed by applying the Splitting of Theorem 1.

In the above example, the system could be decoupled thoroughly into three scalar equations. However, we might encounter different scenarios as well, as illustrated in the following example.

Example 7. [40, p. 57] We recall the system of Example 9 (which is equivalent to (12)) with ξ A = ε:

[A σ A ] ξ 2 A ∂F =           0 1 0 0 0 1 ξ A 0 x           F.
As mentioned before, we have a turning point at x = 0. A 0,0 is nilpotent and the eigenvalues of A 0 (x) are 0 and x, whence s = 1. Let T = diag(1, x, x 2 ) then F = TG yields:

ξ 2 A ∂G = x{           0 1 0 0 0 1 0 0 1           +           0 0 0 0 0 0 1 0 0           x -3 ξ A +           0 0 0 0 -1 0 0 0 -2           x -2 ξ 2 A } G.
Setting ξ Ã = x -3 ε, the former can be rewritten equivalently as:

[ Ãσ Ã ] ξ 2 Ã x 5 ∂G = {           0 1 0 0 0 1 0 0 1           +           0 0 0 0 0 0 1 0 0           ξ Ã +           0 0 0 0 -x 4 0 0 0 -2x 4           ξ 2 Ã } G.
The leading constant matrix Ã0,0 is no longer nilpotent. Hence the system can be decoupled into two subsystems upon setting G = T W where

T =           1 0 1 0 1 1 0 0 1           +           -1 -1 0 -1 -1 -2 -1 -1 0           ξ Ã + O(ξ 2 Ã).
The resulting equivalent system then consists of the two decoupled lower dimension systems where

W = [W 1 , W 2 ] T , σ B = σ C = σ Ã = -3: [B σ B ] ξ 2 B x 5 ∂W 1 = { 0 1 0 0 + -1 0 -1 0 ξ B + 1 -1 1 -1 + x 4 ξ 2 B + O(ξ 3 B ) } W 1 . [C σ C ] ξ 2 C x 5 ∂W 2 = {1 + ξ C + (1 + 2x 4 )ξ 2 C + O(ξ 3 C ) } W 2 .
The exponential part (and hence formal solution) of the second subsystem is clearly ξ

-2 C x -5 (1+ O(ξ C )) dx = 1 2 ε -2 x 2 (1 + O(ε -2 x 2 )
). One can observe that the eigenvalue x of the leading matrix coefficient of the input system [A σ A ] is recovered as expected. This is in accordance with the exponential parts obtained for [START_REF] Barkatou | Formal Solutions of Linear Differential Systems with Essential Singularities in their Coefficients[END_REF]. As for the first subsystem, B 0 (x) and B 0,0 are simultaneously nilpotent. Due to this dual nilpotency, in order to proceed in the formal reduction of [B σ B ], we will make use of the ε-rank reduction in the following subsection.

ε-Rank reduction

We consider again system [A σ A ] given by [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF]. We assume without any loss of generality that A 0 (x) and A 0,0 are simultaneously nilpotent. In this section, we investigate the ε-rank reduction of the system, i.e. we seek to determine the minimal integer ε-rank among all systems equivalent to [A σ A ]. If this minimal integer rank turns out to be non-positive then we continue the reduction by treating the system as in Section 2.1.1. Otherwise, the minimal integer rank gives an upper bound for the ε-exponential order, which allows us to proceed to Section 4.

In analogy to its unperturbed counterpart, we define the ε-Moser rank and ε-Moser invariant of system [A σ A ] to be the following rational numbers respectively:

m ε (A) = max (0, h + rank (A 0 (x)) n ) and µ ε (A) = min {m ε (T [A σ A ]) for all T in GL n (K)}. Definition 4. System [A σ A ] (the matrix A(x, ξ A ) respectively) is called ε-reducible if m ε (A) > µ ε (A).
Otherwise it is said to be ε-irreducible.

If m ε (A) ≤ 1, then h = 0. Hence, we restrict our attention to the case of m ε (A) > 1. This definition is not to be mixed neither with the usual sense of Moser-irreducible unperturbed system [A]; nor with the sense of reduced system already employed in the literature, i.e. the system whose matrix is the leading coefficient matrix of [A σ A ]: x p ∂F = A 0 (x)F. In fact, the motivation behind Moser's work in [START_REF] Moser | The Order of a Singularity in Fuchs' Theory[END_REF] was to determine the nature of the singularity (regular or irregular) of an unperturbed system

x p ∂F = A(x)F, A(x) ∈ M n (C[[x]]
). Consequently, the reduction of the so-called Poincaré rank p to its minimal integer value was investigated. However, given our perturbed system [A σ A ], it seems more plausible to reduce h, rather than p, to its minimal integer value, since a system with h ≤ 0 can be treated as in Subsection 2.1.

Since µ ε (A) cannot be computed from the outset, we aim in this section to generalize Moser's criterion to detect whether the ε-rank of a given perturbed system is minimal among all equivalent systems. We remark that Moser's criterion has been generalized as well to linear functional matrix equations in [6], and borrowed from the theory of differential systems in [START_REF] Jeannerod | A Reduction Algorithm for Matrices Depending on a Parameter[END_REF] to investigate efficient algorithmic resolution of the perturbed algebraic eigenvalue-eigenvector problem. However, despite their utility and efficiency for such univariate systems, algorithms based on this criterion are not considered so far over bivariate fields.

In [START_REF] Abbas | On the Reduction of Singularly-Perturbed Linear Differential Systems[END_REF], we generalize the work of [START_REF] Barkatou | A Rational Version of Moser's Algorithm[END_REF][START_REF] Moser | The Order of a Singularity in Fuchs' Theory[END_REF] for the case σ A = 0. Herein, we establish the following theorem without any restriction on σ A (see [28, Section 5.5, p. 101] and Appendix C for the proof). This yields an algorithm which takes an input system [A σ A ], and outputs a transformation T ∈ GL n (K) and an equivalent system

[ Ãσ Ã ] := T [A σ A ] which has a minimal ε-rank among all systems equivalent to [A σ A ].
Theorem 2. Consider system [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF] given by

[A σ A ] ξ h A x p ∂F = A(x, ξ A )F = ∞ k=0 A k (x)ξ A k F,
and suppose that h > 1 and m ε (A) > 1. The polynomial

θ A (λ) := ξ A rank (A 0 (x)) det(λI + A 0 (x) ξ A + A 1 (x))| ξ A =0 (15)
vanishes identically in λ, if and only if there exists a transformation F = TG, T ∈ GL n (K), such that the equivalent system [ Ãσ Ã ] := T [A σ A ] has either a strictly lower ε-rank or a leading coefficient matrix with a strictly lower algebraic rank. Moreover, σ T , σ Ã ≥ σ where

σ =        σ A + p-1 h if p < 1 σ A otherwise . Example 8. Consider the system [A σ A ] ξ 4 A ∂F = A(x, ξ A )F with σ A = 0 and A(x, ξ A ) =               2xξ 3 A 3x 2 ξ 4 2xξ 2 A (2x + 1)ξ 5 A 0 ξ 4 A 0 0 0 0 ξ 2 A 0 -2x 0 0 0               .
Our algorithm computes the transformation F = TG, where T =

              0 ξ 3 T 0 0 0 0 ξ T 0 ξ 3 T 0 0 0 0 0 0 1               , with σ T = σ A = 0 which results in an equivalent ε-irreducible system [A σ à ] ξ 2 Ã∂G = Ã(x, ξ Ã)G where Ã(x, ξ Ã) =                1 0 0 0 2x 2xξ à 3x 2 2x + 1 0 0 ξ 2 à 0 0 -2xξ à 0 0               
, and σ Ã = 0. One can observe that the ε-rank is diminished by two.

Example 9. Consider the first subsystem resulting from the reduction of Example 7:

[B σ B ] ξ 2 B x 5 ∂W 1 = { 0 1 0 0 + -1 0 -1 0 ξ B + 1 -1 1 -1 + x 4 ξ 2 B + O(ξ 3 B ) } W 1 .
It is easy to see that the leading coefficient and constant matrices coincide and that θ B (λ) = 1.

Hence, the ε-rank which is equal to 2 is minimal.

To proceed in the formal reduction of Example 9, one needs to introduce a ramification in ε (resp. ξ B ). The computation of the necessary ramification will be discussed in Section 4. To establish the results therein, the ε-Moser invariant of equivalent scalar equations is needed. We thus introduce it in Subsection 3.3. Before proceeding however, we discuss the case h = 1 which is not tackled by Theorem 2.

The case h = 1

The limitation of Theorem 2 to h > 1 is of technical nature (see the proof of necessity). In this subsection, we treat the case h = 1. We consider again the system given by ( 8):

[A σ A ] ξ A x p ∂F = A(x, ξ A )F = ∞ k=0 A k (x)ξ k A F.
Let h true denote the minimal integer ε-rank which can be attained upon applying a transformation in GL n (K) and ω ε (A)s denote the ε-exponential order of [A σ A ]. In this subsection, we give a method to decide whether h true = 1 or h true = 0. We first apply the substitution ε = εn+1 to [A σ A ]. We denote its ε-rank by h. Since h = n + 1 > 1, we can then apply the ε-rank reduction of Theorem 2 and we denote by htrue and ω ε the true ε-rank of the resulting system and its ε-exponential order respectively. Under these notations, we have:

Lemma 3. If htrue = 1 then h true = 0.
Proof. Let ω ε = d where and d are coprime natural numbers. Then, on the one hand, d, ≥ 0, d < n and ω ε ≤ h true . On the other hand, ω ε = (n + 1)ω ε and ω ε ≤ htrue (the minimal ε-rank bounds the ε-exponential order, see Corollary 2). Hence, we have,

ω ε = (n + 1)ω ε = (n + 1) d ≥ (n + 1) n ≥ (1 + 1 n ) . 20
Hence, if 0 then

1 < (1 + 1 n ) ≤ ωε ≤ htrue .
Lemma 3 moves the problem from the case h = 1 to the case h > 1. It can be restated as follows: Within the formal reduction, whenever the case h = 1 is encountered, a ramification ε = εn+1 is applied so that the reduction of Theorem 2 can be applied. After this reduction, one of the following two cases arises:

• If the ε-rank of the resulting system is less than or equal to one then h true of the original system is zero, and so we stop our reduction (the ε-exponential part have been computed completely). If one wishes to continue reduction, then the classical Arnold-Wasow approach can be employed: First, put the system in Arnold's form [START_REF] Arnold | On Matrices Depending on Parameters[END_REF] and then use the Iwano-Sibuya's polygon to determine the proper polynomial transformation required to arrive at an equivalent system whose ε-rank is zero [START_REF] Iwano | Reduction of the Order of a Linear Ordinary Differential Equation Containing a Small Parameter[END_REF]. The reduction then continues as explained in Subsection 2.1.1.

• Otherwise, we proceed with the reduction: If the leading coefficient matrix of the ramified system has two distinct eigenvalues, then we treat turning points (if any) and apply splitting lemma. Otherwise, we proceed to Section 4.

We remark that in the implementation, we first try to find whether there exist constant vectors in the left null space of G A (λ) (see Appendix C). If such vectors do exist, we use them to construct a transformation which might reduce the rank. If the ε-rank remains one after all the constant vectors are exhausted, we compute a candidate for the ramification in ξ A from the characteristic polynomial. If this candidate, after applying again the ε-rank reduction results in a system with non-nilpotent leading coefficient, we keep it. Otherwise, we use Arnold-Wasow approach.

ε-Moser invariant of a scalar equation

Due to the equivalence between a scalar differential equation [a σ a ] and its companion system [A σ A ], it is natural to define and question the ε-Moser invariant of the former, in the hope to gain more insight into the problem. This is the goal of this subsection which is fulfilled by generalizing the analogous notion discussed for unperturbed scalar linear differential equations in [START_REF] Moser | The Order of a Singularity in Fuchs' Theory[END_REF]Part IV].

We consider again the singularly-perturbed linear differential equations given by ( 9)

[a σ a ] ∂ n f + a n-1 (x, ξ a ) ∂ n-1 f + . . . + a 1 (x, ξ a ) ∂ f + a 0 (x, ξ a ) f = 0,
where a n (x, ξ a ) = 1. We prove the following proposition:

Proposition 3. Given the differential equation [a σ a ]. Let τ, ν be the smallest integers such that val ε (a i (x, ξ a )) ≥ (i -n)(τ -1) -ν.
In other words, let

κ = min { ∈ N| val ε (a i ) + (n -i) ≥ 0, 0 ≤ i ≤ n} (16) ν = max {(i -n)(τ -1) -val ε (a i ), 0 ≤ i ≤ n} (17)
Then the ε-Moser invariant of the system given by the associated companion matrix is Example 10. Let σ a = -1 and ξ a = x σ a ε. We consider the scalar differential equation

µ ε (a) = κ + ν n .
[a σ a ] ∂ 5 f + (xξ -3 a + x)∂ 4 f + 2x 2 ξ -1 a ∂ 3 f + (ξ -3 a + 1)∂ 2 f + (-3ξ -4 a + x 2 ξ -2 a )∂ f -ξ -2 a f = 0.
First, we plot the points (i, val ε (a i )), 0 ≤ i ≤ n = 5 and the ε-polygon (see Figure 4). Next, we plot the straight line passing through the point (n, 0) = (5, 0) with smallest integer slope such that it stays below all these point. This yields V = 3U -15 = 3(U -5) and τ = 3. Then, we plot the straight line of slope τ -1 = 2 which stays below these points: V = 2U -11 = 2(U -5) -1. This yields ν = 1. Finally, we consider γ = (γ 0 , . . . , γ 4 ) = (-11, -9, -7, -5, -3). We observe that i 0 = n -ν = 4 and so the equivalent system is given by:

[A σ A ] xξ 3 A ∂W =                    11ξ 3 A xξ A 0 0 0 0 9ξ 3 A xξ A 0 0 0 0 7ξ 3 A xξ A 0 0 0 0 5ξ 3 A xξ A xξ 9 A 3xξ 5 A -x 3 ξ 7 A -xξ 4 A -xξ 7 A -2x 3 ξ 4 A -x 2 -(x 2 -3)ξ 3 A                    W, with σ A = σ a = -1. One can verify that the system [A σ A ] is ε-irreducible since θ A (λ) = xλ 4 .
Proof. (Proposition 3) For clarlity, we set ξ := ξ a . We define

γ i = max {κ(i -n), (κ -1)(i -n) -ν}, 0 ≤ i ≤ n -1.
Then, by construction, γ i ≤ val ε (a i ) and equality is attained for at least one i ≥ n -ν (otherwise ν and κ can be minimized). Let i 0 represent the smallest integer such that γ i 0 = κ(i 0n), then i 0 = n -ν. Geometrically, the γ i 's represent a broken line dominated by the (i, val ε (a i )). They will aid in the construction of an ε-irreducible system which is equivalent to [START_REF] Barkatou | A Rational Version of Moser's Algorithm[END_REF]. In fact, let

w i+1 = ξ γ i ∂ i f, 0 ≤ i ≤ n -1.
then we have

                     ∂w i+1 = ξ -κ [        ξw i+2 0 ≤ i < i 0 w i+2 i 0 ≤ i ≤ n -2 + ξ κ σ a γ i x w i+1 ], w n = ξ γ n-1 ∂ n-1 f = ξ -κ ∂ n-1 f, ∂w n = σ a γ n-1 x w n + ξ -κ ∂ n f = ξ -κ [ξ κ σ a γ n-1 x w n + n-1 i=0 α i (x, ε)w i+1 ].
where

α i (x, ξ) = -a i (x, ε)ξ -γ i = ∞ k∈Z a i,k (x)ξ k-γ i . Let W = (w 1 , . . . w n ) T and A(x, ξ) =                                            σ a γ 0 ξ κ xξ . . . σ a γ i 0 -1 ξ κ xξ σ a γ i 0 ξ κ x σ a γ i 0 +1 ξ κ x xα 0 xα 1 . . . xα i 0 xα i 0 +1 . . . xα n-1 + σ a γ n-1 ξ κ                                           
.

Then one can verify that,

[A σ a ] xξ κ ∂W = A(x, ξ)W.
It remains to prove that this system is ε-irreducible. We first remark that A(x, 0) has rank ν = n-i 0 due to the linear independence of its last ν rows. In fact, it's clear that it has ν -1 linearly independent rows (the rows with the x's). Moreover,

α i (x, ξ) = -a i (x, ξ)ξ -γ i = -a i (x, ξ)ξ -val ε (a i ) ,
and so α i (x, 0) 0 for at least one i ∈ {0, . . . , i 0 }. Thus the last ν rows are linearly independent.

Setting ξ = x σ a ε we have:

θ A (λ) = ξ ν det(λI + A 0 (x) ξ + A 1 (x))| ξ=0 = ξ ν λ + δ(σ a γ 0 ) x . . . λ + δ(σ a γ i 0 -1 ) x λ + δ(σ a γ i 0 ) x ξ λ + δ(σ a γ i 0 +1 ) x ξ xα 0 ξ xα 1 ξ
. . .

xα i 0 ξ xα i 0 +1 ξ . . . λ + δ(σ a γ n-1 ) + xα n-1 ξ | ξ=0 = ξ ν-n+i 0 λ + δ(σ a γ 0 ) x . . . λ + δ(σ a γ i 0 -1 ) x xα 0 xα 1 . . . xα i 0 ,
where δ = 1 if κ = 1 and is zero otherwise. Hence, since α i (x, 0) 0 for at least one i ∈ {0, . . . , i 0 }, θ A (λ) does not vanish identically in λ and its highest possible degree is i 0 . It follows from Theorem 2 that the system is ε-irreducible and µ ε (a

) = µ ε (A) = κ + ν n . Corollary 1. Given an ε-irreducible system [A σ A ] ξ h A x p ∂F = A(x, ξ A )F. Let [a σ a ]
be the scalar equation corresponding to a companion system equivalent to [A σ A ]. Then for κ and ν defined in [START_REF] Bender | Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory[END_REF] and [START_REF] Béringer | Multi-variate polynomials and Newton-Puiseux expansions[END_REF] we have: κ = h and ν = rank(A 0 ).

Proof. By the ε-irreducibility of [A σ

A ] and its equivalence to the companion system, we have:

κ + ν n = µ ε (a) = µ ε (A) = m ε (A) = h + r n
Consequently, κ = h and ν = rank(A 0 ).

Corollary 2. Under the notations of Corollary 1, we have:

κ -1 + ν n ≤ ω ε (A) ≤ κ and h -1 + rank(A 0 ) n ≤ ω ε (A) ≤ h.
Proof. Follows from ( 16), [START_REF] Béringer | Multi-variate polynomials and Newton-Puiseux expansions[END_REF], and (13).

Computing the ε-exponential order

We consider again system [A σ A ] which is given by ( 8):

[A σ A ] x p ξ h A ∂F = A(x, ξ A )F = ∞ k=0 A k (x)ξ k A F. (18)
In this section, we assume without loss of generality that h > 0, [A σ A ] is ε-irreducible, and A 0 (x), A 0,0 are both nilpotent. The goal of this section is to compute the ε-exponential order ω ε (A) and ε-polynomials of the system, which give indispensable information within formal reduction. In particular, the former determines the ramification in ε which can lead to a system whose leading matrix is non-nilpotent, so that the process of formal reduction can be resumed. This leads to the recursive Algorithm 1. Eventually, we can construct fundamental matrices of formal solutions in any given subdomain. We illustrate this algorithm by examples and motivate items for further investigation in Appendix E. Let det(λI -A(x, ε)

x p ξ h A ) = λ n + α n-1 (x, ξ A )λ n-1 + • • • + α 0 (x, ξ A ). ( 19
)
such that α n = 1 and

α i (x, ξ A ) = ∞ j=val ε (α i ) α i, j (x)ξ j A ∈ K for i ∈ {0, . . . , n}. We define the ε- polygon N ε (A) of [A σ A ]
by taking P ε (A) of Section 2.1.2 to be the union of P(i, val ε (α i (x, ξ A )) for i ∈ {0, . . . , n}. We can then prove the following:

Theorem 3. Consider the ε-irreducible system [A σ A ]
given by [START_REF] Chen | Solutions Formelles de Systemes d'Equations Differentielles Lineaires Ordinaires Homogenes[END_REF] with h > 0 and [START_REF] Churchill | Cyclic vectors[END_REF].

If h > n -rank(A 0 (x)) then the ε-exponential order of [A σ A ] is given by ω ε (A) = max 0≤i<n ( -val ε (α i ) n -i ),
and its corresponding ε-polynomial is given by the algebraic equation

E ε (X) = k=0 x σ A •val ε (α i k ) α i k ,val ε (α i k ) X (i k -i 0 )
where 0 ≤ i 0 < i 1 < • • • < i = n denote the integers i for which ω ε (ni) = -val ε (α i ) (i.e. lie on the edge of slope ω ε of the ε-polygon

N ε (A) of [A σ A ]); and α i,val ε (α i ) (x) = ξ -val ε (α i ) A α i (x, ξ A )| ξ A =0 .
The proof can be established by adapting to the parametrized setting the proofs of [8, Lemma 3, Lemma 4, Proposition 1, Theorem 1] (see Appendix D). Not only does this theorem compute these invariants of the system, but it also allows a further reduction of the system as follows: Suppose that ω ε (A) = d where , d are relatively prime positive integers. One can then set ε = ε 1/d (or equivalently 3ξA = ξ 1/d A x σ(d-1)/d ) in [A σ A ] and perform the ε-rank reduction (the minimal εrank is ). This will lead to an equivalent system whose ε-rank is and whose leading coefficient matrix has at least d distinct eigenvalues. The system can thus be decoupled into subsystems of lower dimensions. By repeating these procedures for each of the resulting subsystems, we can decouple the initial system into subsystem(s) of lower dimension(s) or zero ε-rank. This leads to the recursive algorithm of Section 5. We remark that the condition h > nrank(A 0 (x)) in Theorem 3 is non-restrictive. We can always arrive to a system satisfying this condition by a generalization of [8, Lemma 5], which is based on applying the affinity (U, V) → (U, d V) for some integer d to the ε-polygon. In fact, we can prove the following lemma:

(A) = ω ε( Ã) d . Hence, ω ε( Ã) ≥ d(h -1 + r n ) ≥ n. But h ≥ ω ε( Ã) and r ≥ 1, which yields h + r ≥ ω ε( Ã) + 1 ≥ d(h -1 + r n ) + 1 ≥ n + 1 > n.

Formal reduction algorithm

With the algorithms of Splitting lemma, turning point resolution, ε-rank reduction, and as can be verified by Algorithm 1 below, we have re-established constructively the well-known general form for a fundamental matrix of formal outer solutions for an input system [A ε ] given by (3):

F = ( ∞ k=0 Φ k (x 1/s ) (x σ ε) k/d ) exp( Q(x 1/s , ε -1/d )), (20)
where s, d are positive integers; σ is a nonpositive rational number; Q is a diagonal matrix whose entries are polynomials in ε -1/d with coefficients in C((x 1/s )). We refer to Q := Q as the εexponential part (logarithms in a root of x might arise as a result of integration); and the entries of the Φ k (x 1/s )'s are root-meromorphic in x (see [START_REF] Schaefke | On the Formal Fundamental Solution of a Differential Equation Depending on a Parameter[END_REF]Introduction] or [START_REF] Schaefke | Formal Fundamental Solutions of Irregular Singular Differential Equations Depending Upon Parameters[END_REF]).

Remark 3. Under the notations and statements of Theorem 3, the leading term of Q is given by

1 ε ω ε (A) diag(X 1 (x), . . . , X deg (E ε ) (x), 0, . . . , 0) dx,
where the X i 's denote the roots of the ε-polynomial E ε (X). A similar statement can be stated in terms of the eigenvalues of the leading matrix coefficient, x p , and ξ h A . We sum up our main results in the formal reduction algorithm, Algorithm 1, which computes the ε-exponential part and consequently a fundamental matrix of formal solutions [START_REF] Fruchard | Composite Asymptotic Expansions[END_REF] in a given subdomain. We recall that:

• If n = 1 then we proceed by integration up to the first h terms;

• If h ≤ 0 then we follow Subsection 2.1.1 (using ISOLDE, miniIsolde, or Lindalg).

Conclusion and further investigations

In this article, we give an algorithm, implemented in Maple 4 , which computes a fundamental matrix of outer formal solutions for singularly-perturbed linear differential systems in a neighborhood of a turning point. The subprocedures discussed are stand-alone algorithms (Splitting, rank reduction, and resolution of turning points). The formal reduction algorithm is based on the generalization of the algorithm given in [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF] for the unperturbed counterparts of such systems.

Our results are presented in the formal setting. However, the growth of the order of poles in x is tracked within the formal reduction. This gives information on the restraining index of the system, and consequently, an adequate stretching transformation can be chosen and performed. This furnishes the first step in resolving Iwano's first problem and computing, for an input system, [ρ] which is given by ( 11):

[ρ] 0 = ρ 0 < ρ 1 < ρ 2 < • • • < ρ m .
In Appendix E, we try to motivate the employment of our proposed formal reduction algorithm to resolve the same problem for a general system: Suppose that we start the reduction with the input system [A σ A ], σ A = 0, given by ( 8) as an input. Then, upon applying Algorithm 1, we can determine the outer solutions and the final restraining index σ f inal which allows full reduction. If the restraining index is nonzero then we have a turning point and ρ m = -1/σ f inal . We can then perform the stretching τ = xε -ρ 1 in the input system [A σ A ], and apply to it Algorithm 1. This determines ρ 2 . We show an example where the same process can be repeated with ρ 2 to determine ρ 3 , and iteratively we reach a system whose final restraining index is zero. In future work, we hope to investigate this approach and its correctness.

Another field of investigation is the relation between the algebraic eigenvalues of the matrix A(x, ε) of such systems and the exponential part of the solution (see Example 21 and [28, Chapter 2]). One can then benefit from the existing work on fractional power series expansions of solutions of bivariate algebraic equations, to compute the restraining index (see, e.g. [START_REF] Béringer | Multi-variate polynomials and Newton-Puiseux expansions[END_REF][START_REF] Aroca | Power Series Solutions for Non-L inear PDEs[END_REF] and references therein). On the other hand, it would be interesting to investigate a differential-like reduction for a two-parameter perturbation of a JCF (see [28, Chapter 2] and references therein for the one-parameter case). In fact, one can observe that the main role in the reduction process is reserved to the similarity term of T [A σ A ], i.e. T -1 A σ A T rather than T -1 ∂T . Hence, for a non-differential operator where T [A σ A ] = T -1 A σ A T , the discussion is not expected to deviate substantially from the discussion presented here for a differential one.

The generalization of other notions and efficient algorithms can be investigated as well (see, e.g. simple system [START_REF] Barkatou | An algorithm computing the regular formal solutions of a system of linear differential equations[END_REF][START_REF] Barkatou | On k-simple forms of first-order linear differential systems and their computation[END_REF]

and [28, Appendix A]).

In Examples 18 and 20 we treat input systems with σ A < 0. The algorithms generalized herein have been generalized to treat systems with an essential singularity [START_REF] Barkatou | Formal Solutions of Linear Differential Systems with Essential Singularities in their Coefficients[END_REF] and difference systems as well(see, e.g. [START_REF] Abbas | Formal Reduction of a class of Pfaffian Systems in Two Variables[END_REF]6]). This motivates investigating the adaptation of our proposed algorithms at least to the difference setting.

There remain as well the questions which fall under the complexity, e.g. studying the complexity of this algorithm; the bounds on number of coefficients needed in computations; developing efficient algorithms for computing a cyclic vector; and comparing it to the results of our direct algorithm.

And finally, to give a full answer on the behavior of the solutions, the related problems of connection, matching, and secondary turning points, are to be studied as well (see, e.g. [START_REF] Fruchard | Composite Asymptotic Expansions[END_REF][START_REF] Hulek | Systémes d' équations différentielles linéaires singuliérement perturbées et développements asymptotiques combinés[END_REF] and references therein).

M i (a): Multiplies the i th column and i th row by a and 1/a respectively. It also adds (-ξ h A x p ∂a a ) to the diagonal entry in the (i, i) th position. C i, j (a): Adds to the j th column the i th column multiplied by a, adds to the i th row the j th row multiplied by -a, and adds (-ξ h A x p ∂a) to the entry in the (i, j) th position. R i, j (-a): Adds to the i th row the j th row multiplied by a, adds to the j th column the i th column multiplied by -a, and adds (ξ h A x p ∂a) to the entry in the (i, j) th position. A careful choice of these transformations will allow us to establish the following result constructively: Theorem 4. Consider the system [A σ A ] given by [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF]. Then there exists a transformation T which is a product of transformations of the forms (A.1) and (A.2) in a root of x, such that the equivalent system [ Ãσ à ] ξ hà x p∂F = ÃF has the block diagonal form:

Ã(x, ξ Ã) = diag( Ã11 (x, ξ Ã), . . . , Ã (x, ξ Ã))
where the Ã11 , . . . , Ã are companion matrices.

The proof of this theorem relies on two lemmas which we establish herein. Before proceeding, we can always assume, without loss of generality, that A 0,0 = A(x = 0, ξ A = 0) is in Frobenius canonical form, i.e.,

A 0,0 = diag(A 11 0,0 , . . . , A 0,0 ) (A.3)
where the block submatrices A 11 0,0 , . . . , A 0,0 are constant matrices in companion form, of dimensions

n 1 ≥ n 2 ≥ • • • ≥ n ≥ 1 (n 1 + n 2 + • • • + n = n).
Lemma 5. Consider the system [A ξ A ] given by (8) whose leading constant coefficient A 0,0 is in Frobenius canonical form. Then, by a sequence of elementary operations of the form (A.2), the equivalent system [ Ãσ à ] ξ h à x p ∂F = ÃF has the following block decomposition:

Ã(x, ξ Ã) = S (x, ξ Ã) U(x, ξ Ã) V(x, ξ Ã) W(x, ξ Ã) , (A.4)
where

• S (x, ξ Ã) is a n 1 -square matrix in companion form: S (x, ξ Ã) =                       0 1 0 0 0 . . . 1 s 1 s 2 . . . s n 1                       ; • U(x, ξ Ã), V(x, ξ Ã) are of dimensions n 1 × (n -n 1 ) and (n -n 1 ) × n 1 respectively, U 0,0 = O n 1 ×(n-n 1 ) , V 0,0 = O (n-n 1 )
×n 1 , and of the forms:

U(x, ξ Ã) =                  0 . . . 0 . . . . . . 0 . . . 0 u 1 . . . u n-n 1                  and V(x, ξ Ã) =             v 1 0 . . . 0 . . . . . . . . . v n-n 1 0 . . . 0             ; • W(x, ξ Ã) is a square matrix of dimension n -n 1 ;
• Furthermore, Ã0,0 = A 0,0 .

Proof. We first assume n 1 > 1. Otherwise, A(x, ξ A ) = [a i j ] 1≤i, j≤n is already in form (A.4). Since A 00 has the form (A.3), a i,i+1 (x = 0, ξ = 0) = 1. Hence, 1/a i,i+1 (x, ξ A ) ∈ R and 1/a i,i+1 (x = 0, ξ A = 0) = 1.

Step 1 For each row number i from 1 to n 1 -1:

(1.1) We use 1/a i,i+1 (x, ξ A ) as a pivot and apply M i+1 (1/a i,i+1 (x, ξ A )) to set the entries in the (i, i + 1) positions to 1. We remark that the term (-ξ h A x p (a)∂(1/a)) is added to the diagonal entry in the (i + 1, i + 1) th position.

(1.2) We now make use of these 1's to set the entries in the (i, j) positions for i ∈ {1, . . . , n 1 -1}, j ∈ {1, . . . , n}, j i + 1 to zeroes: For each row number j i + 1 from 1 to n, we apply C i+1, j (-a i j ). We remark that the term (+ξ h A x p ∂a) is added to the entry in the (i + 1, j) position. Clearly, the rows 1 to i -1 are not altered by these elementary operations.

Step 2 We have attained the anticipated form for the first n 1 rows, i.e. we have shown how to construct an equivalent system whose first n 1 rows have the described properties of S (x, ξ A ) U(x, ξ A ) . If n 1 < n then we can proceed to work on the last nn 1 rows.

(2.1) Again, we make use of these 1's created in Step [(1.1)] to set the entries in the (i, j) positions for i ∈ {n 1 + 1, . . . , n} and j ∈ {2, . . . , n 1 } to zeroes: For each i, j in these ranges, we apply R i, j-1 (-a i j ) successively for each i from n 1 + 1 to n with an inner loop of j from n 1 to 2. Evidently, each such operation does not alter the form created the inner loop shld start from n 1 because the term coming from the derivation settles in the (i, j -1) position and so will be deleted afterwards. I should specify how σ Ã changes w.r.t. σ A .

by Step 1 for the first n 1 rows, and by its preceeding operations within Step 2.

Finally, since each of the operations performed is of the form (A.2), for x = 0 and ξ A = 0, such an operation reduces to I n . Hence, the form of A 0,0 is preserved. Remark 4. One can observe that by a repetitve application of Lemma 5 on the diagonal blocks of dimensions n 1 ≥ n 2 ≥ • • • ≥ n , we arrive at the aformentioned Arnold form. Lemma 6. Consider the system [A σ A ] given by [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF]. Suppose that its leading constant coefficient A 0,0 is in Frobenius canonical form and A(x, ξ A ) is in form (A.4). Then, by a sequence of constant transformations, elementary operations of the form (A.2), and shearings of the form (A.1), the equivalent system [ Ãσ à ] ξ h à x p∂F = ÃF is such that has the following block diagonal form:

Ã(x, ξ Ã) = S (x, ξ Ã) O O W(x, ξ Ã) , (A.5)
where S (x, ξ Ã) is a r-square companion matrix where r ≥ n 1 and W(x, ξ Ã) is a nr square matrix.

Proof. If U(x, ξ A ) and V(x, ξ A ) of (A.4) are smiltaneously zero matrices or n 1 = n then A(x, ξ A ) is already in the form (A.5). Otherwise we suppose that U(x, ξ A ) is nonzero and limit our discussion to this case. If V(x, ξ A ) is nonzero as well, then we can proceed analogously. Let u j = ∞ k=0 u j,k (x)ξ k A for j ∈ {1, . . . , nn 1 } then we define β = min 1≤ j≤n-n 1 ( val ε (u j )) and α = min 1≤ j≤n-n 1 ( val x (u j,β )).

We remark that β ≥ 0, and if β = 0 then α > 0 since U 0,0 is a zero matrix. Let

T = diag(x α ξ β A I n 1 , I n-n 1 ) and [B σ B ] = T [A σ A ].
Then we have:

B(x, ξ B ) = S (x, ξ A ) -x p+α-1 ξ β+h A (α + βσ A )I n 1 x -α ξ -β A U(x, ξ A ) x α ξ β A V(x, ξ A ) W(x, ξ A ) . B 0 (x) = S 0 (x) x -α ξ -β A U 0 (x) O W 0 (x)
and

B 0,0 = S 0,0 U 0,0 O W 0,0
where S 0,0 is a companion matrix of dimension n 1 , U 0,0 is a companion matrix whose last row has at least one nonzero element, and W 0,0 is a block diagonal square matrix whose blocks are in companion form. Let e i denote the i th row of the identity matrix I n and B i 0,0 denote the i th power of B 0,0 . W then have:

e 1 B i 0,0 = i j=1 c i, j e j + e i+1 for i = 1, . . . , n 1 -1 e 1 B n 1 0,0 = n 1 j=1 c n 1 , j e j + n-n 1 j=1 u j, βe n 1 + j
where c i, j ∈ C. Since there exists at least one j ∈ {1, . . . , nn 1 } such that u j,β 0, then the vectors e 1 , e 1 B 0,0 , e 1 B 2 0,0 , e 1 B n 1 0,0 are linearly independent. Hence, denoting by r the dimension of the first block of the Frobenius form of B 0,0 , i.e. the degree of the minimal polynomial of B 0,0 , we have r > n 1 . Hence, by putting the leading constant matrix in Frobenius form, applying Lemma (5), and the procedure described herein, we either arrive at the desired form in a finite number of steps.

Proof. (Theorem 4) The proof can be attained by a recursive application of Lemmas 5 and 6.

Remark 5. We remark that this method can give a cyclic vector and that [ Ãσ Ã ] is not uniquely determined by

[A σ A ].
is a polynomial in (q, ∂q, . . . , ∂ i-1 q). Then, f (x, ξ) verifies ( 9) if and only if q(x, ξ) is a solution of n i=0 a i (x, ξ)P i (q). We denote by e the polar order of q in ξ (ε-order), and write q(x, ξ) = ξ e (s(x) + O(ξ)), q 0. Then, ∀i ∈ {0, . . . , n}, val ε (P i ) = -ie and in fact P i (q(x, ξ)) = ξ -id (s(x) i + O(ξ)).

Thus,

n i=0 a i (x, ξ)P i (q) = n i=0 ξ ν i -id W i (q), where val ε (W i ) ≥ 0. Let v(e) = min 0≤i≤n (ν i -ie) and I(e) = {i ∈ {0, . . . , n} | val ε (a i -ie) = v(e)}. Hence, n i=0 a i (x, ξ)P i (q) = ξ v(e) ( i∈I(e) a i,ν i (x)s i (x) + O(ξ)).
Thus, q(x, ξ) = ξ e (s(x) + O(ξ)), q 0 satisfies n i=0 a i (x, ξ)P i (q) ≡ 0 if and only if s(x) is a non zero solution of the equation

i∈I(e) a i,ν i (x)X i (x) = 0.
This equation has non trivial solutions if and only if |I(e)| > 1, i.e. if and only if e ∈ {e 1 , . . . , e }. And then, X(x) is one of the nonzero roots of the associated determinant equation.

Proof. (Theorem 1) We proceed in two steps:

Step 1: We first block-diagonalize the leading matrix coefficient: Suppose that there exists

T 0 (x) ∈ Gl n (C[[x]]) s.t. the substitution F = T 0 (x)G in [A σ A ] yields an equivalent system [ Ãσ A ] ξ h A x p ∂G = Ã(x, ξ A )G where Ã(x, ξ A ) = ∞ k=0 Ãk (x)ξ k A , Ãk (x) ∈ M n (C[[x]]
) for all k ≥ 0, and the leading coefficient matrix Ã0 (x) is non-zero and has the following block-form in accordance with A 0,0 :

Ã0 (x) = Ã11 0 (x) O O Ã22 0 (x)
.

Then we have:

A T 0 (x) -T 0 (x) Ã = ξ h A x p ∂T 0 (x) which yields A k (x) T 0 (x) -T 0 (x) Ãk (x) = δ kh ∂T 0 (x), where δ kh = 1 if k = h and δ kh = 0 otherwise. (B.1)
In particular (since h > 0):

A 0 (x) T 0 (x) -T 0 (x) Ã0 (x) = 0. (B.2)
We further assume that

T 0 (x) ∈ M n (C[[x]]
) is in the following form

T 0 (x) = I T 12 0 (x) T 21 0 (x) I .
Then, for 1 ≤ ς ≤ 2, we have:

       A 0 (x) -Ã 0 (x) + A ς 0 (x) T ς 0 (x) = O A ς 0 (x) + A ςς 0 (x) T ς 0 (x) -T ς 0 (x) Ã 0 (x) = O (B.3)
Inserting the series expansions A 0 (x) = ∞ i=0 A i,0 x i and Ã0 (x) = ∞ i=0 Ãi,0 x i in (B.3), and equating the power-like coefficients, we get for i = 0:

       A 0,0 = Ã 0,0 A ςς 0,0 T ς 0,0 -T ς 0,0 Ã 0,0 = O
which are satisfied by setting Ãςς 0,0 = A ςς 00 and T ςς 0,0 = I, T ς 0,0 = O. And for i ≥ 1 we get:

A ςς 0,0 T ς i,0 -T ς i,0 A 0,0 = -A ς i,0 - i-1 j=1 (A ςς i-j,0 T ς j,0 -T ς j,0 Ã i-j,0 ) (B.4) Ã i,0 = A i,0 + i-1 j=1 A ς i-j,0 T ς j,0 . (B.5)
It's clear that (B.4) is a set of Sylvester matrix equations that possess a unique solution due to the assumption on the disjoint spectra of A 11 0,0 and A 22 0,0 (see, e.g., [START_REF] Balser | Formal Power Series and Linear Systems of Meromorphic Ordinary Differential Equations[END_REF]Appendix A.1,). Remarking that the right hand side depends solely on the A j,0 , T j,0 , à j,0 with j < i, equations (B.4) and (B.5) are successively soluble. Hence, such T 0 (x) ∈ Gl n (C[[x]]) can be constructed. Clearly, p T = ν T = 0 and σ T ≥ σ A . Moreover, it follows from (B.1) that ν à = 0, p à = 0, and σ à ≥ σ A .

Step 2: Due to the first step, we can assume without loss of generality that the leading coefficient matrix A 0 (x) of the system [A σ A ] has a block-diagonal form in accordance with that of A 0,0 . It suffices to seek T and à in the form T = ∞ k=0 T k (x)ξ and à = ∞ k=0 Ãk (x)ξ respectively where T k (x), Ãk (x) ∈ GL n (C[[x]]) for all k > 0. We set T 0 (x) = I n , Ã0 (x) = A 0 (x), and we rewrite

A = ∞ k=0 A k (x)ξ A = ∞ k=0 B k (x)ξ, where B k (x) = x (σ A -σ)k A k (x). Since σ A ≥ σ, it follows that B k (x) ∈ M n (C[[x]]
) for all k > 0. With the expansions in ξ, we get the following set of ansatz for k ≥ 0:

k i=0 (B k-i (x)T i (x) -T i (x) Ãk-i (x)) = x p (∂ + σ(k -h) x )T k-h , where T k-h = 0 for k < h. (B.6)
It then follows from (B.6) that for k > 1 we have:

B 0 (x)T k (x) -T k (x) Ã0 (x) = k-1 i=1 B k-i (x)T i (x) -T i (x) Ãk-i (x) + x p (∂ + σ(k -h) x )T k-h (x) -Ãk (x) + B k (x).
By inserting the claimed block forms of T k (x) and Ãk (x) for all k > 1 in the above equation, with 1 ≤ ς ≤ 2, we obtain:

à k (x) = B k (x) B ςς 0 (x) T ς k (x) -T ς k (x) à 0 (x) = k-1 i=1 B ςς k-i (x) T ς i (x) -T ς k-i (x) à i (x) + x p (∂ + σ(k -h) x )T ς k-h (x) + B ς k (x).
Hence, for every k > 1, the computation of T k (x) requires only lower order terms of T, B, and Ã. Hence, by setting T k (x = 0) = 0 and expanding each of the equations for a fixed k > 1 w.r.t. x, we arrive at constant Sylvester matrix equations, which can be resolved successively due to the disjoint spectrum of A 11 0,0 and A 22 0,0 , to compute T k (x) up to any desired order in x.

Appendix C. Proofs for Section 3.2: ε-rank reduction

In this section, we establish our constructive proof of the necessity of Theorem 2. For sufficiency, one can consult [28, Section 5.5, p. 101]. We prove the following: Theorem 5. Consider the system (8) given by:

[A σ A ] x p ξ h A ∂F = A(x, ξ A )F = ∞ k=0 A k (x)ξ A k F
with h > 1 and m ε (A) > 1 and set rank (A 0 (x)) = r . Suppose that the polynomial

θ A (λ) := ξ A r det(λI + A 0 (x) ξ + A 1 (x))| ξ A =0
vanishes identically in λ. Then there exists a transformation F = R(x, ξ R )G such that the equivalent system

[ Ãσ à ] ξ à h x p∂G = Ã(x, ξ Ã) G = ∞ k=0 Ãk (x)ξ k à G
satisfies either h < h or h = h and rank( Ã0 (x)) < r. Moreover, such a R(x, ξ R ) can be constructed by [START_REF] Abbas | On the Reduction of Singularly-Perturbed Linear Differential Systems[END_REF]Algorithm 1] where it is always chosen to be a product of unimodular transformations in GL n (C[[x]]) and polynomial transformations in ξ A .

The proof of the theorem will be given after a set of intermediate results.

Lemma 7. Given the system

[A σ A ] x p ξ h A ∂F = A(x, ξ A )F = ∞ k=0 A k (x)ξ k A F,
with r = rank (A 0 (x)). There exists a unimodular transformation U(x) in GL n (C[[x]]) such that the leading coefficient matrix Ã0 (x) of the equivalent system (10) given by

[ Ãσ à ] ξ h à x p ∂G = Ã(x, ξ Ã)G, σ à ≥ min( p h + σ A , σ A )
has the form

Ã0 (x) = Ã11 0 (x) O Ã21 0 (x) O
where Ã11 0 (x) is a square matrix of dimension r and

Ã11 0 (x) Ã21 0 (x)
is a n × r matrix of full column rank r.

Remark 6. In practice, U(x) can be obtained by performing Gaussian elimination on the columns of A 0 (x) taking as pivots the elements of minimum valuation (order) in x.

Proof. By Remark 7, Ã0 (x) = U -1 (x)A 0 (x)U(x). Hence it suffices to search a similarity transformation U(x). Since C

[[x]] is a principal ideal domain (the ideals of C[[x]] are of the form x k C[[x]]), it is well known that one can construct unimodular transformations Q(x), U(x) lying in GL n (C[[x]]) such that the matrix Q(x)A 0 (x)U(x) has the Smith normal form Q(x) A 0 (x) U(x) = diag(x β 1 , . . . , x β r , 0, . . . , 0))
where det(U(0)) 0, det(Q(0)) 0 , and

β 1 , . . . , β r in Z with 0 ≤ β 1 ≤ β 2 ≤ • • • ≤ β r . It follows that we can compute a unimodular matrix U(x) in GL n (C[[x]]) so that its last n -r columns form a C[[x]]-basis of ker (A 0 (x)).
As for finding a lower bound of σ Ã, one can see that since U(x = 0) is invertible, we have:

val x ( Ãk ) ≥        val x ( Ãk ) if k h val x ( Ãh ) + p otherwise. Thus, if p < 0 then σ Ã ≥ σ A + p h . Remark 7. Consider system [A σ A ] given by (8). Let T (x) = (I n + Q(x)) ∈ GL n (C[[x]]) such that Q(x) = [q i j ] 1≤i, j≤n , s.t.        q n j ∈ C[[x]
] , for r + 1 ≤ j < n q i j = 0 elsewhere. One can then verify that the resulting equivalent system (10) is

[ Ãσ à ] ξ h à x p ∂G = Ã(x, ξ Ã)G, σ à ≥ σ A where              Ã0 (x) = (I n + Q(x)) -1 A 0 (x) (I n + Q(x)) Ã1 (x) = (I n + Q(x)) -1 A 1 (x) (I n + Q(x)) -γx p ∂Q(x); Ã2 (x) = . . . (C.1)
where γ = 1 if h = 1 and γ = 0 otherwise. But, due to (15), we can limit our interest to A 0 (x), A 1 (x), Ã0 (x), and Ã1 (x) solely. Hence, if h > 1 then it suffices to investigate T -1 A(x, ξ σ A ) T . We remark that if h = 1 then the term T -1 ∂T (x) should be taken into account. This is the reason the case h = 1 is treated separately. This distinction is purely technical as explained in the following remark and unlike for the Poincaré rank p of an unperturbed system, it does not lead to a classification of regularity.

By Lemma [START_REF] Barkatou | An Algorithm for Computing a Companion Block Diagonal Form for a system of Linear Differential Equations[END_REF], we can suppose without loss of generality that A 0 (x) is of the form

A 0 (x) = A 11 0 (x) O A 21 0 (x) O (C.2)
with r independent columns and (nr) zero columns. We partition A 1 (x) in accordance with

A 0 (x), i.e. A 1 (x) = A 11 1 (x) A 12 1 (x) A 21 1 (x) A 22 1 (x)
, and consider

G A (λ) = A 11 0 (x) A 12 1 (x) A 21 0 (x) A 22 1 (x) + λI n-r . (C.3)
This consideration of G A (λ) gives an ε-reduction criterion equivalent to θ A (λ). In fact, let D(ξ A ) = diag(ξ A I r , I n-r ) where r = rank (A 0 (x)). Then we can write ξ -1

A A(x, ξ A ) = ND -1 where N := N(x, ξ A ) ∈ M n (R). Set D 0 = D(0) and N 0 = N(x, 0). Then we have det(G A (λ)) = det(N 0 + λD 0 ) = det(N + λD)| ξ A =0 = (det( A(x, ξ A ) ξ A + λI n ) det(D))| ξ A =0 = (det( A 0 (x) ξ A + A 1 (x) + λI n ) ξ r A )| ξ A =0 = θ A (λ).
Thus, det(G A (λ)) ≡ 0 vanishes identically in λ if and only if θ A (λ) does. We illustrate our progress with the following simple example.

Example 11. [2, Example 2] Consider [A σ A ] ξ h A ∂F = A(x, ξ)F with σ A = 0, h > 1, and
A(x, ξ A ) =               ξ A -x 3 ξ A (1 + x)ξ A 0 x 2 xξ A 0 -2xξ A -x 0 0 2ξ A 0 2 0 ξ 2 A               . Clearly, A 0 (x) is nilpotent of rank 2 and G A (λ) =               0 0 x + 1 0 x 2 0 0 -2x -x 0 λ 2 0 2 0 λ               . (C.4)
We then have the following proposition:

Proposition 4. Given the system

[A σ A ] ξ h A x p ∂F = A(x, ξ A )F = ∞ k=0 A k (x)ξ k A F,
with r = rank (A 0 (x)). If m ε (A) > 1 and det(G A (λ)) ≡ 0 then there exists a finite product of triangular matrices T (x) = P(I n + Q(x)) where P is a permutation and det T (x) = ±1, such that for the equivalent system

[ Ãσ à ] ξ h à x p ∂G = Ã(x, ξ Ã) G = ∞ k=0 Ãk (x)ξ k à G
we have:

G Ã(λ) =           A 11 0 (x) U 1 (x) U 2 (x) V 1 (x) W 1 (x) + λI n-r- W 2 (x) M 1 (x) M 2 (x) W 3 (x) + λI           , (C.5)
where 0 ≤ ≤ nr, W 1 , W 3 are square matrices of order (nr -) and respectively, and

rank ( A 11 0 (x) U 1 (x) M 1 (x) M 2 (x) ) = rank ( A 11 0 (x) U 1 (x) ), (C.6) rank ( A 11 0 (x) U 1 (x) ) < r. (C.7)
Moreover, σ Ã ≥ min( p h + σ A , σ A ). We shall need the following Remark in the proof of the Proposition. Remark 8. [2, Remark 6] Suppose that G A (λ) has the form (C.5) and there exists a transformation T (x)

∈ GL n (C((x))) such that G T [A σ A ] (λ) has the form G T [A σ A ] (λ) =           A 11 0 U 1 U 2 V 1 W 1 + λI n-r- W 2 O O W3 + λI          
, where 0 ≤ ≤ nr and W3 is upper triangular with zero diagonal. Then,

det(G T [A σ A ] (λ)) = λ det( A 11 0 U 1 V 1 W 1 + λI n-r-
).

If det(G T [A σ

A ] (λ)) ≡ 0 then we have det(G T [A σ A ] (λ)) ≡ 0 as well (rank of leading coefficient matrix is unchanged). Hence,

det( A 11 0 U 1 V 1 W 1 + λI n-r- ) ≡ 0. (C.8)
For a fixed ∈ {0, . . . , n -r} we shall denote by G ( ) 0 the matrix

A 11 0 U 1 V 1 W 1 of (C.5).
Proof. (Proposition 4) Since det(G A (λ)) ≡ 0 then in particular, the matrix G A (λ = 0) is singular.

Let E 1 (respectively E 2 ) be the vector space spanned by the first r (resp. last nr) rows of G A (λ = 0). We have

dim(E 1 + E 2 ) = rank (G A (λ = 0)) < n.
If dim(E 1 ) < r then setting = 0 suffices to fulfill our claim. Otherwise, since

dim(E 1 + E 2 ) = dim(E 1 ) + dim(E 2 ) -dim(E 1 ∩ E 2 ) < n, it follows that either dim(E 2 ) < n-r or dim(E 1 ∩ E 2 ) > 0.
In both cases, there exists at least a row vector (1) (x) = ( (1) 1 (x), . . . , (1) n (x)) with entries in C((x)) in the left null space of G A (λ = 0), such that (1) j (x) 0 for some r + 1 ≤ j ≤ n. We can assume without loss of generality that (1) (x) has its entries in C [[x]]. Indeed, this assumption can be guaranteed by a construction as in Remark 6. Let the constant matrix P (1) denote the product of permutation matrices which exchange the rows of A, so that val x ( (1) n (x)) < val x ( (1) j (x)) , r + 1 ≤ j ≤ n -1, where val x denotes the x-adic valuation (order in x). Let (1) i j (x)] 1≤i, j≤n , s.t.

Q (1) (x) = [q
         q (1) n j (x) = - (1) j (x) (1)
n (x) , for r + 1 ≤ j < n q (1) i j = 0 elsewhere Thus, P (1)

(I n + Q (1) (x)) is unimodular in GL n (C[[x]]).
Set F = F 0 , A = A (0) , Ã = A ( ) and let A (1) be the matrix of the equivalent system ξ h A (1) x p ∂F (1) = A (1) (x, ξ A (1) )F (1) obtained by the transformation (1) . Thus, by Remark [START_REF] Barkatou | An Algorithm for Computing a Companion Block Diagonal Form for a system of Linear Differential Equations[END_REF], G A (1) (λ) has the form (C.5) with (C.6) and = 1. By Remark 8, the matrix G A (1) (λ = 0) is singular and the condition (C.7) does not occur, then one can find, by the same argument as above a permutation matrix and a nozero vector (2) (x) in the left null space of G (1)

F (0) = P (1) (I n + Q (1) (x)) F
(λ = 0). Let Q (2) (x) = [q (2) i j (x)] 1≤i, j≤n , s.t.          q (2) n-1, j (x) = - (2) j (x) (2)
n-1 (x) , for r + 1 ≤ j < n -1 q (2) i j = 0 elsewhere

The matrix G A (2) (λ) is then of the form (C.5) with (C.6) and = 2. Consider the finite sequence of equivalent systems obtained by the transformation

F (s-1) = P (s) (I n + Q (s) (x)) F (s)
where 1 ≤ s ≤ and

Q (s) (x) = [q (s) i j (x)] 1≤i, j≤n , s.t.          q (s) n-s+1, j (x) = - (s) j (x) (s)
n-s+1 (x) , for r + 1 ≤ j < ns + 1 q (s) i j = 0 elsewhere.

Then this process yields an equivalent matrix Ã(x, ξ Ã) := A ( ) (x, ξ A ( ) ) with (C.6) for which either (C.7) occurs or = nr. But in the latter case one has, again by Remark 8, that det(A 11 0 (x)) = 0, and so (C.7) occurs.

Example 12 (Continue Example 11).

A simple calculation shows that det(G A (λ)) ≡ 0 hence A is ε-reducible. From (C.4), for λ = 0, we have the singular matrix

G A (λ = 0) =               0 0 x + 1 0 x 2 0 0 -2x -x 0 0 2 0 2 0 0               . Let T =               1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0              
then the transformation F = TG yields the equivalent system ξξ h Ã∂G = Ã(x, ξ Ã)G where σ à = σ A and λ) has the form (C.5) with = 1 and r = 2. In fact,

Ã(x, ξ Ã) =                ξ Ã -x 3 ξ Ã 0 (1 + x)ξ Ã x 2 xξ Ã -2xξ Ã 0 0 2 ξ 2 Ã 0 -x 0 2ξ Ã 0                . G Ã(
G Ã(λ) =               0 0 0 (1 + x) x 2 0 0 0 0 2 λ 0 -x 0 2 λ               . Lemma 8. Given the system [A σ A ] ξ h A x p ∂F = A(x, ξ A )F = ∞ k=0 A k (x)ξ k A F.
Set r = rank(A 0 (x)) and suppose that m ε (A) > 1 and G A (λ) has the form (C.5) with conditions (C.6) and (C.7) satisfied. Consider the shearing transformation S (ξ

A ) = diag(ξ A I r , I n-r-, ξ A I ) if 0 and S (ξ A ) = diag(ξ A I r , I n-r ) otherwise. Then F = S (ξ A )G yields the equivalent system [ Ãσ à ] ξ h à x p ∂G = Ã(x, ξ Ã)G = ∞ k=0 Ãk (x)ξ k à G, for which rank( Ã0 (x)) < r and σ à ≥        σ A if p > 1 σ A + p-1 h otherwise.
Proof. We partition A(x, ξ A ) as follows ( we drop (x, ξ A ) for clarity)

A(x, ξ A ) =           A 11 A 12 A 13 A 21 A 22 A 23 A 31 A 32 A 33          
where A 11 , A 22 , A 33 are of dimensions r, nr -, and respectively. It is easy to verify then that

Ã(x, ξ A ) = S -1 AS -S -1 ∂S = S -1 AS -ξ h A x p-1 diag(σI r , O n-r-, σI ) =           A 11 -ξ h A x p-1 σI r ξ -1 A A 12 A 13 ξ A A 21 A 22 ξ A A 23 A 31 ξ -1 A A 32 A 33 -ξ h A x p-1 σ A I           .
Hence, the new leading coefficient matrix is

Ã0 (x) =           A 11 0 U 1 O O O O M 1 M 2 O           and rank( Ã0 (x)) = rank(A 11 0 U 1 ) < r. Moreover, due to the entries of the form ξ h A x p-1 σ A in Ã, σ A should be adjusted as claimed if p < 1. Example 13 (Continue Example 12). Let S (ξ Ã) = diag(ξ Ã, ξ Ã, 1, ξ Ã) then G = S (ξ Ã)U yields ξ h Ã∂U = Ã(x, ξ)U where σ à = σ à and Ã(x, ξ Ã) =                ξ à -x 3 ξ à 0 (1 + x)ε x 2 xξ à -2x 0 0 2ξ à ξ 2 à 0 -x 0 2 0               
.

It is clear that the leading term Ã0 (x) has rank 1 < 2 = rank(A 0 (x)).

Proof. (Theorem 5) By Lemma 7, we can assume that A 0 (x) is in the form (C.2). Then, G A (λ) is constructed as in (C.3). Since det(G A (λ)) ≡ 0, it suffices to take the change of basis F = RG = T S G, where T and S are as in Propositions 4 and Lemma 8 respectively.

Remark 9. The ε-reducibility of [A σ A ] implies that the rank of the leading coefficient matrix can be reduced (and consequently the ε-Moser rank) without necessarily reducing the ε-rank h of the system. If the ε-reduction criterion holds for a sufficient number of equivalent systems then a repetitive application of such a transformation results in an equivalent system whose leading coefficient matrix has a zero rank, hence h can be reduced at least by one (e.g. Example 8).

Appendix D. Proofs for Section 4

We give the proof of Theorem 3 after establishing a series of useful lemmas. The following proofs are an adaptation to the parametrized setting of the proofs in [8, Lemma 3, Lemma 4, Proposition 1, Theorem 1] respectively. For clarity within these intermediate proofs, we will express systems in the equivalent notation [A σ A ] ∂F = A(x, ε)F where A(x, ε) ∈ M n (K), rather than

[A σ A ] ξ h A x p ∂F = Ã(x, ξ A )F. Let det(λI -A(x, ε)) = λ n + α n-1 (x, ε)λ n-1 + • • • + α 0 (x, ε). (D.1)
such that α n = 1 and α i (x, ε) = ∞ j=val ε (α i ) α i, j (x)ε j ∈ K for i ∈ {0, . . . , n}. We define the εpolygon N ε (A) of [A σ A ] as in Section 2.1.2, by taking P ε (A) to be the union of P(i, val ε (α i (x, ε)) for i ∈ {1, . . . , n}. We thus prove the following theorem, of which Theorem 3 is a straightforward corollary. Theorem 6. Consider the ε-irreducible system [A σ A ] ∂F = A(x, ε)F where A(x, ε) ∈ M n (K), h > 0, and (D.1). If h > nrank(A 0 (x)) then the ε-formal exponential order is given by

ω ε (A) = max 0≤i<n ( -val ε (α i ) n -i ).
Additionally, the corresponding ε-polynomial is given by the algebraic equation

E ε (X) = k=0 α i k ,val ε (α i k ) X (i k -i 0 )
where 0 ≤ i 0 < i 1 < • • • < i = n denote the integers i for which ω ε (ni) = -val ε (α i ) (i.e. lie on the edge of slope ω ε of the ε-polygon

N ε (A) of [A σ A ]); and α i,val ε (α i ) (x) = ε -val ε (α i ) α i (x, ε)| ε=0 .
Lemma 9. Let A(x, ε), W(x, ε) be matrices in M n (K) and M n (R) respectively. Put h = max(0, -val ε (A(x, ε))) and let

det(λI -A(x, ε)) -det(λI -A(x, ε) + W(x, ε)) = α n-1 (x, ε)λ n-1 + α n-2 (x, ε)λ n-2 + • • • + α 0 (x, ε),
where α n-1 (x, ε), α n-2 (x, ε), . . . , α 0 (x, ε) lie in K. Then

val ε (α n-i ) ≥ (1 -i)h 1 ≤ i ≤ n.
Proof. For any i ∈ {1, . . . , n}, it follows from Cramer's rule that

α n-i = γ l=0 (w l i-1 s=1 a l,s )
where for 0 ≤ l ≤ γ, 1 ≤ s ≤ i -1, w l are entries in W(x, ε) and a l,s are entries in W(x, ε) or A(x, ε). Consequently,

val ε (α n-i ) ≥ (i -1)val ε (A(x, ε)) ≥ h. Lemma 10. Let A(x, ε), B(x, ε) be two matrices in M n (K) such that val ε (A(x, ε)) ≤ val ε (B(x, ε)). Consider the two systems [A σ A ] ∂F = A(x, ε)F and [B σ B ] ∂G = B(x, ε)G. Suppose that there exists T ∈ GL n (K) such that [B σ B ] = T [A σ A ]. Put h = max(0, -val ε (A(x, ε))) and let        det(λI -A(x, ε)) = λ n + α n-1 λ n-1 + α n-2 λ n-2 + • • • + α 0 det(λI -B(x, ε)) = λ n + β n-1 λ n-1 + β n-2 λ n-2 + • • • + β 0 .
Then, we have:

val ε (α n-i (x, ε) -β n-i (x, ε)) ≥ (1 -i)h, 1 ≤ i ≤ n.
Proof. By [31, Lemma 1] we can write T (x, ε) = P(x, ε) ε γ Q(x, ε) where P, Q ∈ GL n (R) and

ε γ = diag(ε γ 1 , . . . , ε γ n ) for some integers (γ 1 ≤ γ 2 ≤ • • • ≤ γ n ). Consider [ Ãσ Ã ] := P[A σ A ] and [ Bσ B ] := Q -1 [B σ B ]. Then we have B = ε -γ Ãε γ , val ε ( Ã) = val ε (A), and val ε ( B) = val ε (B).
It follows that,

             det(λI -Ã) = det(λI -P -1 AP + P -1 ∂P) = det(λI -A + (∂P)P -1 ) det(λI -B) = det(λI -QBQ -1 + Q∂Q -1 ) = det(λI -B + (∂Q -1 )Q) det(λI -B) = det(λI -ε -γ Ãε γ ) = det(λI -Ã).
Since P, Q, P -1 , Q -1 , are units of M n (R), it follows that (∂P)P -1 and (∂Q -1 )Q inherit this property as well. The rest of the proof follows as a consequence of Lemma 9.

Proposition 5. Consider the system [A σ A ] ∂F = A(x, ε)F where A ∈ M n (K) and let:

det(λI -A(x, ε)) = λ n + α n-1 (x, ε)λ n-1 + α n-2 (x, ε)λ n-2 + • • • + α 0 (x, ε). Let [C σ C ] ∂G = C(x, ε)G where C(x, ε) = Companion (c i (x, ε)) 0≤i≤n-1 =                       0 1 0 . . . 0 0 0 1 . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . 1 c 0 c 1 c 2 . . . c n-1                       , (D.
2) be a companion system which is equivalent to [A σ A ] over GL n (K). Then we have,

val ε (α i -c i ) ≥ (1 -(n -i))h 0 ≤ i ≤ n -1. Proof. Let [c σ c ] with c n = 1 be the scalar equation representing [C σ C ] (σ c = σ C ). Consider κ of [c σ c
] as defined in Section 3.3. We define

β i = (n -i)κ, for i ∈ {0, . . . , n -1} and [D σ D ] = ε β [C σ C ]. Then D(x, ε) = ε -β C(x, ε)ε β where ε β = diag(ε β 0 , . . . , ε β n-1 ). It follows that D(x, ε) = ε -κ Companion (ε β i c i (x, ε)) 0≤i≤n-1 . We have, val ε (D(x, ε)) ≥ -κ since val ε (ε β i c i ) = (n -i)κ + val(c i ) ≥ 0. By the equivalence between [A σ A ] and [C σ C ] (resp. [D σ D ]) we have h ≥ κ (resp. val ε (D(x, ε)) ≥ -κ ≥ -h). Let det(λI -D(x, ε)) = λ n + d n-1 (x, ε)λ n-1 + • • • + d 0 (x, ε).
Hence, by Lemma 10, we have

val ε (α i -d i ) ≥ (1 -(n -i))h, 0 ≤ i ≤ n.
Moreover, by Lemma 9

val ε (d i -c i ) ≥ (1 -(n -i)) max (0, -val ε (D(x, ε))) ≥ (1 -(n -i))κ ≥ (1 -(n -i))h.
It follows that

val ε (α i -c i ) ≥ min (val ε (α i -d i ), val ε (d i -c i )) ≥ (1 -(n -i))h. Proof. ( Theorem
ω ε := ω ε (A) = ω ε (C) = max 0≤i<n ( -val ε (c i ) n -i ) and E ε (X) = k=0 c i k ,val ε (c i k ) X (i k -i 0 )
where 0

≤ i 0 < i 1 < • • • < i = n denote the integer i for which ω ε (n -i) = -val ε (c i )
. By Corollary 2 and Proposition 5 one has

val ε (α i -c i ) ≥ (i -n + 1)h = (i -n)ω ε + (i -n)(h -ω ε ) + h ≥ ω ε (i -n) + (-n)(1 - r n ) + h ≥ ω ε (i -n) + r + h -n. It follows from h + r -n > 0 that val ε (α i -c i ) > (i -n)ω ε , 0 ≤ i ≤ n. • If i ∈ {i 0 , . . . i } then ω ε (i -n) = val ε (c i ) and val ε (α i -c i ) > val ε (c i ). Hence, val ε (α i ) = val ε (c i ) and α i, val ε (α i ) = c i, val (c i ) . • Else val ε (α i -c i ) > ω ε (i -n) and val ε (c i ) > ω ε (i -n). Thus, val ε (α i ) ≥ min (val ε (α i -c i ), val ε (c i )) > ω ε (i -n).
This completes the proof.

We also illustrate Lemma 4 with the following example:

Example 14 (Lemma 4). Consider the system

[A σ A ] ξ 2 A ∂F = A(x, ξ A )F with σ A = 0 and A(x, ξ A ) =                                               0 0 ξ A 2 0 ξ A x 0 ξ A 2 ξ A 3 xξ A 0 0 0 0 0 0 0 x 0 0 4 3 ξ A x 2 0 0 0 ξ A 2 xξ A 0 0 xξ A 0 0 0 ξ A 2 0 0 0 xξ A ξ A 2 ξ A 2 0 0 0 ξ A 3 0 0 0 0 xξ A                                               . A(x, ξ A ) is ε-irreducible,
(x, ξC ) =                                                0 0 0 x ξC ξ4 C ξ4 C 0 1/4 ξ2 C x 2 x ξ2 C ξ7 C 0 1/4 ξC ξ3 C -x 2 1/4 ξ4 C x 3/4 ξC x 3 0 0 ξ8 C x ξ2 C 0 0 ξ5 C 0 4 0 0 0 0 0 x ξ3 C 0 0 0 0 ξ5 C 3 ξ2 C x 2 ξ6 C 0 x ξ2 C 0 0 0 0 x 0 0 ξ5 C 0 ξ2 C 0                                               
, for which, h = 5 and r = 2. Thus, the condition h > nr is still not necessarily verified after a random ramification. To guarantee that we will arrive at a system verifying this condition, we make use of Lemma 4 and choose an integer d such that d ≥ 7/(1

+ 2/7) = 49/9. Let d = 7 > 49/9, perform ξ A = ξ7 A in [A σ A ],
and then ε-rank reduction. This yields the system

[B σ B ] ξ11 B ∂G = B(x, ξB )G with σ B = 0 and B(x, ξB ) =                                                ξ4 B x 1/4 ξ4 B x 2 ξ15 B 0 1/4 ξB ξ7 B -x 2 1/4 ξ8 B x 3/4 ξB x 3 0 0 0 ξB x ξ8 B ξ8 B 0 0 0 ξ18 B ξ4 B x 0 0 ξ11 B 4 0 0 0 0 0 0 0 ξ7 B x 0 0 0 ξ11 B 3 ξ4 B x 2 0 ξ14 B ξ4 B x 0 0 0 0 0 x 0 ξ11 B 0 ξ4 B 0                                               
.

For system [B σ B ], we have: h = 11, r = 2, and so the condition h = 11 > nr = 7 -2 holds, and so, Theorem 3 can now be applied. From det(λI -B(x, ξB ) ξ11 B

), we can compute ω varepsilon (M) = 42 4 . And so, we introduce the ramification ξB = ξ4

B and apply ε-rank reduction which gives a εsystem of order 42 and whose leading coefficient has 4 nonzero eigenvalues. In fact, we compute

[ Mξ M ] ξ42 M ∂W = M(x, ξM )W with σ M = 0 and M(x, ξM ) =                                                        ξM 14
x 1/4 ξM

                                                      
.

We remark that one can observe that this procedure was not necessary for such system since gcd(14, 42, 28) = 14. Thus a simplification would leave us with a system of order 3. Thus, in the implementation we first try to ignore the condition h > nr and apply Theorem 3 directly to the system [A σ A ]. For this particular example we obtain ω ε (A) = 3/2. Then by performing ξ A = ξ2

A and the ε-rank reduction, we arrive at a ε-irreducible system of order 3 whose leading coefficient has 4 nonzero eigenvalues. In fact, we compute [ Ãσ à ] ξ3 Ã∂G = Ã(x, ξ)G where σ à = 0 and

Ã(x, ξ Ã) =                                                 0 0 0 x ξ2 Ã ξ2 Ã 0 1/4 ξ Ã x 2 ξ Ã x ξ4 Ã 0 -1/4 x 2 + 1/4 ξ2 Ã 1/4 x ξ2 Ã 3/4 x 3 0 0 ξ5 Ã ξ Ã x 0 0 ξ3 Ã 0 4 0 0 0 0 0 x ξ2 Ã 0 0 0 0 ξ3 Ã 3 ξ Ã x 2 ξ4 Ã 0 ξ Ã x 0 0 0 0 x 0 0 ξ3 Ã 0 ξ Ã 0                                                 .
If Theorem 3 does not lead to the desired result, we then resort Lemma 4. In the case of unperturbed system, it is also a matter of discussion whether this condition is necessary. M. Miyake recently claimed to have such an example (where this condition is necessary).

B(x, ξ B) = 0 1 -1 0 + -x 3/2 0 0 0 ξ B + 0 0 x 3 0 ξ2 B + x 9/2 0 0 -x 9/2 + x 17/2 ξ3 B + O( ξ4 B).
The leading matrix has two distinct eigenvalues. Applying Splitting Lemma we get

[ Bσ B ] ξ3 B x 19/2 ∂R = B(x, ξ B)R where σ B = -3 and B(x, ξ B) = -i 0 0 i + -1 2 x 3/2 0 0 -1 2 x 3/2 ξ B + O( ξ2 B).
Thus this system can be decoupled to any desired precision. The solution follows by straightforward integration. Remark that, as expected, the system has the same ε-formal exponential order and ε-polynomials, as its equivalent scalar equation given in Example 2. The leading term of the ε-exponential part Q is then given by: exp(

       -2i √ x ε 3/2 0 0 2i √ x ε 3/2        ).
Moreover, since σ B = -3, we set ρ 1 = 1/3. Further computations yield the full exponential part of outer solutions of the input system [A σ A ] of Example 9:

x 2 2ε 2 - 1 xε and - 1 20t 5 ε 1/2 - 1 2t 2 ε + 2t ε 3/2 where x = -t 2 .
Moreover, since σ B = -3, we set ρ 1 = 1/3. For the initial system of Example 9 given by:

ε 2 ∂F =           0 1 0 0 0 1 ε 0 x           F,
we have obtained so far outer formal solutions. To compute inner (and probably intermediate) solutions, we set τ = xε -ρ A = xε -1/3 which yields:

ε 5 3 ∂ τ F =           0 1 0 0 0 1 ε 0 τε 1 3           F,
which can be expressed equivalently as:

[E σ E ] ξ5 E ∂ τ F = E(τ, ξE )F =           0 1 0 0 0 1 ξ3 E 0 τ ξE           F, with σ E = 0 and ξ E = ε = ε3 = ξ3 E .
By Algorithm 1, one can compute the transformation F = TG where

T =           0 0 1 0 1 0 1 0 0                     ξ2 E 0 0 0 ξE 0 0 0 1           . The resulting system is [ Ẽσ Ẽ ] ξ4 Ẽ ∂ τ G = Ẽ(x, ξ Ẽ )G =           τ 0 1 1 0 0 0 1 0           G, with σ Ẽ = 0.
One can verify that Ẽ0,0 has three distinct roots and so the system can be decoupled into three scalar equations. The leading term of the ε-exponential part in the inner subdomain is then given by: exp(

               τ dz+... ε 4/3 0 0 0 τ -1+i/ √ 3+...dz 2ε 4/3 0 0 0 τ -1-i/ √ 3+...dz 2ε 4/3                ).
And further computations show that the diagonal entries of the ε-exponential part of a fundamental matrix of formal solutions in the inner subdomain are given by:

τ + (1/6)τ 2 + (1/27)τ 3 + O(τ 4 ) ε 4/3 and τ RootOf(z 2 + z + 1) + (1/6)τ 2 -(1/27)(RootOf(z 2 + z + 1) + 1)τ 3 + O(τ 4 ) ε 4/3 .
Since σ Ẽ = 0, the reduction stops here.

Example 16 (Iwano-Sibuya polygon). Consider the following scalar equation whose σ a = 0:

[a σ a ] ∂ 4 f -ε -4 (2x + ε 3 )∂ 2 f + ε -8 x 3 f = 0, (E.1)
One can verify that ω ε = 2. The Iwano-Sibuya polygon is given by the following set of points: P ω ε (a) = (2, -1), P 00 = (0, 3 4 ), P 20 = (0, 1 2 ), and P 23 = ( 3 2 , 0) leading to only one slope given by σ = -3 4 . Then the behavior can be investigated with the help of [ρ] 0 < ρ 1 = -1/σ = 4/3. Hence, there is only one stretching transformation to consider: τ = xε -4/3 . We now wish to compute formal solutions and the sequence [ρ] with the matricial representation using our algorithm:

Let F = [ f, ε 2 ∂ f, ε 4 ∂ 2 f, ε 6 ∂ 3 f ] T then (E.1
) can be expressed as the following first order differential system (σ A = 0)

[A σ A ] ξ 2 A ∂F = A(x, ξ A )F =               0 1 0 0 0 0 1 0 0 0 0 1 -x 3 0 ξ 3 A + 2x 0               F, where ξ A = ε. (E.2)
• For the turning point treatment, we set x = t 2 and then apply the transformation

T =                       0 0 0 1 0 0 t 2 0 t 3 0 0 0 0 t 4 0 0                       which yields the system [ Ãσ Ã ] ξ 2 Ãt∂ t G = Ã(t, ξ Ã)G =                        -3 ξ 2 Ã 2 0 0 4 + 2 ξ 3 Ãt 5/2 -4 ξ 2 Ã 0 -2 t 2 0 -2 ξ 2 Ã 0 0 0 2 t 0                        G. (E.3)
where σ Ã = -3/2, i.e. ξ Ã = t -3/2 ε.

Applying the Splitting lemma decouples (E.3) into two subsystems (truncated at order 4 in both t and ε)

[S 1 σ S 1 ] ξ 2 S 1 t∂ t U 1 = S 1 (t, ξ S 1 )U 1 and [S 2 σ S 2 ] ξ 2 S 2 t∂ t U 2 = S 2 (t, ξ S 2 )U 2 where σ S 1 = σ S 2 = -3/2, ξ S 1 = ξ S 2 = t -3/2 ε, and
S 1 (t, ξ S 1 ) =          -1/2 ξ 2 S 1 ξ 4 S 1 t 2 + t 2 + 6 -ξ 4 S 1 t 2 + 2 2 ξ 3 S 1 t 5/2 -ξ 4 S 1 t 2 -t 2 + 4 -1/2 ξ 2 S 1 3 ξ 4 S 1 t 2 + t 2 + 8         
whose leading coefficient has two constant nonzero eigenvalues (hence the dimension of the system /order of the equivalent equation can be reduced). And,

S 2 (t, ξ S 2 ) = 1/2 ξ 2 S 2 (3 ξ 4 S 2 t 2 + t 2 -4) γ 1 1/2 t(2 ξ 4 S 2 t 2 + t 2 + 4) γ 2 , (E.4) where        γ 1 = -1/8 t(4 ξ 3 S 2 t 5/2 -9 ξ 4 S 2 t 2 -8 ξ 4 S 2 -8) γ 2 = -1/16 t 2 ξ 2 S 2 (8 ξ 3 S 2 t 5/2 -105 ξ 4 S 2 t 2 -40 ξ 4 S 2 -14 t 2 -16. Since εt -3/2 = εx -3/4 , we have ρ 1 = 4/3.
If we wish to continue the reduction for the second subsystem, then, as usual, we first consider its leading coefficient:

S 2 0 (t) =        0 t 2t + 1/2 t 3 0        .
Hence, by the turning point resolution (here it suffices to factorize t), the system [S 2 σ S 2 ] can be rewritten as:

[ S 2 σ S 2 ] ξ 2 S 2 t 2 ∂ t R = S 2 (t, ε)R where σ S 2 = -2, ξ S 2 = t -2 ε, and S 2 (t, ξ S 2 ) =          1/2 ξ 2 S 2 3 t 4 ξ 4 S 2 + t 2 -4 γ1 t 4 ξ 4 S 2 + 1/2 t 2 + 2 γ2          where          γ1 = -1/2 ξ 3 S 2 t 4 + 9 t 4 ξ 4 S 2 8 + ξ 4 S 2 t 2 + 1 γ2 = -1/16 t 2 ξ 2 S 2 -105 t 4 ξ 4 S 2 + 8 ξ 3 S 2 t 4 -40 ξ 4 S 2 t 2 -14 t 2 -16 .
The leading constant coefficient has two constant eigenvalues and hence can be decoupled.

• We consider again system [A σ A ] which is given by (E.2), and we apply the stretching x = ε 4/3 τ which yields the system [E σ E ]. One can verify that σ E = 0 and so with ξ E = ε = ξ3 E = ε3 , we have:

[E σ E ] ξ2 E ∂ τ H =               0 1 0 0 0 0 1 0 0 0 0 1 -2τ 3 ξ12 E 0 ξ9 E + 2τ ξ4 E 0               H. (E.5)
Then the ε-rank reduction with

T =                       0 0 0 1 0 0 ξ2 E 0 0 ξ4 E 0 0 ξ6 E 0 0 0                       yields [ Ẽσ Ẽ ] ∂ τ V = E(τ, ξ Ẽ )V =                        0 ξ15 Ẽ + 2 τ 0 -2 τ 3 ξ12 Ẽ 1 0 0 0 0 1 0 0 0 0 1 0                        V, σ Ẽ = 0, ξ Ẽ = ε.
Evidently, The ε-exponential part of [ Ẽσ Ẽ ] is zero. A formal fundamental matrix of solutions can be constructed following Subsection 2.1.1.

Example 17 (Roo's equation [START_REF] Nakano | On the complex WKB analysis for a Schrdinger equation with a general three-segment characteristic polygon[END_REF]). Consider the singularly-perturbed scalar differential equation with σ a = 0 (ξ a = ε) :

[a σ a ] ξ 4 a ∂ 2 f = (x 5 + ξ a x 2 + ξ 2 a ) f. (E.6)
Iwano-Sibuya's polygon of [a σ a ] consists of three segments connecting four points [33, p. 607], P 00 = (0, 5/2), P 01 = (1/2, 1), P 02 = (1, 0), and P o ε = (2, -1), and thus having the three slopes, -3, -2, and -1. Consequently, Iwano-Sibuya's sequence [ρ] of positive rational numbers is:

[ρ] 0 < ρ 1 = 1/3 < ρ 2 = 1/2 < ρ 3 = 1.
We now wish to recompute the sequence [ρ] with the matrix representation of [a σ a ] using our algorithm. Let F = [ f, ξ 2 a ∂ f ] T then [a σ a ] can be rewritten as the following first order differential system (σ A = 0)

[A σ A ] ξ 2 A ∂F = A(x, ξ A )F = 0 1 x 5 + ξ A x 2 + ξ 2 A 0 F. (E.7)
• For the resolution of the turning point, let T = 1 0 0 x 5/2 (resolution of turning point) then F = TG yields

[ Ãσ Ã ] ξ 2 Ã x 8 ∂G = Ã(x, ξ Ã) = 0 1 1 + ξ Ã + ξ Ã x -5 2 ξ 2 Ã x 5/2 G (E.8)
with σ Ã = -3 and ξ Ã = x -3 ε. The system [ Ãσ Ã ] can be decoupled into two scalar equations using the Splitting lemma and we have ρ 1 = 1/3.

• We consider again the system [A σ A ] and we perform the stretching τ = xε -1/3 and ramification ξ E = ε = ε3 = ξ3 E . This yields: The leading coefficient matrix is nilpotent. The computation of the ε-exponential order suggests introducing an additional ramification in ε of order 2. Since we already have a ramification of order 3, to avoid a proliferation of notations, we reset ε = ε6 . This yields:

[E σ E ] ξ5 E ∂ τ U = E(τ, ξE )U = 0 1 τ 5 ξ5 E + ξ5 E τ 2 + ξ6
[ Ẽσ Ẽ ] ξ6 Ẽ ∂ τ G = Ẽ(τ, ξ Ẽ )G =        0 ξ2 Ẽ τ 5 + τ 2 + ξ2 Ẽ 0        G, with σ Ẽ = 0, ε = ξ E = ξ6 Ẽ = ε6 . Let T 2 =
ξ Ẽ 0 0 1 then G = T 2 W yields the following ε-irreducible system:

[J σ J ] ξ5 J ∂ τ W = 0 1 τ 5 + τ 2 + ξ2 J 0
W, with σ J = 0, and ξ6 J = ε6 = ε.

Let T 3 = 1 0 0 τ (resolution of turning point) then W = T 3 V yields

[ Jσ J ] ξ5 J τ 4 ∂ τ V = B(x, ξ J )V = 0 1 τ 3 + 1 + ξ2 J -ξ5 J τ 3 V
with σ J = -1, ξ J = τ -1 ε = τ -1 ε 1/6 . Clearly, system [ Jσ J ] can be decoupled into two scalar equations. To find the restraining index, we rewrite τ and ε in terms of x and ε: τ -1 ε = τ -1 ε 1/6 = (x -1 ε 1/3 )ε 1/6 = x -1 ε 1/2 = (x -2 ε) 1/2 , and so ρ 2 = 1/2.

• We consider again system [A σ A ] and perform the stretching τ = xε -1/2 , and the ramification ε = ε2 . This yields:

[C σ C ] ξ3 C ∂ τ Z = C(τ, ξC )Z = 0 1 τ 5 ξ5 C + τ 2 ξ4 C + ε 2 0 Z, with σ C = 0, ξ C = ε = ε2 = ξ2 C . Let T 1 = 1 0 0 ξ2 C ( ε-rank reduction) then Z = T 1 G yields [ Cσ C ] ε 1/2 ∂ τ G = 0 1 τ 5 ε 1/2 + τ 2 + 1 0
G, with σ C = 0 and ξ C = ξC , which can be decoupled into two scalar equations by Splitting lemma. Since σ C = 0, the restraining index is infinity.

And so, by our techniques, we have obtained the two slopes ρ 1 = 1/3 and ρ 2 = 1/2.

Example 18 (( [START_REF] Bender | Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory[END_REF], Example 1, Section 9.5, p. 446)). Consider the following singularlyperturbed scalar differential equation

ε∂ 3 f -∂ f + x f = 0
which can be rewritten as the following differential first order system

[A σ A ] ξ A ∂F = A(x, ξ A )F =           0 ξ A 0 0 0 ξ A -x 1 0           F where F = [ f, ∂ f, ∂ 2 f
] T , σ A = 0, and ξ A = ε. Since A 0 (x) is nilpotent, Splitting lemma cannot be applied but no treatment of turning points is required. The transformation

F =               0 1 0 1 0 0 0 0 1                             1 x 0 0 1 0 0 0 1              
G yields the following ε-irreducible equivalent system:

[ Ãσ à ] ξ Ã∂G = Ã(x, ξ Ã)G =               -ξ à x -x 2 ξ à -ξ à ξ à ξ à ξ à x 0 1 0 0               G,
with σ Ã = 0 and θ Ã(λ) = -λx. Since the leading matrix coefficient is still nilpotent, we have to compute the ε-exponential oder. The characteristic polynomial of

Ã(x,ξ Ã ) ξ Ã
is given by:

λ 3 + ξ Ã -1 ξ Ã λ + x ξ Ã .
Consequently, the ε-exponential oder is given by: ω ε ( Ã) = max { 1 2 , 1 3 } = 1 2 . Let ε = ε2 and ξ = ξ2 then we have:

[ Ãσ Ã ] ξ2 ∂G = Ã(x, ξ)G =               -ξ2 x -x 2 ξ2 -ξ2 ξ2 ξ2 ξ2 x 0 1 0 0               G, with σ Ã = 0. The transformation G =           ξ 0 0 0 ξ 0 0 0 1                         0 1/2 1/2 1 0 0 0 -1/2 -1/2              
U yields the following ε-irreducible system:

[B σ B ] ξB ∂U = B(x, ξB )U =                 ξB x ξB /2 ξB /2 -x 2 -1 ξB -1/2 ξB x -1 -1/2 ξB x -x 2 -1 ξB -1/2 ξB x -1/2 ξB x + 1                 U,
with σ B = 0, ξB = ε. The leading matrix coefficient B 0,0 has three distinct roots and consequently, the system can be decoupled into three first order scalar differential equations by applying the Splitting lemma. A fundamental matrix of formal solutions is then given by: Φ(x, ε 1/2 ) exp(

          0 0 0 0 -x ε 1/2 0 0 0 x ε 1/2           ).
where Φ(x, ε 1/2 ) is the product of all the transformations performed including that of the Splitting lemma. Since σ B = 0, the origin is not a turning point for this system and we do not need stretchings.

Example 19 ([? ]

). We consider the scalar differential equation

[a σ a ] ε 3 ∂ 2 f -(x 3 + ε) f = 0,
with σ a = 0. Setting F = [ f, ∂ f ] T , we get the following differential first order system

[A σ A ] ξ 3 A ∂F = A(x, ξ A )F = 0 ξ 3 A x 3 + ξ A 0 F
with σ A = 0 and ξ A = ε.

• Since A 0,0 is nilpotent and A 0 (x) is not, we first start with the treatment of the turning point at x = 0. It suffices however to factorize x 3 from A 0 (x) which results in σ A = -3 and ξ A = x -3 ε. We then apply the transformation F = ξ A 0 0 1 G to get the ε-irreducible system:

[ Ãσ Ã ] ξ 2 Ã x 6 ∂G = Ã(x, ξ Ã)G = 3x 5 ξ 2 Ã ξ Ã x 6 1 + ξ Ã 0 G
with σ Ã = -3 and ξ Ã = x -3 ε. Since Ã0 (x) is still nilpotent, we proceed to computing the ε-exponential order from the characteristic polynomial of Ã(x,ε) ξ 2 Ã x 6 which is given by:

λ 2 - 3 x λ - 1 + ξ Ã ξ 3 Ã x 6 .
Hence, ω ε ( Ã) = 3 2 and the ε-polynomial is given by: X 2 -x 3 = 0. Let ε = ε2 and ξ Ã = x -3 ε, then with ξ Ã = ξ2 Ã x 3 we have:

[ Ãσ à ] x 12 ξ4 Ã∂G = Ã(x, ξ Ã)G = 3x 11 ξ4 à ξ2 à x 9 1 + x 3 ξ2 à 0 G
with σ Ã = -3. We then perform G = ξ Ã 0 0 1 U to get the ε-irreducible system: Since B00 has two distinct eigenvalues, ±1, one can proceed to the Splitting lemma which decouples the system thoroughly. From the the eigenvalues of B00 or the ε-polynomial, one can read the ε-exponential part of an outer solution of this system:

exp(        
x z 3/2 dz ε 3/2 0 0 -

x z 3/2 dz ε 3/2         ),
and deduce that ρ 1 = 1/3.

• We consider again system [A σ A ] and perform the stretching τ = xε -1/3 . After introducing the ramification ε = ε3 , we have: V, with σ S = 0, ξ S = ε. Due to the nature of the eigenvalues of S 0,0 , the system [ S σ S ] can be decoupled into two scalar equations. σ S = 0 suggests that there are no more stretchings to perform.

• However, if we consider again system [A σ A ] and experiment with the stretching τ = xε -1 , we get: So far, this matricial approach can determine the stretching τ = xε -1/3 but not the stretching τ = xε -1 . However, both stretchings can be obtained from Iwano-Sibuya's polygon by the treatment of this system as a scalar equation [? , Last section]. 56

Example 20 ([16], Example 1, Section 9.5, p. 453). We consider the scalar equation

ε 3 x∂ 2 f + x 2 ∂ f -(x 3 + ε) f = 0,
which can be rewritten as the following first order differential system

[A σ A ] xξ 3 A ∂F = A(x, ξ A )F = 0 ξ 3 A x x 3 + ξ A -x 2 F
where σ A = 0, ξ A = ε, and F = [ f, ∂ f ] T . We first start with the treatment of the turning point at x = 0 since A 0,0 is nilpotent while A 0 (x) is not. It suffices here to factorize x 2 which yields

[ Ãσ à ] x 5 ξ 3 Ã∂F = Ã(x, ξ Ã)F = 0 ξ 3 à x 5 x + ξ à -1 F
where σ Ã = -2 and ξ Ã = x -2 ε. Now that Ã0,0 has two distinct eigenvalues, we can proceed to the Splitting lemma. The transformation F = 0 1 1 x G block-diagonalizes Ã0 (x) which yields

[ Ãσ à ] x 5 ξ 3 Ã∂G = Ã(x, ξ Ã)G =        -ξ 3 à x 6 -1 -ξ 3 à x 7 -ξ 3 à x 5 + ξ à ξ 3 à x 5 ξ 3 à x 6        G,
where σ Ã = -2 and ξ Ã = x -2 ε and from which we can deduce that the ε-exponential part is given by: exp( -x 2 2ε 3 + O(ln(x)ε -2 ) 0 0 O(ln(x)ε -2 ) ).

Since the restraining index is nonzero, we can deduce that ρ 1 = 1/2 and apply the stretching τ = xε -1/2 in the input system [A σ A ] and repeat the formal reduction procedure.

Example 21. In this example, we use our algorithm to explain some of the computations in the introductory example in the light of the eigenvalues of the matrix of the system. Thus, we consider again the system [A σ A ] given by (1):

[A σ A ] ε ∂F = A(x, ε)F = 0 1 x 3 -ε 0 F.
We remark that the eigenvalues of A(x, ε) are given by λ 1,2 = ±(x 3 -ε) 1/2 . We thus have the two possible expansions which correspond to the two subdomains above:

• λ 1,2 = ±(x 3 -ε) 1/2 = ±iε 1/2 (1 -1 2 ε -1 x 3 -1 8 ε -2 x 6 + . . . ), valid for |ε -1 x 3 | < 1.

• λ 1,2 = ±(x 3 -ε) 1/2 = ±x 3/2 (1 -1 2 εx -3 -1 8 ε 2 x -6 + . . . ) , valid for |εx -3 | < 1.

  Let ∂ = d dx and consider the decisively simple-looking differential equation [34, Introduction, Subsection 1.5 ]

Figure 1 :

 1 Figure 1:

Figure 2 :

 2 Figure 2: Geometric interpretation of Example 1: (S a ), (S g ), (S f ), (S s ), are traced in black, green, blue, and red respectively.

Figure 3 :

 3 Figure 3: The ε-polygon associated to (12) in Example 2

Figure 4 :

 4 Figure 4: Geometric Interpretation

Lemma 4 .

 4 Consider an ε-irreducible system [A σ A ] given by[START_REF] Chen | Solutions Formelles de Systemes d'Equations Differentielles Lineaires Ordinaires Homogenes[END_REF] with h > 0 and r = rank (A 0 (x)). Let d be an integer such that d ≥ n h-1+ r n and letξA h x p∂G = Ã(x, ξA )G.be the ε-irreducible differential system obtained by the ramification ε = εd (or equivalently ξA = ξ 1/d A x σ A (d-1)/d and σ Ã = σ A (d -1)) and performing ε-rank reduction. Then µ ε( Ã) = h + r n , where r = rank ( Ã0 (x)) and h + r > n.Proof. By Corollary 2 we have, h-1+ r n ≤ ω ε (A) ≤ h. And due to ε = εd , we have ω ε

6 )

 6 Let [C σ C ] ∂G = C(x, ε)G be as in Proposition 5. Due to their equivalence, [A σ A ] and [C σ C ] have the same ε-exponential order and ε-polynomial. Hence, we have

  and we have n = 7, h = 2, and r := rank(A 0 (x)) = 2. Thus, the condition h > nr of Theorem 3 is not verified. Let us first consider a random ramification in ξ A regardless of Lemma 4. For instance, let us consider ξ A = ξ3A and apply ε-rank reduction. This yields system [C σ C ] ξ5C ∂G = C(x, ξC )G with σ C = 0 and C

  rank reduction) then U = T 1 G yields [ Ẽσ Ẽ ] ξ3 Ẽ ∂ τ G = Ẽ(τ, ξ Ẽ )G = 0 ξ Ẽ τ 5 + τ 2 + ξ Ẽ 0G, with σ Ẽ = 0.

[B σ B ] x 12 9

 129 σ B = -3 and ξB = x -3 ε. 55Now that B 0 (x) has two distinct eigenvalues and B 0,0 is nilpotent, we perform U =x σ B = -3 and ξ B = x -3 ε.

[Sτ 3

 3 S σ S ] ξ5 S ∂ τ U = S (τ, ξS ) = 0 ξ6 + 1 0 U with σ S = 0, ξS = ε, ξ S = ξ3 S .The transformation U =

[ 2 S τ 3 ξ 2 S + 1 0 U

 220 S σ S ] ξ S ∂ τ U = S (τ, ξ S )U = 0 ξwith σ S = 0 and ξ S = ε. By the transformation U = σ S = 0 and ξ S = ε. Due to the nature of the eigenvalues of S 0,0 , the system [ S σ S ] can be decoupled into two scalar equations. Its ε-exponential part is zero.

The turning points other than τ = 0, namely the roots of this equation do not explicitly correspond to the original turning point. They are referred to as secondary turning points[START_REF] Nakano | On a Secondary Turning Point Problem[END_REF].

In fact, we want to set εd = ε. And so we choose an integer s such that ξA d = ξ A x s which is equivalent tox σd εd = x σ+s ε.

The package is available at: http : //www.spec f un.inria. f r/ smaddah/Research .html

Acknowledgement

We would like to thank Dr. Reinhard Schaefke for pointing out inaccuracies in the previous version. Our discussions with him highly contributed to improving the quality of this manuscript.

Algorithm 1 Exp param(h, p, σ A , A(x, ξ A )) : Computes the ε-exponential part of system [A σ A ] (performs formal reduction) . Input: h, p, σ A , A(x, ξ A ) of system (8) (p = σ A = 0 for system (3)) Output: Q(x 1/s , ε -1/d ) Q ← diag(0, . . . 0); while h > 0 and n 1 do if A 0,0 has at least two distinct eigenvalues then Apply the ε-block diagonalization of Section 2.3; Exp param(h, p, σ A , Ã11 (x, ε)); Update Q and σ A ; Exp param(h, p, σ A , Ã22 (x, ε)); Update Q and σ A ; else if A 0,0 has one non-zero eigenvalue then Update Q from the eigenvalues of A 0,0 ; A(x, ξ A ) ← perform eigenvalue shifting ( 14)

Appendix A. Computing a companion block diagonal form

It is well-known that a scalar equation [a σ a ] given by ( 9) can be expressed as a first-order system by setting

However, the opposite direction of this transition is nontrivial although possible. The theoretical possibility stems from the work of [START_REF] Churchill | Cyclic vectors[END_REF] with the so called cyclic vectors. Consider again system [START_REF] Barkatou | An Algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System[END_REF] given by:

We also recall that with the definition of ξ A , A(x,

) for all k > 0. In this appendix, we compute a companion block diagonal form for [A σ A ], in a finite number of steps, in a neighborhood of a turning point. This algorithm is an adaptation of that of [START_REF] Barkatou | An Algorithm for Computing a Companion Block Diagonal Form for a system of Linear Differential Equations[END_REF], developed for the unperturbed counterpart of [A σ A ]. It relies on a sequence of polynomial (shearing) transformations in x, ξ A and elementary row/column operations.

• Shearing transformations:

where the α's and β's are respectively rational numbers and integers. We remark that shearing transformations may alter σ A because of the α's.

• Elementary transformations: We consider transformations of the form

where P i, j is a zero matrix except for entry in the i th row and j th column for i, j ∈ {1, . . . , n} and a ∈ R such that a(x = 0, ξ A = 0) 0. Obviously, T ∈ GL n (R). The equivalent system [ Ãσ Ã ] = T [A σ A ] is such that h = h, p = p, and σ Ã ≥ min( p h + σ A , σ A ). Moreover, it can be easily verified that the effect of this transformaiton on the system [A σ A ] is as follows:

Appendix E. Examples

We treat with our algorithm examples from literature. For the computation of the full the ε-exponential parts, we refer to our Maple package ParamInt [START_REF] Maddah | ParamInt : Maple Package for Symbolic Resolution of Singularly-Perturbed Linear Differential Systems, Prototype Implementation[END_REF].

Example 15 (Continue Example 9). We resume the computation of the outer solution of system [B σ B ] given by:

with σ B = -3. The leading matrix coefficient 0 1 0 0 is nilpotent and ε-irreducible. We compute the characteristic polynomial of B(x,ξ B )

.

Then the ε-formal exponential order is ω ε = 3 2 and the ε-polynomial is X 2 + 1 x = 0. Thus, by setting ε = ε2 , ξ B = ξ2 B x 3 ( ξB = x -3 ε) and applying ε-rank reduction of Section 3.2 via diag(1, ξB ), we get an ε-irreducible system whose leading matrix coefficient is 0 1 -x -3 0 . And so, by the turning point algorithm of Section 3.1, we apply the transformation

which yields (note that we can also express the former using a change of independent variable x = t 2 )

[ Bσ B ] ξ3 B x 19/2 ∂U = B(x, ξ B)U where σ B = -3 and