J. Shelton, S. Crew, P. H. TrinhThe higher-order Stokes phenomenon

# Exponential asymptotics and higher-order Stokes phenomenon in singularly
perturbed ODEs ††thanks: Submitted to the editors 01/03/2023. JS and SC are
joint lead authors. PHT is supported by the Engineering and Physical Sciences
Research Council EP/V012479/1. This work was supported by EPSRC grant
EP/R014604/1.

Josh Shelton Department of Mathematical Sciences, University of Bath, BA2 7AY,
UK<EMAIL_ADDRESS><EMAIL_ADDRESS>Samuel Crew22footnotemark: 2
Faculty of Computer Science, Ruhr University, Bochum, Germany. (),
Max Planck Institute for Security and Privacy, Bochum, Germany.
<EMAIL_ADDRESS><EMAIL_ADDRESS>Philippe H.
Trinh22footnotemark: 2

###### Abstract

Outside the area of exponential asymptotics, the concept of the higher-order
Stokes phenomenon remains somewhat esoteric. The intention of this work is to
provide several examples of relatively simple ordinary differential equations
where the phenomenon arises, and to develop additional practical methodologies
for its study. In particular, we show how the higher-order Stokes phenomenon
may be derived through the hyperterminant representation of the late-term
divergence of an asymptotic expansion developed by Olde Daalhuis (J. Comp.
Appl. Math. vol. 76, 1996, pp. 255-264). Lower-order components of the
factorial-over-power divergence are considered, for which late-late-term
divergence arises. Borel resummation of these lower-order components reveals
the higher-order Stokes phenomenon, in which new components of the late-terms
of the expansion are smoothly switched on with an error function dependence.
The techniques are firstly demonstrated with a second-order linear
inhomogeneous ODE that exemplifies the simplest example of higher-order Stokes
phenomena. Further examples studied include higher-order equations and
eigenvalue problems.

###### keywords:

Stokes phenomenon, exponential asymptotics, beyond-all-orders analysis

##  1 Introduction

Outside the exponential asymptotics community, the concept of the _higher-
order Stokes phenomenon_ (HOSP) is still considered to be relatively obscure.
Yet, as we shall explain in this work, it occurs generically and routinely in
divergent asymptotic expansions. Deep understanding of the phenomenon seems to
be necessary in order to perform beyond-all-orders analysis of multi-
dimensional problems (such as in partial differential equations) arising in
many physical applications (Howls et al., 2004; Chapman & Mortimer, 2005; Body
et al., 2005). The emphasis of this work is to provide some relatively simple
examples of the HOSP in the context of singularly perturbed ordinary
differential equations. For example, we show that it occurs in the asymptotic
analysis of the seemingly innocuous equation,

(1) $\epsilon^{2}y^{\prime\prime}+\epsilon(1+z)y^{\prime}+zy=1$

with $\epsilon\to 0$ and $y=y(z)$. We shall explain why the above is likely
the simplest example of HOSP, and we develop applicable techniques for its
study. Finally, our work serves to survey perspectives of the phenomenon
across disciplines in asymptotic analysis.

We first review the regular Stokes phenomenon. In applications to singularly
perturbed ordinary differential equations, the solution is typically expanded
as an asymptotic expansion. For example, we approximate the solution as

(2) $y(z;\epsilon)\sim y^{(0)}(z;\epsilon)\equiv y_{0}(z)+\epsilon
y_{1}(z)+\cdots+\epsilon^{n}y_{n}(z)+\cdots$

with $\epsilon\to 0$ and $z\in\mathbb{C}$, and refer to such an expansion as
the base series. The late terms, $y_{n}$, with $n\to\infty$ diverge and this
divergence is coupled to the existence of Stokes lines in $z\in\mathbb{C}$.
Across such lines, exponentially-small terms, say
$\sigma(z,\epsilon)\mathrm{e}^{-\chi(z)/\epsilon}$, are switched-on within a
small boundary layer via the (regular) Stokes phenomenon. The quantity
$\sigma$ is called the _Stokes multiplier_. By now, such ideas are well
established [cf. reviews in e.g. Dingle (1973) and Berry (1989)], and there
exists a plethora of tools in asymptotics beyond-all-orders to study these
effects in the context of ordinary and partial differential equations,
difference equations, integral equations, and beyond.

Typically, the divergence of $y_{n}$ is approximated by a factorial-over-power
growth as $n\to\infty$. However, this is generally only the leading order term
in an asymptotic expansion in decreasing powers of $n$. If $\delta=1/n$ is now
considered to be the perturbative parameter, then it is interesting to
consider if $y_{n}$ can, itself, exhibit Stokes phenomena when analytically
continued in $z$. Indeed this is the higher-order Stokes phenomenon associated
with higher-order Stokes lines (HOSL). The standard Stokes phenomenon
described above is linked to the divergence in $y_{n}$, but now, since these
approximations may be non-uniform in $\mathbb{C}_{z}$, the regular Stokes line
structure is consequently altered.

For the practitioner interested in the prediction of exponentially-small
terms, two significant changes can occur as a consequence of HOSP:

  1. (i)

Stokes lines can be truncated or there can be new Stokes lines; and

  2. (ii)

there can be atypical values of the Stokes multiplier(s).

We shall explore both of these consequences in this work.

###  1.1 Prior literature and ongoing work

There are a number of seminal works that study HOSP, though not necessarily in
the same context of our work here. The origins of the phenomenon are generally
attributed to the work of Berk et al. (1982). In their paper, titled “ _New
Stokes’ line in WKB theory_ ”, the altered structure of Stokes lines was
discovered in the context of an integral equation (Pearcey’s integral). The
theory of such “new Stokes lines” was then extensively studied in the exact
WKBJ (Liouville-Green) communities (cf. Aoki 1994 and references therein).
Further WKBJ examples with exact integral solutions were studied by Aoki et
al. (1998, 2001); a singularly perturbed Painlevé equation appears in the work
of Honda (2007); and HOSP appears as well in the analysis of the Henon map in
quantum chaos by Shudo (2008). We refer readers to reviews by _e.g._ Aoki et
al. (2001), Honda et al. (2015), and Takei (2017) and references therein for a
historical survey.

We note an important line of research in the early 2000s, initiated by a
triplet of works by Howls et al. (2004), Chapman & Mortimer (2005), and Body
et al. (2005), which characterise some important connections between beyond-
all-orders analysis of partial differential equations and the HOSP. These
three works apply a number of different methodologies to beyond-all-order
asymptotics, and present complementary views of similar problems. In
particular, the work of Howls et al. (2004) returned to analyse the Pearcey
integral, discussed above, and established a connection between HOSP and the
prior theory of hyperterminants developed by Olde Daalhuis (1996, 1998). A
particularly simple example of the HOSP was presented by Olde Daalhuis (2004)
in the context of an inhomogeneous linear second-order ODE.

There have been a number of recent exciting works concerning HOSP in both
applications and theoretical studies. For example, we highlight the work of
Nemes (2022) who studied atypical Stokes multipliers arising in the asymptotic
representation of the Gamma function. Further applications have been
discovered in fluid dynamics by Lustri et al. (2019) for three dimensional
water waves, and Shelton et al. (2023a, b) for geophysical instabilities. Most
recently, we note unpublished investigations by Olde Daalhuis (2022) into
smoothing of HOSP, and King (2022) on the occurrence of HOSP in linear PDEs.

###  1.2 Outline

In section 2, we present a unified view of three different perspectives for
understanding HOSP: the Borel plane considered by e.g. Howls et al. (2004),
hyperasymptotics and analysis of the transseries (Olde Daalhuis, 1996, 1998),
and the application of a complete late-term ansatz to study lower-order
components (as $n\to\infty$) of the late terms. Thus for example, when our
minimal example (1) is studied in the Borel plane, the transformed solution
exhibits a square root branch point and a simple pole; we shall explain why
these are essentially the minimal ingredients for developing HOSP.

In the context of (1), we shall consider corrections to the late-term
divergence in the form of a $1/n$ transseries expansion (this is analogous to
the hyperterminant representation by Olde Daalhuis 1996, 1998). In section 3,
it is shown that a new divergent series emerges within the late-terms
themselves, which we denote the late-late-terms of the asymptotic expansion.
Borel summation of these late-late-terms reveals the higher-order Stokes
lines, across which new components of the late-terms are seen to smoothly
switch on with an error function dependence.

Finally, in section 4, we apply our general techniques to detect the higher-
order Stokes lines in three previously encountered ODEs. These problems had
previously been studied through methods that cannot be easily generalised to
arbitrary nonlinear differential equations—such as steepest descent analysis
on an integral representation, or evaluation of a known recurrence relation.
Together, the methods studied in this work will provide additional practical
tools for the derivation of HOSP and HOSL.

##  2 Transseries, divergence, and the Borel transform

###  2.1 Transseries relations

We begin with a review of the relationship between exponential asymptotic
expansions (the transseries) and the late terms of an asymptotic series. Note
that in general, the base series, written in (2), represents only one of many
possible contributions to the formal parametric transseries of
$y(z;\epsilon)$. The full transseries can be written as

(3) $\begin{split}y(z;\epsilon)&=\biggl{[}y_{0}(z)+\epsilon
y_{1}(z)+\cdots+\epsilon^{n}y_{n}(z)+\cdots\biggr{]}\\\
&\quad+\sum_{j=1}\sigma_{j}(z;\epsilon)\biggl{[}y^{(j)}_{0}(z)+\epsilon
y^{(j)}_{1}(z)+\cdots+\epsilon^{n}y^{(j)}_{n}(z)+\cdots\biggr{]}\mathrm{e}^{-\chi_{j}(z)/\epsilon},\end{split}$

where the components $y_{n}^{(j)}$ and $\sigma_{j}$ are to be determined, and
following Dingle, the functions $\chi_{j}(z)$ are known as singulants. We
focus on factorially divergent asymptotic series; this structure frequently
arises in solutions to singularly perturbed ordinary differential equations
(Dingle, 1973).

As we have introduced previously, it is possible to re-sum the divergent
series (2) using Borel resummation, or to equivalently study the remainder
after optimal truncation. This procedure gives rise to the Stokes phenomenon,
in which exponentially subdominant terms switch on across Stokes lines. For
the case above, the associated Stokes lines of the asymptotic expansion are
curves in the $\mathbb{C}_{z}$ plane specified by the loci

(4)
$l_{\chi_{j}}:\quad\text{Im}[\chi_{j}(z)]=0\qquad\text{and}\qquad\text{Re}[\chi_{j}(z)]\geq
0.$

Across such curves, terms proportional to $\mathrm{e}^{-\chi_{j}(z)/\epsilon}$
are switched on. Indeed, these additional components of the asymptotic
solution are contained within the transseries (3). Analysis of the Stokes
phenomenon also determines the multipliers, $\sigma_{j}$, which rapidly and
smoothly transition across Stokes lines (Berry, 1989).

Let us assume that the components, $y_{n}(z)$, of the base series are
meromorphic functions with singular behaviour at a point, say $z=z_{1}$. The
leading-order divergence associated with this point is given by

(5) $y_{n}(z)\sim
y^{(1)}_{0}(z)\frac{\Gamma(n+\alpha_{1})}{\chi_{1}(z)^{n+\alpha_{1}}}\qquad\text{as
$n\to\infty$},$

where $\chi_{1}(z_{1})=0$ and typically $\alpha_{1}$ is constant. We denote
(5) as the _late-term divergence_. Note that the singulant, $\chi_{1}$, and
the amplitude function, $y_{0}^{(1)}$, in (5) above are those same functions
that appear in the transseries (3), about the exponential
$\mathrm{e}^{-\chi_{1}/\epsilon}$. This connection between late terms and
exponentials is fundamental in exponential asymptotics.

Surprisingly, the relationship between the late-terms, $y_{n}$, of the base
expansion, and each of the amplitudes $y_{p}^{(j)}$ from the
$\mathrm{e}^{-\chi_{j}/\epsilon}$ transseries component can be found through
the consideration of the lower-order components of $y_{n}$ as $n\to\infty$
(typically in inverse powers of $n$). More specifically, we may include
additional $1/n$ corrections to the late-terms of (5) and organise $y_{n}$ as
follows:

(6)
$y_{n}(z)\sim\sum_{j=1}\left(\sum_{p=0}^{\infty}y^{(j)}_{p}(z)\frac{\Gamma(n-p+\alpha_{j})}{\chi_{j}(z)^{n-p+\alpha_{j}}}\right)\qquad\text{as
$n\to\infty$}.$

Above, the outer sum over $j$ is over the singulants, $\chi=\chi_{j}$. The
inner sum, indexed by $p$, includes the lower-order correction terms.

We note that lower-order corrections to the late-term divergence have
previously been calculated by e.g. King & Chapman (2001) and Shelton et al.
(2023a), but these prior studies considered a series in powers of $n^{-1}$,
and thus the connection with the transseries (3) was not clear. We refer the
reader to e.g. Crew & Trinh (2023) for a more thorough review and explanation
of this connection between the $\epsilon$-transseries expansion (3) and the
$1/n$ expansion of the late terms (6).

###  2.2 The Borel plane

Borel resummation encodes the parametric transseries (3) as a single function
$y_{B}(w,z)$, which is holomorphic in $w$. The location of singularities of
$y_{B}$ encode the singulants, $\chi_{j}(z)$, and the power of the singularity
encodes the constants $\alpha_{j}$. To see this, consider the formal base
series (2). The Borel transform with respect to $\epsilon$ is a power series
on the Borel plane, $w\in\mathbb{C}_{w}$, defined by111The $\delta$ term here
is a formal symbol that allows us to consider the Borel transform of
perturbative series beginning at $O(1)$. It is defined such that
$\mathcal{B}^{-1}[\delta]=1$. In the following we study analytic properties of
the series $y_{1}(z)+wy_{2}(z)+\cdots$ where the leading term plays little
role.

(7) $y_{B}(w,z):=\delta
y_{0}(z)+\sum_{n=0}^{\infty}\frac{y_{n+1}(z)}{n!}w^{n}=y_{0}(z)\delta+y_{1}(z)+y_{2}(z)w+\frac{y_{3}(z)}{2!}w^{2}+\cdots.$

We assume throughout this work that $y_{B}$ is endlessly analytically
continuable (cf. Mitschi et al. 2016) except at a finite number of points
$z\in\mathbb{C}_{z}$, and extends to a function with at most exponential
growth along any ray emanating from the origin. Under these conditions we may
then define the inverse Borel transform of $y_{B}$ via

(8)
$y(z;\epsilon)=\mathcal{B}^{-1}[y_{B}]=\int_{0}^{\infty}y_{B}(w,z)\mathrm{e}^{-w/\epsilon}\operatorname{d\\!}{}{w}=y_{0}(z)+\epsilon
y_{1}(z)+\cdots.$

We now explain how Stokes phenomena can be understood via the Borel plane.
Consider the asymptotic evaluation of the integral (8) for small $\epsilon$.
Provided singularities of $y_{B}$ do not lie on the integration contour, we
may demonstrate using the integral formula of the Gamma function that
$\mathcal{B}^{-1}[y_{B}]=y_{0}(z)+\epsilon y_{1}(z)+\cdots$. However, if we
suppose now that $y_{B}(w,z)$ has an algebraic singularity in the analytic
continuation, of the form

(9)
$y_{B}(w,z)=\frac{1}{(w-\chi_{1}(z))^{\alpha_{1}}}\Big{[}y_{0}^{(1)}(z)+\big{(}w-\chi_{1}(z)\big{)}y_{1}^{(1)}+\cdots\Big{]},$

then whenever the singularity at $w=\chi_{1}(z)$ crosses the integration
contour $(0,\infty)$, we pick up an extra contribution to the asymptotics of
$\mathcal{B}^{-1}[y_{B}]$ from the Hankel contour integral, denoted
$\mathcal{H}_{\chi}$,

(10)
$\int_{\mathcal{H}_{\chi}}\frac{y_{0}^{(1)}(z)}{\big{(}w-\chi_{1}(z)\big{)}^{\alpha_{1}}}\mathrm{e}^{-w/\epsilon}\,\operatorname{d\\!}{}{w}\sim\sigma_{1}(z,\epsilon)y_{0}^{(1)}(z)\mathrm{e}^{-\chi_{1}(z)/\epsilon},$

where $\sigma_{1}$ is the Stokes multiplier and can be derived via an inner-
matching procedure. The above transition occurs across the locus
$l_{\chi_{1}}$ in the physical $\mathbb{C}_{z}$ plane, which corresponds to
the Stokes lines in (4); thus switching-on an exponentially small contribution
in the transseries of $y(z,\epsilon)$. Further details of the above can be
found in Crew & Trinh (2023).

###  2.3 Higher-order Stokes phenomena

Consider a scenario in which the Borel transform, $y_{B}$, has two
singularities, say at $w=\chi_{1}(z)$ and $w=\chi_{2}(z)$, in the Borel plane
$\mathbb{C}_{w}$. In the following, we present a geometrical explanation of
HOSP.

  1. (i)

One singularity, say $w=\chi_{1}$, is assumed to be a branch point in
$\mathbb{C}_{w}$, with a zero $\chi_{1}(z_{1})=0$. Then, according to the
correspondence elucidated in §2.2, this singularity generates divergence in
the base series, with $y^{(0)}_{n}=y_{n}$ diverging as $n\to\infty$. The
result is a standard $\text{B}>1$ Stokes line at $l_{\chi_{1}}$, where
$\chi_{1}$ is real and positive.

  2. (ii)

Note that there may also exist a point, $z=z_{2}$, where $\chi_{2}(z_{2})=0$.
While this location is also associated with a singularity in the analytic
continuation of the Borel transform of the transseries (3), it is not
necessarily present in the base series $y^{(0)}(z,\epsilon)$.

  3. (iii)

Consider a Stokes crossing point in $\mathbb{C}_{z}$, $z=z_{*}$, where
$\chi_{1}=\chi_{2}$. In the Borel plane this is where the two singularities in
$\mathbb{C}_{w}$ coalesce and therefore the Borel power series $y_{B}^{(1)}$,
defined by the Borel transform of the transseries component around
$\mathrm{e}^{-\chi_{1}/\epsilon}$, is singular at $z=z_{2}$. In the
perturbative language it means that the higher transseries components,
$y^{(1)}_{n}(z)$, are singular at $z=z_{2}$. This yields a $1>2$ Stokes line
at $l_{\chi_{2}}$, where $\chi_{2}$ is real and positive.

  4. (iv)

However, since crucially our Borel plane is ramified, $w=\chi_{2}$ may not
always be on the same Riemann sheet as the integration contour $(0,\infty)$
emanating from the distinguished origin (relative to $y_{n}^{(0)}$); then even
if the Stokes line $l_{\chi_{2}}$ is crossed in $\mathbb{C}_{z}$, the
inversion contour of the Borel transform may not produce a Hankel
contribution. The critical locus where $\chi_{2}$ is visible from rays
emanating from the distinguished origin, $w=0$, is the condition that the two
singularities are co-linear on a half line:

(11)
$\text{Im}\bigg{[}\frac{\chi_{2}-\chi_{1}}{\chi_{1}}\bigg{]}=0\qquad\text{and}\qquad\text{Re}\bigg{[}\frac{\chi_{2}-\chi_{1}}{\chi_{1}}\bigg{]}\geq
0.$

The above is a geometrical explanation of HOSP from the perspective of the
Borel plane, which produced condition (11) for a higher-order Stokes line.
Thus, it is only in the correct region of the $\mathbb{C}_{z}$ plane, on one
or the other side of (11), that the Stokes line $l_{\chi_{2}}$ is active
(which side depends on the problem setup). The Stokes line and higher-order
Stokes line is illustrated in figure 1. Our explanation, here, is similar and
complementary to those given in Howls et al. (2004) and Chapman & Mortimer
(2005). We shall provide a concrete example of the above geometry in section
3.

Figure 1: In order for the Borel singularity at $w=\chi_{2}$ to induce the
Stokes phenomenon by crossing the integration contour $(0,\infty)$ of the
inverse Borel transform, it must be present on the principle Riemann sheet
when viewed from the origin, $O$. In ($b$), as the singularity at $w=\chi_{2}$
moves towards the real axis, it is on the principal sheet and results in a
$B>2$ Stokes phenomenon. However, in ($c$), the singularity at $w=\chi_{2}$
crosses the real-axis on a higher Riemann sheet and thus the $B>2$ Stokes line
is truncated. The geometrical colinearity condition (11) divides these two
illustrative scenarios.

One of the basic principles of resurgence is that the singularity structure of
$y_{B}$ is encoded in the asymptotics (in $n$) of the late-terms of the power
series coefficients, i.e. $y_{n}(z)$ as $n\to\infty$. This follows by
Darboux’s lemma and extensions thereof. Throughout the rest of the work we
shall study the HOSP from the perspective of the large-$n$ asymptotics of the
late terms, $y_{n}$. Then, the visibility of $\chi_{2}$ from the origin may be
realised as a Stokes phenomena (in $1/n$) in the expansion of the late terms
near $w=\chi_{1}$, which occurs on the colinear locus (11). In other words,
the geometric condition (11) is equivalent to which divergent components are
present in the late-term asymptotics of $y_{n}$.

In fact, from the late terms perspective, we typically only learn about the
presence of an additional singularity, $w=\chi_{2}(z)$, by noting that the
terms $y_{n}^{(1)}$ contain an extra physical singularity at some $z=z_{2}$
(that is not necessarily a zero of $\chi_{1}$). In this sense, it is
appropriate to refer to the truncated Stokes line associated to $\chi_{2}$ as
a new Stokes line. We now turn to this late terms perspective in more detail.

###  2.4 Late-late-term divergence and optimal truncation

The HOSP may also be detected purely through the study of the late-terms of
the asymptotic expansion. The presence of an additional singularity in the
Borel plane indicates that the component of the transseries around
$\mathrm{e}^{-\chi_{1}/\epsilon}$,

(12) $\mathrm{e}^{-\chi_{1}(z)/\epsilon}\left(y^{(1)}_{0}+\epsilon
y^{(1)}_{1}+\cdots\right),$

forms a divergent asymptotic expansion. By the correspondence discussed in
section 2.1, we have that the $1/n$ late-term representation,

(13)
$y_{n}(z)\sim\sum_{p=0}^{\infty}y^{(1)}_{p}(z)\frac{\Gamma(n-p+\alpha)}{\chi_{1}(z)^{n-p+\alpha}}\qquad\text{as
$n\to\infty$},$

is also divergent as $p\to\infty$. We express this divergence of $y_{p}^{(1)}$
in the factorial-over-power form of

(14) $y_{p}^{(1)}(z)\sim
y^{(2)}_{0}(z)\frac{\Gamma(p+\beta)}{\tilde{\chi}(z)^{p+\beta}}\qquad\text{as
$p\to\infty$}.$

Here, $\beta$ is a constant, and $\tilde{\chi}$ is the late-late-term
singulant for which $\tilde{\chi}(z_{2})=0$. Consistent with the simplicity of
the Borel plane assumed in section 2.3, we assume that $y_{0}^{(2)}$ contains
no further singular point.

Due to the divergence of $y_{p}^{(1)}$ as $p\to\infty$ appearing in the
infinite series (6), optimal truncation is required at $p=P-1$, where
consecutive terms are of the same order. The remainder to this optimally
truncated series then displays the Stokes phenomenon across a HOSL, which are
contours satisfying the Dingle conditions

(15)
$\text{Im}\left[\frac{\tilde{\chi}}{\chi_{1}}\right]=0\qquad\text{and}\qquad\text{Re}\left[\frac{\tilde{\chi}}{\chi_{1}}\right]\geq
0.$

Thus the geometrical condition (11) is now realised precisely as a Stokes
line, but in the $1/n$ asymptotic series. This transition across the HOSL is
smooth, and follows error function dependence over a boundary layer of width
$O(n^{-1/2})$. We demonstrate this in section 3.4 through Borel resumming the
remainder, which is analogous to the method used by Berry (1989) to smooth the
regular Stokes phenomenon. Furthermore, the term switched on across the
higher-order Stokes line (15) is another factorial-over-power component of the
late-term divergence, with a singulant given by $\chi+\tilde{\chi}$. This is
the HOSP, which results in a modified late-term divergence of the base series,
given by

(16)
$y_{n}(z)\sim\tilde{\sigma}(z,\epsilon)y_{0}^{(2)}(z)\frac{\Gamma(n+\alpha+\beta)}{(\chi_{1}+\tilde{\chi})^{n+\alpha+\beta}}+\sum_{p=0}^{P-1}y^{(1)}_{p}(z)\frac{\Gamma(n-p+\alpha)}{\chi_{1}(z)^{n-p+\alpha}}.$

In the above, $\tilde{\sigma}$ is the higher-order Stokes multiplier, which
rapidly transitions across the HOSL (15). Thus, in addition to the Stokes
phenomenon induced by the $\chi_{1}$ divergence, there will be an additional
Stokes phenomenon induced by $\chi_{1}+\tilde{\chi}$, which only occurs where
$\tilde{\sigma}\neq 0$.

This yields the $l_{\chi_{1}}$ and $l_{\chi_{2}}$ Stokes lines previously
specified in section 2.3 with $\chi_{2}=\chi_{1}+\tilde{\chi}$. Following Berk
et al. (1982), contours that satisfy $\text{Im}[\chi_{2}]=0$,
$\text{Re}[\chi_{2}]\geq 0$, and $\tilde{\sigma}\neq 0$ are often referred to
as new Stokes lines. Thus, we find the following transseries components with
singulants of $\chi_{1}$ and $\chi_{1}+\tilde{\chi}$,

(17) $\begin{split}y(z,\epsilon)&\sim y_{0}(z)+\epsilon
y_{1}(z)+\epsilon^{2}y_{2}(z)+\cdots+\sigma_{1}(z,\epsilon)y^{(1)}_{0}(z)\mathrm{e}^{-\chi_{1}/\epsilon}\\\
&\quad+\sigma(z,\epsilon)y^{(2)}_{0}(z)\mathrm{e}^{-(\chi_{1}+\tilde{\chi})/\epsilon}.\end{split}$

The additional exponential present in (17), generated by the new Stokes line,
is often necessary to avoid contradictions arising from the intersection of
regular Stokes lines.

###  2.5 Summary

The HOSP may be viewed from the two perspectives outlined in this section.
These are:

  1. (i)

The late-late-term divergence. The leading-order component of the late-terms
exhibits factorial-over-power divergence of
$O(\Gamma(n+\alpha)/\chi^{n+\alpha})$. Lower order components of this
divergence, of $O(\Gamma(n+\alpha-p)/\chi^{n+\alpha-p})$, were considered in
section 2.4. The amplitude functions at each of these orders diverge in the
limit of $p\to\infty$; this is referred to as late-late-term divergence. Borel
resummation of the optimally truncated late-late-term expansion, performed in
section 3.4, reveals the HOSP.

  2. (ii)

The Borel transform. In the Borel plane, $\mathbb{C}_{w}$, HOSP occurs when
there are at least two singularities, one of which is a branch point with a
non-trivial Riemann sheet structure. The HOSL is then a condition for a
particular singularity to be ‘visible’ from the origin and therefore cause
Stokes phenomena. This leads to a truncated ‘new Stokes line’.

##  3 A model problem

We now introduce a new model example whose Stokes line structure contains the
simplest possible realisation of the HOSP, due to the presence of one branch
point and a simple pole in the Borel plane. Let us consider the following
singularly perturbed second order linear inhomogeneous differential equation:

(18) $\epsilon^{2}y^{\prime\prime}+\epsilon(1+z)y^{\prime}+zy=1,$

where $y(z,\epsilon)$ is the unknown function and $0<\epsilon\ll 1$ is a small
parameter.

###  3.1 The transseries solution

First, we seek the formal transseries expansion by writing

(19)
$y(z,\epsilon)=\sum_{i=0}\mathrm{e}^{-\chi_{i}(z)/\epsilon}y^{(i)}(z,\epsilon),$

where $\chi_{i}(z)$ is the singulant function, and $y^{(i)}(z,\epsilon)$ is
the solution corresponding to the singulant $\chi_{i}$. Substitution of
transseries ansatz (19) into the differential equation (18) yields the three
possible solutions for the singulant, which are given by

(20)
$\chi_{0}(z)=0,\qquad\chi_{1}^{\prime}(z)=z,\qquad\chi^{\prime}_{2}(z)=1.$

The first of these, $\chi_{0}$, corresponds to the base solution, whose
asymptotic expansion will be a standard power series in $\epsilon$.
Furthermore, equations for each of the amplitude solutions,
$y^{(i)}(z,\epsilon)$, are obtained as

(21a)
$\displaystyle\epsilon^{2}y^{(0)\prime\prime}+\epsilon(1+z)y^{(0)\prime}+zy^{(0)}=1,$
(21b) $\displaystyle\epsilon
y^{(1)\prime\prime}+(1-z)y^{(1)\prime}-y^{(1)}=0,$ (21c)
$\displaystyle\epsilon y^{(2)\prime\prime}+(z-1)y^{(2)\prime}=0.$

We begin by expanding the base solution, $y^{(0)}(z,\epsilon)$, in powers of
$\epsilon$ as

(22) $y^{(0)}(z,\epsilon)=\sum_{n=0}^{\infty}\epsilon^{n}y_{n}(z),$

for which substitution into the base equation (21a) yields the early orders of
the solution expansion as

(23)
$y^{(0)}(z,\epsilon)=\frac{1}{z}+\epsilon\frac{(1+z)}{z^{3}}+\epsilon^{2}\frac{(2z^{2}+3z+3)}{z^{5}}+\cdots.$

Note that there is a singularity at $z=0$, near to which the asymptotic series
(23) reorders. The resultant inner analysis in this region is performed in
Appendix A. At $O(\epsilon^{n})$ in the base equation (21a), we find the
equation

(24) $y_{n-2}^{\prime\prime}+(1+z)y_{n-1}^{\prime}+zy_{n}=0,$

which holds for $n\geq 2$. As successive orders of the asymptotic expansion
are determined by the differentiation of previous orders, the solution $y_{n}$
will diverge as $n\to\infty$. This divergence is studied in the next section.

###  3.2 Late-term divergence

To capture the factorial-over-power divergence of $y_{n}$, we consider a late-
term ansatz of the form

(25)
$y_{n}^{(0)}\sim\sum_{p=0}^{\infty}B_{p}(z)\frac{\Gamma(n-p+\alpha)}{[\chi(z)]^{n-p+\alpha}}\qquad\text{as
$n\to\infty$},$

where $\chi$ and $B_{p}$ are the late-term singulant and amplitude functions
respectively, and $\alpha$ is a constant. In (25), in addition to the leading-
order divergence specified by $p=0$, we have also included lower-order
components. In fact, there will be additional components of the late-term
divergence not captured by the ansatz (25). These will occur due to the HOSP,
in which they are smoothly switched on across a boundary-layer, of diminishing
width as $n\to\infty$, that surrounds a HOSL. The key result of this paper is
that this phenomenon may be derived through the study of the late-late-terms,
$B_{p}$, as $p\to\infty$, which is performed in section 3.4.

Substitution of the factorial-over-power ansatz (25) into the
$O(\epsilon^{n})$ equation (24) yields at each order of
$\Gamma(n-p+\alpha)/\chi^{n-p+\alpha}$ the equations

(26a) $\displaystyle p=0:$
$\displaystyle\qquad(\chi^{\prime}-1)(\chi^{\prime}-z)=0,$ (26b)
$\displaystyle p=1:$
$\displaystyle\qquad(1+z-2\chi^{\prime})B_{0}^{\prime}-\chi^{\prime\prime}B_{0}=0,$
(26c) $\displaystyle p\geq 2:$ $\displaystyle\qquad
B_{p-2}^{\prime\prime}+(1+z-2\chi^{\prime})B_{p-1}^{\prime}-\chi^{\prime\prime}B_{p-1}=0.$

Integration from the singularity at $z=0$, in order to satisfy the matching
criteria of $\chi(0)=0$, yields the two solutions

(27) $\chi_{1}(z)=z^{2}/2\qquad\text{and}\qquad\chi_{2}(z)=z.$

Since the latter of these two singulants yields $B_{p}(z)=\text{const}$ for
$p\geq 0$, it will not generate any late-term Stokes switching. We therefore
focus on the first late-term singulant, $\chi=z^{2}/2$. The first of the
amplitude equations, (26b), may then be solved to find

(28) $B_{0}(z)=\frac{\Lambda_{0}}{1-z},$

where $\Lambda_{0}$ is a constant of integration. Note that solution (28) is
singular at $z=1$, and that equation (26b) for $B_{1}$ requires
differentiation of $B_{0}$. Hence, $B_{1}$ will have stronger singular
behaviour at $z=1$. This pattern continues and leads to factorial-over-power
divergence of $B_{p}$ as $p\to\infty$, which we denote the late-late-terms of
the asymptotic expansion. This is considered next in section 3.3.

###  3.3 Late-late-term divergence

In this section, we determine the divergence of $B_{p}$ as $p\to\infty$, which
arises as a consequence of the singularity at $z=1$ in solution (28) for
$B_{0}$. The power of this singularity will grow as we proceed into the
asymptotic series. For instance, equation (26b) may be solved for $B_{1}$,
yielding

(29) $B_{1}(z)=-\frac{\Lambda_{0}}{(1-z)^{3}}+\frac{\Lambda_{1}}{1-z},$

and this growing singular behaviour generates factorial-over-power divergence
of $B_{p}$ as $p\to\infty$. We begin by considering the factorial-over-power
ansatz

(30) $B_{p}(z)\sim
C(z)\frac{\Gamma(p+\beta)}{[\tilde{\chi}(z)]^{p+\beta}}\qquad\text{as
$p\to\infty$},$

where $\tilde{\chi}$ and $C$ are the late-late-term singulant and amplitude
functions, and $\beta$ is a constant. In order to match with an inner solution
for $y^{(0)}_{n}$ near $z=1$, we require that $\tilde{\chi}(1)=0$.
Substitution of ansatz (30) into equation (26c) yields the equations
$\tilde{\chi}^{\prime}=1-z$ and $C^{\prime}=0$, which we integrate to find

(31) $\tilde{\chi}=-\frac{(z-1)^{2}}{2}\qquad\text{and}\qquad
C(z)=\tilde{\Lambda}.$

In (31) above, we have used the matching condition $\tilde{\chi}(1)=0$, and
$\tilde{\Lambda}$ is a constant, determined by an inner matching procedure as
$\tilde{\Lambda}=\mathrm{i}/(2\pi)$ in (75) of Appendix A.

###  3.4 Smoothing the higher-order Stokes phenomenon

In this section, we demonstrate that optimal truncation of the late-term
ansatz (25) at $p=P-1$ yields an exponentially-small (in $n$) remainder that
displays the Stokes phenomenon. This remainder corresponds to an additional
factorial-over-power component of the late-terms with a singulant given by
$\chi+\tilde{\chi}$. Furthermore, we show how this change occurs smoothly
across a HOSL, which is located where $\text{Im}[\tilde{\chi}/\chi]=0$ and
$\text{Re}[\tilde{\chi}/\chi]\geq 0$. This analysis closely follows the method
used by Berry (1989), in which the tail of the divergent expansion is re-
summed, and examination of the resultant integral reveals error function
behaviour localised across the Stokes line.

We begin by truncating the late-term ansatz (25) at $p=P-1$ by writing

(32)
$y_{n}^{(0)}\sim\sum_{p=0}^{P-1}B_{p}(z)\frac{\Gamma(n-p+\alpha)}{[\chi(z)]^{n-p+\alpha}}+\underbrace{\sum_{p=P}^{\infty}B_{p}(z)\frac{\Gamma(n-p+\alpha)}{[\chi(z)]^{n-p+\alpha}}}_{R_{P}},$

and consider $P$ to be large. Note that later in equation (38), we specify $P$
in terms of $n$ by minimising $R_{p}$. In substituting for both the late-late-
term divergent ansatz of $B_{p}$ from (30) and the integral definition of the
Gamma function (since $n-P$ is also large), we also swap the order of
summation and integration. This yields a geometric series, which may be
evaluated exactly to find

(33)
$\begin{split}R_{p}&=\frac{C}{\tilde{\chi}^{\beta}\chi^{n+\alpha}}\int_{0}^{\infty}\int_{0}^{\infty}w^{\beta-1}t^{n+\alpha-1}\mathrm{e}^{-w-t}\sum_{p=P}^{\infty}\bigg{(}\frac{\chi
w}{\tilde{\chi}t}\bigg{)}^{P}\operatorname{d\\!}{}{w}\operatorname{d\\!}{}{t},\\\
&=\frac{C}{\tilde{\chi}^{\beta}\chi^{n+\alpha}}\frac{\chi^{P}}{\tilde{\chi}^{P}}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{\infty}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{\infty}w^{P+\beta-1}t^{n-P+\alpha-1}\mathrm{e}^{-w-t}\frac{1}{1-\frac{\chi
w}{\tilde{\chi}t}}\operatorname{d\\!}{}{w}\operatorname{d\\!}{}{t}.\end{split}$

In the above, we have used the principle value integral due to the simple pole
at $\chi w=\tilde{\chi}t$ introduced by evaluating the infinite geometric
series. Next, we change integration variables from $(w,t)$ to $(s,f)$ where
$s=\chi w/(\tilde{\chi}t)$ and $f=t+w$. This substitution allows for the
independent evaluation of the two integrals, one of which is recognised as the
definition of the Gamma function. Performing this substitution, and
simplifying, yields

(34)
$\begin{split}R_{p}&=\frac{C}{\chi^{n+\alpha+\beta}}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{\infty}\frac{s^{P+\beta-1}}{\big{(}1+\frac{s\tilde{\chi}}{\chi}\big{)}^{n+\alpha+\beta}(1-s)}\operatorname{d\\!}{}{s}\int_{0}^{\infty}f^{n+\alpha+\beta-1}\mathrm{e}^{-f}\operatorname{d\\!}{}{f},\\\
&=C\frac{\Gamma(n+\alpha+\beta)}{\chi^{n+\alpha+\beta}}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{\infty}\frac{s^{P+\beta-1}}{\big{(}1+\frac{s\tilde{\chi}}{\chi}\big{)}^{n+\alpha+\beta}(1-s)}\operatorname{d\\!}{}{s}.\end{split}$

Next we substitute for $x=s-1$, which shifts the location of the pole from
$s=1$ to $x=0$. As the dominant component of the integrand at $x=0$ lies
inside the new range of integration, $-1\leq x<\infty$, we may also extend the
lower range to $-\infty$ without changing the dominant asymptotic behaviour.
This gives

(35)
$\begin{split}R_{p}&=-C\frac{\Gamma(n+\alpha+\beta)}{\chi^{n+\alpha+\beta}}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{-1}^{\infty}\frac{{(1+x)}^{P+\beta-1}}{x\big{[}1+(1+x){\tilde{\chi}}/{\chi}\big{]}^{n+\alpha+\beta}}\operatorname{d\\!}{}{x},\\\
&\sim-C\frac{\Gamma(n+\alpha+\beta)}{(\chi+\tilde{\chi})^{n+\alpha+\beta}}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{-\infty}^{\infty}\frac{1}{x}\mathrm{e}^{(P+\beta-1)\log{(1+x)}-(n+\alpha+\beta)\log{[1+x\tilde{\chi}/(\chi+\tilde{\chi})]}}\operatorname{d\\!}{}{x}.\end{split}$

We now seek to approximate the integral in (35) above, both near the HOSL and
close to the pole at $x=0$. In writing

(36)
$\chi=r\mathrm{e}^{\mathrm{i}\vartheta}\qquad\text{and}\qquad\tilde{\chi}=\tilde{r}\mathrm{e}^{\mathrm{i}\tilde{\vartheta}},$

the HOSL, for which $\text{Im}[\tilde{\chi}/\chi]=0$ is located at
$\tilde{\vartheta}-\vartheta=0$. Thus, in the vicinity of this, and the pole
at $x=0$, both $x$ and $\tilde{\vartheta}-\vartheta$ are small. We may
therefore Taylor expand these components in (35) to find

(37)
$\begin{split}\mathrm{e}^{(P+\beta-1)\log{(1+x)}}&\sim\mathrm{e}^{(P+\beta-1)\big{(}x-\tfrac{x^{2}}{2}+O(x^{3})\big{)}},\\\
\mathrm{e}^{-(n+\alpha+\beta)\log{\big{[}1+\tfrac{x\tilde{\chi}}{(\chi+\tilde{\chi})}\big{]}}}&\sim\mathrm{e}^{-(n+\alpha+\beta)\big{[}\tfrac{\tilde{r}x}{r+\tilde{r}}-\tfrac{\tilde{r}^{2}x^{2}}{2(r+\tilde{r})^{2}}+\mathrm{i}\tfrac{r\tilde{r}x(\tilde{\vartheta}-\vartheta)}{(r+\tilde{r})^{2}}+O(x(\tilde{\vartheta}-\vartheta)^{2};x^{3})\big{]}}.\end{split}$

The remainder, $R_{P}$, is minimised when the largest components of (37),
$\exp{(Px)}$ and $\exp{(-\frac{\tilde{r}xn}{r+\tilde{r}})}$, cancel upon
multiplication. This motivates the specification of our optimal truncation
point by

(38) $P=\frac{\tilde{r}n}{r+\tilde{r}}+\tilde{\rho},$

where $0\leq\tilde{\rho}<1$ ensures that $P$ takes integer values. We may now
verify a priori that our prior assumption of $n-P$ being large is self
consistent, since $n-P\sim rn/(r+\tilde{r})\gg 0$. Substitution of the optimal
value of $P$ from (38) into the integral (35), and using expansions (37),
yields real and imaginary components in the integrand. Analogous to the
methodology used by Berry (1989), it may be verified through direct evaluation
that the real component of this integral is subdominant to the imaginary
component. The leading order imaginary component then simplifies to give

(39)
$\begin{split}R_{P}&\sim\mathrm{i}C\frac{\Gamma(n+\alpha+\beta)}{(\chi+\tilde{\chi})^{n+\alpha+\beta}}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{-\infty}^{\infty}\frac{1}{x}\mathrm{e}^{-\tfrac{r\tilde{r}nx^{2}}{2(r+\tilde{r})^{2}}}\sin{\bigg{(}\frac{r\tilde{r}(\tilde{\vartheta}-\vartheta)n}{(r+\tilde{r})^{2}}x\bigg{)}}\operatorname{d\\!}{}{x},\\\
&\sim\sqrt{2\pi}\mathrm{i}C\frac{\Gamma(n+\alpha+\beta)}{(\chi+\tilde{\chi})^{n+\alpha+\beta}}\int_{0}^{\tfrac{(r\tilde{r})^{1/2}(\tilde{\vartheta}-\vartheta)n^{1/2}}{(r+\tilde{r})}}\mathrm{e}^{-t^{2}/2}\operatorname{d\\!}{}{t}.\end{split}$

Equation (35) is the main result of this section: an explicit solution for the
smooth transition of $R_{P}$ across a boundary layer of width $O(n^{-1/2})$
surrounding the HOSL. Thus, we have that

(40) $R_{P}\sim
2\pi\mathrm{i}C(z)\frac{\Gamma(n+\alpha+\beta)}{(\chi+\tilde{\chi})^{n+\alpha+\beta}}$

switches on in the late-terms of the asymptotic solution across the HOSL,
which is given by

(41)
$\text{Im}\left[\frac{\tilde{\chi}}{\chi_{1}}\right]=0\qquad\text{and}\qquad\text{Re}\left[\frac{\tilde{\chi}}{\chi_{1}}\right]\geq
0.$

We have focused in the above on a concrete calculation that applies to the
simplest form of HOSP occurring in practice in singularly perturbed
differential equations. In fact, it is possible to recast the above
calculation in the hyperterminant language of Olde Daalhuis (1996, 1998) where
the integral (39) may be understood as a hyperterminant. We refer the reader
to the recent conference talk by Olde Daalhuis (2022) for a treatment of the
higher-order Stokes smoothing in this hyperterminant setting.

###  3.5 Stokes and higher-order Stokes line structure

The late-terms of the base expansion contain three distinct divergent
components. Firstly, there are the two factorial-over-power terms with
singulants $\chi_{1}=z^{2}/2$ and $\chi_{2}=z$ from (27). Secondly, there is
the component from (40) switched on across the HOSL, with
$\chi_{3}=\chi_{1}+\tilde{\chi}=z-1/2$. The Stokes phenomenon associated with
these three components produce Stokes lines, which we customarily denote as
$B>1$, $B>2$, and $B>3$, respectively, in order to express the exponential
dominance of one component over another as switching occurs. Furthermore,
examination of the $O(\mathrm{e}^{-\chi_{1}/\epsilon})$ transseries equation
(21b) for $y^{(1)}(z,\epsilon)$ reveals a divergent series whose solutions at
each order of $\epsilon$ are the same as the $B_{p}$ functions found in the
late-late-term divergence of the base expansion from section 3.3. This
relationship was discussed in section 2.1, and this switching results in a
term of $O(\mathrm{e}^{-(\chi_{1}+\tilde{\chi})/\epsilon})$. Thus, there will
also be a $1>3$ Stokes line generated by the divergence of $y^{(1)}$.

Combined, the Stokes lines for this problem are given by

(42a) $\displaystyle(B>1)$
$\displaystyle:\qquad\text{Im}[\chi_{1}]=0,\qquad\text{Re}[\chi_{1}]\geq 0,$
(42b) $\displaystyle(B>3)$
$\displaystyle:\qquad\text{Im}[\chi_{3}]=0,\qquad\text{Re}[\chi_{3}]\geq 0,$
(42c) $\displaystyle(1>3)$
$\displaystyle:\qquad\text{Im}[\tilde{\chi}]=0,\qquad~{}\text{Re}[\tilde{\chi}]\geq
0,$

where the $(B>3)$ Stokes line is active only when the higher-order Stokes
multiplier is non-zero. Note that we have not included the $B>2$ Stokes line
here, as it can be shown through a more detailed matching procedure (to all
orders in $\epsilon$) than Appendix A.1 that the constant prefactor of the
$\chi_{2}=z$ late-term component must equal zero. These are shown in figure 2.
For this example, the new $B>3$ Stokes line generated by the HOSP is necessary
to avoid a contradiction occurring from rotating about the simple pole at
$z=1$.

Figure 2: The Stokes line structure for our model problem (18) is shown in
$z\in\mathbb{C}$. HOSL are shown with thin lines, and regular Stokes lines
with Bold lines.

###  3.6 Borel-Padé and the higher-order Stokes phenomenon

The above calculation may also be interpreted from the Borel plane
perspective. As discussed in section 2.2, singularities in the Borel plane,
say at $w=\chi(z)$, correspond to the singulant functions that characterise
the divergence of late terms in the transseries. Moreover, the local behaviour
of $y_{B}$ near such singularities is determined from the constant, $\alpha$,
that appears in divergent representation (5). For instance, from the late-term
analysis in §3.2, the divergence associated with $\chi_{1}$ has $\alpha=1/2$,
and that associated with $\chi_{3}$ in (40) has $\alpha+\beta=1$. These
constants were determined in Appendix A. Thus, the parametric Borel transform
has a branch point at $w=\chi_{1}(z)$ and a pole at $w=\chi_{3}(z)$. The HOSL
is then the condition for the pole at $\chi_{3}(z)$ to be hidden (with respect
to rays from the distinguished origin) behind the branch cut from
$\chi_{1}(z)$. This geometric condition (11) was realised precisely as a
Stokes line in the late-terms from (41).

The singularity structure of $y_{B}$ can also be explored via numerical
analytic continuation techniques. The Padé approximant is a powerful method to
analytically continue truncated power series, and we refer the reader to the
work by Costin & Dunne (2022) (and references therein) for the state of the
art of using Padé approximants as analytic continuation tool. Here, we shall
demonstrate that Padé approximants can be used as a useful heuristic to
understand the parametric Borel plane singularity structure for a given
parametric series.

In our present example, we compute $2N$ perturbative terms and consider the
following truncated power series, for which $z$ is a parameter,

(43)
$y^{2N}_{B}(w;z)=y_{0}(z)+y_{1}(z)w+\cdots+\frac{y_{2N}(z)}{(2N)!}w^{2N}.$

We may then compute the off-diagonal $[N-1:N]$ Padé approximant for the fixed
values of $z\in\mathbb{C}_{z}$ of interest. That is, we write

(44) $P_{[N-1:N]}(w;z)=\frac{P_{N-1}(w;z)}{Q_{N}(w;z)},$

where $P_{N-1}(w;z)$ and $Q_{N}(w;z)$ are polynomials, whose coefficients
depend on $z$, chosen to agree with the truncated power series at a given
order. The interesting region to explore is near the HOSL where one may see
poles disappearing behind branch cuts, realising the HOSP.

Figure 3: The Borel plane, $\mathbb{C}_{w}$, is shown in ($a$) for
$z=1/2+\mathrm{i}$, and in ($b$) for $z=1/2$. The singularities (nodes)
correspond to $w=\chi_{1}(z)=z^{2}/2$ and $w=\chi_{3}(z)=z-1/2$, and are
computed from a Padé approximant with $N=250$. Note that the situation $(a)$
corresponds to the region ‘above’ the HOSL in figure 2 while $(b)$ corresponds
to ‘below’ the HOSL.

We demonstrate this in figure 3 for our example equation. In this figure we
note the characteristic coalescing poles that the rational approximation
assigns to branch points. We plot the Borel plane for two different values of
$z$: inside the HOSL for $z=1/2+\mathrm{i}$ in $(a)$, and outside for $z=1/2$
in $(b)$. One interesting fact to note is that in $(b)$ when $z=1/2$, the
Borel singularity $w=\chi_{3}(1/2)=0$; however the series $y_{B}$ is
convergent here as the singularity $w=0$ lies on a higher sheet above the
distinguished origin. The full parametric Alien calculus structure of this
example will be explored in future work Crew & Sauzin (2023).

##  4 Further examples of the higher-order Stokes phenomenon

We now consider three previously studied differential equations whose
asymptotic solutions display the HOSP. We demonstrate how this phenomenon may
be resolved in these problems through the application of the methods developed
in this paper, in which the late-late-term divergence of the asymptotic
solution is derived.

The singularly perturbed equations considered in this section are:

  1. (i)

The Pearcey equation in §4.1. This is a linear third-order homogeneous
equation, previously studied by Howls et al. (2004) through consideration of
an integral representation of the solution and a hyperterminant analysis. Our
results are summarised in 4.1 and full details appear in Appendix B.

  2. (ii)

A linear second-order inhomogeneous differential equation in §4.2. This
equation was proposed by Trinh & Chapman (2013), who considered an exact
series representation for the late-terms of the asymptotic solution.

  3. (iii)

A linear second-order differential equation with an eigenvalue in §4.3, which
arises in geophysical fluid dynamics. This equation, and the reusltant
eigenvalue divergence, was studied by Shelton et al. (2023a, b).

The first two examples contain new Stokes lines which emanate from
intersecting Stokes lines. In the third example the higher-order and regular
Stokes lines coincide, which changes the value of the regular Stokes
multiplier.

###  4.1 The Pearcey equation

The Pearcey function is perhaps the most recognisable subject in connection to
the HOSP; it is discussed in the work of Howls et al. (2004), primarily with
respect to the integral representation

(45)
$I(z)=\int_{C}\mathrm{e}^{-f(s;z)/\epsilon}\,\operatorname{d\\!}{}{s}\quad\text{with}\quad
f(s;z)=-\mathrm{i}\left(\frac{1}{4}s^{4}+\frac{1}{2}s^{2}+sz\right).$

Here, $\epsilon\to 0$ and $C$ is a contour from
$\infty\exp(-3\pi\mathrm{i}/8)$ to $\infty\exp(\pi\mathrm{i}/8)$.
Alternatively, manipulation of (45) yields a linear, third-order, homogeneous
differential equation:

(46) $\epsilon^{3}I^{\prime\prime\prime}(z)-\epsilon
I^{\prime}(z)-\mathrm{i}zI(z)=0,$

with decay conditions in the same limits as specified above. The relationship
between these two systems is discussed in section 7.2 of Takei (2017), who
considers the BNR equation from Berk et al. (1982), which is obtained by
rescaling (46). The general Pearcey integral (45) contains a second parameter,
assumed to be fixed in the above formulation (Takei, 2017).

Aided by the integral representation (45), the HOSP can be understood more
easily via a steepest descent analysis and this is presented by Howls et al.
(2004). Interestingly, we could not find a direct asymptotic analysis of the
HOSP – as connected to the late terms – of the Pearcey function in its
singularly-perturbed differential equation form. This we do now using the
techniques from the previous section.

Three different WKBJ solutions are permitted to (46), each indexed by a
leading exponential, $A_{i}(z,\epsilon)\mathrm{e}^{-S_{i}(z)/\epsilon}$ for
$i=1,2,3$, which yields divergent expansions of the form

(47) $\displaystyle I(z)$
$\displaystyle=\mathrm{e}^{-S_{1}/\epsilon}\bigg{[}A_{1}^{(0)}+\cdots+\epsilon^{n}\bigg{(}B_{1}^{(0)}\frac{\Gamma(n+\alpha)}{(S_{2}-S_{1})^{n+\alpha}}+\tilde{\sigma}_{1}C_{1}\frac{\Gamma(n+\alpha)}{(S_{3}-S_{1})^{n+\alpha}}\bigg{)}+\cdots\bigg{]}$
$\displaystyle\phantom{=}+\mathrm{e}^{-S_{2}/\epsilon}\bigg{[}A_{2}^{(0)}+\cdots+\epsilon^{n}\bigg{(}B_{2a}^{(0)}\frac{\Gamma(n+\alpha)}{(S_{3}-S_{2})^{n+\alpha}}+B^{(0)}_{2b}\frac{\Gamma(n+\alpha)}{(S_{1}-S_{2})^{n+\alpha}}\bigg{)}+\cdots\bigg{]}$
$\displaystyle\phantom{=}+\mathrm{e}^{-{S_{3}}/{\epsilon}}\bigg{[}A_{3}^{(0)}+\cdots+\epsilon^{n}\bigg{(}B_{3}^{(0)}\frac{\Gamma(n+\alpha)}{(S_{2}-S_{3})^{n+\alpha}}+\tilde{\sigma}_{3}C_{3}\frac{\Gamma(n+\alpha)}{(S_{1}-S_{3})^{n+\alpha}}\bigg{)}+\cdots\bigg{]}.$

The derivation of each of these components is given in Appendix B. Each of the
six factorial-over-power components of (47) induces the Stokes phenomenon on
another exponential, when the relevant singulant is real and positive.
However, two of these divergences are present only when the corresponding
higher-order Stokes multiplier, $\tilde{\sigma}_{1}$ or $\tilde{\sigma}_{3}$,
is non-zero. The higher-order Stokes lines, across which these higher-order
multipliers change, are also shown in figure 4.

Figure 4: The Stokes line structure is shown in the $z$-plane for the Pearcey
equation (46). Higher-order Stokes lines (HOSL), across which the late-terms
switch-on a new divergent component, are shown by thin lines. Stokes lines,
across which the regular Stokes phenomenon occurs on an exponentially-small
component of the asymptotic expansion, are shown as bold lines. On account of
the HOSP, two Stokes lines are born (or truncated) at the Stokes crossing
point. The notation of $i>j$ corresponds to exponential
$\mathrm{e}^{S_{i}/\epsilon}$ switching on $\mathrm{e}^{S_{j}/\epsilon}$.

The resultant Stokes line structure is shown in figure 4. Note that there are
two singular point at $z_{-}$ and $z_{+}$, which correspond both to
singularities in the leading-order asymptotic solution $A_{i}^{(0)}$, and
zeros of the singulant functions,
$S_{1}(z_{-})=S_{2}(z_{-})=S_{2}(z_{+})=S_{3}(z_{+})=0$. In total, four Stokes
lines emanate from $z_{-}$ and $z_{+}$ (shown thick). The Stokes lines
intersect at two Stokes crossing points on the real axis. These crossing
points also lie along the respective higher-order Stokes lines (shown thin),
across which the respective multiplier $\tilde{\sigma}$ switches on to a non-
zero value. Thus, beyond the Stokes crossing point, an additional divergences
are present in the late-terms of transseries (47), which leads to two new
Stokes lines. These new Stokes lines are marked $1>3$ and $3>1$ in the figure.

###  4.2 An inhomogeneous second-order equation

We now study the following second-order singularly perturbed differential
equation for $A(z)$:

(48) $\epsilon^{2}A^{\prime\prime}+2\epsilon A^{\prime}+(1-z)A=\frac{1}{z-a}.$

Equation (48) was considered as a model problem by Trinh & Chapman (2013) due
to the turning point at $z=1$ and the singularity at $z=a$. In fact, each
order of the asymptotic solution to (48) may be written exactly as a
polynomial in powers of $(1-z)^{-1}$ and $(z-a)^{-1}$, for which the constant
coefficients are determined as the solution to a recurrence relation. Trinh &
Chapman (2013) conjectured that the analysis of the late-orders behaviour of
the recurrence relations could yield details of the HOSP, but this was left as
an open challenge.

We now apply our general techniques to detect HOSP in this problem. In
expanding the solution to equation (48) as $A(z,\epsilon)=A_{0}(z)+\epsilon
A_{1}(z)+\cdots$, the first two orders of this asymptotic expansion may be
determined as

(49) $A_{0}=\frac{1}{(1-z)(z-a)}\quad\text{and}\quad
A_{1}=\frac{2}{(1-z)^{2}(z-a)^{2}}-\frac{2}{(1-z)^{3}(z-a)}.$

At $O(\epsilon^{n})$ in equation (48), we find
$A_{n-2}^{\prime\prime}+2A_{n-1}^{\prime}+(1-z)A_{n}=0$, which holds for
$n\geq 2$. The late-terms, $A_{n}$, of the asymptotic expansion will diverge
as $n\to\infty$ on account of the singularities at $z=1$ and $z=a$ in (49).

#### Late-term divergence

In order to study the divergence of $A_{n}$, we consider a factorial-over-
power ansatz of the form

(50)
$A_{n}(z)\sim\sum_{p=0}^{\infty}B_{p}(z)\frac{\Gamma(n-p+\alpha)}{[\chi(z)]^{n-p+\alpha}}\qquad\text{as
$n\to\infty$},$

where $\alpha$ is a constant, and $\chi$ and $B_{p}$ are the singulant and
amplitude functions, respectively. Substitution of ansatz (50) into the
$O(\epsilon^{n})$ equation yields at each order of
$\Gamma(n-p+\alpha)/\chi^{n-p+\alpha}$ the equations

(51a) $\displaystyle(\chi^{\prime})^{2}-2\chi^{\prime}+(1-z)$
$\displaystyle=0,$ (51b) $\displaystyle
2(1-\chi^{\prime})B_{0}^{\prime}-\chi^{\prime\prime}B_{0}$ $\displaystyle=0,$
(51c) $\displaystyle
B_{p-2}^{\prime\prime}+2(1-\chi^{\prime})B_{p-1}^{\prime}-\chi^{\prime\prime}B_{p-1}$
$\displaystyle=0.$

Here, (51c) holds for $p\geq 2$. The quadratic (51a) for $\chi^{\prime}$, and
(51b) for $B_{0}$ may be solved to find

(52) $\chi_{\pm}^{\prime}(z)=1\pm z^{1/2}\qquad\text{and}\qquad
B_{0}(z)=\frac{\Lambda_{0}}{(1-\chi_{\pm}^{\prime})^{1/2}},$

where $\pm$ denotes each of the two solutions for $\chi^{\prime}$, and
$\Lambda_{0}$ is a constant of integration. From the singulant solution in
(52) above we have $1-\chi^{\prime}_{\pm}=\mp z^{1/2}$, and thus the solution
for $B_{0}$ is singular at $z=0$. This new singularity at $z=0$ will result in
the divergence of $B_{p}$ as $p\to\infty$, which is now studied.

#### Late-late-term divergence

To capture the divergence of $B_{p}$ as $p\to\infty$, we consider a typical
factorial-over-power ansatz [cf. (30)]. Substitution of this into equation
(51c) yields $\tilde{\chi}^{\prime}_{\pm}(z)=2(1-\chi_{\pm}^{\prime})$, and
hence $\tilde{\chi}_{\pm}(z)=\mp 4z^{3/2}/3$ upon imposing the condition
$\tilde{\chi}_{\pm}(0)=0$. Combined with the solution for the prefactor
equation, we find a late-late-term divergence of the form

(53)
$B_{p}(z)\sim\frac{\tilde{\Lambda}}{(\chi^{\prime}_{\pm}-1)^{1/2}}\frac{\Gamma(p+\beta)}{[\mp
4z^{3/2}/3]^{p+\beta}}\qquad\text{as $p\to\infty$},$

where $\tilde{\Lambda}$ is a constant of integration. This will induce the
HOSP in which a new late-term component with singulant $\chi_{\pm}\mp
4z^{3/2}/3$ switches on when $\mp 4z^{3/2}/(3\chi_{\pm})$ is real and
positive.

#### Stokes line structure

The naive Stokes lines may be determined by integrating the solution
$\chi_{\pm}^{\prime}(z)=1\pm z^{1/2}$ from (52). As there are two singular
points, $z=a$ and $z=1$, that generate divergence, we will have four initial
singulants. Only the two generated by the boundary condition $\chi(a)=0$ are
now considered to rectify the issue of crossing Stokes lines found by Trinh &
Chapman (2013). This yields

(54)
$\chi_{1}(z)=z+\frac{2}{3}z^{3/2}-a-\frac{2}{3}a^{3/2}\quad\text{and}\quad\chi_{2}(z)=z-\frac{2}{3}z^{3/2}-a+\frac{2}{3}a^{3/2}.$

Furthermore, the expansion in $n$ for the amplitude functions of the above
four singulants diverges in the manner specified by (53), on account of a
singularity at $z=0$. This new divergent series induces the HOSP in which a
new factorial-over-power contribution to $A_{n}$, with a singulant given by
$\chi+\tilde{\chi}$, switches on across the HOSL
$\text{Im}[\tilde{\chi}/\chi]=0$ and $\text{Re}[\tilde{\chi}/\chi]>0$. The two
new singulants switched on in $A_{n}$ are therefore

(55)
$\chi_{3}(z)=z-\frac{2}{3}z^{3/2}-a-\frac{2}{3}a^{3/2}\quad\text{and}\quad\chi_{4}(z)=z+\frac{2}{3}z^{3/2}-a+\frac{2}{3}a^{3/2}.$

The new singulants in (55), generated through the HOSP, are equal to zero at
locations that correspond to regular points of the differential equation (48).
This was commented upon by (Trinh & Chapman, 2013, p. 418):

> _The curiosity, however, is that these singularities do not appear anywhere
> in the base series, and so they do not seem to be associated with any
> eventual divergence._

As we know now, the divergent contributions of $A_{n}$ with singulants (55)
will switch-off across the HOSL; consequently, they do not generate unexpected
singular behaviour at the apparent singularity, itself. However, the
corresponding singulants still generate divergence in the asymptotic
expansions, within other regions where they are indeed present.

Figure 5: The Stokes line and HOSL structure is shown near the singularity
(circle) at $z=a$ for $a=-1+\mathrm{i}$. The turning point at $z=0$ (circle),
generates a $1>3$ Stokes line as it is a singularity in the
$O(\mathrm{e}^{-\chi_{1}/\epsilon})$ transseries equation. Thick lines
represent Stokes lines, across which exponentially-small terms in the
asymptotic solution switch on. The thin line is a HOSL, across which a new
divergent contribution switches on in the late-terms, $A_{n}$.

In figure 5, we plot the Stokes lines and HOSL associated with the singularity
at $z=a$. Two naive Stokes lines, which we denote by $\text{B}>1$ and
$\text{B}>2$, are generated by (54). The Airy Stokes line found by evaluating
Dingles conditions on $\tilde{\chi}=\mp 4z^{3/2}/3$, for which $1>3$, is also
shown. This results in two Stokes crossing points. The first, between the
$1>3$ and $\text{B}>2$ Stokes lines, is self consistent. However the second,
between $1>3$ and $\text{B}>1$, results in a contradiction. The new divergent
component with $\chi_{3}$ generates a $\text{B}>3$ new Stokes line emanating
from the SCP, which resolves this inconsistency.

###  4.3 The equatorial Kelvin wave

Our last example of the HOSP concerns the study of travelling equatorial
Kelvin waves with small latitudinal shear. Here, the system is governed by a
second-order singularly perturbed differential equation:

(56)
$\displaystyle\epsilon^{2}\bigg{[}(Y-c-1)(Y-c+1)^{2}+(Y-c-1)^{2}(Y-c+1)\bigg{]}V^{\prime\prime}$
$\displaystyle-2Y\bigg{[}(Y-c-1)(Y-c+1)^{2}+(Y-c+1)(Y-c-1)^{2}+\frac{2\epsilon^{2}(Y-c)^{2}}{Y}\bigg{]}V^{\prime}$
$\displaystyle+\bigg{[}\big{[}2Y(Y-c-1)+\epsilon^{2}\big{]}(Y-c+1)^{2}+(2c-2-\epsilon^{2})(Y-c-1)^{2}\bigg{]}V=0.$

The unknown is $V(Y)$, representing the eigenfunction amplitude, corresponding
to the vertical velocity, and $\epsilon$ is a non-dimensional shear. The
wavespeed is $c$, an eigenvalue whose value is obtained from enforcing the
boundary condition $V(0)=1$. For further details of this model, see Shelton et
al. (2023b) and Griffiths (2008). We consider a solution expansion of the form
$V(Y)=V_{0}+\epsilon^{2}V_{1}+\epsilon^{4}V_{2}+\cdots$. Enforcing the
boundary condition at each order of $\epsilon$ is only possible if the
eigenvalue is also expanded as
$c=c_{0}+\epsilon^{2}c_{1}+\epsilon^{4}c_{2}+\cdots$.

At leading order, once the boundary condition $V_{0}(0)=1$ is satisfied, the
solutions are given by

(57) $V_{0}(Y)=(1-Y)^{\frac{1}{2}}\mathrm{e}^{Y/2}\qquad\text{and}\qquad
c_{0}=1.$

This singularity at $Y=1$ results in a divergent series expansion.
Furthermore, the eigenvalue expansion
$c=1+\frac{1}{2}\epsilon^{2}-\frac{1}{8}\epsilon^{4}+\cdots$ will also
diverge.

####  4.3.1 Late-term divergence

Due to the nonlinearity between the solution, $V$, and eigenvalue, $c$, in
equation (56), the $O(\epsilon^{2n})$ equation will contain a countably
infinite number of terms as $n\to\infty$. However, only those appearing at
leading order as $n\to\infty$ will be required to derive the HOSP. We do not
provide this here, and refer to reader to Shelton et al. (2023b) for this.

The $O(\epsilon^{2n})$ equation contains both homogeneous terms, for instance
$V_{n}$ and $V_{n-1}^{\prime}$ and inhomogenous terms such as $c_{n}$. Both of
these contribute to the divergent solution, which contains one term associated
with the divergence of $c_{n}$, and another generated by the singular
behaviour at $Y=1$. This yields

(58) $V_{n}(Y)\sim
B_{0}^{(c_{n})}(Y)\Gamma(n+\alpha)+B_{0}(Y)\frac{\Gamma(n+\gamma)}{(1-Y^{2})^{n+\gamma}},$

where $\alpha$ and $\gamma$ are constants, and the singulant function
$\chi=1-Y^{2}$ satisfies the condition $\chi(1)=0$. Similar to the example
from §4.2, there is an apparent paradox in that a singularity is seemingly
present at $Y=-1$ in the factorial-over-power representation (58), which is
absent from the early orders of the expansion. This may be resolved by
considering lower-order components of $B_{0}^{(c_{n})}$, for which late-late-
term divergence and the HOSP emerges.

####  4.3.2 Late-late-term divergence

In the study by Shelton et al. (2023b), it was assumed that the HOSP occurred,
and that expression (58) would switch off across a HOSL across the imaginary
axis. We now demonstrate how this may be determined through the application of
the methods developed in this present paper.

The late-terms of the eigenvalue expansion will diverge in the factorial-over-
power manner of $c_{n}\sim\delta\Gamma(n+\alpha)$, where $\delta$ and $\alpha$
are constants. The particular solution generated by this eigenvalue divergence
can then be determined with a late-term ansatz of the form

(59)
$V_{n}(Y)\sim\sum_{p=0}^{\infty}B^{(c_{n})}_{p}(Y){\Gamma(n-p+\alpha)}\qquad\text{as
$n\to\infty$}.$

In substituting ansatz (59) into the $O(\epsilon^{2n})$ equation and solving
the resultant equations, we have $B_{0}^{(c_{n})}\sim 1/Y$ as $Y\to 0$. This
singularity at $Y=0$ will result in late-late-term divergence as $p\to\infty$,
for which we consider an ansatz of the form

(60) $B_{p}^{(c_{n})}(Y)\sim
C(Y)\frac{\Gamma(p+\beta)}{[\tilde{\chi}(Y)]^{p+\beta}}\qquad\text{as
$p\to\infty$}.$

The requirement that $\tilde{\chi}(0)=0$ then yields $\tilde{\chi}(Y)=-Y^{2}$.
Thus, across the HOSL $\text{Im}[-Y^{2}/1]=0$ and $\text{Re}[-Y^{2}/1]>0$, a
new factorial-over-power component of $V_{n}$ with a singulant given by
$\chi=1-Y^{2}$ will switch on. This is the same singulant as that generated by
the singularity at $Y=1$. Thus, the apparent contradiction of $\chi=1-Y^{2}$
containing an undesired singularity at $Y=-1$ is resolved, since this
expression switches off as we cross the HOSL.

Figure 6: The Stokes line structure is shown for the Kelvin wave equation
(56). The regular Stokes line for $\chi=1-Y^{2}$, shown bold, lies along the
imaginary axis and between $0$ and $1$ on the real-axis. There is also a HOSL
along the imaginary axis, which causes the regular Stokes line on the negative
real axis to be irrelevant.

##  5 Discussion

Although somewhat obscure, the higher-order Stokes phenomena is generic in the
analysis of singular perturbation problems. As we have noted, it accompanies
situations where the late-term approximation of a $\epsilon$-expansion, say
$y_{n}$, is itself divergent in inverse powers of $n$. Consequently Stokes
phenomena occurs (in $n$), switching on or off late terms and hence Stokes
lines themselves. Alternatively, the genericity of HOSP can be interpreted as
a consequence of the Borel plane possessing a branched structure. The HOSP
often arises in situations where: (i) Stokes lines cross; (ii) additional
singularities, not present in early orders, are predicted in the late terms;
or (iii) a naive divergent series is unable to satisfy necessary boundary
conditions.

The methods developed in this paper have been presented for linear ODEs.
However, they are also applicable to more general nonlinear problems and PDEs.
Our approach is similar, and in some cases more general, than those of Chapman
& Mortimer (2005); Body et al. (2005); Howls et al. (2004); in ours, we have
focused on the analysis of the $1/n$ divergent series of the late terms and
its subsequent Stokes phenomena. Additionally, we have shown that the HOSP and
second-generation Stokes smoothing are both caused by the same type of
divergence, either in the late-late-terms, or found in the subsequent
transseries expansion. Thus, any HOSP that _e.g._ generates a $B>3$ Stokes
line with $\chi_{3}=\tilde{\chi}+\chi_{1}$ must necessarily accompany a
second-generation $1>3$ Stokes line that switches on $\chi_{3}$.

One of the main complications with nonlinearity is the determination of the
constants of integration, $\Lambda$ and $\tilde{\Lambda}$, of the amplitude
functions from the late- and late-late-term ansatzes. For our linear examples,
the series expansion for the outer limit of the inner solution could be
determined explicitly. This will not be the case, in general, for nonlinear
problems, where the associated recurrence relation is more complicated. In the
study of PDEs, integration of the singulant and amplitude equations from the
late- and late-late-term ansatzes can be challenging (Chapman & Mortimer,
2005; Body et al., 2005). The singularity structure in the Borel plane will
also be more complicated and one expects a more intricate Riemann sheet
structure.

In addition to the development of methodology, another contribution of our
work has been to highlight the beautifully simple example (1). Similar to the
example given in Olde Daalhuis (2004), such cases are reminders of the subtle
and complicated nature of divergence and transseries, even in apparently
benign problems. There are a number of exciting ongoing works related to HOSP
that have emerged from a recent Isaac Newton Institute programme, notably by
Olde Daalhuis (2022) on smoothing of HOSP, King (2022) on coincident Stokes
lines and HOSP, and Nemes (2022) on atypical Stokes multipliers due to HOSP.
This is an exciting and advanced area of research with still more to discover.

Acknowledgements. We thank the Isaac Newton Institute for Mathematical
Sciences, Cambridge, for the programme “Applicable resurgent asymptotics:
towards a universal theory”, during which many insightful discussions took
place. In particular, we are thankful for discussions with S. J. Chapman
(Oxford), G. Nemes (Tokyo Metropolitan University), A. Olde Daalhuis
(Edinburgh), and many other participants of the programme.

## References

  * Aoki (1994) Aoki, T. 1994 New turning points in the exact WKB analysis for higher-order ordinary differential equations. Analyse algébrique des perturbations singulières. I .
  * Aoki et al. (1998) Aoki, T., Kawai, T. & Takei, Y. 1998 On the exact WKB analysis for the third order ordinary differential equations with a large parameter. Asian J. Math. 2 (4), 625–640.
  * Aoki et al. (2001) Aoki, T., Kawai, T. & Takei, Y. 2001 On the exact steepest descent method: A new method for the description of Stokes curves. J. Math. Phys. 42 (8), 3691–3713.
  * Berk et al. (1982) Berk, H. L., Nevins, W. M. & Roberts, K. V. 1982 New Stokes’ line in WKB theory. J. Math. Physics 23 (6), 988–1002.
  * Berry (1989) Berry, M. V. 1989 Uniform asymptotic smoothing of Stokes discontinuities. Proc. R. Soc. Lond. A 422, 7–21.
  * Body et al. (2005) Body, G. L., King, J. R. & Tew, R. H. 2005 Exponential asymptotics of a fifth-order partial differential equation. Eur. J. Appl. Math. 16 (5), 647–681.
  * Chapman & Mortimer (2005) Chapman, S. J. & Mortimer, D. B. 2005 Exponential asymptotics and Stokes lines in a partial differential equation. Proc. R. Soc. Lond. A 461 (2060), 2385–2421.
  * Costin & Dunne (2022) Costin, O. & Dunne, G. V. 2022 Uniformization and constructive analytic continuation of Taylor series. Commun. Math. Phys. 392 (3), 863–906.
  * Crew & Sauzin (2023) Crew, S. & Sauzin, D. 2023 Parametric resurgence and piecewise stokes constants. In preparation .
  * Crew & Trinh (2023) Crew, S. & Trinh, P. H. 2023 Resurgent aspects of applied exponential asymptotics. arXiv preprint arXiv:2208.07290 .
  * Dingle (1973) Dingle, R. B. 1973 Asymptotic Expansions: Their Derivation and Interpretation. Academic Press, London.
  * Griffiths (2008) Griffiths, S. D. 2008 The limiting form of inertial instability in geophysical flows. J. Fluid Mech. 605, 115–143.
  * Honda (2007) Honda, N. 2007 On the Stokes geometry of the Noumi-Yamada system. Algebraic, Analytic and Geometric Aspects of Complex Differential Equations and their Deformations. Painleve Hierarchies 2, 45–72.
  * Honda et al. (2015) Honda, N., Kawai, T. & Takei, Y. 2015 Virtual turning points, , vol. 4. Springer.
  * Howls et al. (2004) Howls, C. J., Langman, P. J. & Daalhuis, A. B. Olde 2004 On the higher-order Stokes phenomenon. Proc. R. Soc. Lond. A 460, 2285–2303.
  * King (2022) King, J. 2022 Exponential asymptotics, the Stokes phenomenon and the higher-order Stokes phenomenon in some linear partial differential equations. Issac Newton Institute Applicable resurgent asymptotics: summary workshop presentation 16th Dec.
  * King & Chapman (2001) King, J. R. & Chapman, S. J. 2001 Asymptotics beyond all orders and Stokes lines in nonlinear differential-difference equations. Eur. J. Appl. Math. 12 (4), 433–463.
  * Lustri et al. (2019) Lustri, C. J., Pethiyagoda, R. & Chapman, S. J. 2019 Three-dimensional capillary waves due to a submerged source with small surface tension. J. Fluid Mech. 863, 670–701.
  * Mitschi et al. (2016) Mitschi, C., Sauzin, D., Loday-Richaud, M. & Delabaere, É. 2016 Divergent series, summability and resurgence. Springer.
  * Nemes (2022) Nemes, G. 2022 Dingle’s final main rule, Berry’s transition, and Howls’ conjecture. J. Phys. A-Math. Theor. 55 (49), 494001.
  * Olde Daalhuis (1996) Olde Daalhuis, A. B. 1996 Hyperterminants I. J. Comp. Appl. Math. 76 (1-2), 255–264.
  * Olde Daalhuis (1998) Olde Daalhuis, A. B. 1998 Hyperterminants II. J. Comp. Appl. Math. 89 (1), 87–95.
  * Olde Daalhuis (2004) Olde Daalhuis, A. B. 2004 On higher-order Stokes phenomena of an inhomogeneous linear ordinary differential equation. J. Comp. Appl. Math. 169 (1), 235–246.
  * Olde Daalhuis (2022) Olde Daalhuis, A. B. 2022 Smoothing for the higher-order Stokes phenomenon. Issac Newton Institute Applicable resurgent asymptotics: summary workshop presentation 16th Dec.
  * Shelton et al. (2023a) Shelton, J., Chapman, S. J. & Trinh, P. H. 2023a Pathological exponential asymptotics for a model problem of an equatorially trapped Rossby wave. In Preperation .
  * Shelton et al. (2023b) Shelton, J., Griffiths, S., Chapman, S. J. & Trinh, P. H. 2023b On the exponentially-small instability of the equatorial Kelvin wave. In preparation .
  * Shudo (2008) Shudo, A. 2008 A role of virtual turning points and new Stokes curves in Stokes geometry of the quantum Hénon map. Algebraic analysis of differential equations. Springer.
  * Takei (2017) Takei, Y. 2017 WKB analysis and Stokes geometry of differential equations. In Analytic, algebraic and geometric aspects of differential equations, pp. 263–304. Springer.
  * Trinh & Chapman (2013) Trinh, P. H. & Chapman, S. J. 2013 New gravity–capillary waves at low speeds. Part 2. Nonlinear geometries. J. Fluid Mech. 724, 392–424.

##  Appendix A Inner analyses near $z=0$ and $z=1$

The purpose of the section is to determine the unknown constants appearing in
the late-terms, $y_{n}$, of the asymptotic solution. We note that, in addition
to $\alpha$ and $\beta$ appearing in the late- and late-late-term ansatzes
(25) for $y_{n}$ as $n\to\infty$ and (30) for $B_{p}$ as $p\to\infty$, each of
$B_{p}$ for $p\geq 0$ contains an unknown constant of integration,
$\Lambda_{p}$. However, explicit knowledge of all of these constants,
$\Lambda_{p}$, is not required. Firstly, we require $\Lambda_{0}$ which
appears in the Stokes switching induced by the $\chi_{1}=z^{2}/2$ late-term
component. Secondly, there is the Stokes switching induced by the
$\chi_{1}+\tilde{\chi}=z-1/2$ late-term component generated through the HOSP;
this requires knowledge of $B_{p}$ as $p\to\infty$, for which the amplitude
solution (31) contains the constant of integration $\tilde{\Lambda}$.

###  A.1 Inner analysis near $z=0$

We begin by deriving the inner equation near $z=0$, near to which the early
orders of expansion (22) reorder. This is due to a turning point in the
original differential equation that generates a singularity in the leading
order solution, $y_{0}^{(0)}(z)=z^{-1}$. Since $y_{0}^{(0)}\sim z^{-1}$ and
$\epsilon y_{1}^{(0)}\sim\epsilon z^{-3}$ as $z\to 0$ from equation (23), the
early orders of the expansion reorder when $z\sim\epsilon^{1/2}$. To study the
solution in this region, we introduce the inner variable, $\hat{z}$, and the
inner solution, $\hat{y}$, by

(61) $z=\epsilon^{1/2}\hat{z}\qquad\text{and}\qquad
y_{\text{outer}}=\frac{1}{(\epsilon^{1/2}\hat{z})}\hat{y}_{\text{inner}}(\hat{z}).$

The inner region is defined by $\hat{z}=O(1)$, and the precise definition of
the inner variable in (61) ensures that $\hat{y}_{\text{inner}}=O(1)$ in the
inner region. Next, we derive the inner equation by substituting equations
(61) into the outer equation (21a). As the inner limit of the outer divergent
solution, $\epsilon^{n}y_{n}$, will be seen in (67) to reorder into the
leading-order inner solution, only the leading-order, in $\epsilon$, inner
equation need be considered. This yields

(62)
$\frac{1}{\hat{z}}\frac{\operatorname{d\\!}{}\hat{y}}{\operatorname{d\\!}{}\hat{z}}+\bigg{(}1-\frac{1}{\hat{z}^{2}}\bigg{)}\hat{y}=1.$

####  A.1.1 Inner solution

To motivate the correct form to take for the inner expansion of
$\hat{y}(\hat{z})$, we take the inner limit of the early orders of the outer
expansion from (23) by substituting for $z=\epsilon^{1/2}\hat{z}$ from (61)
and expanding as $\epsilon\to 0$. This yields a 1-term inner limit of 3-terms
of the outer-solution, given by $\hat{y}\sim
1+{\hat{z}^{-2}}+3{\hat{z}^{-4}}+\cdots$. Furthermore, as knowledge of the
inner solution of equation (62) is required only in the limit of
$\hat{z}\to\infty$ in order to match with the outer solution, we will consider
a series solution for $\hat{y}(\hat{z})$ of the form

(63) $\hat{y}(\hat{z})=\sum_{n=0}^{\infty}\frac{\hat{a}_{n}}{\hat{z}^{2n}}.$

Substitution of series solution (63) into the inner equation (62) yields at
each order of $\hat{z}$ the equations

(64) $\hat{a}_{0}=1\qquad\text{and}\qquad\hat{a}_{n}=(2n-1)\hat{a}_{n-1},$

for $n\geq 1$. Typically, and usually for nonlinear problems, this inner
recurrence relation must be solved numerically to a reasonably large value of
$n\approx 100$, for which matching with the outer solution yields an
approximate value for the constants of integration in the amplitude function
of the factorial-over-power ansatz. However, recurrence relation (64) has the
exact solution

(65) $\hat{a}_{n}=\frac{\Gamma(2n)}{2^{n-1}\Gamma(n)}.$

####  A.1.2 Inner-outer matching

We now match with the outer divergent solution by writing the inner solution
(63) in terms of the outer variable, $z$, and outer solution, $y$, from
equation (61). This outer limit of the inner solution is given by

(66)
$y=\frac{1}{\epsilon^{1/2}\hat{z}}\sum_{n=0}^{\infty}\frac{\hat{a}_{n}}{\hat{z}^{2n}}=\sum_{n=0}^{\infty}\epsilon^{n}\frac{\hat{a}_{n}}{z^{2n+1}}.$

We match the $O(\epsilon^{n})$ component of this outer limit with the
$O(\epsilon^{n})$ component of the inner limit (as $z\to 0$) of the outer
divergent solution, which is given by

(67)
$y_{n}=\frac{\epsilon^{n}\Lambda_{0}}{1-z}\frac{\Gamma(n+\alpha)}{(z^{2}/2)^{n+\alpha}}\sim\Lambda_{0}\frac{\Gamma(n+\alpha)}{(z^{2}/2)^{n+\alpha}}.$

Matching the $O(\epsilon^{n})$ component of (66) with (67) requires
$\alpha=1/2$, which then gives an expression for the constant, $\Lambda_{0}$.
Thus, we have found

(68)
$\alpha=1/2\qquad\text{and}\qquad\Lambda_{0}=\lim_{n\to\infty}\bigg{(}\frac{\Gamma(2n)}{2^{2n-1/2}\Gamma(n)\Gamma(n+1/2)}\bigg{)}=\frac{1}{\sqrt{2\pi}}.$

###  A.2 Inner analysis near $z=1$

The constants $\beta$ and $\tilde{\Lambda}$ are now determined by matching
with an inner solution at $z=1$. These constants appear in the divergent form
of $B_{p}$ as $p\to\infty$. One method is to consider a boundary layer of
diminishing width as $n\to\infty$ on account of the late-term reordering
$B_{0}\Gamma(n+\alpha)/\chi^{n+\alpha}\sim
B_{1}\Gamma(n+\alpha-1)/\chi^{n+\alpha-1}$. However, we note that a power
series expansion,

(69) $y^{(1)}(z,\epsilon)=\sum_{p=0}^{\infty}\epsilon^{p}B_{p}(z),$

in the $\chi=z^{2}/2$ transseries equation (21b) yields at each order of
$\epsilon$ the same equations as those satisfied by $B_{p}$ in (26b) and
(26c). Thus, we may determine the unknown constants in the divergent form of
$B_{p}$ through the consideration of a boundary layer, of diminishing width as
$\epsilon\to 0$, at $z=1$ in the governing equation (21b) for $y^{(1)}$.

Since the early orders of expansion (69) reorder when
$\Lambda_{0}(1-z)^{-1}\sim-\epsilon\Lambda_{0}(1-z)^{-3}$, and
$y^{(1)}\sim\Lambda_{0}/(1-z)$, we introduce the inner variable, $\hat{z}$ and
inner solution, $\hat{y}$, via

(70) $(1-z)=\epsilon^{1/2}\hat{z}\qquad\text{and}\qquad
y^{(1)}=\frac{\Lambda_{0}}{\epsilon^{1/2}\hat{z}}\hat{y}(\hat{z}).$

The leading-order inner equation for $\hat{y}$ may be found by substituting
relations (70) into the transseries equation (21b) for $y^{(1)}$, and
retaining the dominant terms as $\epsilon\to 0$. This yields

(71)
$\frac{1}{\hat{z}}\frac{\operatorname{d\\!}{}^{2}\hat{y}}{\operatorname{d\\!}{}\hat{z}^{2}}-\bigg{(}1+\frac{2}{\hat{z}^{2}}\bigg{)}\frac{\operatorname{d\\!}{}\hat{y}}{\operatorname{d\\!}{}\hat{z}}+\frac{2}{\hat{z}^{3}}\hat{y}=0,$

for which a series solution as $\hat{z}\to\infty$ may be considered. The
coefficients of this series solution may be determined exactly, yielding

(72)
$\hat{y}(\hat{z})=\sum_{p=0}^{\infty}\frac{(-1)^{p}\Gamma(2p)}{2^{p-1}\Gamma(p)}\hat{z}^{-{2p}}.$

####  A.2.1 Inner-outer matching

Now that the outer limit of the inner solution near $z=1$ is known from (72)
we may match, in terms of the outer variable $z$, with the inner limit of the
outer solution. Writing (72) in terms of the outer variables $z$ and $y$ from
(70) gives the outer limit of the inner solution as

(73)
$y=\sum_{p=0}^{\infty}\frac{\Lambda_{0}(-1)^{p}\Gamma(2p)}{2^{p-1}\Gamma(p)}\frac{\epsilon^{p}}{(1-z)^{2p+1}}.$

We now match the $O(\epsilon^{p})$ component of (73) with that of the outer
solution (69),

(74)
$B_{p}\sim\tilde{\Lambda}(-1)^{p+\beta}2^{p+\beta}\frac{\Gamma(p+\beta)}{(1-z)^{2p+2\beta}},$

which yields the constants $\beta$ and $\tilde{\Lambda}$ as

(75) $\beta=1/2\qquad\text{and}\qquad\tilde{\Lambda}=\frac{\mathrm{i}}{2\pi}.$

##  Appendix B Details for the Pearcey equation

Equation (46) permits three WKBJ solutions of the form
$A_{m}(z,\epsilon)\exp{(-S_{m}(z)/\epsilon)}$, where $m=1$, $2$, and $3$.
Substitution of this ansatz into (46) yields the following equations for the
singulants, $S_{m}^{\prime}$, and amplitude function, $A_{m}$:

(76a)
$\displaystyle\big{[}S_{m}^{\prime}(z)\big{]}^{3}-S_{m}^{\prime}(z)+\mathrm{i}z=0,$
(76b) $\displaystyle\epsilon^{2}A_{m}^{\prime\prime\prime}-3\epsilon
S_{m}^{\prime}A_{m}^{\prime\prime}+\Big{(}3(S_{m}^{\prime})^{2}-3\epsilon
S_{m}^{\prime\prime}-1\Big{)}A_{m}^{\prime}+\Big{(}3S_{m}^{\prime}S_{m}^{\prime\prime}-\epsilon
S_{m}^{\prime\prime\prime}\Big{)}A_{m}=0.$

It will be seen that the asymptotic expansions for the amplitude functions,
$A_{m}$, diverge. The Stokes phenomenon that occurs as a consequence of this
divergence will switch on one of the other initial singulants. For instance,
across a $1>2$ Stokes line, we will have $A_{1}\exp{(-S_{1}/\epsilon)}\mapsto
A_{1}\exp{(-S_{1}/\epsilon)}+A_{2}\exp{(-S_{2}/\epsilon)}$. However, this
$1>2$ Stokes line will intersect a $2>3$ Stokes line at a Stokes crossing
point, for which a new Stokes line is required in order to avoid a
contradiction. We now demonstrate how this new Stokes line, and the HOSL that
generates it, may be derived through determination of the late-late
divergence.

In expanding $A_{m}(z,\epsilon)=A_{m}^{(0)}+\epsilon A_{m}^{(1)}+\cdots$ in
equation (76b), we find the following equations at each order of $\epsilon$:

(77a) $\displaystyle O(1):$
$\displaystyle\qquad\mathcal{D}\big{[}A_{m}^{(0)}\big{]}\equiv\Big{[}3(S_{m}^{\prime})^{2}-1\Big{]}A^{(0)\prime}_{m}+3S_{m}^{\prime}S_{m}^{\prime\prime}A^{(0)}_{m}=0,$
(77b) $\displaystyle O(\epsilon):$
$\displaystyle\qquad\mathcal{D}\big{[}A_{m}^{(1)}\big{]}=3S_{m}^{\prime}A_{m}^{(0)\prime\prime}+3S_{m}^{\prime\prime}A_{m}^{(0)\prime}+S_{m}^{\prime\prime\prime}A_{m}^{(0)},$
(77c) $\displaystyle O(\epsilon^{n}):$
$\displaystyle\qquad\mathcal{D}\big{[}A_{m}^{(n)}\big{]}=3S_{m}^{\prime}A_{m}^{(n-1)\prime\prime}+3S_{m}^{\prime\prime}A_{m}^{(n-1)\prime}+S_{m}^{\prime\prime\prime}A_{m}^{(n-1)}-A_{m}^{(n-2)\prime\prime\prime},$

where $\mathcal{D}$ is the linear operator defined in the leading-order
equation (77a). The first of these, equation (77a), yields the solution

(78) $A^{(0)}_{m}(z)=\frac{C^{(0)}_{m}}{(3[S_{m}^{\prime}(z)]^{2}-1)^{1/2}},$

where $C_{m}^{(0)}$ is a constant of integration. The singular behaviour of
(78) may be determined by solving the cubic equation (76a) for
$S_{m}^{\prime}$, which yields the three solutions

(79) $S^{\prime}_{1}(z)=\frac{12+p^{2}(z)}{6p(z)}\quad\text{and}\quad
S^{\prime}_{2,3}(z)=-\frac{(1+a\mathrm{i}\sqrt{3})p^{2}(z)+12(1-a\mathrm{i}\sqrt{3})}{12p(z)}.$

Here, $a=1$ for $S_{2}^{\prime}$ and $a=-1$ for $S_{3}^{\prime}$, and
$p(z)=\big{[}12\mathrm{i}\big{(}-9z+\sqrt{81z^{2}+12}\big{)}\big{]}^{1/3}$.

We note that $A_{m}^{(0)}$ is singular whenever $(S_{m}^{\prime})^{2}=1/3$,
for which substitution into (76a) yields two singular points, given by

(80) $z_{+}=\frac{2\sqrt{3}\mathrm{i}}{9}\qquad\text{and}\qquad
z_{-}=-\frac{2\sqrt{3}\mathrm{i}}{9}.$

The singularities in $A_{m}^{(0)}$ lie at $z_{+}$ for $m=2$ and $m=3$, and at
$z_{-}$ for $m=1$ and $m=2$. All of these singularities are of order $1/4$. As
an example, we have that $A_{1}^{(0)}(z)\sim
6C_{1}^{(0)}/(-36\sqrt{3}\mathrm{i}[z-z_{-}])^{1/4}$.

###  B.1 Integrating the singulant solutions

In order to determine the Stokes line structure in the next section, each of
(79) must be integrated to find $S_{m}$. We now determine these explicitly
through examination of the Hamiltonian system for the Borel operator of the
Pearcy equation, much like that performed by Honda et al. (2015).

We consider the Borel operator, $\mathscr{P}_{B}$, associated with the Pearcy
equation (46). The symbol, $\sigma$, of the operator $\mathscr{P}_{B}$ is a
complex Hamiltonian on the phase space
$T^{*}(\mathbb{C}_{w}\times\mathbb{C}_{z})$, where $\mathbb{C}_{z}$ is the
physical plane and $\mathbb{C}_{w}$ is the Borel plane. These are given by

(81)
$\mathscr{P}_{B}=\partial^{3}_{z}-\partial_{w}^{2}\partial_{z}-\mathrm{i}z\partial_{w}^{3}\qquad\text{and}\qquad\sigma=\xi^{3}-\eta^{2}\xi-\mathrm{i}z\eta^{3},$

where $\xi$ and $\eta$ are local coordinates on the cotangent directions. The
Hamiltonian equations arising from the natural Kähler structure on
$T^{*}(\mathbb{C}_{w}\times\mathbb{C}_{z})$,

(82)
$\dot{z}(\tau)=\partial_{\xi}\sigma,\quad\dot{w}(\tau)=\partial_{\eta}\sigma,\quad\dot{\xi}(\tau)=-\partial_{z}\sigma,\quad\dot{\eta}(\tau)=-\partial_{w}\sigma,\quad\sigma(\tau)=0,$

may be solved to find the solutions

(83a) $\displaystyle{\eta}(\tau)$ $\displaystyle=\eta_{0},$ (83b)
$\displaystyle{\xi}(\tau)$ $\displaystyle=\mathrm{i}\eta_{0}^{3}\tau+\xi_{0},$
(83c) $\displaystyle{z}(\tau)$
$\displaystyle=-\eta_{0}^{6}\tau^{3}+3\mathrm{i}\xi_{0}\eta_{0}^{3}\tau^{2}+(3\xi_{0}^{2}-\eta_{0}^{2})\tau+z_{0},$
(83d) $\displaystyle{w}(\tau)$
$\displaystyle=\frac{3\mathrm{i}\eta_{0}^{8}}{4}\tau^{4}+3\eta_{0}^{5}\xi_{0}\tau^{3}+\frac{\mathrm{i}\eta_{0}^{4}-9\mathrm{i}\eta_{0}^{2}\xi_{0}^{2}}{2}\tau^{2}-(2\eta_{0}\xi_{0}+3\mathrm{i}\eta_{0}^{2}z_{0})\tau+w_{0}.$

Here, $\eta_{0}$, $\xi_{0}$, $z_{0}$, and $w_{0}$ are constants of
integration. The desired solution, $S_{m}$, equals $w(\tau)$ from (83d), which
may be expressed as a function of $z$ through knowledge of $z(\tau)$ from
(83c). The four constants of integration in (83) are determined through
imposition of the conditions

(84) $\eta(0)=1,\qquad w(0)=0,\qquad
z(0)=\pm\frac{2\sqrt{3}\mathrm{i}}{9},\qquad\sigma(0)=0.$

The first three of these yield $\eta_{0}=1$, $w_{0}=0$, and $z_{0}=\pm
2\sqrt{3}\mathrm{i}/9$. The last of conditions (84) gives the relationship
$\xi_{0}^{3}-\xi_{0}-\mathrm{i}z_{0}=0$. Since $z_{0}$ is known, this equation
has the four solutions of $\xi_{0}=\pm\sqrt{3}/3$ and $\xi_{0}=\pm
2\sqrt{3}/3$. The last two of these solutions give $S_{m}\neq 0$ at the
singular point of the leading-order outer solution $A_{m}^{(0)}$ from (78).
Thus, we consider only $\xi_{0}=\pm\sqrt{3}/3$ for which the plus sign may be
taken without loss of generality as the resultant expression for $w$ involves
$\xi_{0}^{2}$ only. In defining
$\hat{\tau}=-\mathrm{i}(\tau-\mathrm{i}\xi_{0})$, we substitute $\eta_{0}=1$,
$w_{0}=0$, $z_{0}=\pm 2\sqrt{3}\mathrm{i}/9$, and $\xi_{0}=\sqrt{3}/3$ into
equations (83c) and (83d) to find

(85)
${z}(\hat{\tau})=\mathrm{i}\hat{\tau}(\hat{\tau}^{2}-1)\qquad\text{and}\qquad{w}(\hat{\tau})=\frac{3\mathrm{i}}{4}\left(\hat{\tau}^{2}-1/3\right)^{2}.$

We recognise the first of equations (85) as the same cubic equation that
governs $S_{m}^{\prime}$ in (76a), and hence $S^{\prime}_{m}=\hat{\tau}$.
Furthermore, we may calculate
$\mathrm{d}w/\mathrm{d}z=(\mathrm{d}w/\mathrm{d}\hat{\tau})(\mathrm{d}\hat{\tau}/\mathrm{d}z)=\hat{\tau}=S_{m}^{\prime}(z)$,
which verifies that $S_{m}^{\prime}=w^{\prime}$. Since we have
$S^{\prime}_{m}=\hat{\tau}$, it is not necessary to invert $z(\hat{\tau})$ to
find $\hat{\tau}(z)$. Instead, we simply write our expression for $S_{m}$ in
terms of the previously determined function $S_{m}^{\prime}$, which yields

(86)
$S_{m}(z)=\frac{3\mathrm{i}}{4}\left(\big{[}S_{m}^{\prime}(z)\big{]}^{2}-\frac{1}{3}\right)^{2}.$

Equation (86) is the main result of this section: an explicit solution for the
singulant, $S_{m}$, specified in terms of the function $S^{\prime}_{m}$
previously determined in equation (79).

###  B.2 Late-term divergence

We consider a late-term ansatz of the form

(87)
$A^{(n)}_{m}(z)\sim\sum_{p=0}^{\infty}B^{(p)}_{m}(z)\frac{\Gamma(n+\alpha-p)}{[{\chi_{m}}(z)]^{n+\alpha-p}},$

where $\alpha$ is a constant, $\chi_{m}$ is the singulant, and $B_{m}^{(p)}$
are amplitude functions. The HOSP may be detected through analysis of the
late-late term divergence of $B_{m}^{(p)}$ as $p\to\infty$. Substitution of
(87) into the $O(\epsilon^{n})$ equation (77c) yields at each order of
$\Gamma(n+\alpha-p)/\chi_{m}^{n+\alpha-p}$ the equations

(88a)
$\displaystyle{\chi}_{m}^{\prime}\Big{(}({\chi}_{m}^{\prime})^{2}+3S_{m}^{\prime}{\chi}_{m}^{\prime}+3(S_{m}^{\prime})^{2}-1\Big{)}=0,$
(88b)
$\displaystyle\mathcal{L}\Big{[}B^{(0)}_{m}\Big{]}\equiv\left(3[S_{m}^{\prime}+\chi_{m}^{\prime}]^{2}-1\right)B^{(0)\prime}_{m}+3(S_{m}^{\prime}+\chi_{m}^{\prime})(S_{m}^{\prime\prime}+\chi_{m}^{\prime\prime})B^{(0)}_{m}=0,$
(88c)
$\displaystyle\mathcal{L}\Big{[}B^{(1)}_{m}\Big{]}=3\left(S_{m}^{\prime}+\chi_{m}^{\prime}\right)B^{(0)\prime\prime}_{m}+3\left(S_{m}^{\prime\prime}+\chi_{m}^{\prime\prime}\right)B^{(0)\prime}_{m}+\left(S_{m}^{\prime\prime\prime}+\chi_{m}^{\prime\prime\prime}\right)B^{(0)}_{m},$
(88d) $\displaystyle\begin{aligned}
\mathcal{L}\Big{[}B^{(p)}_{m}\Big{]}&=3\left(S_{m}^{\prime}+\chi_{m}^{\prime}\right)B^{(p-1)\prime\prime}_{m}+3\left(S_{m}^{\prime\prime}+\chi_{m}^{\prime\prime}\right)B^{(p-1)\prime}_{m}\\\
&\quad+\left(S_{m}^{\prime\prime\prime}+\chi_{m}^{\prime\prime\prime}\right)B^{(p-1)}_{m}-B^{(p-2)\prime\prime\prime}_{m},\end{aligned}$

where (88d) holds for $p\geq 2$, and $\mathcal{L}$ is the differential
operator defined in (88b).

We begin by solving equation (88a) to find $\chi_{m}$. Six solutions,
${\chi}_{m}^{\prime}=-3S_{m}^{\prime}/2\pm\sqrt{4-3(S_{m}^{\prime})^{2}}/2$,
are found in total, but two of these are unable to match with inner solutions
near the associated singular points. It may be verified that these solutions
for $\chi_{m}^{\prime}$ yield the difference between another $S_{j}^{\prime}$
and $S_{m}^{\prime}$ (for instance
$\chi_{1}^{\prime}=S_{2}^{\prime}-S_{1}^{\prime}$). Integration of the
expressions for $\chi_{m}^{\prime}$ then yields the four initial singulants

(89) $\left\\{\quad\begin{aligned}
\chi_{1}(z)&=S_{2}(z)-S_{1}(z),\qquad\chi_{3}(z)=S_{2}(z)-S_{3}(z),\\\
\chi_{2a}(z)&=S_{1}(z)-S_{2}(z),\qquad\chi_{2b}(z)=S_{3}(z)-S_{2}(z),\end{aligned}\right.$

which satisfy $\chi_{1}(z_{-})=0$, $\chi_{3}(z_{+})=0$, $\chi_{2a}(z_{-})=0$,
and $\chi_{2b}(z_{+})=0$.

###  B.3 Late-late-term divergence

We now consider equations (88) for the late-term amplitude functions,
$B_{m}^{(p)}$. The first of these may be solved to find the solution

(90)
$B_{m}^{(0)}(z)=\frac{\Lambda_{m}^{(0)}}{\Big{(}3(S_{m}^{\prime}+\chi_{m}^{\prime})^{2}-1\Big{)}^{1/2}},$

where $\Lambda_{m}^{(0)}$ is a constant of integration. For $m=2$,
$B_{2}^{(0)}$ in (90) is singular at the same locations as that predicted by
the factorial-over-power form with $\chi_{2}$. However, for $m=1$ and $m=3$,
(90) is singular at additional locations beyond that where the factorial-over-
power form has $\chi_{m}(z)=0$. For example, when $m=1$, while $\chi_{1}$
predicts a singularity at $z_{-}$ only, the amplitude function (90) is also
singular at $z_{+}$. This will lead to a new divergent series as $p\to\infty$,
which we now study for $m=1$ and $m=3$.

In substituting a late-late-term factorial-over-power ansatz, of the form

(91) $B_{m}^{(p)}(z)\sim
C_{m}(z)\frac{\Gamma(p+\beta)}{[\tilde{\chi}_{m}(z)]^{p+\beta}},$

into equation (88d), we find the equation
$(\tilde{\chi}_{m}^{\prime})^{2}+3(S_{m}^{\prime}+\chi_{m}^{\prime})\tilde{\chi}^{\prime}_{m}+3(S_{m}^{\prime}+\chi_{m}^{\prime})^{2}-1=0$.
For each value of $m$, only one solution of this quadratic equation has the
correct singular behaviour after the boundary condition
$\tilde{\chi}_{1}(z_{+})=0$ or $\tilde{\chi}_{3}(z_{-})=0$ is applied. This
yields our late-late-term singulant as

(92) $\tilde{\chi}_{m}(z)=\left\\{\begin{aligned}
S_{3}(z)-S_{2}(z)\qquad\text{if $m=1$},\\\ S_{1}(z)-S_{2}(z)\qquad\text{if
$m=3$}.\end{aligned}\right.$

From the method developed in section 3.4, these will induce the HOSP, in which
a new factorial-over-power component of $A_{m}^{(n)}$ with a singulant of
$\chi_{m}+\tilde{\chi}_{m}$ will switch on across a HOSL given by
$\text{Im}[\tilde{\chi}_{m}/\chi_{m}]=0$ and
$\text{Re}[\tilde{\chi}_{m}/\chi_{m}]\geq 0$. The dominant components of the
late-terms of the initial amplitude functions, $A_{m}^{(n)}$, are therefore
given by

(93a) $\displaystyle A_{1}^{(n)}(z)$ $\displaystyle\sim
B_{1}^{(0)}(z)\frac{\Gamma(n+\alpha)}{(S_{2}-S_{1})^{n+\alpha}}+\tilde{\sigma}_{1}(z,n)C_{1}(z)\frac{\Gamma(n+\alpha)}{(S_{3}-S_{1})^{n+\alpha}},$
(93b) $\displaystyle A_{2}^{(n)}(z)$ $\displaystyle\sim
B_{2a}^{(0)}(z)\frac{\Gamma(n+\alpha)}{(S_{3}-S_{2})^{n+\alpha}}+B_{2b}^{(0)}(z)\frac{\Gamma(n+\alpha)}{(S_{1}-S_{2})^{n+\alpha}},$
(93c) $\displaystyle A_{3}^{(n)}(z)$ $\displaystyle\sim
B_{3}^{(0)}(z)\frac{\Gamma(n+\alpha)}{(S_{2}-S_{3})^{n+\alpha}}+\tilde{\sigma}_{3}(z,n)C_{3}(z)\frac{\Gamma(n+\alpha)}{(S_{1}-S_{3})^{n+\alpha}},$

where $\tilde{\sigma}_{1}$ and $\tilde{\sigma}_{3}$ are the higher-order
Stokes multipliers. The resultant Stokes lines structure was shown in figure
4.