# Rigorous numerical computation of the Stokes multipliers for linear differential equations with single level one

Michèle Loday-Richaud, Marc Mezzarobba, Pascal Remy

► To cite this version:

Michèle Loday-Richaud, Marc Mezzarobba, Pascal Remy. Rigorous numerical computation of the Stokes multipliers for linear differential equations with single level one. 2026. hal-05446686

HAL Id: hal-05446686

https://hal.science/hal-05446686v1

Preprint submitted on 7 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Rigorous numerical computation of the Stokes multipliers for linear differential equations with single level one

MICHÈLE LODAY-RICHAUD (1), MARC MEZZAROBBA (2) AND PASCAL REMY (3)

**ABSTRACT.** — We describe a practical algorithm for computing the Stokes multipliers of a linear differential equation with polynomial coefficients at an irregular singular point of single level one. The algorithm follows a classical approach based on Borel summation and numerical ODE solving, but avoids a large amount of redundant work compared to a direct implementation. It applies to differential equations of arbitrary order, with no genericity assumption, and is suited to high-precision computations. In addition, we present an open-source implementation of this algorithm in the SageMath computer algebra system and illustrate its use with several examples. Our implementation supports arbitrary-precision computations and automatically provides rigorous error bounds. The article assumes minimal prior knowledge of the asymptotic theory of meromorphic differential equations and provides an elementary introduction to the linear Stokes phenomenon that may be of independent interest.

## Introduction

The aim of the present article is to describe a practical algorithm for the numerical computation of the *Stokes matrices* of a linear differential

Keywords: linear differential equations, Stokes multipliers, summability, algorithms, software, rigorous computing.

2020 Mathematics Subject Classification: 34M03, 34M30, 34M35, 34M40, 65G20.

(1) Université Paris-Saclay, CNRS, Laboratoire de mathématiques d'Orsay, 91405, Orsay, France. — michele.loday@universite-paris-saclay.fr

(2) CNRS-LIX 1 rue Honoré d'Estienne d'Orves - Bâtiment Alan Turing CS35003 - 91120 Palaiseau, France — marc@mezzarobba.net

MM's work was supported in part by ANR grants ANR-20-CE48-0014-02 (NuSCAP) and ANR-22-CE48-0016 (NODE).

(3) Laboratoire de Mathématiques de Versailles, UVSQ (Paris-Saclay) & CNRS (UMR 8100), 45 avenue des Etats-Unis, 78035 Versailles cedex, France — pascal.remy@uvsq.fr
equation $Dy(x) = 0$ with coefficients in $\mathbb{C}[x]$, under the assumption that $D$
has a single level equal to one (Definition 1.1, p. 9) at the singular point of
interest.

The Stokes phenomenon is a discontinuity phenomenon, long observed
by physicists interested in numerical computations, in the asymptotic be-
haviour of analytic functions as the variable approaches a singularity from
a varying direction. In the case of solutions of linear differential equations,
the phenomenon originates in the occurrence in the asymptotic expansions
of exponential factors that alternate between being flat in some directions
and exponentially large in others. At the boundaries of a sector of flat-
ness (*Stokes directions*), the asymptotic expansion of a given solution can
suddenly change albeit the solution itself shows no discontinuity. For the
analysis of the Stokes phenomenon, though, it is more appropriate to look
at the bisectors of the sectors of flatness, called *anti-Stokes* or *singular* di-
rections. Understanding the Stokes phenomenon amounts to characterising
the gaps between analytic solutions having the same asymptotic expansion
on sectors that overlap around an anti-Stokes direction. Once appropriately
normalised<sup>(1)</sup>, these gaps are described by a finite number of constant ma-
trices called Stokes matrices.

Starting with Stokes’ own work on the Airy equation and until the 1970’s,
the study of the Stokes phenomenon focused on the calculation of Stokes mul-
tipliers (nontrivial entries of the Stokes matrices) for differential equations
with explicit solutions, typically given by integral formulae. This mainly
concerns the family of hypergeometric equations, including the generalised
hypergeometric equations, whose solutions are expressible in term of Barnes-
Mellin integrals [14]. Such an approach is possible thanks to the fact that
the Mellin transformation translates these equations into first-order differ-
ence equations.

From a theoretical perspective, the linear Stokes phenomenon was fully
understood with the cohomological approach initiated by Sibuya and Mal-
grange [42], which resulted in its characterisation by the Stokes cocycle [37,
39, 1]. In this line of work, the Stokes matrices are viewed primarily as
the *transcendental local invariants* in the classification of linear differen-
tial equations by meromorphic gauge transformations. The introduction of
the Stokes cocycle makes it clear how to choose a non-redundant family of
Stokes matrices—one for each anti-Stokes direction—in full generality. The
cohomological framework also led to a new approach to the determination of

<footnote><sup>(1)</sup> In older literature, the name *Stokes matrix* is used for matrices connecting arbitrary pairs of analytic fundamental solutions asymptotic to the same formal fundamental solution. Here, in accordance with modern usage, we reserve it for those matrices representing the components of the Stokes cocycle mentioned below.</footnote>
Stokes multipliers, resulting essentially from the comparison between the iso-
morphism theorem of Malgrange-Sibuya and its infinitesimal version based
on the Cauchy-Heine integral. The Stokes multipliers are then obtained as
limits of rapidly converging linearly recurrent sequences. This method was
only developed for linear equations with a single level, either of small order
or satisfying genericity assumptions [35, 44, 41, 36]. (See [41] for additional
references.)

In parallel, the development in the 1980’s of several theories of summa-
bility and of the theory of resurgence [15, 16, 17] provided another access to
Stokes multipliers. In this approach each Stokes matrix results from the com-
parison of the sums of a given formal solution on both sides of an anti-Stokes
direction. At a singular point of a single level one, the classic Borel-Laplace
summation method suffices to obtain integral formulae for the sums that can
be used to compute the Stokes multipliers. An explicit formula for the Stokes
multipliers of a linear differential system with the single level 1 is given in [40,
Thm. 4.3] (see also [22]). In the case of multiple levels, a more powerful the-
ory of summation is required. Écalle’s acceleration method [17, 7, 18] was
the first such theory amenable to effective computations. A variant that has
also been used in practice is Balser’s method of iterated Laplace integrals [3].

Starting in the early 1990’s, Thomann and several coauthors explored the
application of summation methods in numerical computations with divergent
solutions of linear differential equations. They showed how the Borel-Laplace
method and its extensions can be implemented by combining techniques from
symbolic computation with classical numerical methods and used to evalu-
ate the sums of divergent solutions [57, 50, 58, 34, 20] or compute the Stokes
matrices [19]. Their work is, to the best of our knowledge, the only one
that went beyond an ad-hoc implementation dedicated to a specific exam-
ple. However, due in part to the choice of numerical methods, it is mostly
limited to machine-precision computations on equations of small order with
sufficiently generic Stokes values and local exponents. In addition, unfortu-
nately, their software is not publicly available.

Van der Hoeven [26, 28] gave a complete algorithm for computing sums of
solutions of linear differential equations with polynomial coefficients, based
on Écalle's acceleration method. His algorithm is optimised for high-precision
computations; in particular, it can compute the Stokes multipliers of a fixed
equation with an error of at most 2⁻¹² in time O(p log³+o(¹) p). In addition,
it yields rigorous bounds on the difference between the numerical result and
the mathematical value(²). This algorithm is the most complete algorithmic

<footnote>(2) This is in contrast with typical numerical methods, which guarantee at best that the computed result converges to the true value at a certain rate when parameters such</footnote>
treatment to date of the summation of divergent solutions of linear differen-
tial equations, but it has not been fully implemented$^{(3)}$.

In the present work, we focus on the computation of Stokes multipliers,
as opposed to the more general problem of evaluating the sums, and we
restrict ourselves to singular points with a single level equal to one. In this
setting, we give a detailed description and a complete implementation of an
algorithm for computing the Stokes matrices of a linear differential equation
with coefficients in $\bar{\mathbb{Q}}[x]$. There is no additional restriction; in particular, the
algorithm applies to equations of arbitrary order and makes no assumption
on the geometry of Stokes values or the shape of indicial polynomials.

The method is very similar to that of Fauvet, Richard-Jung and Tho-
mann [19], but our algorithm is more complete in that it can handle arbitrary
degeneracies with no user intervention, and produces rigorous error bounds.
It differs from van der Hoeven's method in that the numerical computation
of integral transforms of holonomic functions is replaced by expressions in
terms of special functions. This leads to a marginally better computational
complexity estimate at high precision for the special case we consider. A
more important difference in practice, compared to all existing approaches,
is that we pay attention to avoiding redundancies that occur when computing
several Stokes multipliers of the same equation.

Like virtually all algorithms for the problem of computing Stokes mul-
tipliers, the one we develop mixes exact and approximate computations.
Roughly speaking, the geometry of Stokes values and the structure of local
exponents at singular points in the Borel plane (see Sec. 2.3 and 4) must
be determined exactly because one needs to decide when Stokes values are
aligned or local exponents differ by integers. The rest of the computation
is done numerically but, as already mentioned, with rigorous error bounds.
While, in principle, one could labor to determine a priori how accurately to
perform each step of the computation in order for the final result to satisfy
a given error tolerance, it is difficult to obtain sharp error bounds this way
while keeping the calculations manageable. Instead, and following standard
practice, we content ourselves with keeping track of the accumulated error
over the course of the computation and returning an error bound along with

<footnote>as a step size or a numerical working precision tend to zero or infinity. Nevertheless, most
of the other algorithms discussed above could in principle be adapted to produce rigorous
error bounds as well. Richard-Jung [54] presented an example of the computation of Stokes
multipliers with rigorous error bounds.</footnote>

<footnote>(3) A prototype implementation limited to the case of singular points of a single level 1 and due to the second author is available as part of the ore_algebra package mentioned below.</footnote>
the result. This also allows us to refer to existing work [33, 46] for bound-
ing the error resulting from the key analytic steps of the algorithm (namely
numerical integration of differential equations and the evaluation of special
functions). In case the error bound is insufficient for the application, it is
always possible to run the algorithm again at higher precision, and doing so
eventually produces arbitrarily precise results.

Returning error bounds increases confidence in the computed results but
also has algorithmic applications based on differential Galois theory, via
the Ramis density theorem. For instance, checking that one of the Stokes
multipliers of one of the variational equations of a Hamiltonian system is
nontrivial suffices to prove that the system is not completely integrable with
meromorphic first integrals [48, 6, 53]. Criteria of this kind are available for
a number of other problems such as certifying that the differential opera-
tor *D* is irreducible or that it is the minimal-order operator annihilating
a certain function. Importantly, though, such local criteria typically trans-
late into sufficient conditions only. Full characterizations—of irreducibility,
for instance—in the same style, and complete algorithms based on these
characterizations [26, 21], are possible given a representation of the *global*
differential Galois group through monodromy and Stokes automorphisms at
all singular points *expressed in the same basis*. Moving the Stokes automor-
phisms computed by our algorithm (which are expressed in a basis that
depends on the singular point) to the same point is a connection problem
that effective summation methods can serve to solve, and thus our algorithm,
even when it applies, does not eliminate the need for such methods.

We have implemented our algorithm using the SageMath computer algebra system. Our implementation is part of the ore_algebra package, available at https://www.github.com/mkauers/ore_algebra/. It heavily relies on the FLINT library<sup>(4)</sup>, and on pre-existing code by the second author for the computation of connection matrices between regular singular points [45].

**Outline.** This article is organised as follows.

In Part 1, we develop a first, ‘theoretical’ algorithm that reduces the
numerical computation of Stokes matrices to the numerical solution of con-
nection problems between regular singular points. The algorithm amounts
to comparing the 1-sums of a formal fundamental solution to on either side
of each singular direction. For the convenience of the non-expert reader, and
also because we could not find a fully adequate reference, we recall in that
part most of the necessary background, including the notion of equation of

<footnote><sup>(4)</sup> See https://flintlib.org/. More specifically, we mainly use FLINT's modules for real and complex ball arithmetic, formerly distributed as a separate library called ARB [31].</footnote>
single level 1, the Borel–Laplace summation method, the definition of Stokes matrices based on lateral sums, and the expression of the Stokes multipliers in terms of the connection constants of the Borel transform of the differential operator D.

The regular singular connection problems we are left with are much easier than the general ones alluded to above. Solving them numerically reduces to computing the analytic continuation of solutions of differential equations given by convergent series expansions at singular points, or in other words to solving generalised differential initial value problems. Even so, it is the most computationally expensive stage of the method in practice. In Part 2 of the article, we present a more practical variant of the algorithm from Part 1, designed to minimize the use of numerical analytic continuation and be relatively easy to implement. We give a detailed description of the algorithm in terms of operations readily available in today’s computer algebra libraries. We also introduce our implementation and illustrate its use on several examples.

**Acknowledgments.** We are grateful to Bruno Salvy for enabling this collaboration by initiating the contact between us. We also thank him, as well as Fredrik Johansson, Jean-Pierre Ramis, and Anne Vaugon, for valuable comments and discussions.

## Part 1. Underlying theory

All along this paper we consider a linear differential equation

$$Dy(x) = 0 \qquad (0.1)$$

with polynomial coefficients in the complex variable $x$, which we study at an irregular singular point that we position at $x = 0$. For reasons that will become clear later (cf. Sec. 2, p. 13), we write the operator $D$ in terms of the derivation $\partial = x^2 d/dx$ and, dividing by a convenient power of $x$ if necessary, we assume that $D$ is of the form

$$D \equiv D\left(\frac{1}{x}, \hat{\partial}\right) = a_n\left(\frac{1}{x}\right)\hat{\partial}^n + a_{n-1}\left(\frac{1}{x}\right)\hat{\partial}^{n-1} + \dots + a_1\left(\frac{1}{x}\right)\hat{\partial} + a_0\left(\frac{1}{x}\right) \quad (0.2)$$

where the coefficients

$$a_{\ell}(X) = \sum_{j=0}^{\nu} a_{\ell,j} X^j \qquad (0.3)$$

for $\ell = 0, \dots, n$ are polynomials of degree $\le \nu$. We assume that $D$ has order $n$, that is, $a_n(X) \ne 0$, and at least one of the coefficients $a_\ell$ has degree $\nu$.
Starting from Sec. 1.3 we additionally assume the restrictive condition that the singular point at the origin is of single level one (cf. Definition 1.1, p. 9, Hypothesis 1.3, p. 11).

# 1. The equation in the Laplace plane

The plane $\mathbb{C}$ with coordinate $x$ where the initial equation (0.1) is considered is commonly called the *Laplace plane*, recalling thus that it is the image by a Laplace transformation of the *Borel plane* to be considered in the next section. It is well-known that the equation $Dy = 0$ admits $n$ linearly independent formal solutions of the form

$$ \left( \sum_{j=0}^{s-1} \tilde{g}_{i,j}(x) \ln^j x \right) x^{\mathcal{L}_i} e^{q_i(1/x)}. \quad (1.1) $$

Here the term *formal* refers to the fact that the series $\tilde{g}_{i,j}(x)$ may be divergent. Such a family is collectively called a *formal fundamental solution* and one can choose it in the form

$$ \left[ \tilde{f}_1(x) \quad \tilde{f}_2(x) \quad \dots \quad \tilde{f}_n(x) \right] x^{\mathcal{L}} e^{Q(1/x)} \quad (1.2) $$

where the $\tilde{f}_j(x)$'s are (usually divergent) power series of $x$, $\mathcal{L}$ is a constant matrix in Jordan form and $Q(1/x) = \text{diag}(q_1(1/x), \dots, q_n(1/x))$ with $q_j(1/x)$'s that are polynomials in $1/x$ or in fractional powers of $1/x$ with no constant term. These polynomials are called the *determining polynomials of D*. The matrix $\mathcal{L}$ is called the *matrix of exponents of formal monodromy*.

## 1.1. Newton polygons and single level 1

The Newton polygon of the operator $D$ at $0$ is defined as follows [39, Sec. 3.3.3.1]: one defines the weight $w$ of a monomial by $w(x^{-j}\partial^\ell) = \ell - j$ and the weight of $D$ as the smallest weight of a non-zero monomial in $D$. With these non-zero monomials one associates marked points of coordinates $(\ell, \ell - j)$ in the half-plane $\mathbb{R}_{\ge 0} \times \mathbb{R}$ and second quadrants based at these marked points. The Newton polygon $\mathcal{N}(D)$ of $D$ is the convex hull of this family of quadrants. The Newton polygon may contain a horizontal side, one or more sides with positive slopes, and always ends with a vertical line of abscissa $n$ which we do not consider as a side. One can prove that the slopes of the sides (including $0$) are the degrees of the determining polynomials (including the trivial polynomial $q \equiv 0$). Moreover, if the side of slope $p$ has
a horizontal length equal to $r$ then there are $r$ determining polynomials of degree $p$.

We now explain how to calculate the exponents (i.e., the diagonal terms of $\mathcal{L}$) and the determining polynomials $q_j$. For simplicity—since we are interested in this case only—we suppose from now that the operator is of level one, which means that the largest slope of its Newton polygon is equal to 1. This is a weaker condition than being of single level 1 (Hypothesis 1.3, p. 11). One sees on the Newton polygon $\mathcal{N}(D)$ that $D$ is of level 1 if and only if there exists $j, 0 \le j < n$ such that $a_{n,\nu} a_{j,\nu} \ne 0$.

Like for differential equations with constant coefficients the exponents of the solutions (1.1) without exponential factor are the roots of an algebraic equation called the *indicial equation* and built as follows. Suppose that the horizontal side of $\mathcal{N}(D)$ has length $k_0$ and lies at ordinate $h_0$; then, the operator restricted to this horizontal side reads $\sum_{\ell=0}^{k_0} a_{\ell, \ell-h_0} (1/x)^{\ell-h_0} \partial^\ell$. The indicial equation is obtained by writing that the monomial $x^\mathcal{L}$ is a solution of this restricted operator and it reads

$$ \mathcal{P}(\mathcal{L}) \equiv \sum_{\ell=0}^{k_0} a_{\ell, \ell-h_0}[\mathcal{L}]_{\ell+} = 0 \qquad (1.3) $$

with the following notation.

**Notation 1.1.** — For $\lambda \in \mathbb{C}$ and $j \in \mathbb{Z}_{>0}$ we denote

$[\lambda]_{j-} = \lambda(\lambda-1)\cdots(\lambda-j+1) \quad (j \text{ factors decreasing by 1 from } \lambda);$

$[\lambda]_{j+} = \lambda(\lambda+1)\cdots(\lambda+j-1) \quad (j \text{ factors increasing by 1 from } \lambda).$

Saying that the operator is of level 1 is equivalent to saying that the highest degree of its determining polynomials is equal to 1. However, even determining polynomials of degree 1 might contain fractional powers of $1/x$. If the horizontal length of the side of slope 1 is $k_1$, there are $k_1$ determining polynomials of the form $q_j(1/x) = -\alpha_j/x + o(1/x)$, the other $n-k_1$ having no term of degree 1 (i.e., $\alpha_j = 0$). The coefficients $\alpha_j$, zero or not, are the (possibly multiple) roots of the *characteristic polynomial* associated with the side of slope 1 and obtained as follows. One considers the operator $D$ restricted to the side of slope 1, that is, the operator $\sum_{j=0}^{k_1} a_{n-j,r} (1/x)^r \partial^{n-j}$ for a suitable $r$. When the side of slope 1 of the Newton polygon $\mathcal{N}(D)$ is on the first bisector one has $r=0$, and in general $\mathcal{N}(x^r D)$ has a side of slope 1 on the first bisector of the axes. The *characteristic equation* (of level 1) then reads

$$ P(X) \equiv \sum_{j=0}^{k_1} a_{n-j,r} X^{n-j} = 0. \qquad (1.4) $$
It admits $X = 0$ as a root of multiplicity $n - k_1$. The roots of $P(X)$ are called *Stokes values* (of level 1) of $D$; the non-zero ones are called *characteristic roots*. In what follows, when we refer to ‘the’ characteristic polynomial, characteristic roots, or Stokes values of $D$, we always mean the characteristic polynomial of level 1 and its roots.

**DEFINITION 1.1.** — The operator $D$ is said to be of single level 1 (or of pure level 1) if all its nonzero determining polynomials $q_j(1/x)$ are monomials of degree 1 in $1/x$.

Equivalently, $D$ is of single level 1 if all non-zero differences $q_\ell - q_j$ among the determining polynomials of $D$ are of degree 1. Indeed, suppose that there is a determining polynomial $q(1/x) = \alpha/x + \beta/x^{r/s} + \cdots$ with $r/s \in \mathbb{Q}$, $0 < r/s < 1$ and $\beta \neq 0$. Since the change of variable $x \leftarrow x e^{2\pi i}$ leaves the operator $D$ invariant, there must also be a determining polynomial of the form $\tilde{q}(1/x) = \alpha/x + \beta e^{-2r\pi i/s} x^{r/s} + \cdots$. The difference $q(1/x) - \tilde{q}(1/x)$ is a polynomial of degree $r/s$ with $0 < r/s < 1$, hence the operator $D$ is not of single level. The converse assertion is straightforward.

As the terminology suggests, if $D$ is of single level 1, its Newton polygon $\mathcal{N}(D)$ has a unique side of slope 1 in addition to a possible horizontal side, however this condition is not sufficient.

## 1.2. The operators $D_{[\alpha]}$

To characterise operators of single level 1 it is useful to introduce the operators $D_{[\alpha]}$ deduced from $D$ by the change of variable $y = e^{-\alpha/x} z$ with $\alpha$ a Stokes value (of level 1). The operator $D_{[\alpha]}$ has the elementary following properties.

(1) One has $D_{[\alpha]}(1/x, \hat{\partial}) = D(1/x, \hat{\partial} + \alpha)$, and hence

$$ 
\begin{aligned}
D_{[\alpha]} &= a_n\left(\frac{1}{x}\right)(\hat{\partial} + \alpha)^n + a_{n-1}\left(\frac{1}{x}\right)(\hat{\partial} + \alpha)^{n-1} + \cdots + a_0\left(\frac{1}{x}\right) \\
&= A_n\left(\frac{1}{x}\right)\hat{\partial}^n + A_{n-1}\left(\frac{1}{x}\right)\hat{\partial}^{n-1} + \cdots + A_1\left(\frac{1}{x}\right)\hat{\partial} + A_0\left(\frac{1}{x}\right)
\end{aligned}
$$

where, for $\ell = 0, 1, \dots, n$ and denoting by $\binom{j}{\ell}$ the binomial coefficients, the new coefficients $A_\ell(\frac{1}{x}) = \sum_{j=0}^\nu A_{\ell,j} \frac{1}{x^j}$ read

$$ 
\begin{aligned}
A_{\ell}\left(\frac{1}{x}\right) &= a_{n}\left(\frac{1}{x}\right)\binom{n}{\ell}\alpha^{n-\ell} + a_{n-1}\left(\frac{1}{x}\right)\binom{n-1}{\ell}\alpha^{n-1-\ell} + \dots \\
&\quad + a_{\ell+1}\left(\frac{1}{x}\right)\binom{\ell+1}{\ell}\alpha^{1} + a_{\ell}\left(\frac{1}{x}\right)\binom{\ell}{\ell}\alpha^{\ell}.
\end{aligned}
\quad (1.5) $$

In particular, $A_n = a_n$.
(2) For $\ell = 0, \dots, n$, one has $A_{\ell,\nu} = \sum_{j=\ell}^{n} a_{j,\nu} \binom{j}{\ell} \alpha^{j-\ell} = P^{(\ell)}(\alpha) / (\ell!)$ where $P^{(\ell)}$ is the $\ell^{th}$ derivative of the characteristic polynomial $P$ of $D$ (cf. equation (1.4), p. 8).

(3) If $\alpha$ has multiplicity $k$ as a Stokes value, then $A_{k,\nu} \neq 0$.

(4) Therefore, $D_{[\alpha]}$ has order $n$ and degree $\nu$.

(5) If $D$ is of pure level 1 so is $D_{[\alpha]}$ for any $\alpha$.

(6) If $P(X)$ is the characteristic polynomial of $D$ then that of $D_{[\alpha]}$ is

$$P_{[\alpha]}(X) = P(\alpha + X).$$

And indeed, if the polynomials $q_j(1/x)$'s are the determining polynomials of $D$, the determining polynomials of $D_{[\alpha]}$ are the polynomials $q_j(1/x) + \alpha/x$.

(7) The exponents $\mathcal{L}$ of the solutions of $D$ with exponential part $e^{-\alpha/x}$ are the solutions of the indicial equation (cf. (1.3)) of $D_{[\alpha]}$.

**PROPOSITION 1.2.** — Let $D$ be an operator of order $n$ and level 1 at 0. We exclude the trivial case$^{(5)}$ where $D$ has a unique characteristic root of multiplicity $n$. The operator $D$ is of pure level 1 if and only if the following equivalent properties hold:

(1) For any Stokes value $\alpha$ (including $\alpha = 0$) of $D$ of multiplicity $k$, the Newton polygon $\mathcal{N}(D_{[\alpha]})$ has a horizontal side of length $k$ and a side of slope 1 with horizontal length $n-k$. After a vertical translation of height $\nu$ it looks as depicted on Figure 1.1 below.

<figure><img src="image_2.png" /><figcaption>Figure 1.1. Newton polygon of $x^{\nu}D_{[\alpha]}$</figcaption></figure>

<footnote>$^{(5)}$ With a unique characteristic root $\alpha$ of multiplicity $n$ the change of variable $y = e^{-\alpha/xz}$ changes $D$ into an operator with a regular singular point at 0 if $D$ is of single level 1 or an operator of level $< 1$ if $D$ is of level 1 but not of single level 1.</footnote>
(2) For any Stokes value $\alpha$ (including $\alpha = 0$) of $D$ of multiplicity $k$, the coefficients of $D_{[\alpha]}$ satisfy the relations

$$A_{\ell,j} = 0 \quad \text{for all } \ell, j \text{ satisfying } \ell - j < k - \nu, \\ A_{n,\nu} A_{k,\nu} \neq 0.$$

*Proof.* — (1) The condition for $\alpha = 0$ means that the Newton polygon $\mathcal{N}(D)$ of $D$ has no slope other than 0 and 1. Hence, the determining polynomials of $D$ are all of degree 1. We have to prove that they are monomials. Reasoning by contradiction, suppose one of these determining polynomials reads $q(1/x) = -\alpha/x + \beta/x^r + \dots$ with $\beta \neq 0$ and $r < 1$ (necessarily rational). Then, the polynomial $q(1/x) + \alpha/x = \beta/x^r + \dots$ is a determining polynomial of $D_{[\alpha]}$ of degree $r$. This implies that $\mathcal{N}(D_{[\alpha]})$ has a side with slope $r < 1$. Hence a contradiction with assumption (1) asserting that the only possible slopes for $\mathcal{N}(D_{[\alpha]})$ are 0 and 1. Thus the operator $D$ is of single level 1. The converse assertion is straightforward.

(2) The conditions translate the form of each Newton polygon $\mathcal{N}(D_{[\alpha]})$ on the coefficients of $D_{[\alpha]}$. $\square$

HYPOTHESIS 1.3. — *From here on, we assume that $D$ is of single level 1.*

1.3. Prepared fundamental solution

In full generality, the matrices $\mathcal{L}$ and $Q$ appearing in the formal fundamental solution (1.2) of $D y = 0$ do not commute [2], [38, Cor. 2.5]. Under Hypothesis 1.3, however, the determining polynomials $q_j$ are unramified, so that the matrices $\mathcal{L}$ and $Q$ commute. In this case one can choose the formal fundamental solution in the form

$$Y(x) = [\tilde{y}_1(x) \quad \tilde{y}_2(x) \quad \dots \quad \tilde{y}_n(x)] \\ = [\tilde{f}_1(x) \quad \tilde{f}_2(x) \quad \dots \quad \tilde{f}_n(x)] x^{\mathcal{L}} e^{Q(1/x)} \quad (1.6)$$

where

• $Q(1/x) = \text{diag}(q_1(1/x), \dots, q_n(1/x))$ with determining polynomials $q_j$ of the form $q_j(1/x) = -\alpha_j/x$;

• the matrix $\mathcal{L}$ of the exponents of formal monodromy is a constant matrix in Jordan form commuting with $Q$;

• the $\tilde{f}_j(x)$'s are power series in $x$;

• the regular part $\tilde{h}(x) = [\tilde{f}_1(x) \quad \tilde{f}_2(x) \quad \dots \quad \tilde{f}_n(x)] x^{\mathcal{L}}$ has valuation $> 0$, in the sense that for all monomials $x^\beta \log^\rho x$ appearing in the entries of $\tilde{h}(x)$ one has $\Re(\beta) > 0$.
The general form is classical and we stick to it for the theoretical discussion that follows. Note, though, that in the second part of the article we will switch to a slightly different choice of fundamental solution (see Rem. 5.3 and Sec. 9) for consistency with the implementation.

The valuation condition can always be satisfied with a change of variable of type $y = x^{-r} z$ for a convenient integer $r$. Such a transformation changes equation (0.1), p. 6, into an equation of the same form and one can check that it preserves the Stokes matrices we want to compute (cf. (3.1), p. 24).

In terms of the classification of differential equations, the matrices $\mathcal{L}$ and $Q$ provide a complete set of formal invariants. These invariants are obtained by solving algebraic equations and they determine the anti-Stokes directions (Definition 3.1, p. 23). In contrast, the Stokes phenomenon (i.e., the meromorphic invariants) is determined from the divergent series $\tilde{f}_j(x)$ and, except in very simple cases, it results from transcendental calculations.

We additionally suppose that the solutions are so arranged that the matrix $Q(1/x)$ splits into diagonal blocks

$$Q(1/x) = q_1(1/x)I_{k_1} \oplus q_2(1/x)I_{k_2} \oplus \cdots \oplus q_N(1/x)I_{k_N}$$

with pairwise distinct $q_\ell$'s and $q_\ell = -\alpha_\ell/x$. We denote by

$$\left\{ \begin{array}{l} \mathcal{I}_1 = \{1, 2, \dots, k_1\} \\ \mathcal{I}_2 = \{k_1 + 1, k_1 + 2, \dots, k_1 + k_2\} \\ \dots \\ \mathcal{I}_N = \{k_1 + \dots + k_{N-1} + 1, \dots, k_1 + k_2 + \dots + k_N = n\} \end{array} \right. \qquad (1.7)$$

the sets of row indices corresponding to each block of $Q(1/x)$: the solutions with index $i \in \mathcal{I}_\ell$ are those with exponential part $e^{q_\ell}$ (one also says that they are associated with $q_\ell$ or with $\alpha_\ell$). We further denote

$$K_0 = 0, \quad K_1 = k_1, \quad K_2 = K_1 + k_2, \quad \dots, \quad K_N = K_{N-1} + k_N = n$$

and split the fundamental solution $Y(x)$ according to the splitting of $Q$, obtaining thus the block decomposition

$$Y(x) = [Y_i(x)]_{1 \le i \le N} \quad \text{with} \quad Y_i(x) = [\tilde{y}_u(x)]_{u \in \mathcal{I}_i} = [\tilde{y}_u(x)]_{K_{i-1}+1 \le u \le K_i}.$$

We can now arrange the solutions with same determining polynomial $q$ so that the matrix $\mathcal{L}$ splits into a sum of Jordan blocks (cf. Notation 1.2 below)

$$\mathcal{L} = (\mathcal{L}_1 I_{m_1} + J_{m_1}) \oplus (\mathcal{L}_2 I_{m_2} + J_{m_2}) \oplus \cdots \oplus (\mathcal{L}_{m_M} I_{m_M} + J_{m_M})$$

which refines the splitting of $Q$ into blocks.

**Notation 1.2.** — We denote by $I_r$ the identity matrix of dimension $r$ and by $J_r$ the upper nilpotent Jordan matrix of dimension $r$ (with a super-diagonal of ones and zeroes elsewhere).
*Remark 1.4 (generic case).* — Quite often the situation is not as general as the one considered above. A common case—classically called the generic case—is that all Stokes values $\alpha$ are simple ($k = 1$). Then, the matrix $Q$ splits into diagonal blocks of size 1:

$$Q(1/x) = [q_1(1/x)] \oplus [q_2(1/x)] \oplus \cdots \oplus [q_n(1/x)],$$

the matrix $\mathcal{L}$ is diagonal, and hence the formal solutions in the expanded form (1.1) involve no logarithms. The characteristic equation of $D_{[\alpha]}$ reads (cf. equation (1.4), p. 8)

$$P_{[\alpha]}(X) = P(\alpha + X) = A_{n,\nu}X^n + A_{n-1,\nu}X^{n-1} + \cdots + A_{1,\nu}X^1$$

and the conditions in Proposition 1.2 reduce to $A_{0,\nu} = 0$ and $A_{1,\nu}A_{n,\nu} \neq 0$. The indicial equation associated with the Stokes value $\alpha$ reads

$$\mathcal{P}(\mathcal{L}) \equiv A_{0,\nu-1} + A_{1,\nu}\mathcal{L} = 0,$$

corresponding to a unique exponent $\mathcal{L} = -A_{0,\nu-1}/A_{1,\nu}$.

A fundamental property for the algorithm to be developed below is that under Hypothesis 1.3 the series $\tilde{f}_j(x)$ appearing in the fundamental solution (1.6) are all 1-summable [39, Thm. 5.2.5]. One obtains their sums by applying the so-called Borel-Laplace method, as we discuss in the next two sections.

## 2. The equation in the Borel plane

The Borel plane is a second copy of $\mathbb{C}$, with coordinate $\xi$, to which the equation (0.1) is sent by means of a Borel transformation. A Laplace transformation sends the Borel plane back to the Laplace plane.

As much as possible we use Latin letters in the Laplace plane and Greek letters in the Borel plane.

### 2.1. Borel transformations

We will consider two kinds of Borel transformations:

— The Borel transformation $\mathcal{B}$ for operators is the morphism from the ring of linear differential operators $\sum_{j,\ell} c_{\ell,j} \frac{1}{x^j} \partial^\ell$ 'in the Laplace plane' to that of differential operators $\sum_{j,\ell} \gamma_{\ell,j} \xi^j \frac{d}{d\xi^\ell}$ 'in the Borel plane' defined by

$$\mathcal{B}\left(\frac{1}{x}\right) = \frac{d}{d\xi} \quad \text{and} \quad \mathcal{B}(\partial) = \xi \quad (\text{operator of multiplication by } \xi) \quad (2.1)$$
Observe that an operator $D = D(\frac{1}{x}, \partial)$ with polynomial coefficients in $1/x$ is changed into an operator $\Delta = D(\frac{d}{d\xi}, \xi)$ with polynomial coefficients in $\xi$. This property is the main reason why we chose to write our operators in terms of $1/x$ and $\partial$ (positive powers of $x$ would lead to integro-differential operators). If $D$ has order $n$ and degree $\nu$ in $1/x$ then $\Delta$ has order $\nu$ and degree $n$ in $\xi$.

— The (functional) Borel transformation $B_{\theta}$ depends on a direction $\theta$ and applies to actual functions. Written at 0 (as opposed to infinity) it reads

$$B_{\theta}(f(x))(\xi) = \frac{1}{2\pi i} \int_{\gamma_{\theta}} f(x) e^{\xi/x} \frac{dx}{x^2} \quad (2.2)$$

where the contour $\gamma_{\theta}$ is a ‘Hankel contour at 0’, cf. Figure 2.1. Equivalently, one has $B_{\theta}(f(x))(\xi) = \frac{1}{2\pi i} \int_{\Gamma(-\theta+\pi)} f(1/u) e^{\xi u} du$ where the integral is now along the classical Hankel contour $\Gamma(-\theta+\pi)$ around the half-line $d(-\theta+\pi)$ from $\arg u = -\theta - \pi$ to $\arg u = -\theta + \pi$.

We require that the integral converge in the Lebesgue sense.

<figure><img src="image_3.png" /><figcaption>Figure 2.1. A Hankel contour at 0</figcaption></figure>

The following elementary lemmas are useful in the theory that follows. For more details on Borel transformations we refer to [47, Chap. 5], [39, Sec. 5.3] and [11] among many references.

**LEMMA 2.1.** — For any $\lambda \in \mathbb{C}$ and $\xi \in \mathbb{C}\setminus\{0\}$, the Borel transform of monomials $x^{\lambda} \log^m x$ is given by the following formulae in which the Borel integrals converge on the open half-plane $\Re(\xi e^{-i\theta}) > 0$ bisected by $d_{\theta}$.

(1) If $\lambda \in \mathbb{C}\setminus\mathbb{Z}_{\le 0}$, then $B_{\theta}(x^{\lambda})(\xi) = \frac{\xi^{\lambda-1}}{\Gamma(\lambda)} \neq 0$ for any direction $\theta \in \mathbb{C}$.

If $\lambda \in \mathbb{Z}_{\le 0}$ then$^{(6)}$ $B_{\theta}(x^{\lambda})(\xi) = \frac{\xi^{\lambda-1}}{\Gamma(\lambda)} \equiv 0$.

(2) For any $m \in \mathbb{Z}_{\ge 1}$ and any direction $\theta$, one has

$$B_{\theta}(x^{\lambda} \ln^m x)(\xi) = \frac{d^m}{d\lambda^m} B_{\theta}(x^{\lambda})(\xi). \quad (2.3)$$

<footnote>$^{(6)}$ However, the definition of the Borel transform is classically extended by introducing Dirac masses in this case, see Remark 2.3 below.</footnote>
(3) The case $m=1$ is worth stating explicitly. For any direction $\theta$, one has

$$ B_{\theta}(x^{\lambda} \ln x)(\xi) = \xi^{\lambda-1} \left( \frac{\ln \xi}{\Gamma(\lambda)} + \left(\frac{1}{\Gamma}\right)'(\lambda) \right). $$

In particular, when $\lambda = p$ is a positive integer the formula reads

$$ B_{\theta}(x^p \ln x)(\xi) = \frac{\xi^{p-1}}{(p-1)!} \left( \ln \xi + \gamma - \mathcal{H}_{p-1} \right) $$

where $\gamma$ is the Euler constant and $\mathcal{H}_n = 1+\cdots+1/n$ is the harmonic sum of order $n$. For nonpositive integer $\lambda = -p$ one has

$$ B_{\theta}\left(\frac{1}{x^p} \ln x\right) = (-1)^p \frac{p!}{\xi^{p+1}}, $$

and in particular $B_{\theta}(\ln x)(\xi) = \frac{1}{\xi}$.

*Proof.* — (1) By definition, $2\pi i B_{\theta}(x^{\lambda})(\xi) = \int_{\gamma_{\theta}} x^{\lambda} e^{\xi/x} \frac{dx}{x^2}$.

Setting $u = 1/x$ and denoting as before by $\Gamma(-\theta+\pi)$ the Hankel contour around the half-line $d(-\theta+\pi)$ run from $\arg u = -\theta - \pi$ to $\arg u = -\theta + \pi$ we obtain

$$ B_{\theta}(x^{\lambda})(\xi) = \frac{1}{2\pi i} \int_{\Gamma(-\theta+\pi)} u^{-\lambda} e^{\xi u} du. $$

This expression is valid for all $\xi$ in the open half-plane $\Re(\xi e^{-i\theta}) > 0$, that is, the open half-plane bisected by $d_{\theta}$.

Suppose first that $\theta = 0$. The contour is then $\Gamma_{\pi}$ run from $\arg(u) = -\pi$ to $\arg(u) = +\pi$. The function $u^{-\lambda}$ is taken with its principal determination $u^{\lambda} = e^{\lambda \ln u}$ on $-\pi < \arg(u) < \pi$ and, for $\xi > 0$ (i.e., $\arg(\xi) = 0$), we can apply [13, (4.8.1), p. 296] to obtain the relation

$$ B_{0}(x^{\lambda})(\xi) = \frac{1}{2\pi i} \int_{\Gamma_{\pi}} u^{-\lambda} e^{\xi u} du = \frac{\xi^{\lambda-1}}{2\pi i} \int_{\Gamma_{\pi}} v^{-\lambda} e^{v} dv = \xi^{\lambda-1} \frac{1}{\Gamma(\lambda)}. $$

If $\arg(\xi) = \theta' \in (-\pi/2, \pi/2)$ then $\int_{\Gamma_{\pi+\theta'}} v^{-\lambda} e^{v} dv = \int_{\Gamma_{\pi}} v^{-\lambda} e^{v} dv$ by Cauchy's residues theorem, hence the formula is valid for all $\xi$ satisfying $\Re(\xi) > 0$.

Suppose now $\theta \neq \pi \mod 2\pi$, so that the half-planes bisected by $d_{\theta}$ and by $d_0$ have a non-empty intersection (if $\theta = \pi$, one has to proceed in two steps). Choose $\xi$ in the intersection of these two half-planes so that both the integral defining $B_{\theta}(x^{\lambda})(\xi)$ with contour $\Gamma(-\theta+\pi)$ and the same integral with contour $\Gamma_{\pi}$ are convergent. Cauchy's residues theorem shows again that they are equal. Hence, $B_{\theta}(x^{\lambda})$ and $B_{\theta}(x^{\lambda})$ are analytic continuations of each other and we can say that $B_{\theta}(x^{\lambda})(\xi) = \frac{1}{\Gamma(\lambda)}\xi^{\lambda-1}$ for all $\theta, -\pi < \arg(\theta) \le +\pi$ and all $\xi$ satisfying $\Re(\xi e^{-i\theta}) > 0$. In particular, $B_{\theta}(x^{\lambda})(\xi)$ is identically 0 if and only if $\lambda \in \mathbb{Z}_{\le 0}$.
(2) The formula is a direct application of Lebesgue's derivation theorem.

(3) In the general formula above, in order to calculate $(1/\Gamma)' = -\Gamma'/\Gamma^2$ we use the identity (cf. [13, (4.5.2), p. 294]):

$$ \frac{\Gamma'(z)}{\Gamma(z)} = -\gamma - \frac{1}{z} + \sum_{s \ge 1} \frac{z}{s(z+s)} \quad \text{for all } z \in \mathbb{C} \setminus \mathbb{Z}_{\le 0}. \quad (2.4) $$

If $\lambda = p \in \mathbb{Z}_{>0}$ the partial fraction decomposition $\frac{z}{s(z+s)} = \frac{1}{s} - \frac{1}{z+s}$ gives the result.

If now $\lambda = -p \in \mathbb{Z}_{\le 0}$, then $1/\Gamma(-p) = 0$ and the first term $\xi^{\lambda-1} \ln \xi / \Gamma(\lambda)$ of the general formula vanishes. To obtain the value of the second term we write, using (2.4),

$$ \frac{\Gamma'(z-p)}{\Gamma(z-p)} = -\gamma - \frac{1}{z-p} + \sum_{s \ge 1} \frac{z-p}{s(z-p+s)} = \frac{1}{z}(-1 + O(z)). $$

One has

$$ \frac{1}{\Gamma(z-p)} = (-1)^p \frac{[p-z]_{p-}}{\Gamma(z)} = (-1)^p \Gamma(1+p)z(1+O(z)), $$

hence $-\frac{\Gamma'(z-p)}{\Gamma^2(z-p)} = (-1)^p \Gamma(1+p) + O(z)$ and the given formula follows. $\square$

As the lemma shows, the image of monomials

$$ B_{\theta}(x^{\lambda} \ln^m x)(\xi) = \frac{d^{\lambda}}{d\lambda^{\lambda}} \frac{\xi^{\lambda-1}}{\Gamma(\lambda)} \quad (\lambda \notin \mathbb{Z}_{\le 0} \text{ or } m \ge 1) \quad (2.5) $$

is independent of the direction $\theta$. This formula extends naturally to polynomials and formal series with no constant or polar part. The resulting map

$$ \mathcal{B} : x^{1+\lambda} \mathbb{C}[[x]][\log x] \to \xi^{\lambda} \mathbb{C}[[\xi]][\log \xi] \quad (\lambda \in \mathbb{C} \setminus \mathbb{Z}_{<0}) $$

is called the *formal Borel transformation*. It is classically extended to constant and polar terms as explained in Remark 2.3 below. Later in the article we will also use the ‘naive’ extension $\tilde{B}$ obtained by dropping the condition $\lambda \notin \mathbb{Z}_{\le 0}$ from (2.5), that is, by putting $\tilde{B}(x^{-p}) = 0$ for $p \in \mathbb{N}$.

**LEMMA 2.2.** — Given a series $f(x)$ with no constant or polar terms, but possibly with non integer powers of $x$, and a linear differential operator $D$ as before, we consider $B_{\theta}$ as defined in formula (2.2). Then,

$$ B_{\theta}(D)(\xi) = B(D) \cdot (B_{\theta}(f)(\xi)). \quad (2.6) $$

Here we denote with a dot the action of the operator $B(D)$ on the function $B_{\theta}(f)(\xi)$. We allow ourselves to omit it further on.
*Proof.* — It suffices to consider the case of $f(x) = x^{\lambda}$, $\lambda \notin \mathbb{Z}_{\le 0}$ and either $D = \partial$ or $D = 1/x$. Suppose $D = \partial$. Then, $\mathcal{B}(f)(\xi) = \xi^{\lambda-1}/\Gamma(\lambda)$ and $\mathcal{B}(D) = \xi$, hence $\mathcal{B}(D) \cdot \mathcal{B}(f)(\xi) = \xi^{\lambda}/\Gamma(\lambda)$ and therefore,

$$ \mathcal{B}(Df) = \lambda \xi^{\lambda} / \Gamma(\lambda + 1) = \xi^{\lambda} / \Gamma(\lambda) = \mathcal{B}(D) \cdot \mathcal{B}(f)(\xi). $$

Now suppose $D = 1/x$. Then, $\mathcal{B}(D) = d/d\xi$ and

$$ \mathcal{B}(D) \cdot \mathcal{B}(f)(\xi) = \xi^{\lambda-2} / \Gamma(\lambda - 1) = \mathcal{B}(x^{\lambda-1})(\xi) = \mathcal{B}(Df)(\xi). \quad \square $$

*Remark 2.3.* — Being concerned with differential equations it is desirable to work with morphisms of differential algebras. Consider the differential algebra $(\mathbb{C}[x], d/dx)$ of polynomials equipped with the standard derivation. The Borel transformation as defined so far maps monomials $x^m$ with $m \in \mathbb{Z}_{>0}$ to $\mathcal{B}(x^m) = \xi^{m-1}/(m-1)!$ and one can check that

$$ \mathcal{B}(x^m x^p) = \mathcal{B}(x^m) * \mathcal{B}(x^p) $$

where $*$ is the convolution product given by $(f * g)(\xi) = \int_0^\xi f(t)g(\xi - t) \, dt$. Moreover, we have $\mathcal{B}((d/dx)(x^m)) = \xi\mathcal{B}(x^m)$. Hence, $\mathcal{B}$ maps $\mathbb{C}[x]\C$ into the nonunital differential algebra $\mathbb{C}[\xi]$ endowed with the convolution product and the multiplication by $\xi$ as a derivation. To extend $\mathcal{B}$ into a morphism of differential algebras, one has to set $\mathcal{B}(1)$ to the Dirac mass $\delta_0$ since $\delta_0$ is the unit for the convolution product. Extending $\mathcal{B}$ to polynomials in $x$ and $1/x$ leads to set $\mathcal{B}(1/x^m) = \delta_0^{(m)}$, the $m$th derivative of $\delta_0$.

*Remark 2.4.* — Beware of the difference between the Borel transform of the constant function 1, which is the Dirac mass $\delta_0$, and that of the operator 1 = id, which is also the identity operator. As an example, the Borel transform of the differential operator $\partial + 1$ is the operator of multiplication by $\xi + 1$ since $\mathcal{B}(\partial\tilde{y} + \tilde{y})(\xi) = (\xi + 1)\hat{y}(\xi)$.

## 2.2. Transformed operators

Recall the notation $\nu = \max_{0 \le l \le n} \deg a_l(X)$. For all Stokes values $\alpha$, the operator

$$ D_{[\alpha]} = A_n\left(\frac{1}{x}\right)\partial^n + A_{n-1}\left(\frac{1}{x}\right)\partial^{n-1} + \dots + A_2\left(\frac{1}{x}\right)\partial^2 + A_1\left(\frac{1}{x}\right)\partial + A_0\left(\frac{1}{x}\right) $$

admits as a Borel transform

$$ \Delta_{[\alpha]} = A_n\left(\frac{d}{d\zeta}\right)\zeta^n + A_{n-1}\left(\frac{d}{d\zeta}\right)\zeta^{n-1} + \dots + A_1\left(\frac{d}{d\zeta}\right)\zeta + A_0\left(\frac{d}{d\zeta}\right) \quad (2.7) $$

$$ \text{where } A_l\left(\frac{d}{d\zeta}\right) = \sum_{j=0}^{\nu} A_{l,j}\left(\frac{d^j}{d\zeta^j}\right) \text{ for } l = 1, \dots, n \text{ (thus, one has } \Delta_{[0]} = \Delta\text{).} $$

As the following lemma shows, studying $\Delta$ at $\xi = \alpha$ is equivalent to studying $\Delta_{[\alpha]}$ at $\zeta = 0$.
LEMMA 2.5. — (1) Given a point $\alpha \in \mathbb{C}$ and an analytic function $h(x)$ with $h(0) = 0$, the Borel transform satisfies

$$ \mathcal{B}(h(x)e^{-\alpha/x})(\xi) = \mathcal{B}(h(x))(\xi - \alpha). $$

(2) Let $\zeta = \xi - \alpha$ be the local variable at $\xi = \alpha$. Then,

$$ \Delta_{[\alpha]}(y)(\zeta) = \Delta(y)(\alpha + \zeta). $$

*Proof*. — (1) One has

$$ \mathcal{B}(h(x)e^{-\alpha/x})(\xi) = \int_{\gamma_0} h(x)e^{(\xi-\alpha)/x} \frac{dx}{x^2} = \mathcal{B}(h(x))(\xi - \alpha). $$

(2) Applying a Borel transformation to equality $D_{[\alpha]}(\frac{1}{x}, \partial) = D(\frac{1}{x}, \partial + \alpha)$ (item 1 p. 9) we obtain $\Delta_{[\alpha]}(\frac{d}{d\zeta}, \zeta) = \Delta(\frac{d}{d\zeta}, \zeta + \alpha) = \Delta(\frac{d}{d\xi}, \xi)$ for $\xi = \zeta + \alpha$. Hence, the result. $\square$

To define the map $X \mapsto \mathcal{B}(X e^{-\alpha/x})$ on the algebra $\mathbb{C}[[x]]$ of formal power series one has to set $\mathcal{B}(1 e^{-\alpha/x})(\xi) = \delta_{\alpha}$, the Dirac mass at $\alpha$. To extend it to formal Laurent series or to the polar part of meromorphic functions, according to Lemma 2.5, one has to set $\mathcal{B}(x^{-p} e^{-\alpha/x}) = \frac{d^p}{dx^p} \mathcal{B}(1 e^{-\alpha/x}) = \delta_{\alpha}^{(p)}$ for all $p \in \mathbb{Z}_{\le 0}$.

PROPOSITION 2.6. — The operator $\Delta$ satisfies the following properties:

(1) The singular points of $\Delta$ are the Stokes values $\alpha$ of $D$ at $0$;

(2) All singular points of $\Delta$ are regular singular points;

(3) Given $\Delta_{[\alpha]}$ as in (2.7), the indicial polynomial of $\Delta_{[\alpha]}$ at $\zeta = 0$ reads (cf. Notation 1.1, p. 8)

$$ \Pi_{[\alpha]}(\lambda) = \sum_{\substack{l-j=k-\nu \\ 0 \le l \le n \\ 0 \le j \le \nu}} A_{l,j} [l+\lambda]_{j-}. $$

*Proof*. — (1) The singular points are the zeros of the leading term of

$$ a_{n,\nu} \frac{d^{\nu}}{d\xi^{\nu}} \xi^n + \dots + a_{1,\nu} \frac{d^{\nu}}{d\xi^{\nu}} \xi + a_{0,\nu} \frac{d^{\nu}}{d\xi^{\nu}} = (a_{n,\nu} \xi^n + \dots + a_{1,\nu} \xi + a_{0,\nu}) \frac{d^{\nu}}{d\xi^{\nu}} + \text{terms of lower order.} $$

Thus, the coefficient of $d^{\nu}/d\xi^{\nu}$ is equal to the characteristic polynomial $P(\xi)$ of $D$ (cf. (1.4)). Its roots, that is, the Stokes values of $D$ at 0, are then also the singular points of $\Delta$.

(2) Given a Stokes value $\alpha$ of multiplicity $k$, let us prove that the Newton polygon of $\Delta_{[\alpha]}$ reduces to a single, horizontal edge (cf. Figure 2.2, p. 21).
We know from Proposition 1.2, p. 10, that $A_{\ell,j} = 0$ for all $(\ell,j)$ satisfying $\ell - j < k - \nu$ and $A_{k,\nu}A_{n,\nu} \neq 0$. Thus, all marked points of the Newton polygon $\mathcal{N}(\Delta_{[\alpha]})$ are located above the horizontal line at height $k - \nu$. The order of $\Delta_{[\alpha]}$ is equal to $\nu$ and the point $(\nu, k - \nu)$ is a marked point since $A_{k,\nu} \neq 0$. Hence the result.

(3) The terms of weight $k - \nu$ in $\Delta_{[\alpha]}$ are exactly the terms $A_{\ell,j} \frac{d^j}{d\zeta^j} \zeta^\ell$ satisfying $\ell - j = k - \nu$. One obtains $\Pi_{[\alpha]}(\lambda)$ by writing that the operator restricted to its terms with weight $k - \nu$ vanishes when applied to $\zeta^\lambda$. $\square$

## 2.3. Solutions of the transformed equation

From the previous results we know that the Borel transform of any formal solution of $Dy = 0$ at $x = 0$ with exponential part $\exp(-\alpha/x)$ is a solution of $\Delta_{[\alpha]}\hat{y} = 0$ at $\zeta = 0$ (or equivalently, of $\Delta\hat{y} = 0$ at $\xi = \alpha$) free of exponentials.

In this section we suppose that $\alpha$ is a Stokes value of multiplicity $k$: there exist $k$ linearly independent solutions of $Dy = 0$ with exponential part $\exp(-\alpha/x)$, and their Borel transforms are solutions of $\Delta_{[\alpha]}\hat{y} = 0$ at $\zeta = 0$. However, the operator $D$ has order $n$ while the operator $\Delta_{[\alpha]}$ has order $\nu$. We know that $k$ satisfies $k \le n$, but $k$ may be smaller than, equal to or larger than $\nu$. We now look at what happens in each case, starting with the easiest case $k = \nu$.

**Case $k = \nu$.** The equation $Dy = 0$ has $k$ linearly independent solutions attached to $\alpha$, say $\tilde{h}_1(x)e^{-\alpha/x}$, ..., $\tilde{h}_k(x)e^{-\alpha/x}$ as in the prepared fundamental solution in Sec. 1.3, p. 11. Equivalently, $\tilde{h}_1(x)$, ..., $\tilde{h}_k(x)$ is a maximal family of linearly independent solutions of $D_{[\alpha]}y = 0$ attached to the Stokes value 0. The exponents of $D_{[\alpha]}$ at 0 are the $k$ roots of the indicial polynomial $P_{[\alpha]}(\mathcal{L}) = \sum_{\ell=0}^k A_{\ell,\ell}[\mathcal{L}]_{\ell^+}$.

The Borel transforms $\tilde{h}_1(\zeta)$, ..., $\tilde{h}_k(\zeta)$ of $\tilde{h}_1(x)$, ..., $\tilde{h}_k(x)$ are linearly independent. Since $k = \nu$, they provide a fundamental solution of $\Delta_{[\alpha]}\hat{y} = 0$ at $\zeta = 0$.

We have seen that the exponents of $\tilde{h}_1, \dots, \tilde{h}_k$ are those of $\tilde{h}_1(x), \dots, \tilde{h}_k(x)$ decreased by 1. And indeed, the indicial polynomials $\Pi_{[\alpha]}(\lambda)$ and $P_{[\alpha]}(\mathcal{L})$ satisfy the relation

$$ \Pi_{[\alpha]}(\lambda) = \sum_{\ell=0}^{\nu} A_{\ell,\ell}[\lambda + \ell]_{\ell^-} = \sum_{\ell=0}^{\nu} A_{\ell,\ell}[\lambda + 1]_{\ell^+} = P_{[\alpha]}(\lambda + 1). \quad (2.8) $$
Moreover, taking into account Lemma 2.1, p. 14 one can observe that the solutions $\tilde{h}_1(x)e^{-\alpha/x}$, ..., $\tilde{h}_\nu(x)e^{-\alpha/x}$ contain logarithms at the same powers as their Borel transforms.

**Case** $k < \nu$. As in the previous case, to the $k$ linearly independent solutions $\tilde{h}_1(x)e^{-\alpha/x}$, ..., $\tilde{h}_k(x)e^{-\alpha/x}$ of the equation $Dy = 0$ attached to $\alpha$, there correspond $k$ linearly independent solutions of $\Delta_{[\alpha]}\hat{y} = 0$ given by

$$ \hat{h}_1 = \mathcal{B}(\tilde{h}_1), \dots, \hat{h}_k = \mathcal{B}(\tilde{h}_k). $$

However, the equation $\Delta_{[\alpha]}\hat{y} = 0$ has $\nu > k$ linearly independent solutions. Where do the $\nu - k$ extra solutions come from?

Alongside $D$, consider the differential operator

$$ D_1 = \partial^{\nu-k} D_{[\alpha]} = \partial^{\nu-k} \sum_{\ell=0}^{n} A_{\ell}(1/x) \partial^{\ell}. \qquad (2.9) $$

The Newton polygon $\mathcal{N}(x^\nu D_1)$ of $x^\nu D_1$ is obtained by translating that of $x^\nu D$ (cf. Figure 1.1, p. 10) by $\nu - k$ units along each axis. It has a horizontal side of length $\nu$ when $\mathcal{N}(x^\nu D)$ has a horizontal side of length $k$. The indicial polynomial of $D_1$ is $P_1(\mathcal{L}) = [\mathcal{L} - 1]_{(\nu-k)-} P_{[\alpha]}(\mathcal{L})$. The exponents of $D_1$ are hence those of $D$ and the extra numbers $1, 2, \dots, \nu - k$. Since these exponents all have positive real parts, the discussion in the case $k = \nu$ above applies to $D_1$.

The space of solutions attached to 0 of the equation $D_1y = 0$ is generated by, on the one hand, the solutions $\tilde{h}_1(x), \dots, \tilde{h}_k(x)$ of $D_{[\alpha]}y = 0$, and, on the other hand, solutions of each of the inhomogeneous equations

$$ D_{[\alpha]} y = 1, \quad D_{[\alpha]} y = 1/x, \quad \dots, \quad D_{[\alpha]} y = 1/x^{\nu-k-1}. \qquad (2.10) $$

These solutions can be viewed as *microsolutions* of $Dy = 0$ in the sense of Malgrange [43]; see in particular Theorem 2.2 and Example (2.4) in that reference for a related problem. The Borel transform $\Delta_1$ of $D_1$ is $\zeta^{\nu-k}\Delta_{[\alpha]}$, hence the equations $\Delta_1\hat{y} = 0$ and $\Delta_{[\alpha]}\hat{y} = 0$ have the same series solutions. The indicial equation of $\Delta_1$ reads

$$ \Pi_1(\lambda) = [\lambda]_{(\nu-k)-} \Pi_{[\alpha]}(\lambda) \quad \text{with} \quad \Pi_{[\alpha]}(\lambda) = P_{[\alpha]}(\lambda+1) $$

and we can check that $\Pi_1(\lambda) = P_1(\lambda + 1)$. Thus, the exponents of $\Delta_1$ are, on the one hand, the roots of $\Pi_{[\alpha]}(\lambda)$, that is, the exponents of $D_{[\alpha]}$ shifted by $-1$ (which we call the *free exponents* since they can in principle take any complex value), and, on the other hand, the integers $0, 1, \dots, \nu - k - 1$ (which we call the *trivial exponents*). The Newton polygon $\mathcal{N}(\Delta_1)$ is the Newton polygon $\mathcal{N}(\Delta_{[\alpha]})$ translated vertically to the horizontal axis, see Figure 2.2.
<figure><img src="image_4.png" /><figcaption>Figure 2.2. The Newton polygons of $\Delta_{[\alpha]}$ and $\Delta_1$ when $k < \nu$</figcaption></figure>

As in the case $k = \nu$, the Borel transforms $\hat{h}_1, \dots, \hat{h}_k$ of $\tilde{h}_1, \dots, \tilde{h}_k$ are linearly independent solutions of $\Delta_{[\alpha]}\hat{y} = 0$. Their exponents are the free exponents of $\Delta_{[\alpha]}$ and they contain logarithms at the same power as $\tilde{h}_1, \dots, \tilde{h}_k$.

What about the solutions with trivial exponents? Since trivial exponents differ from each other (and possibly from some of the free exponents) by integers, one may expect these solutions to involve logarithms. More precisely, consider the exponent $\Lambda = \nu - k - 1 - s \ge 0$, so that there are $s$ trivial exponents larger than or equal to $\Lambda$. Suppose that there are $r \ge 0$ integer roots of $\Pi_{[\alpha]}(\lambda)$ (counted with multiplicity) larger than or equal to $\Lambda$. In this situation, the classical Frobenius method for constructing solutions involving logarithms yields expressions that may contain logarithms at powers up to $r + s - 1$, and generically one of them reaches this bound. This turns out not to be the case here: in fact, solutions attached to$^{(7)}$ the exponent $\Lambda$ behave 'as if the other trivial exponents did not exist'.

**PROPOSITION 2.7.** — With $\Lambda, r$, and $s$ as above, a solution attached to the trivial exponent $\Lambda$ of equation $\Delta_{[\alpha]}\hat{y} = 0$ may contain logarithms at a power at most $r$. In particular, if $\Pi_{[\alpha]}(\lambda)$ has no integer root larger than $\Lambda$, a solution attached to the trivial exponent $\Lambda$ contains no logarithm.

*Proof.* — As we have just seen, the space of solutions of $\Delta_{[\alpha]}\hat{y} = 0$ attached to the exponent $\Lambda$ is generated by the Borel transforms of the solutions attached to $\Lambda+1$ of the homogeneous equation $D_{[\alpha]}y = 0$ and the Borel transforms of one solution attached to $\Lambda+1$ of each of the inhomogeneous equations (2.10) that has one.

The method of Frobenius [52, 29] shows that the solutions of $D_{[\alpha]}y = 0$ attached to the exponent $\Lambda+1$ involve logarithms at powers at most $r-1$. In homogeneous form the equation $D_{[\alpha]}y = 1/x^p$ for $p \in \mathbb{Z}_{\ge 0}$ reads $\frac{d}{dx}x^p D_{[\alpha]}y = 0$

<footnote>$^{(7)}$ We mean here solutions in which one can factor $\xi^\Lambda$, the other factor containing only series in non negative integer powers of $\xi$ and logarithms.</footnote>
and has at most one integer exponent larger than or equal to $\Lambda + 1$ in addition to those of $D_{[\alpha]}$. Its solutions attached to $\Lambda + 1$ thus involve logarithms at powers at most $r$. (More precisely, the indicial equation of $\frac{d}{dx}x^p D_{[\alpha]} y = 0$ is $(\mathcal{L} + k - \nu + p)P_{[\alpha]}(\mathcal{L}) = 0$, so that the additional exponent is $\Lambda + 1 + s - p$ and the inhomogeneous equation has a solution attached to $\Lambda + 1$ if and only if $0 \le p \le s$.)

By Lemma 2.1, p. 14, the Borel transforms of all these solutions involve logarithms at powers at most $r$. $\square$

Let us now consider a *generic* operator $D$ of degree $\nu > 1$ in the sense of Remark 1.4, p. 13. From that remark we know that, for any Stokes value $\alpha$, the operator $D_{[\alpha]}$ has a unique exponent $\mathcal{L}_{[\alpha]} = -A_{0,\nu-1}/A_{1,\nu}$. The operator $\Delta_{[\alpha]}$ admits the exponent $\Lambda = \mathcal{L}_{[\alpha]} - 1$ and the $\nu - 1$ trivial exponents $0, 1, 2, \dots, \nu - 2$.

**COROLLARY 2.8.** — For generic $D$ of degree $\nu > 1$, with $\mathcal{L}_{[\alpha]}$ defined as above one can set:

* if $\mathcal{L}_{[\alpha]} \notin \mathbb{Z}_{>0}$, the solutions of $\Delta_{[\alpha]}\hat{y} = 0$ contain no logarithm;

* if $\mathcal{L}_{[\alpha]} \in \mathbb{Z}_{>0}$, then the solutions of $\Delta_{[\alpha]}\hat{y} = 0$ with exponent $\lambda$ satisfying $\lambda \in \{0, 1, 2, \dots, \nu - 2\}$ and $\lambda < \Lambda = \mathcal{L}_{[\alpha]} - 1$ may involve logarithms, but only at the first power.

**Case** $\nu < k$. In this case, the order $\nu$ of the Borel operator $\Delta_{[\alpha]}$ is smaller than the number $k$ of solutions of $D_{[\alpha]}y = 0$ attached to 0.

The conditions for $D_{[\alpha]}$ to have the single level 1 (Proposition 1.2, p. 10) say, in particular, that $A_{\ell,j} = 0$ as soon as $\ell - j < k - \nu$. When $\ell < k - \nu$ these conditions are satisfied for all $j$, which implies $A_\ell(1/x) \equiv 0$. Thus, one can write $D_{[\alpha]} = D_2 \partial^{k-\nu}$ for some differential operator $D_2$, and the rational functions of nonpositive valuation 1, $1/x$, $1/x^2$, ..., $1/x^{k-\nu-1}$ are solutions of $D_{[\alpha]}y = 0$. This cannot happen with an equation prepared as in Sec. 1.3, p. 11. To make $D_{[\alpha]}$ fit the conditions of Sec. 1.3 one is led to make the change of variable $y = x^{-(k-\nu)} z$, after which the equation falls in case $k = \nu$.

### 3. The Stokes phenomenon in the Laplace plane

We now recall some key facts about the *Stokes phenomenon* of the linear differential equation $D y = 0$ (equation (0.1), p. 6), including the definition of Stokes matrices. We still assume that the equation is prepared as in Sec. 1.3,
p. 11 and that we have fixed a formal fundamental solution

$$[\tilde{y}_1(x) \ \tilde{y}_2(x) \ \dots \ \tilde{y}_n(x)] = [\tilde{f}_1(x) \ \tilde{f}_2(x) \ \dots \ \tilde{f}_n(x)] x^{\mathcal{L}} e^{Q(1/x)}.$$

As already mentioned, describing the Stokes phenomenon amounts to measuring the gaps between solutions with same asymptotics on a common sector based at 0. However, there are infinitely many ways of presenting these gaps and one has to organise them in order to understand what happens and to be able to compare various results. From a theoretical perspective, the Stokes phenomenon of equation (0.1) is characterised by a non-abelian 1-cocycle called the Stokes cocycle that can be described as a finite collection of Stokes automorphisms ([37, Th. II.2.1], [1], [39, Sec. 3.5.3]). There is a Stokes automorphism associated with each anti-Stokes direction, that is, each direction $\omega \in \mathbb{R}/(2\pi\mathbb{Z})$ in which there exist determining polynomials $q_j(1/x)$ and $q_\ell(1/x)$ such that $\arg(x) = \omega \in \mathbb{R}/(2\pi\mathbb{Z})$ is a direction of maximal decay for $\exp(q_j - q_\ell)(1/x)$. In the case under consideration where $q_j(1/x) = -\alpha_j/x$ and $q_\ell(1/x) = -\alpha_\ell/x$ this means that $-(\alpha_j - \alpha_\ell)e^{-i\omega}$ is real negative [39, Th. 3.3.5 (iv), (v)].

**DEFINITION 3.1.** — Let $\alpha_1, \dots, \alpha_N$ be the Stokes values of D. The anti-Stokes directions are the directions $\omega$ given by

$$\omega = \arg(\alpha_j - \alpha_\ell) \quad \text{for all } j \neq \ell.$$

If $\tilde{y}$ is a solution attached to the Stokes value $\alpha_\ell$, the anti-Stokes directions associated with the solution $\tilde{y}$ are the directions $\omega$ such that there exists $j \neq \ell$ with $\omega = \arg(\alpha_j - \alpha_\ell)$.

All details of the abstract definition of Stokes automorphisms can be found in the references cited above. For the purposes of this paper, let us simply note that the 0-cochain from which the Stokes cocycle is built consists of analytic functions asymptotic to $[\tilde{f}_1(x) \ \tilde{f}_2(x) \ \dots \ \tilde{f}_n(x)]$. One can prove that these functions coincide with the 1-sums of the corresponding asymptotic series on each side of the anti-Stokes directions.

More precisely, the formal series $\tilde{f}_j(x)$ are 1-summable [39, Def. 5.1.6, Th. 5.2.5], with 1-sums in all directions except, possibly, the finitely many anti-Stokes directions. In an anti-Stokes direction $\omega$ one defines two lateral 'sums', $f_j^-(x)$ to the right of $\omega$ and $f_j^+(x)$ to the left of $\omega$: the sum to the right (resp. the left) is obtained by gluing together the 1-sums of $\tilde{f}_j$ in the directions $\omega - \varepsilon$ (resp. $\omega + \varepsilon$) for all small enough$^{(8)}$ $\varepsilon > 0$. On a neighbourhood of 0 in the open half-plane bisected by $\omega$, the lateral sums

<footnote>$^{(8)}$ It is enough that the sector $(\omega - \varepsilon, \omega + \varepsilon)$ does not contain any anti-Stokes direction other than $\omega$. From now, this condition is assumed.</footnote>
are both defined and both 1-Gevrey asymptotic to $\tilde{f}_j(x)$. However, in general, they are not 1-sums of $\tilde{f}_j(x)$ in the direction $\omega$, as for this they would have to be 1-Gevrey asymptotic to $\tilde{f}_j(x)$ on a wider sector of opening $\pi/2 - \varepsilon$. There is a nontrivial Stokes automorphism in the direction $\omega$ if $f_j^-(x)$ and $f_j^+(x)$ do not agree on this common domain.

### 3.1. Stokes matrices

To obtain a matrix representation of a Stokes automorphism one has to look at the gap between fundamental solutions associated to the lateral sums. To define these analytic fundamental solutions, the formal factors $x^\omega$ and $\exp Q(1/x)$ are changed into actual functions in a natural way, with the following choices:

* the determination of the argument for $\omega$ and neighbouring directions: we use the principal determination $-\pi < \omega \le \pi$;

* the value $e^0 = 1$ for the exponential function at 0.

The products $[f_1^\pm \ f_2^\pm \ \dots \ f_n^\pm] x^\omega \exp Q(1/x)$ then provide two analytic fundamental solutions $[y_1^\pm \ y_2^\pm \ \dots \ y_n^\pm]$ of the equation on a neighbourhood of 0 in the same sector $|\arg x - \omega| < \pi/2$. They are hence connected by a constant invertible matrix and we can state the following definition.

**DEFINITION 3.2.** — The Stokes matrix in the direction $\omega$ is the unique constant invertible matrix $I + C_\omega$ satisfying

$$ [y_1^-(x) \ \dots \ y_n^-(x)] = [y_1^+(x) \ \dots \ y_n^+(x)] (I + C_\omega), \quad (3.1) $$

on the open half plane bisected by $\omega$ near 0. We denote $C_\omega = [C_\omega^{[j;\ell]}]_{1\le j, \ell \le n}$.

One proves that the Stokes matrix $I + C_\omega$ is unipotent (i.e., $C_\omega$ is nilpotent): the diagonal entries are equal to 1; for $j \ne \ell$ the entry $C_\omega^{[j;\ell]}$ is zero unless $\exp(q_j - q_\ell)(1/x) = \exp(-\alpha_j + \alpha_\ell)/x$ has maximal decay, in other words, unless $(-\alpha_j + \alpha_\ell) e^{-i\omega} < 0$. Thus, for a given $\ell$, if there is no $\alpha_j$ such that $\arg(\alpha_j - \alpha_\ell) = \omega$, then the $\ell$th column of $C_\omega$ is zero. For more details we refer to [37].

Note that the Stokes matrix depends not only on the choice of a formal fundamental solution but also on that of the determination $\omega$, whereas the Stokes automorphism it represents depends only on the direction $\omega$.
**Example 3.3.** — Consider the case when the matrix $\mathcal{L}$ contains a Jordan block of size 2 and a block of size 1: $\mathcal{L} = (\mathcal{L}_1 I_2 + J_2) \oplus (\mathcal{L}_3 I_1)$. necessarily $q_1(1/x) = q_2(1/x) := -\alpha_1/x$ and we assume that $q_1(1/x) \neq q_3(1/x)$ with $q_3(1/x) := -\alpha_3/x$. A formal fundamental solution reads

$$ 
\left\{ 
\begin{aligned}
\tilde{y}_1(x) &= \tilde{f}_1(x) x^{\mathcal{L}_1} e^{-\alpha_1/x}, \\
\tilde{y}_2(x) &= (\tilde{f}_2(x) + \tilde{f}_1(x) \ln x) x^{\mathcal{L}_1} e^{-\alpha_1/x}, \\
\tilde{y}_3(x) &= \tilde{f}_3(x) x^{\mathcal{L}_3} e^{-\alpha_3/x}.
\end{aligned}
\right. 
$$

The unique anti-Stokes direction for $y_1(x)$ and $y_2(x)$ is $\omega_1 = \arg(\alpha_3 - \alpha_1)$. The unique anti-Stokes direction for $y_3(x)$ is $\omega_3 = \arg(\alpha_1 - \alpha_3)$.

The Stokes matrices $I + C_{\omega_1}$ and $I + C_{\omega_3}$ are of the form

$$ I + C_{\omega_1} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ C_{\omega_1}^{[3;1]} & C_{\omega_1}^{[3;2]} & 1 \end{bmatrix} \quad \text{and} \quad I + C_{\omega_3} = \begin{bmatrix} 1 & 0 & C_{\omega_3}^{[1;3]} \\ 0 & 1 & C_{\omega_3}^{[2;3]} \\ 0 & 0 & 1 \end{bmatrix}. $$

We observe that both matrices are unipotent and the following relations hold:

$$ 
\left\{ 
\begin{aligned}
y_1^-(x) - y_1^+(x) &= y_3^+(x) C_{\omega_1}^{[3;1]} \\
y_2^-(x) - y_2^+(x) &= y_3^+(x) C_{\omega_1}^{[3;2]}
\end{aligned}
\right. 
$$

and

$$ y_3^-(x) - y_3^+(x) = y_1^+(x) C_{\omega_3}^{[1;3]} + y_2^+(x) C_{\omega_3}^{[3;2]}. $$

Note that $\omega_1$ (resp. $\omega_3$) is not an anti-Stokes direction for $y_3$ (resp. $y_1$ and $y_2$); hence the zeroes in the third column of $C_{\omega_1}$ (resp. in the third row of $C_{\omega_3}$).

## 3.2. Case of the prepared fundamental solution

Let us now focus, in the notation of Sec. 1.3, p. 11, on the first block $\mathcal{L}_1$ of size $m_1$ in the matrix of exponents $\mathcal{L}$. Recall that $\mathcal{L}_1$ is associated with a sub-block of the first block $Q_1 = q_1(1/x) I_{k_1}$ of the matrix $Q$ of determining polynomials. Considering the first $m_1$ formal solutions will allow us to determine the first $m_1$ columns of the Stokes matrix $I + C_{\omega}$. There is no loss of generality in restricting ourselves to $\mathcal{L}_1$ since any block of the decomposition can be brought in first position by a permutation of the solutions.

The first $m_1$ formal solutions $\tilde{Y}_1(x) = [\tilde{y}_1(x) \ \tilde{y}_2(x) \ \dots \ \tilde{y}_{m_1}(x)]$ are of the form

$$ 
\left\{ 
\begin{aligned}
\tilde{y}_1(x) &= \tilde{f}_1(x) x^{\mathcal{L}_1} e^{q_1(1/x)} \\
\tilde{y}_2(x) &= (\tilde{f}_2(x) + \tilde{f}_1(x) \ln x) x^{\mathcal{L}_1} e^{q_1(1/x)} \\
\dots \\
\tilde{y}_{m_1}(x) &= (\tilde{f}_{m_1}(x) + \dots + \tilde{f}_1(x) \ln^{m_1-1} x) x^{\mathcal{L}_1} e^{q_1(1/x)}
\end{aligned}
\right. 
\quad (3.2) $$
They are associated with the anti-Stokes directions

$$ \arg(\alpha_2 - \alpha_1), \dots, \arg(\alpha_N - \alpha_1). $$

Given one such direction $\omega$, with principal determination $\underline{\omega}$ we set

$$ \beta_1 = \alpha_{i_1} - \alpha_1, \quad \beta_2 = \alpha_{i_2} - \alpha_1, \quad \dots, \quad \beta_r = \alpha_{i_r} - \alpha_1 \quad (3.3) $$

where $\alpha_{i_1}, \dots, \alpha_{i_r}$ are the Stokes values lying on the ray $d_{\omega}$ originating from 0 in the direction $\omega$, ordered by increasing distance from 0.

In restriction to the first $m_1$ columns, the first $k_1$ rows in $C_{\underline{\omega}}(x)$ are zero and only the rows $i \in \mathcal{I}_{i_1} \cup \mathcal{I}_{i_2} \cup \dots \cup \mathcal{I}_{i_r}$ corresponding to the blocks

$$ Q_{i_{\ell}}(1/x) = q_{i_{\ell}}(1/x) I_{k_{i_{\ell}}} = -\alpha_{i_{\ell}}/x I_{k_{i_{\ell}}} \quad \text{for } \ell = 1, 2, \dots, r $$

may be nonzero. From the definition (3.1) of the Stokes matrix we derive

$$ h_{\ell}^{-}(x) - h_{\ell}^{+}(x) = \sum_{s=1}^{r} e^{-\beta_{s}/x} \left( \sum_{j \in \mathcal{I}_{i_s}} h_{j}^{+}(x) C_{\underline{\omega}}^{[j:\ell]} \right) \quad \text{for } \ell = 1, 2, \dots, m_1 \quad (3.4) $$

where $h_j(x)$ is defined by $y_j(x) = h_j(x)e^{-\alpha_i/x}$ when $j \in \mathcal{I}_i$.

## 4. The Stokes phenomenon viewed from the Borel plane

The idea is now to rewrite the definition (3.1) seen from the Borel plane in order to get a system of linear equations for the Stokes multipliers. As in the previous section, we focus on the block

$$ [\tilde{y}_1 \quad \dots \quad \tilde{y}_{m_1}] = [\tilde{h}_1 \quad \dots \quad \tilde{h}_{m_1}] e^{-\alpha_1/x} $$

with $[\tilde{h}_1 \quad \dots \quad \tilde{h}_{m_1}] = [\tilde{f}_1 \quad \dots \quad \tilde{f}_{m_1}] x^{\mathcal{L}_1}$ and $\mathcal{L}_1 = \mathcal{L}_1 I_{m_1} + J_{m_1}$.

The theory says that the power series $\tilde{f}_{\ell}(x)$ are 1-summable in all direction but the anti-Stokes directions and that one can obtain these sums by applying a Borel transformation followed by a Laplace transformation. Given the anti-Stokes direction $\omega$, we are interested in the sums in the directions $\omega \pm \varepsilon$ for $\varepsilon$ small enough. Thanks to the condition on the valuation required for prepared solutions we can apply these transformations directly to the $\tilde{h}_{\ell}(x)$. Indeed, since $\tilde{h}_{\ell}$ has positive valuation, the valuation of $\hat{h}_{\ell} = \mathcal{B}(\tilde{h}_{\ell})$ is $> -1$ and the Laplace integrals converge at 0. As a solution of a linear differential equation with no singular point on the rays $d_{\omega \pm \varepsilon}$, $\hat{h}_{\ell}$ is continuous up to infinity along these rays. The resulting functions have at most exponential growth at infinity because the operator $\Delta_{[\alpha_1]}$ is of level $\le 1$ at infinity, and hence the Laplace integrals converge.
To obtain the left hand-side of formula (3.1) one should a priori first analytically continue $\hat{h}_{\ell}(\zeta)$ to the right and to the left of $d_{\omega}$, then take a Laplace integral of $\hat{h}_{\ell}(\zeta)$ from $\zeta = 0$ to infinity in direction $\omega - \varepsilon$, and then a Laplace integral from infinity to $\zeta = 0$ in direction $\omega + \varepsilon$. It is more convenient to split the integration path $\gamma$ into $r$ paths $\gamma_j$ going around each singular point $\beta_1, \beta_2, \dots, \beta_r$ in the $\zeta$-plane (i.e., the points $\alpha_{i_1}, \alpha_{i_2}, \dots, \alpha_{i_r}$ in the $\xi$-plane) as depicted on Figure 4.1.

<figure><img src="image_5.png" /><figcaption>Figure 4.1. The path $\gamma_{\omega}$ split into finitely many paths $\gamma_j$; the big black dots represent the singular points of $\Delta$</figcaption></figure>

The solutions $\hat{h}_{\ell}(\zeta)$, for $\ell = 1, \dots, m_1$ of the linear differential equation $\Delta_{[\alpha_1]}\hat{y} = 0$ may be analytically continued along any path that avoids the singular points of $\Delta_{[\alpha_1]}$. The analytic continuation of the $\hat{h}_{\ell}(\zeta)$'s along the path of integration $\gamma_{\omega}$ should be done as follows: first, one defines $\hat{h}_{\ell}(\zeta)$ to the right of $d_{\omega}$ by analytic continuation along a path $\gamma_{0\to\beta_r}$ as shown on Figure 4.2. Then, once the value of $\hat{h}_{\ell}(\zeta)$ is known on an arc of each $\gamma_j$, one considers the analytic continuations along both branches of the $\gamma_j$'s from these values. We denote by $\hat{h}_{\ell}^-$ the resulting function. Note that in practice, only the analytic continuation along the path $\gamma_{0\to\beta_r}$ will need to be computed explicitly, and this can be done in a finite number of steps by means of the Cauchy-Weierstrass method (cf. sec. 6, p. 31).

<figure><img src="image_6.png" /><figcaption>Figure 4.2. The path $\gamma_{0\to\beta_r}$ in the $\zeta$-plane ($\zeta = \xi - \alpha_1$)</figcaption></figure>
Integrating along the split path $\gamma$ we obtain, for $\ell = 1, \dots, m_1$,

$$h_{\ell}^{-}(x) - h_{\ell}^{+}(x) = \sum_{s=1}^{r} e^{-\beta_{s}/x} \int_{\gamma_{0}} \hat{h}_{\ell}^{-}(\beta_{s} + \zeta) e^{-\zeta/x} d\zeta \quad (4.1)$$

where $\gamma_0$ is the Hankel type path around $\zeta = 0$ obtained from $\gamma_{i_s}$ by the translation $-\beta_s$.

## 5. Comparing the two approaches

Comparing equations (3.4), p. 26, and (4.1) we obtain, for $\ell = 1, 2, \dots, m_1$,

$$\sum_{s=1}^{r} e^{-\beta_s/x} \left( \sum_{j \in \mathcal{I}_{i_s}} h_j^+(x) C_{\underline{\omega}}^{[j; \ell]} \right) = \sum_{s=1}^{r} e^{-\beta_s/x} \int_{\gamma_0} \hat{h}_{\ell}^{-} (\beta_s + \zeta) e^{-\zeta/x} d\zeta.$$

Perhaps the most subtle point in the approach is that one can identify the exponential terms on each side despite the fact that they have functions as coefficients. The reason for this is that these functions are summable-resurgent [40, Lemma 4.2]. We deduce for all column indices $\ell = 1, \dots, m_1$ and row block indices $s = 1, \dots, r$ (recall that $\mathcal{I}_{i_s}$ is a set of $k_{i_s}$ consecutive rows indices) the relation

$$\sum_{j \in \mathcal{I}_{i_s}} h_j^+(x) C_{\underline{\omega}}^{[j; \ell]} = \int_{\gamma_0} \hat{h}_{\ell}^{-} (\beta_s + \zeta) e^{-\zeta/x} d\zeta. \quad (5.1)$$

This formula is to be understood as follows. Given $s$, the functions $h_j^+(x)$ for $j \in \mathcal{I}_{i_s}$ form a basis of a linear space of dimension $k_{i_s}$ and we must find the coordinates $C_{\underline{\omega}}^{[j; \ell]}$ of the function of the right-hand side in this basis. While the functions $h_j^+$ and $\hat{y}_{\ell}^-$ are theoretically well-defined, it is convenient for the effective calculation of the coordinates to replace both sides of the equation (5.1) by their asymptotic expansions as $x \to 0$. In other words, (5.1) is advantageously replaced by

$$\sum_{j \in \mathcal{I}_{i_s}} \tilde{h}_j(x) C_{\underline{\omega}}^{[j; \ell]} = T_{(x=0)} \int_{\gamma_0} \hat{h}_{\ell}^{-} (\beta_s + \zeta) e^{-\zeta/x} d\zeta \quad (5.2)$$

where the asymptotic expansion of the integral on the right hand-side is easily calculated with the help of the following lemma.

**LEMMA 5.1.** — Let $\Gamma_{\underline{\omega}}$ be a Hankel contour around the half-line $d_{\omega}$, oriented positively from $\arg \zeta = \underline{\omega} - 2\pi$ to $\arg \zeta = \underline{\omega}$, and let $x \in \mathbb{C}$ be a point of the half-plane $|\arg(x) - \underline{\omega}| < \pi/2$.
(1) Suppose that $\varphi(\zeta)e^{-\zeta/x}$ is integrable at infinity in the direction $\omega$ and that $\varphi(\zeta)$ is analytic on and inside $\Gamma_{\omega}$. Then,

$$ \int_{\Gamma_{\omega}} \varphi(\zeta) e^{-\zeta/x} d\zeta = 0. $$

(2) For any $\lambda \in \mathbb{C}$, one has

$$ I_{\lambda}(x) := \int_{\Gamma_{\omega}} \zeta^{\lambda} e^{-\zeta/x} d\zeta = \frac{2\pi i}{\Gamma(-\lambda)} (x e^{-\pi i})^{1+\lambda}. $$

In particular, $I_{\lambda} = 0$ for all $\lambda \in \mathbb{N}$ in accordance with item (1).

(3) For any $\lambda \in \mathbb{C}$ and $p \in \mathbb{N}$, one has

$$ J_{\lambda;p}(x) := \int_{\Gamma_{\omega}} \zeta^{\lambda} \ln^p(\zeta) e^{-\zeta/x} d\zeta = \frac{d^p}{d\lambda^p} I_{\lambda}(x) $$

and $J_{\lambda;p}(x)$ is of the form $x^{1+\lambda}R(\ln(x))$ with $R(\ln(x)) \in \mathbb{C}[\ln(x)]$.
In particular, for $\lambda \in \mathbb{N}$,

$$ J_{\lambda;1} = 2\pi i \Gamma(1 + \lambda) x^{1+\lambda} $$

$$ J_{\lambda;2} = 4\pi i \Gamma(1 + \lambda) (-\gamma + \mathcal{H}_{\lambda} - \pi i + \ln x) x^{1+\lambda}. $$

where $\mathcal{H}_{\lambda} = \sum_{j=1}^{\lambda} \frac{1}{j}$ is the harmonic sum of order $\lambda$.

**Remark 5.2.** — The integrands in the lemma are functions with no singular point other than 0 on a sector neighbouring $d_{\omega}$. Thus, the contour of integration may be deformed so as to look like the contours $\gamma_0$ in formula (4.1) or $\gamma_s$'s on figure 4.1, p. 27.

*Proof.* — Item (1) results from the residue theorem of Cauchy.

In the next two items the integrals $I$ and $J$ converge at infinity for any $x$ such that $-\pi/2 < \omega - \arg(x) < \pi/2$.

To prove (2), suppose first that $\arg(x) = \omega$ (the case we are interested in). The change of variable $u = \zeta e^{i\pi/x}$ leads to

$$ \int_{\Gamma_{\omega}} \zeta^{\lambda} e^{-\zeta/x} d\zeta = (x e^{-i\pi})^{\lambda+1} \int_{\Gamma_{\pi}} u^{\lambda} e^{u} du $$

The integral to the right is the one appearing in the integral formula for $1/\Gamma(-\lambda)$ [13, (4.8.1) p. 296] valid for all $\lambda \in \mathbb{C}$. The case $\arg(x) = \omega + \varphi$ with $|\varphi| < \pi/2$ reduces to the previous one by rotation.

The general formula in point (3) results from (2) and the theorem of derivation of Lebesgue. The proof of the special cases is similar to that of Lemma 2.1 (3), p. 14, except that we now use the second-order expansion

$$ \frac{1}{\Gamma(-\lambda - z)} = (-1)^{\lambda} \Gamma(1 + \lambda) z (1 + (\mathcal{H}_{\lambda} - \gamma)z + O(z^2)) $$
(which can be deduced from the series expansion $1/\Gamma(z) = z + \gamma z^2 + O(z^3)$
[51, 5.7.i] by applying the relation $\Gamma(z+1) = z\Gamma(z)$ repeatedly). Substituting
this expansion into the expression of $I_\lambda$, we get

$$
\begin{align*}
I_{\lambda+z}(x) &= 2\pi i (-1)^\lambda \Gamma(1+\lambda) z (1 + (\mathcal{H}_\lambda - \gamma)z + O(z^2)) \\
&\quad \cdot (xe^{-\pi i})^{1+\lambda} (1 + (\ln x - \pi i)z + O(z^2)) \\
&= 2\pi i x^{1+\lambda} \Gamma(1+\lambda) (z + (\mathcal{H}_\lambda - \gamma + \ln x - \pi i)z^2 + O(z^3))
\end{align*}
$$

and the result follows.

<figure><img src="image_7.png" /></figure>

Remark 5.3. — In Sec. 1.3, p. 11, while defining the prepared fundamen-
tal solution, we multiplied the solutions of the original equation by a power
of x to reduce to solutions of positive valuation. Doing so ensures that the
Laplace transforms of the Borel transforms of the solutions are well defined.
However, it also increases the order of the transformed operator: suppose that
we had to make the change of variable z = x⁻¹y in the original equation
D₁z = 0 to get the prepared equation Dy = 0. Then, the Borel transform Δ
of D has order u more than the Borel transform Δ₁ of D₁. Thus, from a
practical perspective, it would be more economical to work with Δ₁, and as
we show below this is indeed possible.

The Stokes multipliers are provided by formula (5.2), which is relative
to the prepared operator *D*. Denote by $\tilde{k}_j$ and $\hat{k}_\ell^-$ the quantities analogous
to $\tilde{h}_j$ and $\hat{h}_\ell^-$ but defined relative to the equation $D_1 z = 0$. Let us first
rewrite the right-hand side of (5.2) in terms of $\hat{k}_\ell^-$. For this, split $\tilde{h}_j(x)$
as $\tilde{P}_j(x) + \tilde{H}_j(x)$ with a polynomial part $\tilde{P}_j(x)$ of degree $u$. Accordingly,
write $\tilde{k}_j(x) = \tilde{R}_j(1/x) + \tilde{K}_j(x)$ where $\tilde{K}_j(x) = (1/x^u)\tilde{H}_j(x)$ and $\tilde{R}_j(1/x)$
is the polar part. After a Borel transformation, we get the decompositions
(compatible with analytic continuation)

$$
\hat{h}_{\ell}^{-}(\zeta) = \hat{P}_{\ell}(\zeta) + \hat{H}_{\ell}^{-}(\zeta), \quad \hat{k}_{\ell}^{-}(\zeta) = \hat{R}_{\ell}(\delta) + \hat{K}_{\ell}^{-}(\zeta)
$$

where $\hat{P}_\ell(\zeta)$ is a polynomial of degree at most $u-1$ and $\hat{R}_\ell(\delta)$ is a linear
combination of derivatives of Dirac masses at $\zeta = \xi - \alpha_1 = 0$. The rela-
tion $\tilde{K}_\ell(x) = (1/x^u)\tilde{H}_\ell(x)$ implies $\hat{K}_\ell(\zeta) = \frac{d^u}{d\zeta^u}\hat{H}_\ell(\zeta)$. Moreover, $\hat{H}_\ell(\zeta)$ and
its derivatives up to the order $u$ have at most exponential growth at infinity,
hence one has

$$
\int_{\gamma_0} \hat{K}_{\ell}^{-}(\beta_s + \zeta) e^{-\zeta/x} d\zeta = \int_{\gamma_0} \frac{d^u}{d\zeta^u} \hat{H}_{\ell}^{-}(\beta_s + \zeta) e^{-\zeta/x} d\zeta
$$
and after $u$ integrations by parts

$$ 
\begin{aligned}
\int_{\gamma_0} \hat{K}_{\ell}^{-}(\beta_s + \zeta) e^{-\zeta/x} d\zeta &= \frac{1}{x^u} \int_{\gamma_0} \hat{H}_{\ell}^{-}(\beta_s + \zeta) e^{-\zeta/x} d\zeta \\
&= \frac{1}{x^u} \int_{\gamma_0} \hat{h}_{\ell}^{-}(\beta_s + \zeta) e^{-\zeta/x} d\zeta.
\end{aligned}
$$

Turning to the left-hand side of (5.2), we have $\tilde{k}_j(x) = (1/x^u)\tilde{h}_j(x)$, hence

$$ \sum_{j \in \mathcal{I}_{i_s}} \tilde{k}_j(x) C_{\tilde{\omega}}^{[j; \ell]} = \frac{1}{x^u} \sum_{j \in \mathcal{I}_{i_s}} \tilde{h}_j(x) C_{\tilde{\omega}}^{[j; \ell]} $$

and equation (5.2) still holds with the $h$'s replaced by $k$'s.

In summary, given $\ell \in \{1, \dots, m_1\}$ (in the notation introduced p. 12) and $s$, the $k_{i_s}$ Stokes multipliers $C_{\tilde{\omega}}^{[j; \ell]}$ appearing in identity (5.1) are located in rows $j \in \mathcal{I}_{i_s}$ and column $\ell$ of the Stokes matrix. Replacing all functions by their asymptotic expansions and identifying the coefficients of the expansions on both sides of the equality, we obtain an infinite inhomogeneous linear system whose solutions are the Stokes multipliers $(C_{\tilde{\omega}}^{[j; \ell]})_{j \in \mathcal{I}_{i_s}}$.

To be able to compute the asymptotic expansion of the integral on the right-hand side, we need to know the asymptotic expansion of $\tilde{y}_{\ell}(\zeta)$ at each of its singular points $\beta_s$. Computing these expansions is the main task of the effective analytic continuation procedure discussed in the next section.

By repeating these computations for all $\ell \in \{1, \dots, m_1\}$ and $s = 1, \dots, r$, we can determine all Stokes multipliers in the first $m_1$ columns that are associated with $\mathcal{L}_1$ and $e^{-\alpha_1/x}$. Next, the same procedure must be applied to the other sub-blocks $\mathcal{L}_i$ associated with $e^{-\alpha_1/x}$, and finally to all exponential parts $e^{-\alpha_j/x}$, $j = 1, \dots, N$.

# 6. Analytic continuation

To conclude this first part, we now recall how the analytic continuation of $\hat{h}_{\ell}$ along the path $\gamma = \gamma_0 \to \beta_r$ depicted on Figure 4.2 can be carried out in practice using the Cauchy-Weierstrass method of series expansion on overlapping discs.

As previously, we consider a solution $\tilde{y}_{\ell}(x) = \tilde{h}_{\ell}(x) e^{-\alpha_1/x}$, $1 \le \ell \le m_1$. The Borel transform $\hat{h}_{\ell}(\zeta)$ of $\tilde{y}_{\ell}(x)$, with $\zeta = \xi - \alpha_1$, is a sum of products of power series, complex powers of $\zeta$, and logarithms. The power series are convergent and we interpret the complex powers and logarithms as analytic
functions by choosing the principal determination of arg ζ. Also recall that
$\hat{h}_{\ell}(\zeta)$ is a solution of $\Delta_{[\alpha_1]}$ and that

$$
\beta_1 = \alpha_{i_1} - \alpha_1, \beta_2 = \alpha_{i_2} - \alpha_1, \dots, \beta_r = \alpha_{i_r} - \alpha_1
$$

are the singular points of $\Delta_{[\alpha_1]}$ on the ray $d_\omega$, ordered from from $\zeta = 0$ to infinity.

Since we need to know the behaviour of $\hat{h}(\zeta)$ at all $\beta_s$, $s = 1, 2, \dots, r$, we choose a covering of $\gamma$ by a chain of discs $D_0, D_1, \dots, D_t$ whose centres include all singular points $\beta_s$. We denote by $\chi_0, \chi_1, \dots, \chi_t$ their centres (in the $\zeta$-plane) and by $1/\rho_0, \dots, 1/\rho_t$ their radii. The first disc $D_0$ is centred at $\chi_0 = 0$; the last one $D_t$ is centred at $\chi_t = \beta_r$; and the discs contain no singular point of $\Delta_{[\alpha_1]}$ except possibly at their centres. We number the discs

<figure><img src="image_8.png" /><figcaption>Figure 6.1. Covering Cov(γ₀→β_r) of the path γ₀→β_r; the big black dots represent the singular points of Δ[α₁] limiting the size of the discs, the blue dots indicate the remaining centres, the points η₁, ..., η_t are chosen in the successive intersections</figcaption></figure>

so that any two consecutive discs overlap and we fix points

$$
\eta_1 \in D_0 \cap D_1, \eta_2 \in D_1 \cap D_2, \dots, \eta_t \in D_{t-1} \cap D_t.
$$

We denote by $\zeta_0, \zeta_1, \dots, \zeta_t$ the local coordinates at $\chi_0, \chi_1, \dots, \chi_t$, obtained
by setting

$$
\zeta_0 = \xi - \alpha_1, \quad \zeta_1 = \zeta_0 - \chi_1, \quad \dots, \quad \zeta_{t-1} = \zeta_0 - \chi_{t-1}, \quad \zeta_t = \zeta_0 - \chi_t
$$

and choosing the principal determinations $-\pi < \arg \zeta_p \le +\pi$ for all $p$.
For each disc $D_p$, we choose a local fundamental solution $\hat{Y}_{[\chi_p]}(\zeta_p)$ of the equation $\Delta_{[\chi_p]}\hat{y} = 0$ at $\zeta_p = 0$ in the classical form

$$ \begin{aligned} \hat{Y}_{[\chi_p]}(\zeta_p) &= [\hat{y}_{[\chi_p],1}(\zeta_p) \quad \hat{y}_{[\chi_p],2}(\zeta_p) \quad \dots \quad \hat{y}_{[\chi_p],\nu}(\zeta_p)] \\ &= [\hat{h}_{[\chi_p],1}(\zeta_p) \quad \hat{h}_{[\chi_p],2}(\zeta_p) \quad \dots \quad \hat{h}_{[\chi_p],\nu}(\zeta_p)] \zeta_p^{\Lambda_p} \end{aligned} \quad (6.1) $$

where the $\hat{h}$'s are power series. The matrix of exponents $\Lambda_p$ is zero when $\chi_p$ is an ordinary point; otherwise, it is in Jordan form with diagonal terms given by proposition 2.6, p. 18.

Denote by $\phi_0(\zeta_0)$ the sum of the series $\hat{y}_\ell(\zeta_0)$ on its disc of convergence $D_0$. By evaluating this sum and its first $\nu-1$ derivatives at $\zeta_0 = \eta_1$, we express it in the basis $\hat{Y}_{[\chi_1]}(\zeta_1)$ at the point $\zeta_1 = \eta_1 - \chi_1$ (choosing the principal determination of the argument $-\pi < \arg(\eta_1 - \chi_1) \le +\pi$). This yields a constant vector $T_1$ such that

$$ \phi_0(\eta_1) = \hat{Y}_{[\chi_1]}(\eta_1 - \chi_1) T_1 \quad (6.2) $$

and a function

$$ \phi_1(\chi_1 + \zeta_1) = \hat{Y}_{[\chi_1]}(\zeta_1) T_1 $$

which is the analytic continuation of $\phi_0(\zeta_0)$ to $D_1$. Iterating the process, we obtain a function

$$ \phi_p(\chi_p + \zeta_p) = \hat{Y}_{[\chi_p]}(\zeta_p) T_p \quad (6.3) $$

defined in the neighbourhood of each $\chi_p$ which is the analytic continuation along $\gamma$ of $\phi_0(\zeta_0)$. Due to the uniqueness of analytic continuation, the result is independent of the choice of the discs $D_0, \dots, D_t$. In practice, for better performance, it is desirable to have large intersections so that we can choose the points $\eta$ not too close to the boundary.

The local expansions of $\hat{h}_\ell^-$ at the singular points $\beta_p$, $p = 1, 2, \dots, r$ that enter into the asymptotic expansion on the right-hand side of (5.2) can then be read off equation (6.3) for the disc $D_p$ with centre $\chi_p = \beta_s$, since the basis functions in $\hat{Y}_{[\chi_p]}$ were defined from their local expansions in the first place.

## Part 2. Algorithm and implementation

We now aim to turn the procedure for computing Stokes matrices described in the previous two sections into a detailed algorithm that can be implemented using existing features of computer algebra systems. The detailed algorithm eliminates numerous redundancies in the method, making it much faster in practice. However, we leave for future work a full complexity analysis and confine ourselves to observing that the algorithm can be implemented so as to compute the Stokes multipliers of any fixed operator $D$ within an error of $2^{-p}$ in $O(p \log(p)^3)$ operations (see Remark 13.8, p. 64).
Compared to the algorithm sketched at the end of Sec. 5, p. 31, the complete version works as much as possible with fundamental matrices, or at least with blocks of solutions attached to a given Stokes values and exponent in the decomposition of Sec. 1.3, p. 11, instead of considering each solution $\tilde{y}_j(x)$ individually. In the same spirit, it computes the Stokes matrices in all directions simultaneously. The main benefit of this modified structure is that the number of analytic continuation steps needed is minimized.

We have implemented this algorithm—with some minor changes with respect to the present description—using the SageMath computer algebra system. Our implementation is freely available under the GNU General Public Licence (version 2 or later) as part of the ore_algebra package, which can be downloaded from https://www.github.com/mkauers/ore_algebra/. The algorithm is implemented in the module ore_algebra.analytic.stokes. Below we give some examples of its use and occasionally comment on implementation details. For more information and additional examples, we refer the reader to the documentation shipped with the package. After installing ore_algebra and starting SageMath, one can access the documentation of the main function of the stokes module with the commands

```
sage: from ore_algebra.analytic.stokes import stokes_dict
```

```
sage: ?stokes_dict
```

In the following sections, we begin by illustrating the usage of the implementation on an example (Sec. 7) that will also serve as a running example when detailing the algorithm. Then we discuss the representation of complex numbers (Sec. 8) and generalised power series (Sec. 9) on which the algorithm operates. We next explain how, in each Stokes matrix, the block associated with a given pair of Stokes values decomposes as a product of three auxiliary matrices that can be computed independently and reused for several blocks (Sec. 10). This allows us to state the main procedure (Sec. 11), initially in terms of subroutines for computing the auxiliary matrices. In the next two sections, we detail these subroutines, the most important one being that corresponding to the analytic continuation step (Sec. 13). We conclude with some additional examples (Sec. 14).

## 7. An example

A. Duval and C. Mitschi [14] compute the Stokes matrices of the confluent generalised hypergeometric operator

$$D_{q,p} = (-1)^{q-p} z \prod_{j=1}^{p} \left( z \frac{d}{dz} + \mu_j \right) - \prod_{j=1}^{q} \left( z \frac{d}{dz} + \nu_j - 1 \right) \quad (7.1)$$
<figure><img src="image_9.png" /><figcaption>Figure 7.1. The Stokes values of the operator $D$ in our running example, with the spanning tree considered in Sec. 13.</figcaption></figure>

in terms of values of the gamma function. As a running example we use the special case

$$p = 1, \quad q = 7, \quad \mu = \left(\frac{1}{2}\right), \quad \nu = \left(1, 0, 0, 0, \frac{2}{3}, \frac{1}{3}\right) \qquad (7.2)$$

considered in Example 2.45 of Mitschi and Sauzin's lecture notes [47]. In Example 14.1 (p. 65), we discuss several instances of the case $q = 2, p = 1$ to illustrate various degenerate situations.

*Example 7.1.* — After the change of variables $x = z^{-1/6}$ and $\partial = x^2 d/dx$ as in Part 1, the operator (7.1) specialised to the parameters (7.2) becomes

$$D = \frac{1}{x^6}\partial^7 + \frac{9}{x^5}\partial^6 + \frac{58}{x^4}\partial^5 + \frac{272}{x^3}\partial^4 + \frac{897}{x^2}\partial^3 + \frac{1875}{x}\partial^2 + \left(-\frac{46656}{x^6} + 1875\right)\partial + \frac{139968}{x^5}.$$

This is an operator of pure level 1, with Stokes values 0 and $6e^{j\pi/3}$, $-2 \le j < 2$ (Figure 7.1).

In order to reproduce Mitschi's example numerically, we load the relevant part of `ore_algebra` and set up an algebra of differential operators with the commands

```
sage: from ore_algebra.analytic.stokes import stokes_dict
```

```
sage: R.<x> = PolynomialRing(QQ)
```

```
sage: A.<Dx> = OreAlgebra(R)
```
These last two lines define x as the indeterminate of the ring $R = \mathbb{Q}[x]$ and $Dx$ as the indeterminate of a skew polynomial ring $A = \mathbb{Q}[x] \langle D_x \rangle$ where $D_x$ acts on expressions involving $x$ as the standard derivation $d/dx$. We can thus define the operator $D$ with

```
sage: d = x^2*Dx
```

```
sage: diffop = (1/x^6*d^7 + 9*1/x^5*d^6 + 58*1/x^4*d^5
```

```
....: + 272*1/x^3*d^4 + 897*1/x^2*d^3 + 1875*1/x*d^2
```

```
....: + (-46656*1/x^6 + 1875)*d + 139968*1/x^5)
```

Then, to compute the Stokes matrices of $D$, we run the implementation of the algorithm described in this section with the command

```
sage: stokes = stokes_dict(diffop, 10^-50)
```

The second argument, $10^{-50}$, is a rough indication of the accuracy we would like to obtain and is used internally to select the working precision for intermediate computations. We will discuss its role in more detail in what follows. The computation on this example takes around 1 s on a mid-range laptop purchased in 2020. The `stokes_dict` function returns a dictionary `stokes` mapping $e^{i\omega}$ where $\omega$ is an anti-Stokes direction to the corresponding Stokes matrices. The choice of $e^{i\omega}$ rather than $\omega$ itself as an indexing key is motivated by technical details of the infrastructure provided by SageMath, which make it more convenient to use algebraic numbers as keys.

In particular, for the direction $\omega = 0$, we obtain a matrix of the form (the output has be reformatted and intermediate digits omitted for readability)

```
sage: stokes[1]
```

$$
\begin{bmatrix}
1.0 & 0 & \cdots \\
0 & 1.0 & \cdots \\
0 & 0 & \cdots \\
[\pm 2.11 \cdot 10^{-47}] + [0.14 \dots 78 \pm 3.86 \cdot 10^{-47}]_i & 0 & \cdots \\
0 & [-2.5 \pm 5.88 \cdot 10^{-35}] + [4.33 \dots 09 \pm 7.62 \cdot 10^{-35}]_i & \cdots \\
0 & 0 & \cdots \\
[16.0 \pm 6.61 \cdot 10^{-42}] + [\pm 6.61 \cdot 10^{-42}]_i & 0 & \cdots \\
\hline
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
1.0 & 0 & 0 & 0 & 0 \\
0 & 1.0 & 0 & 0 & 0 \\
0 & 0 & 1.0 & 0 & 0 \\
[2.5 \pm 2.36 \cdot 10^{-23}] + [4.33 \dots 86 \pm 3.95 \cdot 10^{-23}]_i & 0 & 0 & 1.0 & 0 \\
0 & [\pm 4.36 \cdot 10^{-44}] + [-221.7 \dots 25 \pm 5.58 \cdot 10^{-44}]_i & 0 & 0 & 1.0 & 0
\end{bmatrix}
$$

An entry $[c \pm r]$ in this output indicates that the corresponding exact value is contained in the interval $[c - r, c + r]$. Observe that some of the radii $r$ are significantly larger than the tolerance $10^{-50}$ specified on input, which is treated as indicative only. However, decreasing the tolerance results in
tighter enclosures, as illustrated here on the nontrivial entry of the third column (note that row and column indices in Sage start from zero):

```
sage: stokes_dict(diffop, 10^-100)[1][5,2]
```

$$[2.5 \pm 4.41 \cdot 10^{-69}] + [4.3301270189 \dots 1744862983 \pm 6.95 \cdot 10^{-69}]i$$

By further decreasing it one can compute, in principle, arbitrarily tight enclosures of the exact Stokes multipliers.

We note in passing that, exceptionally, this Stokes multiplier is an algebraic number, so that it is easy to guess its minimal polynomial from the approximation, here using the `algdep` function from PARI [56], and from there the exact value $5e^{i\pi/3}$:

```
sage: algdep(stokes_dict(diffop, 10^-100)[1][5,2].mid(), 2)
x^2 - 5*x + 25
```

Similarly, for $\omega = \pi/6$, we have

```
sage: stokes[exp(I*pi/6)]
```

$$\begin{bmatrix}
1.0 & 0 & 0 \\
0 & 1.0 & 0 \\
0 & 0 & 1.0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
[-5.5 \pm 5.33 \cdot 10^{-31}] + [9.52 \dots 82 \pm 8.31 \cdot 10^{-31}]_i & 0 & 0 \\
0 & [5.5 \pm 1.06 \cdot 10^{-32}] + [9.52 \dots 23 \pm 1.26 \cdot 10^{-32}]_i & 0
\end{bmatrix}$$

We omit the remaining ten Stokes matrices with $-\pi < \omega \le \pi$, which, due to the symmetry resulting from the deramification step, are easily expressed in terms of the previous two.

Mitschi's results are relative to a formal fundamental solution $\mathcal{Z}(x)$ specified in Propositions 1.3 and 2.1 of [14] and related to our formal fundamental solution $\mathcal{Y}(x)$ (detailed in Example 9.4 below) by

$$\mathcal{Z}(x) = \mathcal{Y}(x) \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & c \\ 0 & c & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & c & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & c & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & c & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & c & 0 & 0 \end{bmatrix} \quad \text{where } c = \frac{4\pi^{5/2}}{\sqrt{3}}.$$

In addition, her Stokes matrix corresponds to the inverse of ours. Taking this into account, our evaluation of the Stokes multipliers agrees with the formulae listed in [47, Example 2.43]$^{(9)}$, except that we find $\eta = 16$ instead of $\eta = -16$.

<footnote>$^{(9)}$ An error appears to have slipped in the specialisation of these formulae for the parameters (7.2) in [47, Example 2.45].</footnote>
8. Representation of complex numbers

**Computable complex numbers and intervals.** The Stokes multipliers are in general transcendental numbers with no known 'exact' representation simpler than their definition. Our algorithm must be able to compute these numbers within an error of no more than some given $\varepsilon$, and, by necessity, does so by operating on complex numbers that are known approximately only. This means that some mechanism is needed for keeping track of the accumulation of approximation errors along the course of the computation.

Perhaps the first such mechanism that comes to mind is to determine ahead of time how accurately each intermediate operation needs to be carried out in order for the final result to satisfy the error tolerance. This strategy is extremely tedious and error-prone, and we do not attempt to give any bounds of this type.

For theoretical purposes, however, we can mostly ignore the issue if we work with computable complex numbers, in the sense of the following classical definition.

**DEFINITION 8.1.** — A complex number $x$ is computable if there exists an algorithm that takes as input a tolerance $\varepsilon > 0$ and returns an approximation $\tilde{x} \in \mathbb{Q}[i]$ of $x$ satisfying $|\tilde{x} - x| \le \varepsilon$. We denote by $C_\mathcal{O}$ the set of computable complex numbers.

If $x$ and $y$ are computable complex numbers, then $x + y$ and $xy$ are clearly computable; and, if $x$ is a nonzero computable number, its inverse $1/x$ is computable. Similarly, the standard functions such as $\log x$ and $1/\Gamma(x)$ that we use in what follows all map computable numbers from some open domain to computable numbers. However, there is no algorithm to test if two computable numbers are equal.

Formally, the algorithms presented in this article output matrices of computable complex numbers. The fact that the entries of the matrices are computable means that it is possible, given any $\varepsilon > 0$, to obtain approximations with an error guaranteed not to exceed $\varepsilon$.

In principle, a computable complex number can be implemented as an object that contains an approximation of a number $x \in \mathbb{C}$ but also keeps track of the definition of $x$ in terms of other computable numbers (e.g., [25, 49]). When asked for an $\varepsilon$-approximation of $x$, the object refines the current approximation by recursively recomputing the quantities on which $x$ depends to higher and higher precision until the error tolerance is met. In this model, once the basic functions on which our algorithm relies are available as functions on $C_\mathcal{O}$, all error propagation required for the whole algorithm
to produce rigorous arbitrary-precision approximations is automatic. No explicit error tolerances appear as input of the algorithms below since the output is, effectively, a program taking as input a tolerance and returning an approximation of the mathematical result.

This model is occasionally used in practice but comes with substantial computational overhead. As already illustrated in Example 7.1, in our implementation, we replace the use of computable complex numbers by interval arithmetic.

**DEFINITION 8.2.** — We denote by $\mathbb{C}_\bullet$ the set of complex intervals $[a,b]+[c,d]i$ where the endpoints $a,b,c,d$ are in a fixed dense set of rational numbers. Arithmetic operations and standard elementary and special functions $f: \mathbb{C}^s \to \mathbb{C}$ are extended to partial functions $\mathbb{C}_\bullet^s \to \mathbb{C}_\bullet$ in such a way that, for $\mathbf{x}_1, \dots, \mathbf{x}_s \in \mathbb{C}_\bullet$, the value $f(\mathbf{x}_1, \dots, \mathbf{x}_s)$, if defined, satisfies

$$f(\mathbf{x}_1, \dots, \mathbf{x}_s) \supseteq \{f(x_1, \dots, x_s) : x_1 \in \mathbf{x}_1, \dots, x_s \in \mathbf{x}_s\}$$

and $f(\mathbf{x}_1, \dots, \mathbf{x}_s)$ is undefined if there is any $(x_1, \dots, x_s) \in \mathbf{x}_1 \times \dots \times \mathbf{x}_s$ where $f(x_1, \dots, x_s)$ is undefined.

For instance, $\pi$ might be represented by $[3.14, 3.15] + [-0.1, 0.1]i$ and $[3.14, 3.15] + [0.1, 0.2]$ might be defined to be $[3.2, 3.4]$. In addition to the parameters of the version working over $\mathbb{C}_\circ$, the interval versions take as input a working precision $p \in \mathbb{N}$ specifying the number of significant digits to be used during the computation$^{(10)}$. All basic operations on computable complex numbers are replaced by their interval counterpart. The working precision determines the accuracy of the output of basic operations when it is not already limited by that of the input: for example, while $1/[0.9, 1.0]$ must be an interval of width at least 0.1 regardless of the working precision, $1/[3, 3]$ might be defined to be $[0.3, 0.4]$ at precision 1 and $[0.33, 0.34]$ at precision 2. Note that interval operations may be undefined if their input contains a point where the function is not defined, or simply if the input interval is too wide. In this case, the whole algorithm raises an error.

In general, the output of an interval operation is only guaranteed to contain the corresponding exact result, not to satisfy any preset error tolerance. However, we are only using intervals as approximations of computable complex numbers. This means that we can assume the basic interval operations to be implemented in such a way that none of the interval operations fails when the working precision is large enough, and that each output tends to a point as the working precision tends to infinity. As illustrated in Example 7.1, to compute the Stokes constants with an error bounded by $\epsilon > 0$,

<footnote>$^{(10)}$ In Example 7.1, the working precision is chosen automatically based on the value of the error tolerance parameter.</footnote>
the user<sup>(11)</sup> of our implementation can first run the whole algorithm with a working precision *p* such that $2^{-p} \approx \varepsilon$, and repeat while the algorithm fails or the accuracy of the result exceeds the tolerance.

**Exact complex numbers.** Not all operations we need can be performed using computable complex numbers. Indeed, finding the Stokes values lying on a given ray or deciding if two local exponents differ by an integer requires *exact* equality tests. These tests operate on quantities derived from the coefficients of the operator *D* by algebraic operations.

For this reason, we assume from now on that *D* has coefficients in $K(x)$ for some subfield $K \subset \mathbb{C}$ where the required operations can be carried out. Specifically, we assume that $K$ is algebraically closed and closed under complex conjugation, and that we have a way of representing its elements with a finite amount of data. We also need algorithms operating on this representation for

* testing if two given elements of $K$ are equal,

* performing the field operations $+, -, \times, /$ in $K$,

* computing the complex conjugate of an element of $K$,

* computing the roots of univariate polynomials with coefficients in $K$,

* given $a \in K$ and $\varepsilon > 0$, computing an approximation $\tilde{a} \in \mathbb{Q}[i]$ of $a$ such that $|\tilde{a} - a| \le \varepsilon$.

These assumptions imply that we can compute $\text{Re}(a)$, $\text{Im}(a)$, $|a|$ given $a \in K$, and that we can decide whether $a < b$ for given $a, b \in K \cap \mathbb{R}$. The last one means that $K \subseteq \mathbb{C}_5$.

The prototypal example of a field with these properties is $K = \overline{\mathbb{Q}}$ with elements represented as roots of univariate polynomials with rational coefficients together with isolating intervals. Johansson [32] describes a practical system for working with such ‘exact’ subfields of $\mathbb{C}$ covering both $\overline{\mathbb{Q}}$ and, under suitable number-theoretic assumptions, extensions by transcendental constants. Our implementation is currently limited to the algebraic case.

<footnote><sup>(11)</sup> It tends to be useful to leave it to the user to repeat the computation until the tolerance is met rather than doing it automatically, since, in many applications, one is actually interested in an output that is accurate enough for checking some property or performing some subsequent computation, with no easy way to tell in advance exactly how accurate it needs to be. See also the Arb documentation [30] for more on this point.</footnote>
9. Local bases of solutions

The symbolic part of our algorithm operates on solutions of D and their Borel transforms. To represent these objects using a finite amount of data, we encode them by coordinates in fixed ‘local’ bases of (formal or analytic) solutions of the operators $D_{[\alpha]}$ and $\Delta_{[\alpha]}$.

**Formal solutions in the Laplace plane.** For each Stokes value $\alpha$ of $D$, we denote by $\text{FSol}_0(D_{[\alpha]})$ the space of formal solutions free of exponentials of $D_{[\alpha]}$. In other words, $\text{FSol}_0(D_{[\alpha_i]})$ is the vector space generated by the block $(\tilde{h}_j)_{j \in \mathcal{I}_i}$ of the regular part $\tilde{h}$ of the prepared formal fundamental solution introduced in Sec. 1.3, with $\mathcal{I}_i$ as in (1.7).

We denote by $\mathcal{H}_{[\alpha]}$ a fixed basis of $\text{FSol}_0(D_{[\alpha]})$. For most of the discussion, any basis arranged by exponent modulo $\mathbb{Z}$ will do, including $(\tilde{h}_j)_{j \in \mathcal{I}_i}$ once normalized as explained in Sec. 1.3. However, for compatibility with preexisting code, our implementation uses a specific choice of $\mathcal{H}_{[\alpha]}$ that deviates slightly from the form considered in that section$^{(12)}$. Some subgroups below are stated to be applicable to this basis, or, more generally, to any basis $\mathcal{H}_{[\alpha]}$ that is *echelonized* in the following sense, and require minor adaptations otherwise.

To obtain an echelonized basis, we start from any basis of $\text{FSol}_0(D_{[\alpha]})$ consisting of series of the form

$$x^{\mathcal{L}} \sum_{r=0}^{s-1} \sum_{m \ge 0} c_{r,m} x^m \ln(x)^r \quad (9.1)$$

where $\mathcal{L}$ is one of the exponents of $D$ attached to the Stokes value $\alpha$ and $s$ is the dimension of the corresponding Jordan block in $\mathcal{L}$. For each class $\mathcal{L} + \mathbb{Z}$ of exponents differing by integers, we write the corresponding basis elements with a common $\mathcal{L}$ of maximal real part (the ‘group leader’ of Definition 9.1 below), and enumerate the coefficients of each of these expansions in the order

$$c_{s-1,0}, c_{s-2,0}, \dots, c_{0,0}, \quad c_{s-1,1}, c_{s-2,1}, \dots, c_{0,1}, \quad \dots \quad (9.2)$$

Then we put these resulting family of sequences in reduced echelon form. One can show (e.g., [52, Chapter V]) that the pivots that appear are the coefficients $c_{\rho,\mu}$ where $(\mathcal{L}, \mu, \rho)$ belongs to the set $\mathcal{E}(D_{[\alpha]})$ defined as follows.

<footnote>$^{(12)}$Specifically, the basis it works with is compatible with the block and sub-block structure of Sec. 1.3, but it cannot necessarily be factored as the product of a matrix of formal power series by $x^{\mathcal{L}}$ with $\mathcal{L}$ in Jordan form. The echelon form property that we require simplifies some aspects of the algorithms but the loss of the factorization complicates others, and it is not clear to us at the moment what choice would be best with no consideration for backward compatibility.</footnote>
**DEFINITION 9.1.** — Given any differential operator $D$, we say that an exponent $\mathcal{L}$ of $D$ at the origin is a group leader when none of the numbers $\mathcal{L}-1, \mathcal{L}-2, \dots$ is also an exponent. We denote by $\mathcal{E}(D)$ the set of triples $(\mathcal{L}, \mu, \rho) \in \mathbb{K} \times \mathbb{N} \times \mathbb{N}$ such that $\mathcal{L}$ is a group leader and $\mathcal{L}+\mu$ is an exponent of multiplicity strictly larger than $\rho$, and we call the set $\mathcal{E}(D)$ the structure of the local solutions free of exponentials of $D$.

An echelonized basis is naturally indexed by $\mathcal{E}(D_{[\alpha]})$ in the sense that the coefficient of $x^{\mathcal{L}+m} \ln(x)^r$ in a given basis element is 1 for exactly one triple $(\mathcal{L}, m, r) \in \mathcal{E}(D_{[\alpha]})$ and 0 for all others. We will make free use of this correspondence and refer, for instance, to rows and columns of matrices using indices taken from $\mathcal{E}(D_{[\alpha]})$.

We assume from now on that the $\mathcal{H}_{[\alpha]}$ are echelonized, and we denote by

$$ \mathcal{Y}(x) = \left( e^{-\alpha_1/x} \mathcal{H}_{[\alpha_1]}, \dots, e^{-\alpha_N/x} \mathcal{H}_{[\alpha_N]} \right) = (y_1(x), \dots, y_n(x)) $$

the formal fundamental solution obtained by collecting the families $\mathcal{H}_{[\alpha]}$ for all Stokes values $\alpha$.

**Remark 9.2.** — When the coefficients $c_{j,m}$ of a solution (9.1) are enumerated in the order (9.2), each coefficient can be computed from the previous ones using simple recurrence relations deduced from the differential operator [52, Chapter V]. This means that the series expansions of the solutions $y_n$ are easily computed to any desired order.

**Remark 9.3.** — The conditions above determine $\mathcal{Y}$ up to permutations of the Stokes values and classes of exponents. In our implementation, the families $\mathcal{H}_{[\alpha]}$ are always sorted in such a way that whenever two solutions are asymptotically comparable, the dominant one as $x \to 0$ comes first. To be completely specific, we sort the solutions first by increasing $\mathrm{Re}(\mathcal{L})$, then, in case of ties, by decreasing degree with respect to $\ln(\zeta)$ of the coefficient of $\zeta^\mathcal{L}$, then by decreasing absolute value of $\mathrm{Im}(\mathcal{L})$ (so that purely real exponents come last for a given $\mathrm{Re}(\mathcal{L})$), and finally by increasing $\mathrm{Im}(\mathcal{L})$.

Similarly, we sort the Stokes values firstly by increasing real part, so that the ordering of $y_1, \dots, y_n$ as a whole also reflects the asymptotic dominance relation as $x \to 0$ with $x > 0$. When several Stokes values have the same real part, we order them by decreasing absolute value of the imaginary part, so that, for each real part, purely real values come last. Finally, complex conjugates are sorted by increasing imaginary part. In particular, our Stokes matrices in the direction $\omega = 0$ are lower triangular.
*Example 9.4 (Continued from Example 7.1).* — Denoting $\theta = e^{i\pi/3}$, the formal fundamental solution of the operator $D$ from Sec. 7 used in our implementation is $\mathcal{Y} = (y_1, \dots, y_7)$ where

$$y_1(x) = e^{+6/x} x^{-2} \left( 1 + \frac{2}{9}x + \frac{49}{648}x^2 + \frac{67}{17496}x^3 + \dots \right)$$

$$y_2(x) = e^{-6\theta^{-2}/x} x^{-2} \left( 1 + \frac{\theta^{-2}+6}{27}x + \frac{49}{3888}\theta^{-2}x^2 + \frac{67}{17496}x^3 + \dots \right)$$

$$y_3(x) = e^{-6\theta^2/x} x^{-2} \left( 1 + \frac{\theta^2+6}{27}x + \frac{49}{3888}\theta^2 x^2 + \frac{67}{17496}x^3 + \dots \right)$$

$$y_4(x) = x^3 + \frac{315}{128}x^6 + \dots$$

$$y_5(x) = e^{-6\theta^{-1}/x} x^{-2} \left( 1 + \frac{\theta^{-1}-6}{27}x + \frac{49}{3888}\theta^{-2}x^2 - \frac{67}{17496}x^3 + \dots \right)$$

$$y_6(x) = e^{-6\theta/x} x^{-2} \left( 1 + \frac{\theta-6}{27}x + \frac{49}{3888}\theta x^2 - \frac{67}{17496}x^3 + \dots \right)$$

$$y_7(x) = e^{-6/x} x^{-2} \left( 1 - \frac{2}{9}x + \frac{49}{648}x^2 - \frac{67}{17496}x^3 + \dots \right).$$

Each of the spaces $\text{FSol}_0(D_{[\alpha]})$ is one-dimensional, with $\text{FSol}_0(D_{[0]}) = \mathbb{C}y_4$ and $\text{FSol}_0(D_{[6]}) = \mathbb{C}e^{6/x}y_7$ for instance.

The Borel transform of $D$ is equal to

$$\Delta = \Delta_{[0]} = (\zeta^7 - 46656\zeta) \frac{d^6}{d\zeta^6} + (51\zeta^6 - 139968) \frac{d^5}{d\zeta^5} + 958\zeta^5 \frac{d^4}{d\zeta^4} \\ + 8332\zeta^4 \frac{d^3}{d\zeta^3} + 34521\zeta^3 \frac{d^2}{d\zeta^2} + 62289\zeta^2 \frac{d}{d\zeta} + 36015\zeta.$$

Its indicial polynomial is $\tilde{\Pi}_{[0]}(\lambda) = [\lambda]_{5-}(\lambda-2)$, and the corresponding basis structure is

$$\mathcal{E}(\Delta_{[0]}) = \{(0, 0, 0), (0, 1, 0), (0, 2, 1), (0, 2, 0), (0, 3, 0), (0, 4, 0)\}.$$

The basis $\hat{\mathcal{Y}}_0$ of the space $\text{Sol}(\Delta_{[0]})$, again with the conventions of our implementation, takes the form

$$\begin{aligned} \hat{\mathcal{Y}}_0(\zeta) = & \left( 1 + \frac{2401}{8957952}\zeta^6 + \mathcal{O}(\zeta^{12}), \right. \\ & \zeta + \frac{64}{382725}\zeta^7 + \mathcal{O}(\zeta^{12}), \\ & \zeta^2 \log(\zeta) + \left( \frac{1}{8192} \log(\zeta) - \frac{19}{589824} \right) \zeta^8 + \mathcal{O}(\zeta^{12}), \\ & \zeta^2 + \frac{1}{8192}\zeta^8 + \mathcal{O}(\zeta^{12}), \\ & \zeta^3 + \frac{125}{1285956}\zeta^9 + \mathcal{O}(\zeta^{12}), \\ & \zeta^4 + \frac{14641}{179159040}\zeta^{10} + \mathcal{O}(\zeta^{12}) \right). \end{aligned}$$

For all other Stokes values $\alpha$, one has $\tilde{\Pi}_{[\alpha]}(\lambda) = [\lambda]_{5-}(\lambda+3)$.
**Local bases in the Borel plane.** Similarly, for any $\alpha \in \mathbb{C}$, we denote by $\text{Sol}_0(\Delta_{[\alpha]})$ the space of germs of solutions of $\Delta_{[\alpha]}$ analytic in a neighborhood of the origin slit along the negative real axis, and we fix a basis

$$ \hat{\mathcal{Y}}_{\alpha} = (\hat{y}_{\alpha,1}, \dots, \hat{y}_{\alpha,\nu}) $$

of $\text{Sol}_0(\Delta_{[\alpha]})$. We choose each basis element $\hat{y}_{\alpha,i}$ in the form

$$ \zeta^{\lambda} \sum_{r=0}^{s-1} \sum_{m \ge 0} c_{r,m} \zeta^m \ln(\zeta)^r \quad (9.3) $$

where $\lambda$ is one of the exponents of $\Delta_{[\alpha]}$ at 0 and $s \in \mathbb{N}$. As above in the case of formal solutions in the Laplace plane, we assume that within families of series (9.3) that have the same exponent modulo $\mathbb{Z}$, the sequences

$$ c_{s-1,0}, c_{s-2,0}, \dots, c_{0,0}, \quad c_{s-1,1}, c_{s-2,1}, \dots, c_{0,1}, \quad \dots \quad (9.4) $$

are in reduced echelon form. This way, the basis $\hat{\mathcal{Y}}_{\alpha}$ of $\text{Sol}_0(\Delta_{[\alpha]})$ is naturally indexed by $\mathcal{E}(\Delta_{[\alpha]})$. In the implementation, the exponent classes $\lambda + \mathbb{Z}$ are ordered as specified in Remark 9.3.

## 10. Factorization of Stokes blocks

As discussed in Sec. 3.1, p. 24, the Stokes matrix in the direction $\omega$ has a (potentially) nontrivial block associated with each pair $(\alpha, \beta)$ of Stokes values such that $\arg(\beta - \alpha) = \omega$. Algorithm 1 below computes each such block as a product of a *Borel transform matrix*, a *connection matrix*, and a *connection-to-Stokes matrix*, which we now define.

**Connection matrices.** Let $\alpha, \beta$ be two Stokes values of $D$. For any path $\gamma : [0, 1] \to \mathbb{C}$ avoiding the singular points of $\Delta$ and such that

* $\gamma(0)$ is close enough to $\alpha$, but not on the local branch cut $\alpha + \mathbb{R}_{\le 0}$,

* $\gamma(1)$ is close enough to $\beta$, but not on the local branch cut $\beta + \mathbb{R}_{\le 0}$,

analytic continuation along $\gamma$ provides an isomorphism from $\text{Sol}_0(\Delta_{[\alpha]})$ to $\text{Sol}_0(\Delta_{[\beta]})$. This isomorphism depends on the choice of $\gamma$. We make the following specific choice: given $\eta > 0$, we consider a straight-line path $\gamma$ such that (see also Figure 13.1, p. 58)

* $\gamma(0) = \alpha + (\eta - i\eta^2)(\beta - \alpha)$, i.e., $\gamma$ starts close to $\alpha$, slightly to the right(13) of the ray pointing to $\beta$,

<footnote>(13) The choice of the perturbation factor $\eta - i\eta^2$ instead of the simpler $\eta - i\eta$ one might be tempted to use ensures that the perturbed point is not right on the branch cut of the local argument.</footnote>
• $\gamma(1) = \beta + (\eta - i\eta^2)(\beta - \alpha)$, i.e., $\gamma$ ends just behind $\beta$ seen from $\alpha$, slightly to the right of the same ray.

In particular, $\gamma$ passes to the right of any singular point of $\Delta$ lying on the open line segment $(\alpha, \beta)$. For small enough $\eta$, the map from $\text{Sol}_0(\Delta_{[\alpha]})$ to $\text{Sol}_0(\Delta_{[\beta]})$ obtained by analytic continuation along $\gamma$ does not depend on the choice of $\eta$.

**DEFINITION 10.1.** — We denote by $\mathbf{T}_{\alpha,\beta} \in \mathbb{C}^{\nu \times \nu}$ the matrix of this map in the bases $\hat{\mathcal{Y}}_{\alpha}$ and $\hat{\mathcal{Y}}_{\beta}$.

In other words, $\mathbf{T}_{\alpha,\beta}$ is the change-of-basis matrix from $\hat{\mathcal{Y}}_{\alpha}$ to $\hat{\mathcal{Y}}_{\beta}$, both viewed as bases of the space of solutions of $\Delta$ on a simply connected neighborhood of $\gamma$.

**Borel transform matrices.** Let $\alpha$ be a Stokes value of $D$ of multiplicity $k$. The formal Borel transformation (2.5), p. 16, extended by mapping $x^{-j}$ for nonnegative integer $j$ to zero as explained just after its definition, defines a map $\tilde{\mathcal{B}}_{\alpha} : \text{FSol}_0(D_{[\alpha]}) \to \text{Sol}_0(\Delta_{\alpha})$ between the local solution spaces defined in the previous section.

**DEFINITION 10.2.** — We denote by $\mathbf{B}_{\alpha} \in \mathbb{C}^{\nu \times k}$ the matrix of this map in the bases $\mathcal{H}_{[\alpha]}$ and $\hat{\mathcal{Y}}_{\alpha}$.

It results from the discussion in Sec. 2.3 that, when the exponents of $D_{[\alpha]}$ all belong to $\mathbb{C}\Z_{\le 0}$, one has $k \le \nu$ and the matrix $\mathbf{B}_{\alpha}$ has full column rank.

**Example 10.3 (Continued from Example 9.4).** — On the problem considered in Example 7.1, the transformed operator $\Delta$ has order $\nu = 6$ and all Stokes values have multiplicity $k = 1$. Since the formal Borel transform of $y_4(x) = x^3 + O(x^6)$ is $\zeta^2/2 + O(\zeta^5)$, the Borel transform matrix associated with the Stokes value 0 is

$$ \mathbf{B}_0 = (0, 0, 0, \frac{1}{2}, 0, 0)^{\mathrm{T}}. $$

Similarly, the Borel transform matrix for $\alpha = 6$ is

$$ \mathbf{B}_6 = \left(0, -\frac{67}{17496}, -\frac{9347}{2519424}, \frac{56135}{22674816}, -\frac{13289119}{29386561536}, -\frac{57551105}{1057916215296}\right)^{\mathrm{T}} $$

where the entries are the coefficients of $\zeta^{\lambda}$ for $\lambda = -3, 0, 1, 2, 3, 4$ (that is, for $\lambda$ ranging over the roots of $\tilde{\Pi}_{[6]}$ enumerated in increasing order) in the expansion of $\tilde{\mathcal{B}}_{\alpha}(y_7)$.

**Connection-to-Stokes matrices.** Fix a Stokes value $\beta$ of $D$ of multiplicity $k'$ and a direction $\omega \in (-\pi, \pi]$. We denote by $\zeta$ the local variable given by $\xi = \beta + \zeta$. Like in Sec. 5, let $\gamma_0$ be a Hankel-type path in the $\zeta$-plane,
going around the origin positively and with both branches going to infinity just to the left of the direction $\omega$.

Consider an arbitrary element $\hat{y}(\zeta)$ of $\text{Sol}_0(\Delta_{[\beta]})$. Initially, $\hat{y}(\zeta)$ is defined in a neighborhood of the origin of the $\zeta$-plane slit along $\mathbb{R}_{\le 0}$. We still denote by $\hat{y}$ its analytic continuation starting from small $\zeta > 0$ and along $\gamma_0$. Since, as a solution of $\Delta_{[\beta]}$, this analytic continuation has at most exponential growth $e^{a\Re\zeta}$ at infinity, we have for small enough $x > 0$

$$ \partial \int_{\gamma_0} \hat{y}(\zeta) e^{-\zeta/x} d\zeta = \int_{\gamma_0} \zeta \hat{y}(\zeta) e^{-\zeta/x} d\zeta, \quad \frac{1}{x} \int_{\gamma_0} \hat{y}(\zeta) e^{-\zeta/x} d\zeta = \int_{\gamma_0} \hat{y}'(\zeta) e^{-\zeta/x} d\zeta. $$

Therefore the function $L_{\beta, \omega} \hat{y}(x) = \int_{\gamma_0} \hat{y}(\zeta) e^{-\zeta/x} d\zeta$ is a solution of $D_{[\beta]}$ defined on the Borel disc with diameter $[0, e^{i\omega}/a]$. Taking its asymptotic expansion as $x \to 0$ with $\arg x = \omega$, we obtain an element $\tilde{L}_{\beta, \omega} \hat{y}$ of $\text{FSol}_0(D_{[\beta]})$.

While the function $L_{\beta, \omega} \hat{y}$ depends on $\omega$, rotating the path $\gamma_0$ by a small angle only changes $L_{\beta, \omega} \hat{y}$ by an exponentially small amount as $x \to 0$. It follows that for any two $\omega, \omega'$, the asymptotic expansion of $L_{\beta, \omega} \hat{y}$ as $x \to 0$ with $\arg x = \omega$ is equal to the asymptotic expansion of $L_{\beta, \omega'} \hat{y}$ as $x \to 0$ with $\arg x = \omega'$. In other words, $\tilde{L}_{\beta, \omega} \hat{y}$ is independent of $\omega$. (Another way to see this is that $\tilde{L}_{\beta, \omega} \hat{y}$ is determined by the initial terms of the local expansion of $\hat{y}$ at $\beta$ using Lemma 2.1.) Dropping the unnecessary $\omega$ from the notation, this construction defines a linear map $\tilde{L}_{\beta} : \text{Sol}_0(\Delta_{[\beta]}) \to \text{FSol}_0(D_{[\beta]}$).

**DEFINITION 10.4.** — We denote by $\mathbf{L}_{\beta} \in \mathbb{C}^{k' \times \nu}$ the matrix of the map $\tilde{L}_{\beta}$ in the bases $\hat{\mathcal{Y}}_{\beta}$ and $\mathcal{H}_{[\beta]}$.

**Example 10.5 (Continued from Example 9.4).** — On the problem considered in Example 7.1, one has

$$ \mathbf{L}_0 = (0, 0, 4\pi i, 0, 0, 0), $$

expressing that the only entry of $\hat{\mathcal{Y}}_0$ for which the integral in Definition 10.4 does not vanish is the one starting with $x^2 \log(x)$, and the value of the integral is then $4\pi i(x^3 + \dots)$. Similarly, one has $\mathbf{L}_6 = (\pi i, 0, 0, 0, 0, 0)$.

See also Example 14.2, p. 67 for a slightly more complicated example. The computation of $\mathbf{B}_{\alpha}$ and $\mathbf{L}_{\beta}$ is discussed in detail in Sec. 12.

**Factorization of Stokes blocks.** With these definitions, the block of the Stokes matrix associated with a given pair of Stokes values decomposes as follows.

**LEMMA 10.1.** — Let $\alpha, \beta$ be two distinct Stokes values of $D$. Let $\mathcal{I}, \mathcal{I}'$ be the corresponding blocks of indices in the decomposition (1.7), p. 12, of
the formal fundamental solution. Then the block $C_{\omega}^{[\mathcal{I}';\mathcal{I}]}$ of the Stokes matrix of $D$ in the direction $\omega = \arg(\beta - \alpha)$ is equal to $\mathbf{L}_{\beta}\mathbf{T}_{\alpha,\beta}\mathbf{B}_{\alpha}$.

*Proof*. — The result is essentially a rephrasing of equation (5.2), p. 28, taking into account that, in virtue of Remark 5.3, it is not necessary to assume the absence of monomials of nonpositive integer degree in the solutions. More precisely, for a given $\ell \in \mathcal{I}$, consider the element $y_\ell$ of $\mathcal{Y}$. Since this $y_\ell$ is a solution of the *unprepared* equation, it corresponds to $\tilde{z}_\ell$ in the notation of Remark 5.3. As in the remark, write $\tilde{z}(x) = \tilde{k}(x)e^{-\alpha/x}$ and let $\tilde{K}(x)$ be the series obtained from $\tilde{k}(x)$ by dropping the terms that belong to $\mathbb{C}[1/x]$.

The matrix $\mathbf{B}_\alpha$ maps the coordinates of $\tilde{k}(x)$ in the basis $\mathcal{Y}_\alpha(x)$ of $\text{FSol}_0(D_{[\alpha]})$ to those of $\tilde{K}_\ell$ in the basis $\hat{\mathcal{Y}}_\alpha(\zeta)$ of $\text{Sol}_0(\Delta_{[\alpha]})$.

The matrix $\mathbf{T}_{\alpha,\beta}$ maps these coordinates to those of $\hat{K}_{\ell}^{-}(\beta + \zeta)$ in the basis $\hat{\mathcal{Y}}_{\beta}(\zeta)$ of $\text{Sol}_0(\Delta_{[\beta]})$.

Finally, the matrix $\mathbf{L}_\beta$ maps the coordinates of $\hat{K}_\ell^-(\beta + \zeta)$ in the basis $\hat{\mathcal{Y}}_\beta(\zeta)$ of $\text{Sol}_0(\Delta_{[\beta]})$ to the coordinates in the basis $\mathcal{H}_{[\beta]}(x)$ of the expansion as $x \to 0$ ($\arg x = \underline{\omega}$) of $\int_{\gamma_0} \hat{K}^-(\alpha + \zeta)e^{-\zeta/x}d\zeta$.

By Remark 5.3, these coordinates are exactly the $C_{\underline{\omega}}^{[j;\ell]}$ for $j \in \mathcal{I}'$. $\square$

*Example 10.6 (Continued from Examples 10.3 and 10.5).* — In the Stokes matrix in the direction $\underline{\omega} = 0$ computed in Example 7.1, the block corresponding to the pair of Stokes values $(\alpha, \beta) = (0, 6)$ reduces to a single entry $c \approx -221.7025i$ appearing in the fourth column (corresponding to $y_4(x) = e^{-\alpha/x}(x^{-2} + \dots)$) and last row (corresponding to $y_7(x) = e^{-\beta/x}(x^{-2} + \dots)$). Lemma 10.1 expresses this $c$ as a product $c = \pi i \times a \times \frac{1}{2}$ where $a \approx 141.1420$ is the entry in the fourth column (corresponding to the element $\zeta^2 \log(\zeta) + \dots$ of $\hat{\mathcal{Y}}_0$) and first row (corresponding to the element $\zeta^{-3} + \dots$ of $\hat{\mathcal{Y}}_6$) of the connection matrix $\mathbf{T}_{0,6}$.

See also Example 14.1 (p. 65) for several other examples of the decomposition $\mathbf{L}_{\beta}\mathbf{T}_{\alpha,\beta}\mathbf{B}_{\alpha}$, including a case with nonpositive integer exponents.

## 11. The main algorithm

With this notation in place, we can now describe an algorithm that takes as input a differential operator and computes interval enclosures of all its Stokes matrices at the origin. The structure of the algorithm is presented in Algorithm 1 below. It first computes all Borel transform matrices $\mathbf{B}_\alpha$,
connection matrices $\mathbf{T}_{\alpha,\beta}$, and connection-to-Stokes matrices $\mathbf{L}_{\beta}$ separately,
and then fills the Stokes matrices with blocks of the form $\mathbf{L}_{\beta}\mathbf{T}_{\alpha,\beta}\mathbf{B}_{\alpha}$.

For readability, Algorithm 1 refers to directions $\omega$ as elements of $(-\pi, \pi]$. However, all directions it uses are of the form $\omega = \arg\theta$ for some $\theta \in \mathbb{K}$, so that they can be manipulated exactly by encoding them by a suitable $\theta$.

ALGORITHM 1. — Stokes matrices

Input: An operator $D = D(x^{-1}, \partial) \in \mathbb{K}[x^{-1}][\partial]$.

Output: A set of pairs $(\omega, \mathbf{C}_{\omega})$ where $\omega = \arg\theta$ for some $\theta \in \mathbb{K}$ and
$\mathbf{C}_{\omega} \in \mathbb{C}_{\zeta}^{n \times n}$.

(1) Compute the Borel transform $\Delta = D(\frac{d}{d\xi}, \xi) \in \mathbb{K}[\xi][\frac{d}{d\xi}]$ of $D$.

(2) Compute the set $\Sigma \subset \mathbb{K}$ of singular points of $\Delta$ and their multiplicities.

(3) Compute the set $\Omega = \{\arg(\beta - \alpha) : \alpha, \beta \in \Sigma, \alpha \neq \beta\}$ and, for each $\omega \in \Omega$, the set $\mathcal{A}(\omega) = \{(\alpha, \beta) : \arg(\beta - \alpha) = \omega\}$.

(4) Call Algorithm 4 to compute the matrices $\mathbf{T}_{\alpha,\beta}$ for all $(\alpha, \beta) \in \Sigma^2$.

(5) For each $\alpha \in \Sigma$:

(a) Compute $\Delta_{[\alpha]} \in \mathbb{K}[\zeta][\frac{d}{d\zeta}]$ by replacing $\xi$ by $\alpha + \zeta$ in $\Delta$.

(b) Call Algorithm 2 to compute the matrix $\mathbf{B}_{\alpha}$.

(c) Call Algorithm 3 to compute the matrix $\mathbf{L}_{\alpha}$.

(6) For each $\omega \in \Omega$:

(a) Initialize $\mathbf{C}_{\omega} \in \mathbb{C}_{\zeta}^{n \times n}$ to the zero matrix.

(b) For each $(\alpha, \beta) \in \mathcal{A}(\omega)$, set the block $\mathbf{C}_{\omega}^{[\mathcal{I}';\mathcal{I}]}$ of the matrix $\mathbf{C}_{\omega}$ to $\mathbf{L}_{\beta}\mathbf{T}_{\alpha,\beta}\mathbf{B}_{\alpha}$, where $\mathcal{I}$ and $\mathcal{I}'$ are the sets of row and column indices corresponding respectively to $\alpha$ and $\beta$ (see Sec. 1.3).

(7) Return $\{(\omega, \mathbf{C}_{\omega}) : \omega \in \Omega\}$.

Remark 11.1. — Step 3 can be implemented by first sorting pairs $(\alpha, \beta)$ according to $\omega = \arg(\beta - \alpha)$, and then grouping together pairs with the same $\omega$. In fact, in view of the computation of connection matrices (see Remark 13.7 below), it is convenient to represent each $\mathcal{A}(\omega)$ as a set of lists each containing all Stokes values lying on a certain line in the direction $\omega$. For instance, in the case of the operator considered in Sec. 7, one has for the direction $\omega = 0$ three alignments

$$[6e^{2i\pi/3}, 6e^{i\pi/3}], \ [-6, 0, 6], \ [-6e^{-2i\pi/3}, 6e^{-i\pi/3}]$$

encoding five pairs $(\alpha, \beta)$ corresponding to the five nontrivial blocks (each reduced to a single entry) of the Stokes matrix. Most of this computation can be done numerically using interval arithmetic, with (slower) exact comparisons only needed in the presence of arguments very close to each other or to $\pi$.
This algorithm correctly computes the Stokes matrices of operators of
single level one. More explicitly, this means that the following result holds.

THEOREM 11.2. — Assume that the operator $D$ given as input to Algorithm 1 is of single level one. The algorithm returns a set $\{ (\omega, C_{\omega}) : \omega \in \Omega \}$ such that

(1) $\Omega$ is the set of anti-Stokes directions of $D$;

(2) each $C_{\omega}$ is a matrix of computable numbers such that $I + C_{\omega}$ is the Stokes matrix of $D$ in the direction $\omega$, relative to the fundamental matrix $\mathcal{Y}$.

*Remark 11.3 (Interval version and error bounds).* — Using interval arithmetic instead of computable complex numbers, either the algorithm raises an error or its output satisfies item (1) above and

(2) for all $\omega \in \Omega$ and $1 \le j, l \le n$, the entry $C_{\omega}^{[j;l]} \in \mathbb{C}_{\bullet}$ of $C_{\omega}$ contains $C_{\omega}^{[j;l]}$, where $I + C_{\omega}$ is the Stokes matrix of $D$ in the direction $\omega$, relative to the fundamental matrix $\mathcal{Y}$.

Additionally, for a fixed input operator $D$, the values of the working precision for which the algorithm can fail are bounded, and the diameters of all $C_{\omega}^{[j;l]}$ tend to zero as the working precision tends to infinity.

*Proof*. — The correctness of Algorithms 2, 3, and 4 is proved respectively in Propositions 12.3, 12.6, and 13.6 below. We have seen in Sec. 3 that $\Omega$ is the set of anti-Stokes directions of $D$. Now fix $\omega \in \Omega$. As discussed in Sec. 3.1, the block $C_{\omega}^{[I';I]}$ is zero when $\alpha = \beta$ and when $\alpha \neq \beta$ with $(\alpha, \beta) \notin A(\omega)$. Finally, Lemma 10.1 states that the blocks with $(\alpha, \beta) \in A(\omega)$ are correctly computed. $\square$

Unlike the number of Borel transform matrices and connection-to-Stokes matrices, the number of connection matrices grows quadratically with the number $N$ of Stokes values. However, many of these matrices can be deduced from one another, so that the number of connection matrices that need to be computed *by numerical analytic continuation* is only $O(N)$, as discussed in Sec. 13 below. In our experience, the computation of these $O(N)$ core transition matrices is the most computationally expensive step in practice.

## 12. Borel transform and connection-to-Stokes matrices

In this section, we discuss the computation of the matrices $B_{\alpha}$ and $L_{\alpha}$ introduced in Sec. 10. These parts of the algorithm respectively implement
the formal Borel transformation of solutions of D and the application of connection-to-Stokes formulae to the analytic continuation of solutions of Δ. The main difference with the sketch of these operations in Part 1 is that we now work explicitly with coordinates in the bases $\mathcal{H}_{[\alpha]}$ and $\hat{\mathcal{Y}}_{\alpha}$.

## 12.1. Borel transform matrices

Recall that we have defined $\mathbf{B}_{\alpha}$ as the matrix of the Borel transformation, with monomials $x^{-j}$ with $j \in \mathbb{Z}_{\ge 0}$ mapped to zero, viewed as a linear map

$$ \tilde{\mathcal{B}}_{\alpha} : \mathrm{FSol}_0(D_{[\alpha]}) \to \mathrm{Sol}(\Delta_{[\alpha]}). $$

To compute this matrix, we need to determine the coordinates in the destina- tion basis $\hat{\mathcal{Y}}_{\alpha}$ of the Borel transform of each element of the source basis $\mathcal{H}_{[\alpha]}$. By our choice of $\hat{\mathcal{Y}}_{\alpha}$, these coordinates are the coefficients of the monomials $\zeta^{\lambda+\mu'} \log(\zeta)^{\rho'}$ for $(\lambda, \mu', \rho') \in \mathcal{E}(\Delta_{[\alpha]})$. We have seen in Sec. 2.1, p. 16, that the naive formal Borel transform of a monomial is given by

$$ \tilde{\mathcal{B}}(x^{\lambda+1} \ln(x)^m) = \frac{d^m}{d\lambda^m} \frac{\zeta^{\lambda}}{\Gamma(\lambda+1)}, \quad (12.1) $$

interpreted to give 0 when $\lambda \in \mathbb{Z}_{<0}$ and $m = 0$. Thus, given $y \in \mathrm{FSol}_0(D_{[\alpha]})$, the coefficient of a monomial $\zeta^{\lambda+\mu'} \log(\zeta)^{\rho'}$ in the series $\tilde{\mathcal{B}}_{\alpha}(y)$ only depends on a finite number of terms of $y$, and we can obtain the coordinates of $\tilde{\mathcal{B}}_{\alpha}(y)$ in the basis $\hat{\mathcal{Y}}_{\alpha}$ by computing those terms of $y$ and applying formula (12.1). The procedure, which systematizes the one we applied in Example 10.3, is detailed in Algorithm 2.

**ALGORITHM 2.** — Computation of $\mathbf{B}_{\alpha}$

**Input:** Two operators $D \in \mathbb{K}[x^{-1}][\partial]$ and $\Delta \in \mathbb{K}[\zeta][\frac{d}{d\zeta}]$.

**Output:** A matrix $\mathbf{B} \in \mathbb{C}_{\mathcal{O}}^{\nu \times k}$ where $\nu = \mathrm{ord}(\Delta)$ and $k = \mathrm{dim} \, \mathrm{FSol}(D)$.

(1) Compute the structure $\mathcal{E}(D)$ of the local solutions of $D$ free of ex- ponentials (cf. Definition 9.1).

(2) Compute the structure $\mathcal{E}(\Delta)$ of the local solutions of $\Delta$.

(3) Initialize $\mathbf{B} \in \mathbb{C}_{\mathcal{O}}^{\nu \times k}$ to the zero matrix.

(4) For each triple $(\mathcal{L}, \mu, \rho) \in \mathcal{E}(D)$:

(a) If $\mathcal{L} \notin \mathbb{Z}$, set $\lambda = \mathcal{L} - 1$ and $\delta = 0$. If $\mathcal{L} \in \mathbb{Z}$ and there is an integer $\lambda$ such that $(\lambda, 0, 0) \in \mathcal{E}(\Delta)$, then set $\lambda$ to the unique integer with that property and set $\delta = \lambda - \mathcal{L} + 1$. Otherwise, continue to the next value of $\mathcal{L}$.
(b) Consider the unique solution

$$y(x) = x^{\mathcal{L}} \sum_{r=0}^{s-1} \sum_{m \ge 0} c_{r,m} x^m \ln(x)^r \in \text{FSol}(D)$$

such that $c_{\rho,\mu} = 1$ and $c_{r,m} = 0$ for $(\mathcal{L}, m, r) \in \mathcal{E}(D) \setminus \{(\mathcal{L}, \mu, \rho)\}$. Compute the coefficients $c_{r,\delta+\mu'} \in \mathbb{C}_\mathcal{O}$ for all $0 \le r < s$ and $\mu' \ge \mu - \delta$ such that $(\lambda, \mu', 0) \in \mathcal{E}(\Delta)$ (cf. Sec. 9).

(c) For each $(\mu', \rho')$ such that $(\lambda, \mu', \rho') \in \mathcal{E}(\Delta)$ and $\mu' \ge \mu - \delta$:
(i) Compute the coefficients $\gamma_0, \dots, \gamma_{s-1-\rho'}$ of the expansion

$$\frac{1}{\Gamma(\lambda + \mu' + 1 + X)} = \gamma_0 + \gamma_1 X + \dots + \gamma_{s-1-\rho'} X^{s-1-\rho'} + O(X^{s-\rho'}). \quad (12.2)$$

(ii) Set the entry of $\mathbf{B}$ at row $(\lambda, \mu', \rho')$ and column $(\mathcal{L}, \mu, \rho)$ to

$$\sum_{r=\rho'}^{s-1} \frac{r!}{\rho'!} \gamma_{r-\rho'} c_{r,\delta+\mu'}.$$

(5) Return $\mathbf{B}$.

**Remark 12.1.** — (1) Steps 1, 2, and 4a involve exact computations on elements of $\mathbb{K}$. We present the remaining steps as computations in $\mathbb{C}_\mathcal{O}$ because they can be performed numerically. However, the entries of $\mathbf{B}_\alpha$ belong to the field extension $\mathbb{K}(\mathbf{\Gamma})$ where

$$\mathbf{\Gamma} = \{(1/\Gamma)^{(i)}(a) : a \in \mathbb{K}, i \in \mathbb{N}\}$$

and could also be computed exactly in terms of the elements of $\mathbf{\Gamma}$.

(2) When $\lambda + \mu'$ is a negative integer and $\rho' = 0$, one has $\gamma_0 = 0$ in (12.2). If, additionally, $\delta + \mu'$ is equal to $\mu$, then by definition of the series $y(x)$ one has $c_{r,\delta+\mu'} = 0$ for all $r > \rho'$, and hence the corresponding iteration of the loop (4c) can be skipped.

**Remark 12.2.** — As presented here, the algorithm contains several easy-to-avoid redundancies.

(1) For simplicity, we repeat the computation for each pair $(\mu', \rho')$. When $\lambda + 1 + \mu'$ is a multiple exponent of $\Delta$, though, one can reuse for all pairs $(\mu', \rho')$ with $\rho' > 0$ the coefficients $\gamma_i$ computed for the pair $(\mu', 0)$.

(2) The series expansion of $\Gamma(\lambda + \mu' + 1 + X)$ for $\mu' \ne 0$ can be deduced from that of $\Gamma(\lambda + 1 + X)$ using the functional equation for the gamma function.

**PROPOSITION 12.3.** — *Algorithm 2, called with $D$ set to $D_{[\alpha]}$ and $\Delta$ set to $\Delta_{[\alpha]}$, computes the Borel transformation matrix $\mathbf{B}_\alpha$ defined in Sec. 10.*
*Proof.* — As noted above, computing $\mathbf{B}_\alpha$ is equivalent to computing the coefficients of $\zeta^{\lambda+\mu'} \log(\zeta)^{\rho'}$ in $\tilde{\mathcal{B}}_\alpha(y)$ for each $y \in \mathcal{Y}_\alpha$ and $(\lambda, \mu', \rho') \in \mathcal{E}(\Delta_{[\alpha]})$. Each iteration of the main loop deals with an element of $\mathcal{Y}_\alpha$ and fills the corresponding column of $\mathbf{B}_\alpha$. Since the indicial polynomial of $\Delta_{[\alpha]}$ has degree $\nu$, the set $\mathcal{E}(\Delta_{[\alpha]})$ has cardinality $\nu$, leading to a $\nu \times k$ matrix.

Fix $(\mathcal{L}, \mu, \rho) \in \mathcal{E}(D_{[\alpha]})$, and let $y \in \mathcal{Y}_\alpha$ be the corresponding basis element. Due to the choice of $\mathcal{Y}_\alpha$, the series $y$ belongs to $x^{\mathcal{L}}\mathbb{C}[x][\ln x]$. Its image $\tilde{\mathcal{B}}_\alpha(y)$ is defined as the sum of the Borel transforms of the terms of $y$ (minus the polar and constant parts). It follows from (12.1) that $\mathcal{B}(x^{\lambda+1} \ln(x)^m)$ is of the form $\zeta^\lambda P(\ln \zeta)$ for some polynomial $P$ of degree at most $m$, and hence $\tilde{\mathcal{B}}_\alpha(y)$ lies in $\zeta^{\mathcal{L}+\mu-1}\mathbb{C}[[\zeta)][\ln \zeta]$. In particular, the only triples $(\lambda, \mu', \rho')$ in $\mathcal{E}(\Delta_{[\alpha]})$ that can correspond to a nonzero coordinate are those with $\lambda - \mathcal{L} \in \mathbb{Z}$. There is at most one$^{(14)}$ $\lambda$ with this property. Step 4a of the algorithm computes this $\lambda$, using the fact that, when $\lambda \notin \mathbb{Z}$, one has $\lambda = \mathcal{L} - 1$ by equation (2.8) or (2.10). After Step 4a, one has $\lambda + 1 = \mathcal{L} + \delta$.

The consequence of equation (12.1) noted above also implies that the coefficient of $\zeta^{\lambda+\mu'} \log(\zeta)^{\rho'}$ in $\tilde{\mathcal{B}}_\alpha(y)$ only depends on the coefficients of $x^{\mathcal{L}+m} \ln(x)^\rho$ in $y$ where $\mathcal{L} + m = \lambda + \mu' + 1$ and $\rho \ge \rho'$. The coordinates of $\tilde{\mathcal{B}}_\alpha(y)$ in the basis $\hat{\mathcal{Y}}_\alpha$ are thus determined by the coefficients of $x^{\mathcal{L}+\delta+\mu'} \ln(x)^r$ in $y$ for $(\lambda, \mu', \rho') \in \mathcal{E}(\Delta_{[\alpha]})$ and $r \ge \rho'$ such that $\delta + \mu' \ge \mu$. These coefficients are among the ones computed at step 4b, since, if $(\lambda, \mu', \rho')$ belongs to $\mathcal{E}(\Delta_{[\alpha]})$, then $(\lambda, \mu', 0)$ too.

Then, the inner loop iterates over the potentially nonzero coordinates of $\tilde{\mathcal{B}}_\alpha(y)$. By (12.1), the contribution of $\mathcal{B}(x^{\mathcal{L}+\delta+\mu'} \ln(x)^r)$ to the coordinate of index $(\lambda, \mu', \rho') \in \mathcal{E}(\Delta_{[\alpha]})$ is the coefficient of $\zeta^{\lambda+\mu'} \ln(\zeta)^{\rho'}$ in

$$ \begin{aligned} \frac{d^r}{d\lambda^r} \frac{\zeta^{\lambda+\mu'}}{\Gamma(\lambda+\mu'+1)} &= \sum_{i=0}^{r} \binom{r}{i} \zeta^{\lambda+\mu'} \ln(\zeta)^i \frac{d^{r-i}}{d\lambda^{r-i}} \frac{1}{\Gamma(\lambda+\mu'+1)} \\ &= \zeta^{\lambda+\mu'} \sum_{i=0}^{r} \frac{r!}{i!} \gamma_{r-i} \ln(\zeta)^i. \end{aligned} $$

The body of the loop computes the sum of these contributions for all $r \ge \rho'$ and stores the result in the corresponding entry of the matrix. $\square$

<footnote>$^{(14)}$ We note in passing that the case where $\mathcal{L} \in \mathbb{Z}$ but there is no suitable $\lambda \in \mathbb{Z}$, where the algorithm directly continues to the next iteration, may only occur when $k > \nu$ and the element $y(x)$ of $\mathcal{Y}_\alpha$ of index $(\mathcal{L}, \mu, \rho)$ considered at step 4b is a trivial solution $x^{-j}$, $j \in \mathbb{Z}_{\ge 0}$. Even when $k > \nu$ and $y(x)$ is a trivial solution, though, it may also happen that $\Delta$ has an integer exponent, in which case the fact that $y(x)$ does not contribute to the Stokes matrix is only detected later in the loop body.</footnote>
12.2. Connection-to-Stokes matrices

The computation of the Connection-to-Stokes matrices is similar, with
the roles of $D_{[\beta]}$ and $\Delta_{[\beta]}$ interchanged, and Lemma 5.1 providing the coef-
ficients. We have seen a simple example in Example 10.5.

ALGORITHM 3. — Computation of L_β

**Input:** Two operators $\Delta \in \mathbb{K}[\zeta][\frac{d}{d\zeta}]$ and $D \in \mathbb{K}[x^{-1}][\tilde{d}]$.

**Output:** A matrix $\mathbf{L} \in \mathbb{C}_{\zeta}^{k \times \nu}$ where $k = \dim \mathrm{FSol}(D)$ and $\nu = \mathrm{ord}(\Delta)$.

(1) Compute the structure $\mathcal{E}(D)$ of the local solutions of $D$ free of ex-
ponentials (cf. Sec. 9.1).

(2) Compute the structure $\mathcal{E}(\Delta)$ of the local solutions of $\Delta$.

(3) Initialize $\mathbf{L} \in \mathbb{C}_{\zeta}^{k \times \nu}$ to the zero matrix.

(4) For each triple $(\lambda, \mu, \rho) \in \mathcal{E}(\Delta)$:

(a) If $\lambda \notin \mathbb{Z}$, then set $\mathcal{L} = \lambda + 1$ and $\delta = 0$. If $\lambda \in \mathbb{Z}$ and there is an integer $\mathcal{L}$ such that $(\mathcal{L}, 0, 0) \in \mathcal{E}(D)$, then set $\mathcal{L}$ to the unique integer with that property and set $\delta = \mathcal{L} - \lambda - 1$. Otherwise, continue to the next value of $\lambda$.

(b) Consider the unique solution

$$
\hat{y}(\zeta) = \zeta^{\lambda} \sum_{r=0}^{s-1} \sum_{m \ge 0} c_{r,m} \zeta^m \ln(\zeta)^r \in \mathrm{Sol}_0(\Delta)
$$

such that $c_{\rho,\mu} = 1$ and $c_{r,m} = 0$ for $(\lambda, m, r) \in \mathcal{E}(\Delta) \setminus \{(\lambda, \mu, \rho)\}$.
Compute the coefficients $c_{r,\delta+\mu'}$, $0 \le r < s$, for all $\mu' \ge \mu - \delta$
such that $(\mathcal{L}, \mu', 0) \in \mathcal{E}(D)$ (cf. Sec. 9).

(c) For each $(\mu', \rho')$ such that $(\mathcal{L}, \mu', \rho') \in \mathcal{E}(D)$ and $\mu' \ge \mu$:

(i) Compute the coefficients $a_0, \dots, a_{s-1-\rho'}$ of the expansion

$$
\frac{e^{-i\pi(\mathcal{L}+\mu'+X)}}{\Gamma(1-\mathcal{L}-\mu'-X)} = a_0 + a_1 X + \cdots + a_{s-1-\rho'} X^{s-1-\rho'} + O(X^{s-\rho'}). \quad (12.3)
$$

(ii) Set the entry of $\mathbf{L}$ at row $(\mathcal{L}+\mu', \rho')$ and column $(\lambda+\mu, \rho)$
to

$$
2\pi i \sum_{r=\rho'}^{s-1} \frac{r!}{\rho'!} a_{r-\rho'} c_{r,\delta+\mu'}.
$$

(5) Return L.

Remark 12.4. — Step 4(c)i may deserve some elaboration. In the for-
mula (12.3), both $e^{-i\pi(1+\lambda+X)}$ and $1/\Gamma(-\lambda-X)$ are analytic functions of $X$
for every $\lambda$. The required series expansion can be computed by comput-
ing the expansions to order $s-1-\rho'$ of these two functions (as truncated
series with coefficients in $\mathbb{C}_{\zeta}$) and multiplying them together. We refer to
Johansson [33] for a thorough discussion of the rigorous arbitrary-precision computation of the reciprocal gamma function and its derivatives.

As pointed out to us by F. Johansson, our approach effectively reduces the computation of Hankel integrals of solutions of linear differential equations to evaluations of the gamma function, but some methods used to evaluate the gamma function are in turn based on the numerical computation of Hankel integrals of solutions of differential equations (e.g., [55])—and others on that of connection constants [33, Sec. 6]. This may hint at a more direct approach to the approximation of Stokes multipliers that would be at least as efficient as the present one. As it is, though, state-of-the-art algorithms for the direct evaluation of the Hankel integrals we are computing have slightly worse complexity than the algorithm going through the gamma function when all parameters except the precision are fixed (see Remark 13.8 below). From a practical perspective, being able to take advantage of existing well-optimised implementations of $\Gamma$ and related functions is highly beneficial.

**Remark 12.5.** — Comments analogous to Remark 12.2 apply here: $\mathcal{E}(\Delta)$ can be deduced from $\mathcal{E}(D)$; the series expansion (12.3) only needs to be computed once for each $(\mathcal{L}, \mu')$, and part of that computation can be shared among all iterations corresponding to the same $\mathcal{L}$ thanks to the functional equation for the gamma function. Additionally, the value of $\Gamma(-\lambda - \mu' - X)$ does not depend on $\omega$.

**PROPOSITION 12.6.** — *Algorithm 3, called with $\Delta$ set to $\Delta_{[\beta]}$ and $D$ set to $D_{[\beta]}$, computes the connection-to-Stokes matrix $\mathbf{L}_{\beta}$ defined in Sec. 10.*

*Proof.* — The proof is similar to that of Proposition 12.3. The algorithm fills each column of $\mathbf{L}$ with the coefficients of $x^{\mathcal{L}+\mu'} \ln(x)^{\rho'}$ where $(\mathcal{L}, \mu', \rho') \in \mathcal{E}(D_{[\beta]})$ in the asymptotic expansion as $x \to 0$ in the direction $\arg x = \omega$, for a given element $\hat{y}$ of $\hat{\mathcal{Y}}_{\beta}$, of the integral

$$ \int_{\gamma_0} \hat{y}(\zeta)e^{-\zeta/x}d\zeta, $$

where $\gamma_0$ is the same Hankel type path (depending on $\omega$) as in equation (4.1).

Fix $(\lambda, \mu, \rho) \in \mathcal{E}(\Delta_{[\beta]})$ and consider the corresponding $\hat{y}$.

When $\lambda$ is an integer and $\lambda + \mu + \mathbb{N}$ contains no exponent of $\Delta_{[\beta]}$ except possibly trivial ones (cf. p. 20), the function $\hat{y}$ is analytic at 0 and the integral vanishes, so that the whole column is zero. This is detected at step 4a and in this case the loop iteration is skipped. Otherwise, step 4a sets $\delta \in \mathbb{Z}$ so that $\mathcal{L} = \lambda + 1 + \delta$ (beware that the relation between $\mathcal{L}$, $\lambda$, and $\delta$ is not the same as in the proof of Proposition 12.3). Integrating the series expansion
of $\hat{y}(\zeta)$ termwise yields the asymptotic expansion

$$ \int_{\gamma_0} \hat{y}(\zeta) e^{-\zeta/x} d\zeta \sim \sum_{r=0}^{s-1} \sum_{m \ge \mu} c_{r,m} \int_{\Gamma_{\omega}} \zeta^{\lambda+m} \ln(\zeta)^r e^{-\zeta/x} d\zeta, \quad (12.4) $$

where, by Lemma 5.1, one has

$$ \begin{aligned} \int_{\Gamma_{\omega}} \zeta^{\lambda+m} \ln(\zeta)^r e^{-\zeta/x} d\zeta &= \frac{d^r}{d\lambda^r} \left( \frac{2\pi i}{\Gamma(-\lambda-m)} e^{-i\pi(\lambda+m+1)} x^{\lambda+m+1} \right) \\ &= 2\pi i \sum_{j=0}^{r} \binom{r}{j} \left( \frac{d^{r-j}}{d\lambda^{r-j}} \frac{e^{-i\pi(\lambda+m+1)}}{\Gamma(-\lambda-m)} \right) x^{\lambda+m+1} \ln(x)^j \end{aligned} $$

for $\arg x = \omega$.

The only terms of the expansion (12.4) contributing to the coefficient of $x^{\mathcal{L}+\mu'} \ln(x)^{\rho'}$ in the result are those with $\lambda + m + 1 = \mathcal{L} + \mu'$, that is, $m = \delta + \mu'$, and $r \ge \rho'$. The corresponding coefficients $c_{r,m}$ are computed at step 4b. Then the inner loop computes the cofactors

$$ a_{r-j} = \frac{1}{(r-j)!} \frac{d^{r-j}}{d\lambda^{r-j}} \frac{e^{-i\pi(\lambda+m+1)}}{\Gamma(-\lambda-m)} $$

and collects the terms

$$ 2\pi i c_{r,m} \binom{r}{j} [(r-j)! a_{r-j}] x^{\lambda+m+1} \ln(x)^j $$

with $m = \delta + \mu'$ and $j = \rho'$. $\square$

<figure><img src="image_10.png" /></figure>

See Example 14.2 (p. 67) for an example of the computation of $\mathbf{L}_\beta$ in the presence of logarithms.

## 13. Computing all connection matrices

We turn to the computation of connection matrices, as required by step 4 of Algorithm 1. This is the main step of the algorithm from several perspectives. In particular, it is the most computationally expensive step in practice, and the only one that may introduce transcendental constants not expressible in terms of values at points in $\mathbb{K}$ of the functions $e^x$, $e^{i\pi x}$, and $(1/\Gamma)^{(i)}(x)$.

Recall that the connection matrix $\mathbf{T}_{\alpha,\beta}$ is defined as the matrix of the analytic continuation map from $\alpha$ to $\beta$, along a path $\gamma$ specified in Sec. 10, of solutions of the operator $\Delta$, expressed in the fixed bases of solutions $\hat{Y}_\alpha$ and $\hat{Y}_\beta$. In principle, $\mathbf{T}_{\alpha,\beta}$ can be computed using the procedure sketched in Sec. 6. This procedure can also be viewed as a way of realising a rigorous ODE solver supporting generalised initial values at regular singular points
(and any other ODE solver satisfying these requirements could also be used in its place). Numerical methods for differential equations that follow this general scheme are called Taylor methods. Efficient algorithms are available for computing the sums of series solutions of linear differential equations with rational coefficients (e.g., [5, 10, 27]), leading to Taylor methods well suited to arbitrary-precision computations. The error bounds needed to make the computation rigorous can be derived, essentially, from classical proofs of the Cauchy existence theorem for solutions of differential equations; see [46] for details and tighter bounds. Our code uses the ODE solver available as part of ore_algebra [45], which implements a method of this family.

**Algorithm 1** requires us to compute $T_{\alpha,\beta}$ for all pairs $(\alpha, \beta)$ of distinct Stokes values. Using the above method separately for each pair leads however to highly redundant computations. Indeed, given three points $p, q, r$, analytic continuation from $p$ to $r$ is the same as analytic continuation from $p$ to $q$ and then from $q$ to $r$, up to correcting factors associated with the singular points crossed when deforming the corresponding paths into one another. When the triangle $pqr$ does not contain any singular point in its interior or the interior of its edges, the correcting factors are simply local monodromy matrices around some of the vertices. We now give details of an algorithm that uses numerical ODE solving to compute the matrices $T_{\alpha,\beta}$ for $(\alpha, \beta)$ ranging over a spanning tree of $\Sigma$, and recovers all remaining connection matrices by composing known ones and applying correcting factors as appropriate.

**Local monodromy matrices.** Let $\alpha$ be a Stokes value of $D$. The monodromy map sending a local solution of $\Delta_{[\alpha]}$ at the origin to its analytic continuation along a simple positive loop around the origin is a linear endomorphism of $\text{Sol}_0(\Delta_{[\alpha]})$.

**Definition 13.1.** — The local monodromy matrix at $\alpha \in \mathbb{C}$ is the matrix $\mathbf{M}_{\alpha} \in \mathbb{C}^{\nu \times \nu}$ in the basis $\hat{\mathcal{Y}}_{\alpha}$ of the monodromy map around $\alpha$.

Since $\alpha$ is a regular singular point of $\Delta_{[\alpha]}$, the local monodromy at $\alpha$ coincides with the formal monodromy. It can be computed from the exponent structure $\mathcal{E}(\Delta_{[\alpha]})$ and selected monomials of the series $\hat{\mathcal{Y}}_{\alpha}$ by a reasoning similar to that of Sec. 12, but simpler.

**Lemma 13.1.** — The entry at row $(\lambda, \mu, \rho) \in \mathcal{E}(\Delta_{[\alpha]})$ of the column of $\mathbf{M}_{\alpha}$ associated with an element

$$ \hat{y}_{\alpha,i}(\zeta) = \zeta^{\lambda} \sum_{r=0}^{s-1} \sum_{m \ge 0} c_{r,m} \zeta^m \ln(\zeta)^r \in \hat{\mathcal{Y}}_{\alpha} \quad (13.1) $$
of the local basis at $\alpha$ is given by

$$e^{2\pi i \lambda} \sum_{\delta=0}^{s-1-\rho} c_{\rho+\delta, \mu} \binom{\rho+\delta}{\delta} (2\pi i)^{\delta}. \qquad (13.2)$$

The entries corresponding to exponents not in $\lambda + \mathbb{Z}$ are zero.

*Proof.* — The image of $\hat{y}_{\alpha,i}$ by the monodromy map is equal to

$$e^{2\pi i \lambda} \zeta^{\lambda} \sum_{r=0}^{s-1} \sum_{m \ge 0} c_{r,m} \zeta^m (\ln(\zeta) + 2\pi i)^r.$$

Equation (13.2) follows by extracting the coefficient of $\zeta^{\lambda+m} \ln(\zeta)^\rho$. $\square$

<figure><img src="image_11.png" /></figure>

**Composition rules.** Let $\Sigma$ denote the set of singular points of $\Delta$.

**DEFINITION 13.2 (Void triangle).** — We say that a triangle $pqr$ is void when its vertices do not lie on a single line and its closure does not intersect $\Sigma\setminus\{p, q, r\}$.

Given a void triangle $pqr$, we need to understand how to compute $T_{p,r}$ from $T_{p,q}$ and $T_{q,r}$. For this it is convenient to introduce the following geo-metric predicates.

**DEFINITION 13.3 (Bottom-to-top order).** — For two points $\alpha, \beta \in \mathbb{C}$, we write $\alpha < \beta$ when $\mathrm{Im}(\alpha) < \mathrm{Im}(\beta)$ or $\mathrm{Im}(\alpha) = \mathrm{Im}(\beta) \wedge \mathrm{Re}(\alpha) > \mathrm{Re}(\beta)$.

Equivalently, we have $\alpha < \beta$ when $\alpha \neq \beta$ and $0 < \arg(\beta - \alpha) \le \pi$. The definition in the case of $\mathrm{Im}(\alpha) = \mathrm{Im}(\beta)$ is motivated by the fact that, with the convention $-\pi < \arg(\xi - \alpha) \le \pi$, a point $\beta$ that lies to the left of $\alpha$ is 'above' the local branch cut. The relation $<\uparrow$ is a strict total order.

**DEFINITION 13.4 (Orientation).** — For any three points $p, q, r \in \mathbb{C}$ not all lying on a same line, set $\mathrm{orient}(p, q, r) = +1$ when $p, q, r$ are ordered counterclockwise around a point of the interior of their convex hull, that is, when $r$ lies to the left of the oriented segment $(p, q)$, and $\mathrm{orient}(p, q, r) = -1$ otherwise.

Transition matrices associated with the edges of a void triangle can be composed according to the following rule.

**LEMMA 13.2.** — Let $t = pqr$ be a void triangle. The connection matrices along the edges of $t$ satisfy

$$T_{p,r} = V_r T_{q,r} V_q T_{p,q} V_p \qquad (13.3)$$
<figure><img src="image_12.png" /><figcaption>Figure 13.1. The paths $\gamma$ (with hooks added for clarity) and $\gamma^*$</figcaption></figure>

where

$$V_p = \begin{cases} \mathbf{M}_p & \text{if } q < \uparrow p < \uparrow r \text{ and } \operatorname{orient}(p, q, r) = -1, \\ \mathbf{M}_p^{-1} & \text{if } r < \uparrow p < \uparrow q \text{ and } \operatorname{orient}(p, q, r) = +1, \\ \mathbf{I} & \text{otherwise,} \end{cases} \quad (13.4)$$

$$V_q = \begin{cases} \mathbf{M}_q^{-1} & \text{if } q = \min_{<\uparrow}(p, q, r) \\ & \text{or } (p < \uparrow q < \uparrow r \text{ or } r < \uparrow q < \uparrow p) \\ & \text{and } \operatorname{orient}(p, q, r) = +1, \\ \mathbf{I} & \text{otherwise,} \end{cases} \quad (13.5)$$

$$V_r = \begin{cases} \mathbf{M}_r & \text{if } q < \uparrow r < \uparrow p \text{ and } \operatorname{orient}(p, q, r) = -1, \\ \mathbf{M}_r^{-1} & \text{if } p < \uparrow r < \uparrow q \text{ and } \operatorname{orient}(p, q, r) = +1, \\ \mathbf{I} & \text{otherwise.} \end{cases} \quad (13.6)$$

*Proof.* — In order to obtain the composition rules (13.3)–(13.6) for the matrices $T_{\alpha,\beta}$, we first define variants $T^*_{\alpha,\beta}$ of these matrices which satisfy simpler relations, and express $T^*_{\alpha,\beta}$ in terms of $T_{\alpha,\beta}$. Then we prove analogues of (13.3)–(13.6) in the case of $T^*_{\alpha,\beta}$, and deduce (13.3)–(13.6) from these analogues.

To define $T^*_{\alpha,\beta}$, let $\alpha, \beta \in \mathbb{C}$ be two points such that $\operatorname{Im}(\alpha) \neq \operatorname{Im}(\beta)$ and the open segment $(\alpha, \beta)$ does not contain any singular point of $\Delta$. Then $T^*_{\alpha,\beta}$ is the connection matrix defined in a similar way to $T_{\alpha,\beta}$ (cf. Sec. 10) with the analytic continuation path $\gamma$ replaced by

$$\gamma^* = [(1 - \eta)\alpha + \eta\beta + i\eta^2, (1 - \eta)\beta + \eta\alpha + i\eta^2],$$

that is, by a segment contained in $(\alpha, \beta)$, shifted slightly up to avoid any ambiguity in the definition when $\operatorname{Im}(\alpha) = \operatorname{Im}(\beta)$. Observe that we have $T^*_{\beta,\alpha} = (T^*_{\alpha,\beta})^{-1}$.

Extending $\gamma^*$ by a small circular arc connecting $\gamma^*(0)$ to $\gamma(0)$ contained in the analyticity domain of the local basis $\hat{Y}_\alpha$ does not change the matrix $T^*_{\alpha,\beta}$.
Similarly, one may connect $\gamma^*(1)$ to $\gamma(1)$ without leaving the analyticity domain of $\hat{Y}_\beta$ and this does not change $\mathbf{T}_{\alpha,\beta}^*$. These extensions result in a path $\gamma_{\text{ext}}^*$ with the same endpoints as $\gamma$. When $\alpha < \uparrow \beta$, the paths $\gamma$ and $\gamma_{\text{ext}}^*$ are homotopic in $C\backslash\Sigma$ (Figure 13.1, left), so that $\mathbf{T}_{\alpha,\beta}^*$ is equal to $\mathbf{T}_{\alpha,\beta}$. In contrast, when $\beta < \uparrow \alpha$, the loop $(\gamma_{\text{ext}}^*)^{-1}\gamma$ goes around $\beta$ once in the positive direction (Figure 13.1, right), and hence one has $(\mathbf{T}_{\alpha,\beta}^*)^{-1}\mathbf{T}_{\alpha,\beta} = \mathbf{M}_\beta$. We conclude that

$$ \mathbf{T}_{\alpha,\beta}^* = \begin{cases} \mathbf{T}_{\alpha,\beta}, & \alpha < \uparrow \beta, \\ \mathbf{M}_{\beta}^{-1}\mathbf{T}_{\alpha,\beta}, & \beta < \uparrow \alpha. \end{cases} \qquad (13.7) $$

Now consider the composition of the matrices $\mathbf{T}_{p,q}^*, \mathbf{T}_{q,r}$, and $\mathbf{T}_{r,p}$. Again, we can connect the corresponding paths $\gamma^*$ by small arcs contained in the analyticity domains of local solution bases at $p, q, r$ to get a closed loop. When the vertex $v \in \{p, q, r\}$ lying vertically between the other two is located to the left of the opposite side of $t$ (i.e., when the local branch cut $v+(-\infty, 0]$ points outside of $t$), since $t \cap \Sigma \subseteq \{p, q, r\}$, the loop is contractible in $C\backslash\Sigma$. The same conclusion holds in the limiting case where the bottom edge of $t$ is horizontal, so that the middle vertex according to $<\uparrow$ is the left endpoint of that edge. Otherwise, that is, when the middle vertex in $<\uparrow$ order lies to the right of the opposite edge, including when the top edge is horizontal, the loop goes around that vertex once, positively or negatively depending on $\text{orient}(p, q, r)$. Letting $u, v, w$ denote the points $p, q, r$ ordered so that $u < \uparrow v < \uparrow w$, we thus have

$$ \mathbf{T}_{u,v}^* \mathbf{T}_{w,u}^* \mathbf{T}_{v,w}^* = \mathbf{M}_v^\tau, \quad \tau := \begin{cases} 0, & \text{orient}(u, v, w) = -1, \\ 1, & \text{orient}(u, v, w) = +1. \end{cases} $$

As $\mathbf{T}_{\beta,\alpha}^* = (\mathbf{T}_{\alpha,\beta}^*)^{-1}$, this implies

$$ \mathbf{T}_{u,v}^* = \mathbf{M}_v^\tau \mathbf{T}_{w,v}^* \mathbf{T}_{u,w}^*, \quad \mathbf{T}_{v,u}^* = \mathbf{T}_{w,u}^* \mathbf{T}_{v,w}^* \mathbf{M}_v^{-\tau}, $$

$$ \mathbf{T}_{w,u}^* = \mathbf{T}_{v,u}^* \mathbf{M}_v^\tau \mathbf{T}_{w,v}^*, \quad \mathbf{T}_{u,w}^* = \mathbf{T}_{v,w}^* \mathbf{M}_v^{-\tau} \mathbf{T}_{u,v}^*, $$

$$ \mathbf{T}_{v,w}^* = \mathbf{T}_{u,w}^* \mathbf{T}_{v,u}^* \mathbf{M}_v^\tau, \quad \mathbf{T}_{w,v}^* = \mathbf{M}_v^{-\tau} \mathbf{T}_{u,v}^* \mathbf{T}_{w,u}^*. $$

and therefore, using (13.7),

$$ \mathbf{T}_{u,v} = \mathbf{M}_v^{-1} \mathbf{T}_{w,v} \mathbf{T}_{u,w}, \quad \mathbf{M}_u^{-1} \mathbf{T}_{v,u} = \mathbf{M}_u^{-1} \mathbf{T}_{w,u} \mathbf{T}_{v,w} \mathbf{M}_v^{-\tau}, $$

$$ \mathbf{M}_u^{-1} \mathbf{T}_{w,u} = \mathbf{M}_u^{-1} \mathbf{T}_{v,u} \mathbf{M}_v^{-1} \mathbf{T}_{w,v}, \quad \mathbf{T}_{u,w} = \mathbf{M}_w^{-1} \mathbf{T}_{v,w} \mathbf{M}_v^{-\tau} \mathbf{T}_{u,v}, $$

$$ \mathbf{T}_{v,w} = \mathbf{T}_{u,w} \mathbf{M}_u^{-1} \mathbf{T}_{v,u} \mathbf{M}_v^\tau, \quad \mathbf{M}_v^{-1} \mathbf{T}_{w,v} = \mathbf{M}_v^{-\tau} \mathbf{T}_{u,v} \mathbf{M}_u^{-1} \mathbf{T}_{w,u}. $$

The expressions (13.4)–(13.6) follow by discussing according to the matching between $(u, v, w)$ and $(p, q, r)$. For example, in the case of $\mathbf{V}_p$:

* when $(p, r) = (v, w)$, that is, $q < \uparrow p < \uparrow r$, one has $\mathbf{V}_p = \mathbf{M}_v^\tau = \mathbf{M}_p^\tau$ and $\text{orient}(u, v, w) = -\text{orient}(p, q, r)$;
• when $(p, r) = (v, u)$, i.e., $r < \uparrow p < \uparrow q$, one has $\mathbf{V}_p = \mathbf{M}_p^{-\tau}$, and in this case $\mathrm{orient}(u, v, w) = \mathrm{orient}(p, q, r)$;

• in all other cases, one has $\mathbf{V}_p = \mathbf{I}$.

<figure><img src="image_13.png" /></figure>

We can deduce $\mathbf{T}_{q,p}$ from $\mathbf{T}_{p,q}$ using similar arguments.

LEMMA 13.3. — For two points $p < \uparrow q$ such that $(p,q) \cap \Sigma = \emptyset$, one has $\mathbf{T}_{q,p} = \mathbf{M}_p \mathbf{T}_{p,q}^{-1}$.

*Proof*. — The loop obtained by connecting the paths $\gamma$ associated with $\mathbf{T}_{p,q}$ and $\mathbf{T}_{q,p}$ (in this order) by small arcs around at $p$ and $q$ not passing to the left of the corresponding points goes around $p$ once, counterclockwise. $\square$

Finally, it will be useful to handle the case where $q$ lies between $p$ and $r$. The corresponding rule is a direct consequence of the definition of $\mathbf{T}_{\alpha,\beta}$.

LEMMA 13.4. — Let $pqr$ be a flat triangle with $q \in (p,r)$. Then one has $\mathbf{T}_{p,r} = \mathbf{T}_{q,r} \mathbf{T}_{p,q}$.

**Algorithm.** The above discussion leads us to Algorithm 4 for computing all connection matrices simultaneously.

ALGORITHM 4. — Connection matrices

**Input:** An operator $\Delta \in \mathbb{K}[\zeta][\frac{d}{d\zeta}]$ with regular singular points.

**Output:** A family $(\mathbf{T}_{\alpha,\beta})_{\alpha,\beta \in \Sigma}$ of matrices $\mathbf{T}_{\alpha,\beta} \in \mathbb{C}_{\mathcal{G}}^{\nu \times \nu}$, where $\Sigma$ is the set of singular points of $\Delta$ and $\nu = \mathrm{ord}(\Delta)$.

(1) Compute a Euclidean minimum spanning tree $\mathcal{S}$ of the singular points of $\Delta$.

(2) Initialize a mapping $(\mathbf{T}_{\alpha,\beta})_{\alpha,\beta \in \Sigma}$ by setting $\mathbf{T}_{\alpha,\beta} = \perp$ for $\alpha \neq \beta$ and $\mathbf{T}_{\alpha,\alpha} = \mathbf{I}$ for all $\alpha$.

(3) For each $\alpha \in \Sigma$, compute the matrix $\mathbf{M}_{\alpha}$ using (13.2).

(4) For each edge $\{\alpha, \beta\} \in \mathcal{S}$, with $\alpha < \uparrow \beta$:

(a) Compute $\mathbf{T}_{\alpha,\beta}$ by solving the equation $\Delta(\hat{y}) = 0$ numerically.

(b) Set $\mathbf{T}_{\beta,\alpha} = \mathbf{M}_{\alpha} \mathbf{T}_{\alpha,\beta}^{-1}$.

(5) Compute the set $\mathcal{V}$ of triples $\{p, q, r\} \subseteq \Sigma$ such that the triangle $pqr$ is void (Definition 13.2) and the set $\mathcal{F}$ of triples $\{p, q, r\} \subseteq \Sigma$ such that $p, q, r$ are aligned and all distinct. Set $\mathcal{T} = \mathcal{V} \cup \mathcal{F}$.

(6) While there exist $\alpha, \beta \in \Sigma$ such that $\mathbf{T}_{\alpha,\beta} = \perp$:

(a) Choose a triple $t = \{p, q, r\} \in \mathcal{T}$ such that $\mathbf{T}_{p,q} \neq \perp$ and $\mathbf{T}_{q,r} \neq \perp$.

(b) If $t \in \mathcal{F}$:

(i) If $q \in (p, r)$, set $\mathbf{T}_{p,r} = \mathbf{T}_{q,r} \mathbf{T}_{p,q}$ and $\mathbf{T}_{r,p} = \mathbf{T}_{q,p} \mathbf{T}_{r,q}$.

(ii) If $p \in (q, r)$, set $\mathbf{T}_{p,r} = \mathbf{T}_{q,r} \mathbf{T}_{q,p}^{-1}$ and $\mathbf{T}_{r,p} = \mathbf{T}_{p,q}^{-1} \mathbf{T}_{r,q}$.

(iii) If $r \in (p, q)$, set $\mathbf{T}_{r,p} = \mathbf{T}_{q,p} \mathbf{T}_{q,r}^{-1}$ and $\mathbf{T}_{p,r} = \mathbf{T}_{r,q}^{-1} \mathbf{T}_{p,q}$.
(c) Else:

(i) Compare $p, q, r$ according to $<\uparrow$.

(ii) Compute orient($p, q, r$).

(iii) Compute $T_{p,r}$ and $T_{r,p}$ using (13.3).

(d) Remove from $\mathcal{T}$ all triples $\{p, q', r\}$ such that $T_{p,q'} \neq \perp$ and $T_{q',r} \neq \perp$.

(7) Return $(T_{\alpha,\beta})_{\alpha,\beta \in \Sigma}$.

**Example 13.5 (Continued from Example 7.1).** — In the running example introduced in Sec. 7, our implementation chooses for $\mathcal{S}$ the spanning tree depicted on Figure 7.1. The six core connection matrices associated with the edges of the spanning tree oriented in the bottom-to-top order are computed by calling the ODE solver from `ore_algebra` with a tolerance slightly lower than the one given on input to `stokes_dict`. For instance, the computation of $T_{6,0}$ can be reproduced separately with the command

```
sage: diffop.borel_transform().numerical_transition_matrix()
....:
[6,0], 1e-50/26)
```

and yields a 6 × 6 matrix whose first two columns are

$$
\begin{align*}
\begin{bmatrix} 0.207 & 79 \pm 9.31 \cdot 10^{-56} \\ -0.455 & 77 \pm 3.68 \cdot 10^{-56} \end{bmatrix} + \begin{bmatrix} -0.007 & 24 \pm 8.03 \cdot 10^{-57} \\ [0.007 & 24 \pm 8.03 \cdot 10^{-57}] \end{bmatrix}i & \begin{bmatrix} 126.0 & 22 \pm 2.50 \cdot 10^{-52} \\ -200.0 & 33 \pm 3.34 \cdot 10^{-52} \end{bmatrix}i \\
\begin{bmatrix} -0.192 & 97 \pm 1.99 \cdot 10^{-55} \\ -0.018 & 49 \pm 2.45 \cdot 10^{-56} \end{bmatrix} + \begin{bmatrix} -0.018 & 49 \pm 2.45 \cdot 10^{-56} \\ [0.018 & 49 \pm 2.45 \cdot 10^{-56}] \end{bmatrix}i & \begin{bmatrix} -62.47 & 78 \pm 2.65 \cdot 10^{-53} \\ -62.47 & 78 \pm 2.65 \cdot 10^{-53} \end{bmatrix}i \\
\begin{bmatrix} 0.243 & 32 \pm 3.94 \cdot 10^{-55} \\ [0.024 & 32 \pm 3.94 \cdot 10^{-55}] \end{bmatrix} + \begin{bmatrix} 0.024 & 86 \pm 3.26 \cdot 10^{-56} \\ [0.024 & 86 \pm 3.26 \cdot 10^{-56}] \end{bmatrix}i & \begin{bmatrix} 96.43 & 32 \pm 4.14 \cdot 10^{-53} \\ [96.43 & 32 \pm 4.14 \cdot 10^{-53}] \end{bmatrix}i \\
\begin{bmatrix} 0.034 & 24 \pm 1.02 \cdot 10^{-56} \\ [0.034 & 24 \pm 1.02 \cdot 10^{-56}] \end{bmatrix} + \begin{bmatrix} \pm 3.38 \cdot 10^{-60} \\ [\pm 3.38 \cdot 10^{-60}] \end{bmatrix}i & \begin{bmatrix} 9.263 & 28 \pm 3.96 \cdot 10^{-54} \\ [9.263 & 28 \pm 3.96 \cdot 10^{-54}] \end{bmatrix}i \\
\begin{bmatrix} -0.001 & 52 \pm 4.45 \cdot 10^{-57} \\ [0.001 & 52 \pm 4.45 \cdot 10^{-57}] \end{bmatrix} + \begin{bmatrix} 4.037 & 94 \cdot 10^{-5} \pm 5.46 \cdot 10^{-60} \\ [4.037 & 94 \cdot 10^{-5} \pm 5.46 \cdot 10^{-60}] \end{bmatrix}i & \begin{bmatrix} -0.336 & 15 \pm 1.74 \cdot 10^{-55} \\ [0.336 & 15 \pm 1.74 \cdot 10^{-55}] \end{bmatrix}i
\end{align*}
$$

We note that, in this special case, we could have used the fact that the operator $\Delta$ is invariant by $\xi \mapsto e^{i\pi/3}\xi$ to obtain a full set of core connection matrices from a single one instead of computing them independently. As a similar but more widely applicable optimisation, since $\Delta$ has coefficients in $\mathbb{R}[\xi]$, we could have deduced the connection matrices associated with the edges of $\mathcal{S}$ leading to Stokes values in the open upper half plane from the others by complex conjugation. Our code currently performs neither of these reductions automatically.

The connection matrices along the edges of $\mathcal{S}$ with the reverse orientation are computed using Lemma 13.3. For instance, $T_{0,6}$ is computed as $M_6 T_{6,0}^{-1}$ where $M_6$ is identity except for the first column which reads $(1, -67i\pi/17496, -9347i\pi/2519424, \dots)^T$, and in the first row of the product we find the coefficient $a = [-141.140 \dots 09 \pm 9.36 \cdot 10^{-41}]$ which made an appearance in Example 10.6. We note in passing a significant loss of accuracy due to the inversion of $T_{6,0}$ in interval arithmetic. Similarly, $T_{-6,0}$ is computed as $M_0 T_{0,-6}^{-1}$ where $M_0$ is identity except for a subdiagonal entry equal to $2\pi i$ in the third column.
The remaining connection matrices are computed in the main loop, using Lemmas 13.2 and 13.4. Let $\theta = e^{i\pi/3}$. At first one can compute, among others,

• $\mathbf{T}_{-6,6} = \mathbf{T}_{0,6}\mathbf{T}_{-6,0}$, and $\mathbf{T}_{6,-6} = \mathbf{T}_{0,-6}\mathbf{T}_{6,0}$ by Lemma 13.4,

• $\mathbf{T}_{0,6\theta} = \mathbf{T}_{6,6\theta-1}\mathbf{T}_{0,6}$ by Lemma 13.2 in the case $r < \uparrow q < \uparrow p$ and $\mathrm{orient}(p, q, r) = -1$,

• $\mathbf{T}_{6\theta-1,0} = \mathbf{T}_{6,0}\mathbf{M}_6^{-1}\mathbf{T}_{6\theta-1,6}$ by Lemma 13.2 with $p < \uparrow q < \uparrow r$ and $\mathrm{orient}(p, q, r) = +1$,

• $\mathbf{T}_{6\theta,6\theta-1} = \mathbf{T}_{6,6\theta-1}\mathbf{T}_{6\theta,6}$ and $\mathbf{T}_{6\theta-1,6\theta} = \mathbf{T}_{6,6\theta}\mathbf{M}_6^{-1}\mathbf{T}_{6\theta-1,6}$ by the same two rules.

Once $\mathbf{T}_{6\theta-1,0}$ is known, one can use it to compute, for instance

$$ \mathbf{T}_{6\theta-1,6\theta-2} = \mathbf{M}_{6\theta-2}^{-1}\mathbf{T}_{0,6\theta-2}\mathbf{T}_{6\theta-1,0} $$

by the case $p < \uparrow r < \uparrow q$ and $\mathrm{orient}(p, q, r) = +1$ of Lemma 13.2. The algorithm continues in this way until all 42 possible connection matrices are known.

**PROPOSITION 13.6.** — Assume that the algorithm used at step 4a computes $\mathbf{T}_{\alpha,\beta} \in \mathbb{C}_{\mathcal{G}}^{\nu \times \nu}$ given $\alpha, \beta$, and $\Delta$. Then Algorithm 4, when called on the Borel transform $\Delta$ of the operator $D$, computes the family $(\mathbf{T}_{\alpha,\beta})_{\alpha,\beta \in \Sigma}$ of connection matrices defined in Sec. 10.

*Proof.* — Step 4a makes sense and can be implemented as sketched above because $\Delta$ has regular singular points. The fact that the matrices $\mathbf{T}_{\alpha,\beta}$ computed at steps 4a, 4b, 6(b)i–6(b)iii, and 6(c)iii are correct follows, respectively, from the assumption, Lemma 13.3, Lemma 13.4, and Lemma 13.2. It remains to prove that the algorithm computes $\mathbf{T}_{\alpha,\beta}$ for all $(\alpha, \beta)$.

Let $\mathcal{T}_0 = \mathcal{V} \cup \mathcal{F}$ be the initial contents of the set $\mathcal{T}$. If $|\Sigma| \le 2$, then $\mathcal{T}_0 = \emptyset$ and the algorithm terminates without going through the main loop. In this case all $\mathbf{T}_{\alpha,\beta}$ have been computed. Assume from now on that $|\Sigma| \ge 3$. To start with, observe that every time the algorithm tests the loop condition at step 6:

• for all $\alpha, \beta$, one has $\mathbf{T}_{\alpha,\beta} = \perp$ if and only if $\mathbf{T}_{\beta,\alpha} = \perp$, and the pairs $\{\alpha, \beta\}$ such that $\mathbf{T}_{\alpha,\beta} \ne \perp$ are the edges of an undirected connected graph $\mathcal{G}$;

• $\mathcal{T}$ is equal to the set of triples $\{p, q, r\} \in \mathcal{T}_0$ such that at least one of $\mathbf{T}_{p,q}, \mathbf{T}_{q,r}$, and $\mathbf{T}_{p,r}$ is $\perp$.

If the loop terminates without error, then $\mathcal{G}$ is the complete graph on $\Sigma$ and the algorithm's output is correct. Thus it suffices to prove that every time step 6a is reached, there exists a triangle $\{p, q, r\} \in \mathcal{T}_0$ with exactly two edges
in $\mathcal{G}$. Since the edges of the Euclidean spanning tree $\mathcal{S}$ do not cross, we can
complete $\mathcal{S}$ into a full triangulation$^{(15)}$ $\mathcal{C}$ of $\Sigma$ [12, Lemma 3.1.2].

If $\mathcal{C}$ contains an edge that is not yet in $\mathcal{G}$, let $\mathcal{U}$ be the graph dual to $\mathcal{C}\setminus\mathcal{S}$,
that is, the graph whose vertices are the faces of the triangulation $\mathcal{C}$ and
where two vertices are adjacent when the corresponding faces are separated
by an edge of the triangulation that is not part of the spanning tree. The
graph $\mathcal{U}$ is connected since $\mathcal{S}$ is acyclic and acyclic since $\mathcal{S}$ is connected, so
$\mathcal{U}$ is a tree. Now consider the faces of $\mathcal{C}$ in a postfix traversal of $\mathcal{U}$ rooted on
the face at infinity. The first face whose boundary is not contained in $\mathcal{G}$ is a
void triangle with exactly one edge not in $\mathcal{G}$.

If now all edges of $\mathcal{C}$ are in $\mathcal{G}$ but there exists an element of $\mathcal{V}$ with an edge
not in $\mathcal{G}$, consider a triangulation $\mathcal{C}'$ obtained from $\mathcal{C}$ by a *flip*, that is, by
replacing the common edge of two adjacent triangles by the other diagonal of
the convex quadrilateral consisting of the reunion of the two triangles. The
new faces created by the flip are elements of $\mathcal{T}_0$ with at least two edges in $\mathcal{G}$.
If the remaining edge is not in $\mathcal{G}$, we are done; otherwise, we can iterate
the argument starting from $\mathcal{C}'$. Since any element of $\mathcal{V}$ can be completed
into a triangulation and any triangulation can be connected to any other by
a sequence of flips [12, Sec. 3.4.1], we will eventually reach a triangle with
exactly two edges in $\mathcal{G}$.

Finally, if the edges of elements of $\mathcal{V}$ are all in $\mathcal{G}$, then the only $\{\alpha, \beta\} \notin \mathcal{G}$
with $\alpha \neq \beta$ in $\Sigma$ are such that the segment $[\alpha, \beta]$ contains a third point
of $\Sigma$. Indeed, consider a segment $[\alpha, \beta]$ on which there is no other point
of $\Sigma$. Either the elements of $\Sigma$ are all colinear, in which case $\{\alpha, \beta\} \in \mathcal{S} \subseteq \mathcal{G}$,
or there exists a point $\beta' \in \Sigma$ closest to the line through $\alpha$ and $\beta$ but not on
this line, and $\{\alpha, \beta, \beta'\} \in \mathcal{T}$. Then for any closest $\alpha, \beta$ with $\{\alpha, \beta\} \notin \mathcal{G}$, there
is a point $\beta' \in [\alpha, \beta]$ such that $\{\alpha, \beta, \beta'\} \in \mathcal{F}$ has exactly two edges in $\mathcal{G}$. $\square$

Remark 13.7. — (1) Any plane spanning tree of $\Sigma$ can be used as the spine $\mathcal{S}$. Taking a Euclidean minimum spanning tree helps limit the cost of the computation, but it would be even better to weight the edges by a more precise estimate of the cost of computing the corresponding $\mathbf{T}_{\alpha,\beta}$, possibly also taking into account its estimated condition number.

(2) At steps 6(b)i and 6(c)iii, we compute $\mathbf{T}_{p,r}$ and $\mathbf{T}_{r,p}$ separately
rather than deducing the latter from the former using Lemma 13.3
in order to limit the loss of numerical accuracy from repeated matrix

<footnote>(15) Here by *full triangulation* of Σ we mean a simplicial complex whose vertices are exactly the elements of Σ and whose reunion is the convex hull of Σ [12, Definition 2.2.1, 1.0.2]. For example, the only full triangulation of three distinct colinear points consists of the two segments joining the middle one to the other two, the three points themselves, and the empty face.</footnote>
inverses. Similary, adapting the algorithm used at step 4a to com-
pute $\mathbf{T}_{\alpha,\beta}^{-1}$ as the same time as $\mathbf{T}_{\alpha,\beta}$ may yield more accurate results.
Processing the triples $\{p, q, r\}$ by decreasing accuracy of ($\mathbf{T}_{p,q}, \mathbf{T}_{q,r}$)
may also help.

(3) At step 5, we can restrict ourselves to triangles that can be checked not to contain any other singular points using interval arithmetic, provided that we verify that the number of matrices computed is equal to $|\Sigma|(|\Sigma| - 1)/2$ before exiting the algorithm and restart the whole computation with a higher working precision if not.

(4) As noted in Remark 11.1, the computation of $\mathcal{F}$ (which involves exact operations in $\mathbb{K}$) can be combined with that of the sets $\mathcal{A}(\omega)$ from Algorithm 1.

(5) While we only have $O(n^2)$ matrices to compute, the cardinality of $\mathcal{V}$ can reach $\Theta(n^3)$, computing it (step 5) in the obvious way can take up to $\Theta(n^4)$ operations, and keeping $\mathcal{T}$ up to date (step 6d) up to $\Theta(n^3)$ operations in total. (It is not too hard to arrange that step 6a executes in constant time by modifying step 6d to maintain a partition of $\mathcal{T}$ according to the number of 'known' edges.) These naïve estimates can likely be improved using of geometric data structures; however, estimating the complexity of the algorithm by assigning a unit cost to operations in $\mathbb{K}$ or in $\mathbb{C}_{\mathcal{O}}$ is highly unrealistic.

*Remark 13.8.* — (Complexity with respect to the precision) When the operator *D* is fixed, the output of Algorithm 1, viewed as a family of matrices over *C*<sub>*C*</sub>, is entirely determined, but it makes sense to ask about the cost of computing 2<sup>-p</sup>-approximations of the entries of these matrices as *p* tends to infinity. In this setting, the algorithm reduces to a fixed sequence of operations on computable complex numbers. In addition, for each intermediate result *x*<sub>*i*</sub>, there exists a constant *c*<sub>*i*</sub} > 0 such that, when *p* is large enough, *x*<sub>*i*</sub> can be evaluated within an error bounded by 2<sup>-p</sup> starting from 2<sup>-p-c<sub>*i*</sub></sup>-approximations of *x*<sub>0</sub>, ..., *x*<sub>*i*−1</sub>.

Assume $\mathbb{K} = \bar{\mathbb{Q}}$, that is, the coefficients of the polynomial coefficients of the differential equation $Dy = 0$ are algebraic numbers. Step 4a of Algorithm 4 can then be performed in $O(M(p \log(p)^2))$ bit operations, where $M(p)$ denotes the complexity of multiplying $p$-bit integers [24]. The remaining steps of Algorithm 1 that depend on $p$ boil down to additions and multiplications of complex numbers, evaluations of exponentials and logarithms, and evaluations of derivatives of $1/\Gamma$ at points in $\mathbb{K}$, all of which can be performed within the same complexity bound [8]. The complexity of Algorithm 1 with respect to the target precision alone is thus $O(M(p \log(p)^2))$, or $O(p \log(p)^3)$ using the bound $M(p) = O(p \log p)$ [23]. This is better by a factor $\asymp \log(p)$
than the best published complexity bound in this case, due to van der Ho-
even [26]. However, van der Hoeven's algorithm applies to singularities of
arbitrary multilevel.

For a more general computable $\mathbb{K}$, classical algorithms for the required
basic operations [4, 9] yield a bound of $O(p^{3/2} \log(p)^{O(1)})$, not counting the
time needed to compute numerical approximations of the coefficients of $D$.
This holds both for our method and for van der Hoeven's.

14. Further examples

We conclude with a few additional examples highlighting various aspects
of the algorithm and implementation: one illustrating some degeneracies that
can happen and the role of exponents on a simple hypergeometric equation,
one where the formal solutions in the Laplace plane already involve log-
arithms, and one where the numerical analytic continuation step is more
challenging. We also provide some data on the performance of our imple-
mentation on some of these examples. This is intended to give the reader
an impression of the problems that are within reach with no serious optimi-
sation effort, not as a thorough analysis. All measurements were performed
on a laptop equipped with an Intel i7-10810U CPU, running SageMath ver-
sion 10.7 and ore_algebra revision 07a1adbb.

*Example 14.1 (Variations on the confluent hypergeometric equation D<sub>2,1</sub>).*
Let us consider again the operator (7.1) (p. 34), now in the simple case q = 2,
p = 1. In terms of x = z<sup>-1</sup> and ∂, the operator reads

$$
D = x^{-2} \partial^2 - ((\nu_1 + \nu_2 - 1)x^{-1} + x^{-2})\partial + ((\nu_1 - 1)(\nu_2 - 1) + \mu x^{-1}).
$$

The Stokes values are 0 and 1 and our standard formal fundamental solution
(Sec. 9) is given by

$$
y_1 = x^{\mu} (1 + (\mu - \nu_1 + 1)(\mu - \nu_2 + 1)x + \dots),
$$

$$
y_2 = e^{-1/x} x^{\nu_1 + \nu_2 - \mu - 1} (1 - (\mu - \nu_1)(\mu - \nu_2)x + \dots).
$$

It is known that the associated Stokes matrix in the direction $\omega = 0$ is equal
to

$$
\begin{pmatrix} 1 & 0 \\ c & 1 \end{pmatrix}, \quad c = \frac{2\pi i}{\Gamma(1 + \mu - \nu_1)\Gamma(1 + \mu - \nu_2)}.
$$

The Borel transform of $D$ for general parameters is

$$
\Delta = \xi(\xi - 1) \frac{d^2}{d\xi^2} + ((-\nu_1 - \nu_2 + 5)\xi + \mu - 2) \frac{d}{d\xi} + (\nu_1 - 2)(\nu_2 - 2), \quad (14.1)
$$

with exponents (0, μ - 1) at 0 and (0, ν₁ + ν₂ - μ - 2) at 1. When ν₁ = 1, one
can replace D by xD while preserving both the solution space and the fact
that coefficients are polynomials in $x^{-1}$. The Borel transform then simplifies to

$$ \Delta = \xi(\xi - 1)\frac{d}{d\xi} + (-\nu_2 + 2)\xi + \mu - 1 \quad (14.2) $$

and only the free exponents $\lambda_{[0]} = \mu - 1$ at 0 and $\lambda_{[1]} = \nu_2 - \mu - 1$ at 1 remain. Our implementation automatically performs this simplification.

First consider the case $\mu = 1$, $\nu_1 = 1$, $\nu_2 = 3/2$. After simplification, one has $\Delta = x(x-1)\frac{d}{d\xi} + \frac{1}{2}x$, $\lambda_{[0]} = 0$, $\lambda_{[1]} = -1/2$. It turns out that $\hat{Y}_1$ coincides with the analytic continuation of $-\hat{Y}_0$, though the algorithm only detects the coincidence numerically, by computing $\mathbf{T}_{0,1} = [-1.00...]$.

$$ \mathbf{L}_1 = [2\pi i e^{i\pi(\lambda_{[1]}+1)} \Gamma(\lambda_{[1]})^{-1}] = [2\sqrt{\pi}], $$

Moreover, $\hat{Y}_0$ is just $[\mathcal{B}(y_0)]$, that is, $\mathbf{B}_0 = \mathbf{I}_1$, and one has so that the Stokes multiplier in the direction $\omega = 0$ is obtained as the product $1 \times (-1.00...)(-2\sqrt{\pi} = -3.54...)$.

In the direction $\omega = \pi$, the Stokes matrix is trivial since all solutions of $\Delta$ are analytic at the origin, leading in the algorithm to $\mathbf{L}_0 = 0$.

$$ \begin{aligned} \hat{Y}_0(\zeta) &= \left( 1 + \frac{2}{9}\zeta \log \zeta + \frac{5}{9}\zeta + \dots, \quad 1 + \frac{2}{9}\zeta + \dots \right), \\ \hat{Y}_1(\zeta) &= \left( 1 - \frac{2}{9}(\zeta - 1) \log(\zeta - 1) - \frac{5}{9}(\zeta - 1) + \dots, \quad 1 - \frac{2}{9}(\zeta - 1) + \dots \right). \end{aligned} $$

The Borel transform matrices both send the relevant solution of $D$ to the element of $\hat{Y}_\alpha$ free of logarithms, that is, $\mathbf{B}_0 = \mathbf{B}_1 = [0, 1]^T$. As discussed in Sec. 2.3, the operator $\Delta$ has ‘extra’ solutions that can be interpreted as the Borel transforms of solutions of $D(y) = a$ for $a \in \mathbb{C}$. The algorithm computes the connection matrix

$$ \mathbf{T}_{0,1} = \begin{pmatrix} 0.90854... + 0.00000...i & 0.27566... + 0.00000...i \\ 0.63318... - 2.85427...i & -0.90854... - 0.86602...i \end{pmatrix} $$

and deduces $\mathbf{T}_{1,0} = \mathbf{T}_{0,1}^{-1} = \mathbf{T}_{0,1}^{-1} + 0.86602...i$. One has $\mathbf{L}_0 = \mathbf{L}_1 = (2\pi i, 0)$, leading to $c = -1.73205...i$, and in this case the Stokes matrix in the direction $\omega = \pi$ is the transpose of that in the direction $\omega = 0$.

Finally set $\mu = 0$, $\nu_1 = 0$, $\nu_2 = 1/2$. This leads to an ‘unprepared’ equation where neither of the exponents $\mu = 0$ of $y_1(x) = 1 + \frac{1}{2}x + \dots$ and $\nu_1 + \nu_2 - \mu - 1 = -1/2$ of $y_2 = e^{-1/x}x^{-1/2}$ lies in the open right-hand plane. The formal Borel transform

$$ \mathcal{B}(y_1)(\zeta) = \delta_0(\zeta) + \frac{1}{2} + \frac{3}{4}\zeta + \dots $$
contains a Dirac term, while

$$B(e^{1/x} y_2)(\zeta) = -\frac{\zeta^{-3/2}}{2\sqrt{\pi}}$$

is not integrable at the origin. One has

$$\Delta = \zeta(\zeta - 1) \frac{d^2}{d\zeta^2} + \left(\frac{9}{2}\zeta - 2\right) \frac{d}{d\zeta} + 3,$$

and

$$\hat{Y}_0(\zeta) = \left( \zeta^{-1} + \frac{1}{2} \log \zeta + \frac{3}{4} \zeta \log \zeta + \dots, 1 + \frac{3}{2} \zeta + \dots \right),$$

$$\hat{Y}_1(\zeta) = \left( (\zeta - 1)^{-3/2}, 1 - \frac{6}{5}(\zeta - 1) + \dots \right).$$

From these truncated expansions one computes $\mathbf{B}_0 = [0, 1/2]^T$, mapping $y_1$ to its Borel transform minus the Dirac terms. The connection matrix is

$$\mathbf{T}_{0,1} = \begin{pmatrix} 0.69314 \dots i & -1.00000 \dots \\ 2.00000 \dots & 0.00000 \dots \end{pmatrix}$$

and one has $\mathbf{L}_0 = (2\pi i e^{i\pi(\lambda[0]+1)}\Gamma(\lambda[0])^{-1}, 0) = (-4\sqrt{\pi}, 0)$, leading to the value $c = -3.54490 \dots i$. In the other direction, one has $\mathbf{B}_1 = [-1/(2\sqrt{\pi}), 0]^T$ and $\mathbf{L}_0 = (2\pi i, 0)$, so the Stokes matrix is trivial.

*Example 14.2 (Pólya walks).* — Denote by $w_\ell^{(d)}$ the number of nearest-neighbor walks on $\mathbb{Z}^d$ starting from the origin and returning to it after $\ell$ steps. For $3 \le d \le 15$, consider ‘the’ differential operator $\Delta^{(d)}$ of minimum order that annihilates the generating series $\sum_{\ell=0}^\infty w_\ell^{(d)}$ and let $D^{(d)}$ be its inverse Borel transform.

In the case $d=3$, one has

$$\begin{aligned} \Delta^{(d)} ={}& \xi^2(4\xi^2 - 1)(36\xi^2 - 1)\frac{d^3}{d\xi^3} + (1296\xi^5 - 240\xi^3 + 3\xi)\frac{d^2}{d\xi^2} \\ & + (2592\xi^4 - 288\xi^2 + 1)\frac{d}{d\xi} + 864\xi^3 - 48\xi, \end{aligned}$$

$$\begin{aligned} D^{(d)} ={}& 144x^{-3}\partial^6 - 1296x^{-2}\partial^5 + (-40x^{-3} + 2592x^{-1})\partial^4 + (240x^{-2} - 864)\partial^3 \\ & + (x^{-3} - 288x^{-1})\partial^2 + (-3x^{-2} + 48)\partial + x^{-1}, \end{aligned}$$
and the components of the standard fundamental solution of $D^{(d)}$ are

$$
\begin{align*}
y_1(x), y_6(x) &= e^{\pm 1/(6x)} x^{3/2} \left(1 \pm \frac{27}{4} x + \dots \right), \\
y_2(x), y_5(x) &= e^{\pm 1/(2x)} x^{3/2} \left(1 \pm \frac{17}{4} x + \dots \right), \\
y_3(x) &= x \ln x + 12x^3 \ln x + 32x^3 + \dots, \\
y_4(x) &= x + 12x^3 + \dots.
\end{align*}
$$

The operator $\Delta^{(d)}$ has exponents $(0, 0, 0)$ at 0 and $(0, 1/2, 1)$ at each of the other four singular points. At 0, our implementation uses the local basis

$$
\begin{align*}
\hat{Y}_0(\zeta) &= \left( \frac{1}{2} \ln^2 \zeta + 3\zeta^2 \ln^2 \zeta + 7\zeta^2 \ln \zeta + \mathcal{O}(\zeta^4 \ln^2 \zeta) \right), \\
&\quad \ln \zeta + 6\zeta^2 \ln \zeta + 7\zeta^2 + \mathcal{O}(\zeta^4 \ln^2 \zeta), \\
&\quad 1 + 6\zeta^2 + \mathcal{O}(\zeta^4).
\end{align*}
$$

Since the solutions $y_3(x)$ and $y_4(x)$ of $D$ attached to the Stokes value 0 only involve logarithms linearly and, by Lemma 2.1, the Borel transform cannot increase the degree in $\log(x)$, we already know that the first coordinate in this basis of their Borel transforms will be zero. Writing

$$
\begin{align*}
\mathcal{B}(y_3)(\zeta) &= \ln \zeta + \gamma + \frac{\zeta^2}{\Gamma(3)} (12 \ln \zeta + 12\gamma - 18 + 32) + \dots, \\
\mathcal{B}(y_4)(\zeta) &= \zeta + \frac{12\zeta^2}{\Gamma(3)} + \dots
\end{align*}
$$

we obtain

$$
\mathbf{B}_0 = \begin{bmatrix} 0 & 0 \\ 1 & 0 \\ \gamma & 1 \end{bmatrix}.
$$

Similarly, the expansion of the integrals

$$
\begin{align*}
\int_{\Gamma_0} (\ln \zeta + 6\zeta^2 \ln \zeta + 7\zeta^2 + \dots) e^{-\zeta/x} d\zeta &= 2\pi i x + 24\pi i x^3 + \dots \\
\int_{\Gamma_0} \left( \frac{1}{2} \ln^2 \zeta + \dots \right) e^{-\zeta/x} d\zeta &= \frac{1}{2} [2\pi i x (2 \ln x - 2\gamma - 2\pi i)] + \dots
\end{align*}
$$

yields

$$
\mathbf{L}_0 = \begin{bmatrix} 2\pi i & 0 & 0 \\ 2\pi^2 - 2\pi i \gamma & 2\pi i & 0 \end{bmatrix}.
$$

The calculation of $\mathbf{B}_\alpha$ and $\mathbf{L}_\alpha$ for the other Stokes values $\alpha$ is closer to the previous examples and we omit it. After computing the connection matrices over each of the four segments delimited by two consecutive Stokes values and forming the appropriate combinations, we eventually obtain the following
Rigorous computation of Stokes multipliers

<table><caption>Table 14.1. Time in seconds to compute all the Stokes matrices of $D^{(d)}$ (which has order $n$ and degree $ν = d$ in $1/x$) with an absolute error on each entry not exceeding $ε$. Missing entries correspond to cases where some intermediate computation failed due to insufficient precision.</caption><thead><tr><th>ν</th><th>n</th><th>ε = 10<sup>-10</sup></th><th>ε = 10<sup>-100</sup></th><th>ε = 10<sup>-1000</sup></th></tr></thead><tbody><tr><td>3</td><td>6</td><td>0.3</td><td>0.9</td><td>4.3</td></tr><tr><td>4</td><td>7</td><td>0.6</td><td>0.5</td><td>6.8</td></tr><tr><td>5</td><td>10</td><td>1.3</td><td>2.0</td><td>7.8</td></tr><tr><td>6</td><td>11</td><td>1.4</td><td>2.2</td><td>9.3</td></tr><tr><td>7</td><td>14</td><td>3.5</td><td>5.6</td><td>23.5</td></tr><tr><td>8</td><td>15</td><td>7.7</td><td>6.9</td><td>26.5</td></tr><tr><td>9</td><td>18</td><td>23.9</td><td>18.8</td><td>62.5</td></tr><tr><td>10</td><td>19</td><td>29.7</td><td>24.0</td><td>81.9</td></tr><tr><td>11</td><td>22</td><td>53.8</td><td>38.0</td><td>146.0</td></tr><tr><td>12</td><td>23</td><td>—</td><td>55.2</td><td>199.3</td></tr><tr><td>13</td><td>26</td><td>129.2</td><td>70.0</td><td>306.2</td></tr><tr><td>14</td><td>27</td><td>—</td><td>96.7</td><td>404.1</td></tr><tr><td>15</td><td>30</td><td>261.0</td><td>115.3</td><td>558.8</td></tr></tbody></table>

Stokes matrix in the direction $ω = 0$:

$$ \begin{bmatrix} 1 & 0 & 0 & \dots \\ 0.000\dots-18.000\dots i & 1 & 0 & \dots \\ 13.540\dots+0.000\dots i & 0.000\dots+1.504\dots i & 1 & \dots \\ -7.815\dots-14.179\dots i & 0.000\dots-0.868\dots i & 0 & \dots \\ 36.000\dots+0.000\dots i & 0.000\dots+8.000\dots i & 0.000\dots-2.930\dots i & \dots \\ 0.000\dots-4.000\dots i & 1.333\dots-0.000\dots i & -0.976\dots+1.772\dots i & \dots \\ & 0 & 0 & 0 \\ & 0 & 0 & 0 \\ & 0 & 0 & 0 \\ & 1 & 0 & 0 \\ & 0.000\dots-5.077\dots i & 1. & 0 \\ & -1.692\dots+0.000\dots i & 0.000\dots-0.666\dots i & 1 \end{bmatrix} $$

(along with a second, upper triangular one in the direction $ω = π$).

More generally, $D^{(d)}$ is an operator of order $n = 2d$ when $d$ is odd and $n = 2d - 1$ when $d$ is even, with coefficients of degree $ν = d$ with our usual conventions. Its Stokes values include 0 with multiplicity $d - 1$ and $[(d+3)/2]$ pairs of opposite rational numbers; in particular, they all lie on the real axis. The Stokes value 0 has multiplicity $d - 1$ and a unique associated exponent equal to 0, so that the corresponding solutions of $D^{(d)}$ involve
logarithms up to the power $d-2$. All other Stokes values have multiplicity one and an associated exponent equal to $d/2 - 1$. Correspondingly, the solutions of $\Delta^{(d)}$ at $0$ involve logarithms up to the power $d-1$, while, at other singular points, there is a basis of solutions consisting of $d-1$ power series and one solution with either half-integer exponents (for odd $d$) or logarithmic terms (for even $d$).

Table 14.1 shows the time taken by our implementation to compute the Stokes matrices of $D^{(d)}$ to various accuracies. In each case, we first call the `stokes_dict` function with the indicated $\epsilon$. If the maximum radius $\delta$ of the computed intervals exceed $\epsilon$, we repeat the computation with $\epsilon$ replaced by $\epsilon^{2.1}/\delta$, and so on until all Stokes multipliers are known with an error of less than the initial $\epsilon$. The reported running time corresponds to the total of all runs. Low-precision computations are slower than intermediate-precision ones for large instances because this procedure can require more iterations when the initial precision is small. More generally, we observe a very slow (less than linear in the number of digits!) growth of running times as the error tolerance decreases, suggesting a significant loss of accuracy over the course of the computation and weaknesses in automatic working precision management.

*Example 14.3 (A tunnel).* — For $\eta > 0$, we consider the operator

$$D = x^{-2}\partial^4 + (1 - 4x^{-2})\partial^3 + (\eta^2 x^{-1} + (5 + \eta^2)x^{-2})\partial + 2(1 + \eta^2)x^{-2}.$$

The Stokes values are 0, 2, and $1 \pm \eta i$. The transformed operator $\Delta$ has exponents $(0,0)$ at the origin, $(-4,1)$ at $2$, and $(1, \mp(2/\eta - \eta/2)i)$ at $1 \pm \eta i$. In the direction $\omega = 0$, there is a single nontrivial Stokes multiplier, corresponding to the pair of Stokes values $(0,2)$. The numerical computation of this Stokes multiplier is challenging for small $\eta$ because the line segment connecting these Stokes values passes between $1 \pm \eta i$, where large imaginary exponents are present.

To test the behaviour of our implementation in this situation, we use the following SageMath code, which computes the Stokes matrices of $D$ and extracts the Stokes multiplier of interest:

```python
sage: d = x^2*Dx
sage: dop = (d^4 + (x^2 - 4)*d^3 + (eta^2*x + 5 + eta^2)*d^2
....:      + (2*x - 2 - 2*eta^2)*d + (2 + 2*eta^2)*x)
sage: stokes = stokes_dict(dop, 1e-50)
sage: stokes[1][3,0]
```

Table 14.2 shows the output for various values of $\eta$ along with the corresponding running times. We observe that the code forgoes computing the
<table><caption>Table 14.2. Computed Stokes constants and timings for Example 14.3.</caption><thead><tr><th>η</th><th>computed Stokes constant</th><th>time (s)</th></tr></thead><tbody><tr><td>10<sup>-1</sup></td><td>[±1.75 · 10<sup>11</sup>] + [−3.34009 ... 269 · 10<sup>52</sup> ± 3.58 · 10<sup>11</sup>]i</td><td>0.4</td></tr><tr><td>10<sup>-2</sup></td><td>[±1.32 · 10<sup>490</sup>] + [−6.39027 ... 429 · 10<sup>543</sup> ± 3.01 · 10<sup>489</sup>]i</td><td>3.5</td></tr><tr><td>10<sup>-3</sup></td><td>[±4.92 · 10<sup>5403</sup>] + [−3.73685 ... 673 · 10<sup>5455</sup> ± 6.22 · 10<sup>5403</sup>]i</td><td>100</td></tr></tbody></table>

result with an absolute error matching the tolerance of $10^{-50}$ passed as input but does return a result with a *relative* error of this order of magnitude. This outcome is somewhat sensitive to the choice of the tolerance. With $\eta = 10^{-4}$, it becomes necessary to decrease the tolerance for the computation to succeed at all.

Bibliography

[1] D. G. BABBITT & V. S. VARADARAJAN, “Formal reduction theory of meromorphic differential equations: a group theoretic view”, *Pac. J. Math.* **109** (1983), p. 1-80.

[2] W. BALSER, W. B. JURKAT & D. A. LUTZ, “A general theory of invariants for meromorphic differential equations. I: Formal invariants”, *Funkc. Ekvacioj, Ser. Int.* **22** (1979), p. 197-221.

[3] W. BALSER, “Summation of Formal Power Series Through Iterated Laplace Integrals”, *Mathematica Scandinavica* **70** (1992), no. 2, p. 161-171.

[4] P. B. BORWEIN, “Reduced complexity evaluation of hypergeometric functions”, *Journal of Approximation Theory* **50** (1987), no. 3, p. 193-199.

[5] A. BOSTAN, F. CHYZAK, F. OLLIVIER, B. SALVY, E. SCHOST & A. SEDOGLAVIC, “Fast computation of power series solutions of systems of differential equations”, in *SODA’07*, SIAM, 2007, p. 1012-1021.

[6] D. BOUCHER & J.-A. WEIL, “Linear Differential Equations, Differential Galois Groups, First Integrals of Differential Systems” (Luminy, France), 2007.

[7] B. L. J. BRAAKSMA, “Multisummability and Stokes multipliers of linear meromorphic differential equations”, *Journal of Differential Equations* **92** (1991), no. 1, p. 45-75.

[8] R. P. BRENT, “A Fortran Multiple-Precision Arithmetic Package”, *ACM Transactions on Mathematical Software* **4** (1978), no. 1, p. 57-70.

[9] D. V. CHUDNOVSKY & G. V. CHUDNOVSKY, “Approximations and complex multiplication according to Ramanujan”, in *Ramanujan revisited*, Academic Press, 1988, p. 375-472.

[10] ——, “Computer algebra in the service of mathematical physics and number theory”, in *Computers in Mathematics* (D. V. Chudnovsky & R. D. Jenks, eds.), Lecture Notes in Pure and Applied Mathematics, vol. 125, Dekker, 1990, p. 109-232.

[11] O. COSTIN, *Asymptotics and Borel summability*, Chapman Hall/CRC Monogr. Surv. Pure Appl. Math., vol. 141, Boca Raton, FL: Chapman & Hall/CRC, 2009.

[12] J. A. DE LOERA, J. RAMBAU & F. SANTOS, *Triangulations*, Algorithms and Computation in Mathematics, vol. 25, Springer, Berlin, Heidelberg, 2010.

[13] J. DIEUDONNÉ, *Calcul infinitésimal*, Hermann, Paris, 1968.

[14] A. DUVAL & C. MITSCHI, “Matrices de Stokes et groupe de Galois des équations hypergéométriques confluentes généralisées”, *Pacific J. Math.* **138** (1989), no. 1, p. 25-56.