# Travelling Waves and Exponential Nonlinearities in the Zeldovich-Frank-
Kamenetskii Equation

S. Jelbart, K. U. Kristiansen & P. Szmolyan

###### Abstract

We prove the existence of a family of travelling wave solutions in a variant
of the Zeldovich-Frank-Kamenetskii (ZFK) equation, a reaction-diffusion
equation which models the propagation of planar laminar premixed flames in
combustion theory. Our results are valid in an asymptotic regime which
corresponds to a reaction with high activation energy, and provide a rigorous
and geometrically informative counterpart to formal asymptotic results that
have been obtained for similar problems using high activation energy
asymptotics. We also go beyond the existing results by (i) proving smoothness
of the minimum wave speed function $\overline{c}(\epsilon)$, where
$0<\epsilon\ll 1$ is the small parameter, and (ii) providing an asymptotic
series for a flat slow manifold which plays a role in the construction of
travelling wave solutions for non-minimal wave speeds
$c>\overline{c}(\epsilon)$. The analysis is complicated by the presence of an
exponential nonlinearity which leads to two different scaling regimes as
$\epsilon\to 0$, which we refer to herein as the convective-diffusive and
diffusive-reactive zones. The main idea of the proof is to use the geometric
blow-up method to identify and characterise a $(c,\epsilon)$-family of
heteroclinic orbits which traverse both of these regimes, and correspond to
travelling waves in the original ZFK equation. More generally, our analysis
contributes to a growing number of studies which demonstrate the utility of
geometric blow-up approaches to the study dynamical systems with singular
exponential nonlinearities.

Keywords: geometric singular perturbation theory, geometric blow-up,
travelling waves, combustion theory, exponential nonlinearity

MSC2020: 34E13, 34E15, 34E10, 35K57, 80A25

##  1 Introduction

Exponential nonlinearities appear differential equation models across the
natural sciences and beyond, ranging from models of earthquake faulting [3,
12, 23] to electrical oscillators [16, 18, 30] and plastic deformation in
metals [4, 13]. It is also common for modelers to introduce exponential
nonlinearities for mathematical reasons, for example in the regularisation of
non-smooth systems using a $\tanh$ function [22], the study of gene regulatory
networks which converge to a non-smooth approximating system at an exponential
rate [15, 17, 32, 34], or more recently, via the introduction of singular
exponential coordinates in order to study polynomial systems using techniques
based on tropical geometry [28, 39].

All of the problems cited above share a common mathematical feature, namely,
that the exponential nonlinearities lead to singular perturbation problems
which have distinct dynamics in different regions of phase space. The basic
reason for this can be seen by considering the nonlinearity
${\mathrm{e}}^{x/\epsilon}$, which is exponentially small, $\mathcal{O}(1)$,
or exponentially large as $\epsilon\to 0$ when $x<0$,
$x=\mathcal{O}(\epsilon)$, or $x>0$ respectively. Thus, if $x$ is a state
variable in a system of differential equations, then one expects to obtain
distinct dynamics for $x<0$, $x=0$ and $x>0$ as $\epsilon\to 0$. Regions of
phase space where the limit is not defined due to unbounded growth, in this
case when $x>0$, can often be controlled after a suitable ‘normalisation’
(positive transformation of time). For example, the authors in [18] applied a
normalisation which amounted to a division of the governing equations by
$1+{\mathrm{e}}^{x/\epsilon}$ in order to control the (otherwise unbounded)
growth of the nonlinearity ${\mathrm{e}}^{x/\epsilon}$. Since

$\frac{{\mathrm{e}}^{x/\epsilon}}{1+{\mathrm{e}}^{x/\epsilon}}\rightarrow\begin{cases}0&x<0\\\
1&x>1\end{cases}\quad\mbox{for}\quad\epsilon\to 0,$

the normalised system is piecewise-smooth in the limit $\epsilon\to 0$. From
here, the loss of smoothness can be ‘resolved’ via a geometric blow-up
procedure which involves the insertion of an inner/switching layer which
describes the dynamics for $x=\mathcal{O}(\epsilon)$, i.e. ‘in between’ the
distinct regimes $x<0$ and $x>0$, using techniques that have in recent years
been developed and applied in a large number of works, e.g. in [2, 8, 19, 20,
24, 25, 26, 27, 29, 31].

In this work we show that similar techniques can also be used to describe
travelling fronts in a reaction-diffusion equation known as the Zeldovich-
Frank-Kamenetskii (ZFK) equation, which arises in combustion theory as a model
for the propagation of planar laminar premixed flames [5, 6, 7, 42, 43]. We
focus on one of the simplest variants of the ZFK equation, taken from [7] (see
Section 2 for details), but we also expect our methods to apply to more
complicated models for planar flame propagation away from so-called
flammability limits which – like the simple ZFK equation considered herein –
can be described by a (system of) reaction-diffusion equation(s) for which the
reaction rate is determined by an Arrhenius law with high activation energy
(this is similar to the electrical oscillator models from [16, 30] considered
in [18]). The small parameter is $\epsilon:=\beta^{-1}$, where $\beta$ is the
so-called Zeldovich number, a dimensionless quantity related to the activation
energy of the Arrhenius reaction, and $\beta\to\infty$ ($\epsilon\to 0$) is
the high activation energy limit. This leads to an exponential nonlinearity
which, despite the relative abundance formal asymptotic studies based on high
activation energy asymptotics (HAEA) (here we refer again to [5, 6, 7, 42,
43]), has often been viewed as an obstacle to rigorous geometric analysis
using modern mathematical approaches like geometric singular perturbation
theory (GSPT) [14, 21, 41]. We note that the problem considered herein, which
is representative of problems with exponential nonlinearities quite generally,
is smooth for each $\epsilon>0$, but non-smooth in the limit $\epsilon\to 0$.
This distinguishes it from the problems considered in e.g. [11, 33, 35, 36,
37, 38], where travelling waves that are identified using GSPT and geometric
blow-up feature a non-smooth cutoff which renders the equations non-smooth for
$\epsilon>0$.

The travelling wave problem can be formulated as a boundary value problem in a
planar system of singularly perturbed ODEs. There are two important scaling
regimes: a convective-diffusive zone where the reaction term is exponentially
small, and a reactive-diffusive zone where the convective term can be
neglected [7, Ch. 2.1.2]. We use the geometric blow-up method to glue these
two zones together, and subsequently identify a heteroclinic orbit which
corresponds to a travelling wave in the original ZFK equation. Specifically,
we show that for each $\epsilon>0$ sufficiently small there exists a
$(c,\epsilon)$-family of travelling waves, where $c>0$ is the wave speed, and
we prove the existence of a minimal wave speed $c=\overline{c}(\epsilon)$.
These results can be seen as a rigorous and geometrically informative
counterpart to existing results that have been derived for similar problems
using formal methods based on HAEA. In addition to providing a rigorous
justification of existing results, we (i) prove the smoothness of the minimum
wave speed function $\overline{c}(\epsilon)$, and provide the asymptotics of
$\overline{c}(\epsilon)$ up to and including $\mathcal{O}(\epsilon)$. This
calculation shows that $\overline{c}(\epsilon)>1$ for all $0<\epsilon\ll 1$.
Moreover, we (ii) identify and provide an asymptotic expansion for a slow
manifold which characterises part of the wave profile for non-minimal wave
speeds $c>\overline{c}(\epsilon)$ in the convective-diffusive zone. We use a
normal form from [9] in order to derive the smoothness results in (i), and we
use an approach inspired by [22] to derive the asymptotics in (ii). The
approach used in (ii) is of independent interest, insofar as we expect that it
can be applied in order to determine the asymptotics of flat slow manifolds
quite generally.

The manuscript is structured as follows: In Section 2 we introduce the ZFK
equation and derive the singularly perturbed system of ODEs that we use in
order to study the travelling wave problem. Section 3 is devoted to the
separate geometric analyses in the convective-diffusive and reactive-diffusive
zones. We state the main result in Section 4, which is proven using geometric
blow-up and normal forms in Section 5. We conclude in Section 6 with a summary
and outlook.

##  2 The ZFK equation

Our starting point is the scalar ZFK equation

$\theta_{t}=\theta_{xx}+\omega(\theta,\beta),$ (1)

where $t\geq 0$, $x\in\mathbb{R}$ and $\theta\in[0,1]$ is the so-called
reduced temperature, a non-dimensional variable which measures the ratio of
burnt to unburnt gas in a simple reaction involving two species (corresponding
to burnt and unburnt premix). In particular, $\theta=0$ corresponds to
completely unburnt premix/fuel, and $\theta=1$ corresponds to a completely
burnt state. The function $\omega(\theta,\beta)$ models the reaction term,
which depends upon the Zeldovich number $\beta>0$. Equation (1) can be
obtained via a reduction from a two component system of reaction-diffusion
equations which models the propagation of planar flames if the diffusivity of
the two reactants are equal, i.e. when the so-called Lewis number $L\equiv 1$.
We refer to [5, 6, 7, 42, 43] and the many references therein for details on
the physical interpretation and derivation of the model. We consider the case
for which the reaction term is given by

$\omega(\theta,\beta)=\lambda\theta(1-\theta){\mathrm{e}}^{-\beta(1-\theta)}.$
(2)

where $\lambda>0$ is a parameter. The asymptotic regime $0<\beta\ll 1$ is
known as the KPP regime, since equation (1) reduces to the KPP equation when
$\beta\to 0$. We shall be interested in the opposite regime $\beta\gg 1$,
known as the ZFK regime, and the corresponding high activation energy limit
$\beta\to\infty$. The parameter $\lambda>0$ can be scaled out using the
parabolic rescaling $(x,t)\mapsto(x/\sqrt{\lambda},t/\lambda)$, however, it is
customary to fix a particular choice for $\lambda$ which ensures that the
minimum wave speed derived in later sections is of order unity. We follow this
convention here, and set

$\lambda=\frac{\beta^{2}}{2}:\quad\omega(\theta,\beta)=\frac{\beta^{2}}{2}\theta(1-\theta){\mathrm{e}}^{-\beta(1-\theta)},$

without loss of generality. This particular choice of $\omega$ and $\lambda$
was used in [7, Ch. 8.3] in order to study the so-called KPP-ZFK transition
from [40, Part. F]. We shall write

$\epsilon:=\frac{1}{\beta}:\quad\omega(\theta,\epsilon^{-1})=\frac{\epsilon^{-2}}{2}\theta(1-\theta){\mathrm{e}}^{-\epsilon^{-1}(1-\theta)}$

and consider the asymptotic behaviour as $\epsilon\to 0$ (which corresponds to
$\beta\to\infty$).

######  Remark 2.1.

We restrict attention to the physically meaningful domain $\theta\in[0,1]$,
but do note that $\omega$ is unbounded for $\theta>1$ for $\epsilon\to 0$. The
graph of $\theta\mapsto\omega(\theta,\epsilon^{-1})$, $\theta\in[0,1]$, for
five different values of $\epsilon$ is shown in Fig. 1. Notice that
$\omega(1,\epsilon^{-1})=0$ for all $\epsilon>0$ and that $\omega$ is
exponentially small with respect to $\epsilon\to 0$ for fixed
$\theta\in[0,1)$. Nevertheless, simple calculations shhow that $\omega$ is
unbounded on $\theta\in[0,1]$ due a local maximum
$\theta=\theta_{c}(\epsilon)=1+\mathcal{O}(\epsilon)$ with
$\omega(\theta_{c}(\epsilon),\epsilon)=\mathcal{O}(\epsilon^{-1})$. As we
shall see, this lack of uniformity with respect to $\epsilon\to 0$ leads to
distinct limiting problems in two distinct regimes.

Figure 1: Graph of $\theta\mapsto\omega(\theta,\epsilon^{-1})$ for five
different values of $\epsilon$, equidistributed between $\epsilon=0.03$ (blue)
and $\epsilon=0.003$ (green).

######  Remark 2.2.

In combustion theory, one typically encounters the form

$\omega(\theta,\beta)=\lambda(1-\theta){\mathrm{e}}^{-\beta(1-\theta)}$ (3)

instead of (2), see e.g. [7, Ch. 2.1 & Ch. 8.2] and [5, 6, 42]. This form can
be derived from the chemistry, however, the travelling wave problem is ill-
posed because $\omega(0,\beta)=\lambda{\mathrm{e}}^{-\beta}>0$ for all
$\beta>0$, which means that the boundary condition at $\theta=0$ (see equation
(5) below) cannot be satisfied. This is known as the cold-boundary difficulty
[7, 42]. In order to circumvent it, one can either (i) adjust the reaction
term so that $\omega(0,\beta)=0$, as in (2), or (ii) introduce a cutoff by
imposing that $\omega(\theta,\beta)\equiv 0$ for all $\theta<\theta_{c}$,
where $\theta_{c}\in(0,1)$ is close to zero. The former approach is typically
preferred in analytical studies such as ours, while the latter tends to be
preferable for applications and simulations. Analytically, we expect that the
ZFK equation (1) with (3) and a cutoff can be treated using a combination of
our approach with the geometric methods developed for problems with cutoffs in
[11, 33, 35, 36, 37, 38].

###  2.1 Travelling waves

We study the existence of travelling wave solutions for $0<\epsilon\ll 1$
within the physically meaningful regime $\theta\in[0,1].$ Rewriting equation
(1) in the travelling wave frame $z=x+ct$, where $c\in\mathbb{R}$ is the
(presently unknown) wave speed, we obtain the following second order ODE for
travelling waves $\theta(x,t)=\widetilde{\theta}(z)$:

$c\frac{{\mathrm{d}}\widetilde{\theta}}{{\mathrm{d}}z}=\frac{{\mathrm{d}}^{2}\widetilde{\theta}}{{\mathrm{d}}z^{2}}+\frac{1}{2}\epsilon^{-2}\widetilde{\theta}(1-\widetilde{\theta}){\mathrm{e}}^{-\epsilon^{-1}(1-\widetilde{\theta})}.$
(4)

We drop the tildes henceforth. As in [7, Ch. 2.1.2 & 8.2], the boundary
conditions are

$\lim_{z\to-\infty}\theta(z)=0,\qquad\lim_{z\to\infty}\theta(z)=1.$ (5)

We now set $\eta:={\mathrm{d}}\theta/{\mathrm{d}}z$ and rewrite equation (4)
as the first order system

$\begin{split}\dot{\theta}&=\eta,\\\
\dot{\eta}&=c\eta-\frac{1}{2}\epsilon^{-2}\theta(1-\theta){\mathrm{e}}^{-\epsilon^{-1}(1-\theta)},\end{split}$
(6)

where the overdot denotes differentiation with respect to $z$. There are two
equilibria:

$p_{-}:(0,0),\qquad p_{+}:(1,0),$ (7)

which exist for all $c\in\mathbb{R}$ and $\epsilon>0$. Direct calculations
show that $p_{-}$ is a proper unstable node for all
$c>2\epsilon^{-2}{\mathrm{e}}^{-\epsilon^{-1}}$, whereas $p_{+}$ is a saddle.
The strong and weak eigenvalues at the node $p_{-}$, which we denote here by
$\lambda_{s}$ and $\lambda_{w}$ respectively, satisfy

$\displaystyle\lambda_{s}:=\frac{c}{2}\left(1+\sqrt{1-\frac{2}{c\epsilon^{2}{\mathrm{e}}^{\epsilon^{-1}}}}\right)>\lambda_{w}:=\frac{c}{2}\left(1-\sqrt{1-\frac{2}{c\epsilon^{2}{\mathrm{e}}^{\epsilon^{-1}}}}\right),$

and the corresponding strong and weak eigenvectors
$(1,v_{s})^{\textnormal{T}}$ and $(1,v_{w})^{\textnormal{T}}$ converge to
$(1,c)^{\textnormal{T}}$ and $(0,1)^{\textnormal{T}}$ respectively as
$\epsilon\to 0$. On the other hand, $p_{+}$ lies right at the interface

$\Sigma=\\{(1,\eta):\eta\in\mathbb{R}\\}$ (8)

between $\theta<1$, where the reaction term is exponentially small, and
$\theta>1$, where the right-hand side of (6) is undefined for $\epsilon\to 0$.

Our aim is to identify and characterise (generally $(c,\epsilon)$-dependent)
heteroclinic connections between $p_{-}$ and $p_{+}$ as $\epsilon\to 0$, which
correspond to travelling wave solutions in the ZFK equation (1). Note that
while heteroclinic orbits that connect to $p_{-}$ tangentially to the weak
eigendirection are robust, heteroclinic orbits which connects to $p_{-}$
tangentially to the strong eigendirection are not. This follows from the
uniqueness of strong stable manifold, which will be clarified in Lemma 3.1
below. In the following, we also refer to the heteroclinic connections
themselves as ‘weak’ or ‘strong’ in order to distinguish these two cases.

######  Remark 2.3.

System (6) has a symmetry

$(\eta,z,c)\leftrightarrow(-\eta,-z,-c).$ (9)

In the following we restrict attention to positive wave speeds $c>0$, keeping
in mind that corresponding statements can be derived for negative wave speeds
using (9).

##  3 Geometric singular perturbation analysis

We turn now to the geometry and dynamics in the singular limit $\epsilon\to
0$. As identified above, there are two important regimes to consider:

  * •

A convective-diffusive zone defined by $1-\theta=\mathcal{O}(1)$, where the
reaction term in (6) is small;

  * •

A reactive-diffusive zone defined by $1-\theta=\mathcal{O}(\epsilon)$, where
the advection term $c\eta$ is small relative to the reaction and diffusion
terms due to the spatial scale.

We consider the geometry and dynamics in each zone separately.

###  3.1 Equations in the convective-diffusive zone

Taking $\epsilon\to 0$ in system (6) with $\theta\in[0,1)$ leads to the
limiting problem in the convective-diffusive zone:

$\begin{split}\dot{\theta}&=\eta,\\\ \dot{\eta}&=c\eta,\end{split}$ (10)

which does not depend on $\omega$. System (10) has a $1$-dimensional critical
manifold

$S:=\left\\{(\theta,0):\theta\in[0,1)\right\\},$

and the Jacobian along $S$ has a single non-trivial eigenvalue $\lambda=c>0$.
Thus, $S$ is normally hyperbolic and repelling. Note that $p_{-}\in S$, but
$p_{+}\notin S$ because $\theta<1$ in the convective-diffusive zone. The fast
fibers are contained within lines of the form

$\eta(\theta)=c\theta+\eta(0),$

see Fig. 2, and the strong stable manifold of $p_{-}$ is described by the
following.

Figure 2: Limiting geometry and dynamics of system (6) as $\epsilon\to 0$ in
the convective-diffusive zone. We show the critical manifold $S$ (in red), and
the set $\Sigma$ (in orange). $S$ is normally hyperbolic and repelling for
$\theta\in[0,1)$. The fast fibers (in green) are contained in lines of slope
$c$, with double arrows to indicate the ‘fast flow’ (with respect to
$z=x+ct$). The ‘slow flow’ on $S$ is increasing between the equilibria $p_{-}$
and $p_{+}$, and indicated by single arrows. The colouring and labelling
conventions adopted here are applied throughout.

######  Lemma 3.1.

Fix any $\kappa\in(0,1)$ and $c>0$. Then there exists an
$\epsilon_{0}=\epsilon_{0}(\kappa,c)>0$ such that for all
$\epsilon\in[0,\epsilon_{0})$ the strong unstable manifold of $p_{-}$ is
contained in a $C^{\infty}$-smooth graph defined by

$\displaystyle\eta=c\theta+R(\theta,\epsilon),$ (11)

over $\theta\in[0,1-\kappa]$, where $R$ is infinitely flat with respect to
$\epsilon\rightarrow 0^{+}$, i.e.

$\displaystyle\frac{\partial^{n}}{\partial\epsilon^{n}}R(\theta,0)=0\quad\forall\,n\in\mathbb{N},$

such that $R(\theta,\epsilon)=\mathcal{O}(\epsilon^{N})$ for all
$N\in\mathbb{N}$.

###### Proof.

This follows directly from the Stable Manifold Theorem when
$\epsilon\in(0,\epsilon_{0})$, i.e. when $p_{-}$ is a proper node [10, Thm.
2.8]. However, it remains to prove the result for $\epsilon=0$. Notice that
Fenichel theory yields a strong unstable manifold given by a $C^{N}$-smooth
graph of the form $\eta=c\theta+\mathcal{O}(\epsilon^{N})$ for any
$N\in\mathbb{N}$, however, it does not yield the $C^{\infty}$ property. In
order to prove the $C^{\infty}$ property we apply a polar blow-up of
${p_{-}}:(\theta,\eta)=(0,0)$. This turns the strong unstable manifold into a
saddle and the result then follows from [10, Thm. 2.7]. In the present case,
we perform the directional blowup:

$\displaystyle(r_{1},\theta_{1})\mapsto\begin{cases}\theta&=r_{1}\theta_{1},\\\
\eta&=r_{1}.\end{cases}$

Plugging this into (6) gives

$\displaystyle\dot{r}_{1}$
$\displaystyle=r_{1}\left(c-\frac{1}{2}\epsilon^{-2}\theta_{1}(1-r_{1}\theta_{1}){\mathrm{e}}^{-\epsilon^{-1}(1-r_{1}\theta_{1})}\right),$
(12) $\displaystyle\dot{\theta}_{1}$
$\displaystyle=1-\theta_{1}\left(c-\frac{1}{2}\epsilon^{-2}\theta_{1}(1-r_{1}\theta_{1}){\mathrm{e}}^{-\epsilon^{-1}(1-r_{1}\theta_{1})}\right).$

Here $(r_{1},\theta_{1})=(0,c^{-1})$ is a hyperbolic saddle for $\epsilon=0$
with the stable manifold contained within the invariant set $r_{1}=0$. In
fact, the system is $C^{\infty}$ smooth in a neighborhood of
$(r_{1},\theta_{1})=(0,c^{-1})$, also with respect to $\epsilon\geq 0$.
Consequently, there is a hyperbolic saddle for any $0\leq\epsilon\ll 1$ with
an unstable manifold that is $C^{\infty}$, also with respect to $c>0$ and
$\epsilon\geq 0$. This invariant manifold corresponds to the strong unstable
manifold upon blowing down. The expansion (11) is a consequence of (12) being
flat with respect to $\epsilon\rightarrow 0$. ∎

In order to justify the orientation of the ‘reduced flow’ indicated by single
arrows on $S$ in Fig. 2, note that Fenichel theory implies the existence of a
locally invariant slow manifold of the form

$S_{\epsilon}=\left\\{(\theta,h(\theta,\epsilon)):\theta\in[0,1-\kappa]\right\\},$
(13)

for an arbitrarily small but fixed constant $\kappa\in(0,1)$. The function $h$
is $C^{N}$-smooth for $\epsilon\in[0,\epsilon_{0}(N))$. Determining the slow
flow on $S_{\epsilon}$ is non-trivial, however, because the perturbation in
system (6) is flat as $\epsilon\to 0$.

######  Proposition 3.2.

The slow manifold $S_{\epsilon}$ has the following asymptotic expansion

$\begin{split}\eta=h(\theta,\epsilon)&=\sum_{k=1}^{\infty}\frac{1}{2^{k}}F_{k}(\theta,\epsilon)\epsilon^{-3k+1}{\mathrm{e}}^{-k\epsilon^{-1}(1-\theta)}\\\
&=\frac{1}{2c}\theta(1-\theta)\epsilon^{-2}{\mathrm{e}}^{-\epsilon^{-1}(1-\theta)}+\mathcal{O}(\epsilon^{-5}{\mathrm{e}}^{-2\epsilon^{-1}(1-\theta)}),\end{split}$
(14)

with $F_{1}(\theta,\epsilon)=\frac{1}{c}\theta(1-\theta)$ and each $F_{k}$,
$k\geq 2$, smooth satisfying the recursion

$\displaystyle
cF_{k}=\sum_{j=1}^{k-1}\left(\frac{\partial}{\partial\theta}F_{j}(\theta,\epsilon)\epsilon+jF_{j}(\theta,\epsilon)\right)F_{k-j}(\theta,\epsilon),\qquad
k\geq 2.$ (15)

###### Proof.

As already noted, the slow manifold is flat with respect to
$\epsilon\rightarrow 0$ and therefore we cannot determine the flow on the slow
manifold in the usual way. Inspired by [22], we proceed instead by looking for
a way to rewrite (6) as an extended system in which the flat term
$\frac{1}{2}\epsilon^{-2}\theta(1-\theta){\mathrm{e}}^{-\epsilon^{-1}(1-\theta)}$
(or some scalings thereof) defines a new variable. In the present case, it
turns out to be useful to introduce $\tilde{\eta}$ and $\zeta$ as follows:

$\displaystyle\tilde{\eta}=\epsilon^{-2}\eta,\quad\zeta=\frac{1}{2}\epsilon^{-4}{\mathrm{e}}^{-\epsilon^{-1}(1-\theta)}.$
(16)

In this way, we can eliminate the flat terms by working in the extended
$(\theta,\tilde{\eta},\zeta)-$space where the equations are given by

$\displaystyle\dot{\theta}$ $\displaystyle=\epsilon^{2}\tilde{\eta},$ (17)
$\displaystyle\dot{\tilde{\eta}}$
$\displaystyle=c\tilde{\eta}-\theta(1-\theta)\zeta,$
$\displaystyle\dot{\zeta}$ $\displaystyle=\epsilon\tilde{\eta}\zeta.$

Notice that the set

$\displaystyle\mathcal{Q}:=\left\\{(\theta,\tilde{\eta},\zeta)\,:\,\zeta=\frac{1}{2}\epsilon^{-4}{\mathrm{e}}^{-\epsilon^{-1}(1-\theta)}\right\\},$
(18)

defines an invariant set for all $\epsilon\geq 0$, upon which (17) reduces to
(6). System (17) has a critical manifold defined by
$\tilde{\eta}=\frac{1}{c}\theta(1-\theta)\zeta$, which is normally hyperbolic
and repelling. Fenichel theory implies the existence of a slow manifold for
system (17) of the form
$\tilde{\eta}=\frac{1}{c}\theta(1-\theta)\zeta+\mathcal{O}(\epsilon)$ on
compact subsets ($\theta\in[0,1-\kappa]$ as in (13)), for all $0<\epsilon\ll
1$.

At the same time, $\tilde{\eta}=\zeta=0$ is also invariant for (17), and the
linearization about any point in this set gives a single nonzero eigenvalue
$c$ for all $\epsilon\geq 0$. The center space at a point $(\theta_{0},0,{0})$
is spanned by the vectors $(1,0,0)^{T}$,
$(0,-\theta_{0}(1-\theta_{0}),c)^{T}$. We therefore have a center manifold
that is a graph $\tilde{\eta}=G(\theta,\zeta;\epsilon)$ over $(\theta,\zeta)$.
In particular, $G$ can be expanded as a formal series in $\zeta$ having
$(\theta,\epsilon)$-dependent coefficients:

$\displaystyle\tilde{\eta}=\sum_{k=1}^{\infty}G_{k}(\theta,\epsilon)\zeta^{k}.$

We find that each $G_{k}$, $k\geq 2$, is uniquely determined by the recursion
relation

$\displaystyle
cG_{k}=\epsilon\sum_{j=1}^{k-1}\left(\frac{\partial}{\partial\theta}G_{j}(\theta,\epsilon)\epsilon+jG_{j}(\theta,\epsilon)\right)G_{k-j}(\theta,\epsilon),\quad
k\geq 2.$

We have $G_{1}(\theta,\epsilon)=\frac{1}{c}\theta(1-\theta)$,
$G_{k}(\theta,0)\equiv 0$ for all $k\geq 2$, and the center manifold is
therefore a slow manifold (i.e. an invariant set that is
$\mathcal{O}(\epsilon)$-close to the critical manifold).

Making the ansatz $G_{k}=F_{k}\epsilon^{k-1}$, $k\geq 1$, we find that
$\epsilon^{k-1}$ cancels on both sides of the recursion relation and obtain
the following equation for $F_{k}=F_{k}(\theta,\epsilon)$:

$\displaystyle
cF_{k}=\sum_{j=1}^{k-1}\left(\frac{\partial}{\partial\theta}F_{j}(\theta,\epsilon)\epsilon+jF_{j}(\theta,\epsilon)\right)F_{k-j}(\theta,\epsilon),\quad
k\geq 2.$

This leads to the following expansion of the slow manifold of (17):

$\displaystyle\tilde{\eta}=G(\theta,\zeta;\epsilon)=\sum_{k=1}^{\infty}F_{k}(\theta;\epsilon)\epsilon^{k-1}\zeta^{k}.$
(19)

We then take the intersection of (19) with the invariant set $Q$, recall (18),
and project the result back onto the $(\theta,\eta)$-plane using (16). This
leads to a slow manifold of (6) of the form

$\eta=\epsilon^{2}G\left(\theta,\frac{1}{2}\epsilon^{-4}{\mathrm{e}}^{-\epsilon^{-1}(1-\theta)};\epsilon\right),$

which together with the expansion in (19) completes the result. ∎

Restricting system (6) to $S_{\epsilon}$, we find that

$\dot{\theta}\big{|}_{S_{\epsilon}}=\frac{1}{2}c^{-1}\epsilon^{-2}\theta(1-\theta){\mathrm{e}}^{-\epsilon^{-1}(1-\theta)}+\mathcal{O}\left(\epsilon^{-5}{\mathrm{e}}^{-2\epsilon^{-1}\theta(1-\theta)}\right)$

as $\epsilon\to 0$. This shows that $\theta$ is increasing between $0$ and
$1-\kappa\in(0,1)$ for all $\epsilon\in(0,\epsilon_{0}(\kappa))$ (i.e. between
$p_{-}$ and $p_{+}$ as $\epsilon\to 0$), as indicated in Fig. 2.

######  Remark 3.3.

The critical manifold $S$ can be extended up to $\theta=1$, i.e. up to
$p_{+}$, but Fenichel theory cannot be applied near $p_{+}$ because the right-
hand side of system (6) is not a $C^{1}$-perturbation on $\theta\in[0,1]$
(only on a compact subset, see Fig. 1). This is also evident from (14), which
shows that the validity of the series expansions for $S_{\epsilon}$ breaks
down near $\theta=1$.

######  Remark 3.4.

A more ad-hoc way of arriving at the leading order asymptotics in the second
line of (14) is to substitute the ansatz

$h(\theta,\epsilon)=\tfrac{1}{2c}\epsilon^{-2}\theta(1-\theta){\mathrm{e}}^{\epsilon^{-1}(1-\theta)}+\tilde{h}(\theta,\epsilon)$
(20)

into (6), and then use formal matching of terms to arrive at

$\tilde{h}(\theta,\epsilon)=\mathcal{O}(\epsilon^{-5}{\mathrm{e}}^{-2\epsilon^{-1}(1-\theta)}).$

We present a more detailed statement in Proposition 3.2 in order to
demonstrate the utility of a more general approach which is (i) of independent
interest as a dynamical systems approach to the derivation of asymptotic
series for flat slow manifolds, and (ii) applicable in situations where the
ansatz (20) is not ‘obvious’.

###  3.2 Equations in the reactive-diffusive zone

In order to understand the dynamics in the reactive-diffusive zone, i.e. close
to $\Sigma$, we introduce the following coordinate translation and rescaling:

$\theta=1+\epsilon\theta_{2}.$

Inserting this into system (6) leads to

$\begin{split}\theta_{2}^{\prime}&=\eta,\\\
\eta^{\prime}&=\frac{1}{2}\theta_{2}{\mathrm{e}}^{\theta_{2}}+\epsilon\left(c\eta+\frac{1}{2}\theta_{2}^{2}{\mathrm{e}}^{\theta_{2}}\right),\end{split}$
(21)

where the dash now denotes differentiation with respect to $z_{2}=z/\epsilon$,
and the subscript ‘2’ is chosen in order for consistency with the usual
notation for dynamics in the ‘rescaling chart’ in geometric blow-up analyses
(our primary tool in Section 5).

Figure 3: Phase portrait for the limiting system (22) in the reactive-
diffusive zone. The point $p_{+}$ is a resonant saddle with global
stable/unstable manifolds $w_{0}^{s/u}(p_{+})$ (shown in blue/red) given by
(24). We use a left-pointing arrow to indicate that the stable manifold
${w_{\epsilon}^{s}}(p_{+})$ is asymptotic to the point $q$ as
$\theta_{2}\to\infty$.

System (21) is a regular perturbation problem, and the limiting system

$\begin{split}\theta_{2}^{\prime}&=\eta,\\\
\eta^{\prime}&=\frac{1}{2}\theta_{2}{\mathrm{e}}^{\theta_{2}},\end{split}$
(22)

is Hamiltonian with solutions contained in level sets
$H(\theta_{2},\eta)=const.$, where

$H(\theta_{2},\eta)=\frac{\eta^{2}}{2}-\frac{\theta_{2}-1}{2}{\mathrm{e}}^{\theta_{2}}.$
(23)

System (22) has a hyperbolic saddle at $p_{+}:(0,0)$, with eigenvalues
$\lambda_{\pm}=\pm 1/\sqrt{2}$. Note that we permit a slight abuse of notation
here; earlier we defined $p_{+}$ in (7) in $(\theta,\eta)$-coordinates. Using
$H(0,0)=1/2$ we find that the the global stable and unstable manifolds of
$p_{+}$ are given by
$w_{0}^{s/u}(p_{+})=\\{(\theta_{2},h^{s/u}(\theta_{2}):\theta_{2}\in\mathbb{R}\\}$,
where

$h^{s}(\theta_{2})=-\operatorname{sign}(\theta_{2})\sqrt{1+(\theta_{2}-1){\mathrm{e}}^{\theta_{2}}},\quad
h^{u}(\theta_{2})=\operatorname{sign}(\theta_{2})\sqrt{1+(\theta_{2}-1){\mathrm{e}}^{\theta_{2}}}.$
(24)

Notice that $h^{s}(\theta_{2})\to 1$ as $\theta_{2}\to-\infty$. This suggests
that we look for solutions in the convective-diffusive zone which connect to
the point $q:=(0,1)\in\Sigma$. The geometry and dynamics described above are
sketched in Fig. 3.

###  3.3 Singular heteroclinic orbits

Using the above, we may construct an entire family of singular heteroclinic
orbits. We define

$\Gamma(c):=\Gamma^{0}(c)\cup\Gamma^{1}(c)\cup\Gamma^{2},\qquad
c\in[1,\infty),$ (25)

where

$\begin{split}\Gamma^{0}(c)&:=\left\\{(\theta,0):\theta\in\left[0,1-c^{-1}\right]\right\\},\\\
\Gamma^{1}(c)&:=\left\\{(\theta,c\theta+1-c):\theta\in\left[1-c^{-1},1\right)\right\\},\\\
\Gamma^{2}&:=\left\\{(0,\eta):\eta\in[0,1]\right\\}.\end{split}$

Notice that we only consider singular orbits $\Gamma(c)$ with $c\geq 1$, since
connections for $c<1$ enter the nonphysical domain $\theta<0$. In this sense,
there is a minimum wave speed ($c=1$) for $\epsilon=0$. The singular orbit
$\Gamma(1)$, corresponding to the minimum wave speed, is distinguished by the
fact that $\Gamma^{0}(1)=p_{-}$ and corresponds to a strong connection; recall
the discussion prior to Remark 2.3. In this case, there are only two
components:

$\Gamma(1)=\Gamma^{1}(1)\cap\Gamma^{2}.$

We sketch singular orbits with $c=1$ and $c>1$ in Fig. 4. The connections
$\Gamma(c)$ with $c>1$ correspond to weak connections, with the additional
segment $\Gamma^{0}(c)$ corresponding to a slow segment along the critical
manifold $S$.

(a) $c=1$

(b) $c>1$

Figure 4: Two representative singular heteroclinic orbits $\Gamma(c)$ (with
components shown in purple), as defined in (25), for two different cases
depending on the wave speed $c$. Left: Minimal wave speed $c=1$, for which
$\Gamma(1)=\Gamma^{1}(c)\cup\Gamma^{2}$ has only two distinct sections.
Center: A fixed (non-minimal) wave speed $c>1$, for which
$\Gamma(c)=\Gamma^{0}(c)\cup\Gamma^{1}(c)\cup\Gamma^{2}$ has three distinct
sections. Right: Schematic representation of the dimensions of the triangle
with sides $\Gamma^{1}(c)$, $\Gamma^{2}$ and $S\setminus\Gamma^{0}(c)$,
showing how $\Gamma^{1}(c)$ and $S\setminus\Gamma^{0}(c)$ (and therefore
$\Gamma^{0}(c)$) change with $c$.

##  4 The main result

We are now in a position to state the main result, namely, the existence of a
family of heteroclinic solutions for all $\epsilon>0$ sufficiently small in
system (6) which converge to the (non-smooth) singular orbits $\Gamma(c)$ as
$\epsilon\to 0$.

######  Theorem 4.1.

Consider the ZFK system (6). Fix $\kappa>0$ small enough and $\sigma>0$ large
enough. Then there exists an $\epsilon_{0}=\epsilon_{0}(\sigma)>0$ and a
$C^{\infty}$-smooth function $\overline{c}:[0,\epsilon_{0})\to\mathbb{R}$,
satisfying

$\overline{c}(0)={1},\quad\overline{c}^{\prime}(0)=\int_{-\infty}^{0}\left(1-h^{s}(x)\right){\mathrm{d}}x-1\approx
0.34405>0,$ (26)

recall the definition of $h^{s}$ in (24), such that the following assertions
hold true for each $\epsilon\in(0,\epsilon_{0})$:

  1. (i)

There are no heteroclinic connections between $p_{-}$ and $p_{+}$ within
$\theta\in[0,1]$ for $c<\overline{c}(\epsilon)$.

  2. (ii)

There is a strong heteroclinic connection when $c=\overline{c}(\epsilon)$,
whose corresponding orbit converges to $\Gamma(1)$ in Hausdorff distance as
$\epsilon\to 0$.

  3. (iii)

There is a weak heteroclinic connection for each
$c\in(\overline{c}(\epsilon),\sigma]$, and for each fixed
$c\in(1+\kappa,\sigma]$, the corresponding orbit converges to $\Gamma(c)$ in
Hausdorff distance as $\epsilon\to 0$.

The heteroclinic connection described by Theorem 4.1 corresponds to a
travelling wave solution of the ZFK equation (1). In essence, Theorem 4.1
provides a rigorous justification of formal asymptotic results that have been
derived for similar problems using HAEA, most notably the existence of
travelling wave solutions for each fixed $c$ greater than a minimum wave speed
satisfying $c=\overline{c}(\epsilon)\sim 1$; we refer again to [7, 42] and the
references therein. Assertion (i) shows that there is minimum wave speed for
$c=\overline{c}(\epsilon)$, and Assertions (ii) and (iii) describe the
asymptotic properties of the wave profile as $\epsilon\to 0$. These two
assertions describe an important qualitative difference between the cases
$c=\overline{c}(\epsilon)$ and $c>\overline{c}(\epsilon)$, namely, that they
correspond to profiles with two vs. three distinct components as $\epsilon\to
0$ respectively; see again Fig. 4. The asymptotics of the ‘extra’ component
close to $\Gamma^{0}(c)$ when $c>\overline{c}(\epsilon)$ are described by the
slow manifold asymptotics in Proposition 3.2. To the best of our knowledge,
the smoothness of the minimum wave speed function $\overline{c}(\epsilon)$ has
also not been proven. Finally, we emphasize that

$\overline{c}(\epsilon)>1\quad\forall\,0<\epsilon\ll 1,$

cf. (26).

In Figs. 5 and 6 we illustrate the results of numerical computations. The
results are in agreement with Theorem 4.1. In particular, in Fig. 5 we show
the stable manifold $W^{s}(p_{+})$ (magenta) and the strong unstable manifold
$W^{uu}(p_{-})$ (green) for $\epsilon=0.01$ and two different $c$-values:
$c=1.5$ in (a) and $c=\overline{c}(\epsilon)$ (using the linear approximation
$\overline{c}(\epsilon)\approx 1+0.34405\times\epsilon$ provided by (26)) in
(b). The results were obtained using the linear approximations offered at
$p_{\pm}$ and Matlab’s ODE45 with low tolerances ($10^{-12}$). In (a), with
$W^{s}(p_{+})$ lying below $W^{uu}(p_{-})$, we have a weak connection. The
resulting profile is shown in Fig. 6 (dashed line). Due to exponentially slow
flow along $S$, recall Proposition 3.2, the decay $\theta(z)\rightarrow 0$ for
$z\rightarrow-\infty$ is slow and it is not visible in Fig. 6(a). Fig. 5(b)
shows that the linear approximation offered by (26) for
$\overline{c}(\epsilon)$ is accurate for this value of $\epsilon$. In fact,
$W^{uu}(p_{-})$ and $W^{s}(p_{+})$ is inseparable on the scale shown (being of
the order $\sim 10^{-4}$). The resulting profile is shown in Fig. 6 (full
line). See figure captions for further details.

(a) $c=1.5$

(b) $c=\overline{c}(\epsilon)$

Figure 5: $W^{s}(1,0)$ for $c=1.5$ (in (a)) and $c=\overline{c}(\epsilon)$ (in
(b)) for $\epsilon=0.01$. We have used the linear approximation
$\overline{c}(\epsilon)\approx 1+0.34405\times\epsilon$ provided by (26) in
(b). The vector-field is illustrated in blue; notice that for illustrative
purposes, the arrows for $\theta>1$ have been moderated.

(a) $c=1.5$

(b) $c=\overline{c}(\epsilon)$

Figure 6: Profile $\theta(z)$ for $c=1.5$ (dashed) and
$c=\overline{c}(\epsilon)$ (full line), with $\epsilon=0.01$. (b) is a zoom of
(a) around the sharp decay towards $\theta=1$. For $c=1.5$ the decay towards
$\theta\rightarrow 0^{+}$ for $z\rightarrow-\infty$ is extremely slow.

######  Remark 4.2.

The travelling wave corresponding to the minimum wave speed
$c=\overline{c}(\epsilon)$ has received the most attention in the HAEA
literature [7, 42], primarily because it is expected to be stable, see e.g.
[1] for recent work on the so-called marginal stability conjecture. We do not
address questions of stability or selection in this work.

##  5 Proof and geometric blow-up

As usual, we begin by considering the extended system

$\begin{split}\dot{\theta}&=\eta,\\\
\dot{\eta}&=c\eta-\frac{1}{2}\epsilon^{-2}\theta(1-\theta){\mathrm{e}}^{-\epsilon^{-1}(1-\theta)},\\\
\dot{\epsilon}&=0,\end{split}$ (27)

which is obtained from system (6) by appending the trivial equation
$\dot{\epsilon}=0$. The extended system (27) is not defined at the interface
$\Sigma=\left\\{(1,\eta,0):\eta\in\mathbb{R}\right\\}$, or more generally on
the set defined by $\theta>1$, $\epsilon=0$. Note that we permit a slight
abuse of notation by using the same notation for $\Sigma$ as in (8) (the two
sets are naturally identified). In order to obtain a well-defined system on
$\theta>1$, $\epsilon=0$, we follow [18] and “normalise the equations” through
division of the right-hand side by

$\epsilon^{-2}\left(1+\frac{1}{2}{\mathrm{e}}^{\epsilon^{-1}(\theta-1)}\right).$
(28)

This corresponds to a reparametrization of time for $\epsilon>0$; we denote
the new time by $s$. In this way, we obtain an equivalent system for
$\epsilon>0$, having a well-defined piecewise smooth limit $\epsilon\to 0$,
with $\Sigma$ being the discontinuity set/switching manifold. In particular,

$\begin{cases}\frac{{\mathrm{d}}\theta}{{\mathrm{d}}s}&=0,\\\
\frac{{\mathrm{d}}\eta}{{\mathrm{d}}s}&=\theta(\theta-1),\end{cases}\quad\mbox{for}\quad\theta>1,$
(29)

for $\epsilon\to 0$. In line with [24], we then gain smoothness along $\Sigma$
by applying the cylindrical blow-up transformation:

$r\geq 0,\ (\bar{\theta},\bar{\epsilon})\in S^{1}\cap\\{\bar{\epsilon}\geq
0\\}\ \mapsto\ \begin{cases}\theta=1+r\bar{\theta},\\\
\epsilon=r\bar{\epsilon},\end{cases}$ (30)

which fixes $\eta$. Notice specifically that
${\mathrm{e}}^{-\epsilon^{-1}(1-\theta)}={\mathrm{e}}^{-\bar{\epsilon}^{-1}\bar{\theta}}$
under (30), which defines a smooth function on
$(\bar{\theta},\bar{\epsilon})\in S^{1}$ with $\bar{\epsilon}\geq 0$ and
$\bar{\theta}\leq 0$.

To analyse the normalised equations under the blowup transformation (30), we
work in local coordinate charts defined by $K_{1}:\bar{\theta}=-1$ and
$K_{2}:\bar{\epsilon}=1$. Local coordinates in charts are related to the
original coordinates $(\theta,\eta,\epsilon)$ by

$K_{1}:(\theta,\eta,\epsilon)=(1-r_{1},\eta,r_{1}\epsilon_{1}),\qquad
K_{2}:(\theta,\eta,\epsilon)=(1+\epsilon\theta_{2},\eta,\epsilon),$

which we insert into (27) (and apply appropriate desingularization to ensure
that the system is well-defined and smooth). The local change of coordinates
formulae are given by

$\kappa_{12}:\
r_{1}=-\epsilon\theta_{2},\quad\epsilon_{1}=-1/\theta_{2},\quad\theta_{2}<0,\qquad\kappa_{21}:\
\theta_{2}=-1/\epsilon_{1},\quad\epsilon=r_{1}\epsilon_{1},\quad\epsilon_{1}>0.$
(31)

The coordinates in chart $K_{2}$ are just the coordinates of the inner
equations used in Section 3.2, except that we now view this system in the
extended $(\theta_{2},\eta,\epsilon)$-space. (In this way, (30) also
compactifies (3.2).) Recall that (3.2) is a smooth regular perturbation
problem on compact subsets. Thus, it remains to consider the (smooth)
equations in chart $K_{1}$, where the matching between the inner and outer
system occurs.

###  5.1 Singular geometry and dynamics in $K_{1}$

The equations in chart $K_{1}$ are given by

$\begin{split}r_{1}^{\prime}&=-r_{1}\eta,\\\
\eta^{\prime}&=cr_{1}\eta-\frac{1}{2}\epsilon_{1}^{-2}(1-r_{1}){\mathrm{e}}^{-\epsilon_{1}^{-1}},\\\
\epsilon_{1}^{\prime}&=\epsilon_{1}\eta,\end{split}$ (32)

after a transformation of time which amounts to multiplication of the right-
hand side by $r_{1}$. We permit a slight abuse of notation by using a dash to
indicate differentiation with respect to the new independent variable (as
opposed to the usage in system (21)). Notice that system (32) is smooth on
$r_{1}\geq 0$, $\epsilon_{1}\geq 0$, as desired, and that the planes defined
by $r_{1}=0$ and $\epsilon_{1}=0$, as well as their intersection; the line
$r_{1}=\epsilon_{1}=0$, are invariant. In particular, the set

$L:=\\{(0,\eta,0):\eta\in\mathbb{R}\\}$

defines a line of resonant saddles for system (32), with eigenvalues $-\eta$,
$0$ and $\eta$ (except for a non-hyperbolic point at $(0,0,0)$, which does not
play a role in the analysis).

The limiting dynamics from the convective-diffusive zone appear within
$\\{\epsilon_{1}=0\\}$, and are governed by

$\begin{split}r_{1}^{\prime}&=-r_{1}\eta,\\\
\eta^{\prime}&=cr_{1}\eta.\end{split}$

Indeed, this system is equivalent to (10) on $r_{1}>0$ upon setting
$r_{1}=1-\theta$. We therefore also re-discover the repelling critical
manifold $S$ as

$S_{1}=\\{(r_{1},0,0):r_{1}\geq 0\\}.$

Considered within $\\{\epsilon_{1}=0\\}$, the Jacobian evaluated along $S_{1}$
has eigenvalues $0$ and $cr_{1}\geq 0$. Trajectories are contained in straight
lines of the form $\eta(r_{1})=-cr_{1}+\eta(0)$. On the other hand, the
limiting dynamics in the reactive-diffusive zone appear within
$\\{r_{1}=0\\}$. We are particularly interested in the extension of
${w_{0}^{s}}(p_{+})$ into chart $K_{1}$. Using (24) and (31), we obtain

${w_{0,1}^{s}(p_{+})}:=\kappa_{21}\left({w_{0}^{s}}(p_{+})\times\\{0\\}\right)=\\{(0,h_{1}^{s}(\epsilon_{1}),\epsilon_{1}):\epsilon_{1}\geq
0\\},$ (33)

where
$h_{1}^{s}(\epsilon_{1}):=\sqrt{1-\epsilon_{1}^{-1}(1+\epsilon_{1}){\mathrm{e}}^{-\epsilon_{1}^{-1}}}$,
the subscript “1” indicates that we are viewing the object in chart $K_{1}$,
and the solution within $w_{0,1}^{s}(p_{+})$ is backward asymptotic to the
point $q:(0,1,0)\in L$. The singular geometry and dynamics is summarised in
Fig. 7.

Figure 7: Singular geometry and dynamics after applying the blow-up
transformation (30) and desingularising in order to resolve the loss of
smoothness along $\Sigma$. In contrast to earlier figures, we also indicate
the dynamics for $\theta>1$, which is well-defined after the normalization of
the system through division of (27) by the quantity (28), see also (29) and
[18, 24] for further details on this procedure. The dynamics in the
convective-diffusive and reactive-diffusive zones extend and ‘connect’ along a
line of hyperbolic saddles $L$, shown here in magenta, which contains the
point $q$ identified as the asymptote of the stable manifold
${w_{0}^{s}}(p_{+})$ in the Subsection 3.2. As in Fig. 4, we show two
different scenarios depending on whether $c=1$ or $c>1$. Left: $c=1$, with the
improved singular heteroclinic $\Gamma^{1}(1)\cup w_{0}^{s}(p_{+})$. Right:
$c>1$, with the improved singular heteroclinic
$\Gamma^{0}(c)\cup\Gamma^{1}(c)\cup w_{0}^{s}(p_{+})$.

The preceding analysis allows for the construction of an improved family of
singular heteroclinics

$\Gamma(c)=\Gamma^{0}(c)\cup\Gamma^{1}(c)\cup\Gamma^{2},$ (34)

considered now in the blown-up space. The definitions for $\Gamma^{0}(c)$ and
$\Gamma^{1}(c)$ remain unchanged (after a natural embedding into the extended
$(\theta,\eta,\epsilon)$-space or representation in chart $K_{1}$
coordinates), but $\Gamma^{2}$ is replaced by

$\Gamma^{2}=q\cup{({w_{0}^{s}(p_{+})}\times\\{0\\})}.$

Improved singular heteroclinic orbits for $c=1$ and $c>1$ are shown in Fig. 7.

###  5.2 Perturbation

It remains to describe the perturbation of the improved singular heteroclinics
in (34). The main task is to understand the local passage near $q\in L$. To
this end, we divide the right-hand side of system (35) by $\eta$ (which is
positive close to $q$) and consider the flow induced by the system

$\begin{split}r_{1}^{\prime}&=-r_{1},\\\
\eta^{\prime}&=cr_{1}-\frac{1}{2}\eta^{-1}\epsilon_{1}^{-2}(1-r_{1}){\mathrm{e}}^{-\epsilon_{1}^{-1}},\\\
\epsilon_{1}^{\prime}&=\epsilon_{1}.\end{split}$ (35)

System (35) can be brought into a kind of local normal form by straightening
the fibers of the stable and unstable manifolds with base points along $L$,
which are contained in the invariant planes $\\{r_{1}=0\\}$ and
$\\{\epsilon_{1}=0\\}$ respectively. This can be achieved with an
$(r_{1},\epsilon_{1})$-fibered diffeomorphism defined by

$(r_{1},y_{1},\epsilon_{1})\mapsto\eta=-cr_{1}+f_{1}(y_{1},\epsilon_{1}),$
(36)

where

$f_{1}(y_{1},\epsilon_{1})=y_{1}\sqrt{1-y_{1}^{-2}\left(\epsilon_{1}^{-1}+1\right){\mathrm{e}}^{-\epsilon_{1}^{-1}}}$

and the inverse is given by

$y_{1}=\sqrt{(\eta+cr_{1})^{2}+\left(\epsilon_{1}^{-1}+1\right){\mathrm{e}}^{-\epsilon_{1}^{-1}}}.$
(37)

Indeed, in these coordinates the sets defined by $(r_{1},y_{1},0),\,r_{1}\geq
0$ and $(0,y_{1},\epsilon_{1}),\,\epsilon_{1}\geq 0$ with $y_{1}>0$ fixed are
the stable and unstable manifolds of the point $(0,y_{1},0)$; see Fig. 8. Here
we have used the conservation of $H(-\epsilon_{1}^{-1},\eta)=\text{const.}$
within $r_{1}=0$ to write the fibers of the unstable manifold with base point
$(r_{1},y_{1},\epsilon_{1})=(0,y_{1},0)$ as

$r_{1}=0,\quad H(-\epsilon_{1}^{-1},\eta)=\frac{y_{1}^{2}}{2}.$ (38)

Solving the last expression for $\eta$ gives $\eta=f_{1}(y_{1},\epsilon_{1})$.
Notice that

$\displaystyle f_{1}(1,\epsilon_{1})=h^{s}(-\epsilon_{1}^{-1}),$ (39)

since $\eta=h^{s}(\theta_{2})$ corresponds to the level set
$H(\eta,\theta_{2})=\frac{1}{2}$ (i.e. $y_{1}=1$). This leads to the following
equations:

$\displaystyle\dot{r}_{1}$ $\displaystyle=-r_{1},$ (40)
$\displaystyle\dot{y}_{1}$
$\displaystyle=r_{1}F_{1}(r_{1},y_{1},\epsilon_{1},c),$
$\displaystyle\dot{\epsilon}_{1}$ $\displaystyle=\epsilon_{1},$

where

$F_{1}(r_{1},y_{1},\epsilon_{1},c):=\frac{1}{2}y_{1}^{-1}\epsilon_{1}^{-2}{\mathrm{e}}^{-\epsilon_{1}^{-1}}\frac{f_{1}(y_{1},\epsilon_{1})-c}{f_{1}(y_{1},\epsilon_{1})-cr_{1}}.$
(41)

Notice that $F_{1}$ is infinitely flat in $\epsilon_{1}\rightarrow 0$, i.e.

$\frac{\partial^{n}}{\partial\epsilon_{1}^{n}}F_{1}(r_{1},y_{1},0,c)=0,$ (42)

for all $n\in\mathbb{N}$. We shall adopt the following notation: for any
$N\in\mathbb{N}\cup\\{\infty\\}$, we let $j_{N}G(\cdot)$ denote the $N$th-
order Taylor-jet of a $C^{N}$-smooth function $x\mapsto G(x)$ defined in a
neighborhood of $x=0$. Then we can write (42) as
$j_{\infty}F_{1}(r_{1},y_{1},\cdot,c)=0$.

Figure 8: Singular geometry and dynamics in chart $K_{1}$ before (left) and
after (right) applying the normal form transformation (36), which rectifies
the stable and unstable manifolds within the invariant planes $\\{r_{1}=0\\}$
and $\\{\epsilon_{1}=0\\}$ locally near $L$.

System (40) is similar (but not identical) to the normal form in [9, Prop.
3.2] (here $F_{1}$, called $k$ in [9, Eq. (2.13)], only depends upon $r_{1}$
through $r_{1}\epsilon_{1}=\epsilon$). Let
$\pi:\Delta_{1}^{in}\to\Delta_{1}^{out}$ denote the transition map induced by
the flow of system (40), where

$\begin{split}\Delta^{in}_{1}&=\left\\{(\rho,y_{1},\epsilon_{1}):y_{1}\in
I,\epsilon_{1}\in[0,\delta]\right\\},\\\
\Delta^{out}_{1}&=\left\\{(r_{1},y_{1},\delta):r_{1}\in[0,\rho],y_{1}\in\tilde{I}\right\\},\end{split}$
(43)

with $I:=[1-\beta,1+\beta]\subset(0,\infty)$ and
$\tilde{I}:=[1-\tilde{\beta},1+\tilde{\beta}]$ for some fixed $\beta\in(0,1)$,
$\tilde{\beta}>\beta$, and small but fixed $\rho,\delta>0$. Due to resonances
along $L$, one would in general expect logarithmic terms in the expansion of
$\pi$ (see e.g. [9]). However, all resonant terms are absent in the present
case due to the flatness (42) with respect to $\epsilon_{1}$. Indeed, we have
the following:

######  Proposition 5.1.

Fix any $N\in\mathbb{N}$ and a compact interval $J\subset(0,\infty)$. Then for
$\rho>0,\delta>0$ sufficiently small, there exists a $C^{N}$-smooth function
$Y:I\times[0,\delta)\times J\rightarrow\mathbb{R}$ such that the map
$\pi:\Delta_{1}^{in}\to\Delta_{1}^{out}$ is well-defined and given by

$\pi:(\rho,y_{1},\epsilon_{1})\mapsto\left(\frac{\rho}{\delta}\epsilon_{1},Y(y_{1},\epsilon_{1},c;\rho,\delta),\delta\right).$
(44)

###### Proof.

The first component of the right-hand side in (44) follows directly from
$\epsilon=\epsilon_{1}r_{1}=\text{const.}$ It therefore remains to show that
the function $Y$ is well-defined and $C^{N}$-smooth. In order to do so, we
show that system (40) can be brought into the local normal form in [9, Prop.
3.2], which, based on additional arguments in [9], is sufficient to prove
$C^{N}$-smoothness of $\pi$. Since $\rho,\delta>0$ are fixed as parameters, we
shall suppress the dependency of $Y$ on these quantities. We also suppress the
dependency on $c$; the smoothness with respect to this parameter should be
obvious from the proof.

First, we claim that for any $N\in\mathbb{N}$ there is an
$(r_{1},\epsilon_{1})$-fibered $C^{\infty}$-diffeomorphism

$\displaystyle y_{1}=y_{N}+\alpha_{N}(r_{1},y_{N},\epsilon_{1}),$

with $j_{\infty}\alpha_{N}(r_{1},y_{N},\cdot)=0$, such that

$\begin{split}\dot{r}_{1}&=-r_{1},\\\
\dot{y}_{N}&=r_{1}^{N}F_{N}(r_{1},y_{N},\epsilon_{1}),\\\
\dot{\epsilon}_{1}&=\epsilon_{1},\end{split}$ (45)

with $j_{\infty}F_{N}(r_{1},y_{N},\cdot)=0$. Following [9], we proceed by
induction. We have already shown the base case $N=1$ in the derivation of
system (40) using (36). We therefore proceed directly to the induction step
and suppose that the statement is true for $N=n$. We write

$y_{n+1}=y_{n}+r_{1}^{n}\beta_{n}(y_{n},\epsilon_{1}),$ (46)

with the (as yet undetermined) $\beta_{n}$ satisfying
$j_{\infty}\beta_{n}(y_{n},\cdot)=0$. Writing

$F_{n}(r_{1},y_{n},\epsilon_{1})=F_{n,0}(y_{n},\epsilon_{1})+r_{1}F_{n,1}(r_{1},y_{n},\epsilon_{1}),$

we obtain

$\dot{y}_{n+1}=r_{1}^{n}\left\\{F_{n,0}(y_{n+1},\epsilon_{1})-n\beta_{n}(y_{n},\epsilon_{1})+\epsilon_{1}\frac{\partial}{\partial\epsilon_{1}}\beta_{n}(y_{n},\epsilon_{1})\right\\}+\mathcal{O}(r_{1}^{n+1}).$
(47)

The $\mathcal{O}(r_{1}^{n})$ terms vanish if we impose the following
solvability condition on $\beta_{n}$:

$\epsilon_{1}\frac{\partial}{\partial\epsilon_{1}}\beta_{n}(y_{n},\epsilon_{1})=n\beta_{n}(y_{n},\epsilon_{1})-F_{n,0}(y_{n},\epsilon_{1}),\quad
j_{\infty}\beta_{n}(y_{n},\cdot)=0,$ (48)

which may be solved to give

$\beta_{n}(y_{n},\epsilon_{1})=-\epsilon_{1}\int_{0}^{\epsilon_{1}}s^{-2}F_{n,0}(y_{n},s){\mathrm{d}}s.$
(49)

Now fix $N+3\in\mathbb{N}$ and use $\epsilon=r_{1}\epsilon_{1}$ to rewrite
system (45) as

$\begin{split}\dot{r}_{1}&=-r_{1},\\\
\dot{y}_{N+3}&=\epsilon^{N+3}R(r_{1},y_{N+3},\epsilon_{1}),\\\
\dot{\epsilon}_{1}&=\epsilon_{1},\end{split}$ (50)

where $R:=\epsilon_{1}^{-(N+3)}F_{N+3}$. Notice that this is smooth since
$j_{\infty}F_{N+3}(r_{1},y_{n},\cdot)=0$. System (50) is – up to a change of
notation – precisely the local normal form appearing in [9, Prop. 3.2] (case
$a=0$). $C^{N}$-smoothness of the function $Y$ then follows from the arguments
in [9, Sec. 3.2]. ∎

Using Proposition 5.1 and a number of elements of the proof thereof, we can
derive the leading order asymptotics for $Y(y_{1},\epsilon_{1},c)$.

######  Lemma 5.2.

The function $Y$ has the following asymptotic expansion in $\epsilon_{1}$:

$Y(y_{1},\epsilon_{1},c)=y_{1}+\epsilon_{1}\rho\int_{0}^{\delta}s^{-2}F_{1}(0,y_{1},s,c){\mathrm{d}}s+\mathcal{O}(\epsilon_{1}^{2}),$
(51)

where the function $F_{1}$ is given in (41).

###### Proof.

From (36) and the proof of Proposition 5.1 we have the following relations:

$y_{N+3}^{out}=y_{N+3}^{in}+\mathcal{O}(\epsilon^{N}),\quad
y_{1}=y_{N+3}+\alpha_{N+3}(r_{1},y_{N+3},\epsilon_{1}),\quad
y_{N+3}=y_{1}+\tilde{\alpha}_{1}(r_{1},y_{1},\epsilon_{1}),$

where $\tilde{\alpha}_{1}(r_{1},y_{1},\epsilon_{1})$ is $C^{N}$-smooth and
$\mathcal{O}(\epsilon_{1})$. Using these relations to evaluate
$y_{1}^{out}=Y(y_{1},\epsilon_{1},c)$, we obtain

$Y(y_{1},\epsilon_{1},c)=y_{1}+\hat{\alpha}_{1}(\rho,y_{1},\epsilon_{1}),$
(52)

where the function $\hat{\alpha}_{1}(\rho,y_{1},\epsilon_{1})$ is
$C^{N}$-smooth and $\mathcal{O}(\epsilon_{1})$ as $\epsilon_{1}\to 0$. In
order to determine $\hat{\alpha}_{1}(\rho,y_{1},\epsilon_{1})$ to leading
order in $\epsilon_{1}\to 0$, we use the transformation
$y_{2}=y_{1}+\alpha_{1}(r_{1},y_{1},\epsilon_{1})=y_{1}+r_{1}\beta_{1}(y_{1},\epsilon_{1})$
and the fact that

$y_{2}^{out}=y_{2}^{in}+\hat{\alpha}_{2}(\rho,y_{2}^{in},\epsilon_{1}^{in}),$

where the function $\hat{\alpha}_{2}(\rho,y_{2}^{in},\epsilon_{1}^{in})$ is
$C^{N-1}$-smooth and $\mathcal{O}((\epsilon_{1}^{in})^{2})$ as
$\epsilon_{1}^{in}\to 0$ (this follows from relations analogous to those which
led to (52)). We obtain

$\begin{split}Y(y_{1},\epsilon_{1},c)&=y_{2}^{out}-\frac{\rho}{\delta}\epsilon_{1}\beta_{1}(Y(y_{1},\epsilon_{1},c),\delta)\\\
&=y_{2}^{in}-\frac{\rho}{\delta}\epsilon_{1}\beta_{1}(y_{1},\delta)+\mathcal{O}(\epsilon_{1}^{2})\\\
&=y_{1}-\frac{\rho}{\delta}\epsilon_{1}\beta_{1}(y_{1},\delta)+\mathcal{O}(\epsilon_{1}^{2}).\end{split}$

Substituting the expression for $\beta_{1}(y_{1},\delta)$ in (49) leads to
(51), as required. ∎

Taken together, Proposition 5.1 and Lemma 5.2 imply that the transition map
$\pi:\Delta_{1}^{in}\to\Delta_{1}^{out}$ is $C^{N}$-smooth of the form (44),
with $Y(y_{1},\epsilon_{1},c)$ given by (51). We now use Lemma 5.2 in order to
derive a bifurcation equation on the section $\Delta_{1}^{out}$. More
precisely, we introduce a distance function on $\Delta_{1}^{out}$ via

$\mathcal{D}(c,\epsilon):={Y_{u}(c,\epsilon)-Y_{s}(c,\epsilon)}.$ (53)

Here $Y_{u}(c,\epsilon)$ denotes the $y_{1}$-coordinate of the intersection
${w_{\epsilon,1}^{u}}(p_{-})\cap\Delta_{1}^{out}$, where
${w_{\epsilon,1}^{u}}(p_{-})$ denotes the (unique) strong unstable manifold of
the point $p_{-}$ in $K_{1}$ coordinates, and $Y_{s}(c,\epsilon)$ denotes the
$y_{1}$-coordinate of the intersection
${w_{\epsilon,1}^{s}}(p_{+})\cap\Delta_{1}^{out}$, where
${w_{\epsilon,1}^{s}}(p_{+})$ denotes the (unique) stable manifold of the
saddle $p_{+}$, which is a perturbation of the limiting stable manifold in
(33). Notice that ${w_{0,1}^{u}}(p_{{-}})=\Gamma^{1}(1)$,
${w_{0,1}^{s}}(p_{+})=\Gamma^{2}{\setminus p_{+}}$, and that zeros of (53)
correspond to heteroclinic solutions of system (6).

######  Lemma 5.3.

For each $N\in\mathbb{N}$ there exists an $\epsilon_{0}=\epsilon_{0}(N)>0$
such that $Y_{u}$ is $C^{N}$-smooth on
$(c,\epsilon)\in(1-\beta,1+\beta)\times[0,\epsilon_{0})$ and satisfies

${Y_{u}(c,\epsilon)}=c+\epsilon\int_{0}^{\delta}s^{-2}F_{1}(0,c,s,c){\mathrm{d}}s+\mathcal{O}(\epsilon^{2})$
(54)

as $\epsilon\to 0$.

###### Proof.

Fixing $c\in(1-\beta,1+\beta)$ ensures that ${w_{0}^{u}}(p_{-})$, as an orbit
of (10), intersects $\Sigma$ in a point $(1,c)$. Regular perturbation theory
on compact subsets of $\theta<1$ implies that the fast fiber corresponding to
this connection lifts to a $C^{\infty}$-smooth curve of the form

$\gamma^{in}(c,\epsilon):=\left\\{(\rho,c+\mathcal{O}(\epsilon^{N}),\rho^{-1}\epsilon):\epsilon\in[0,\epsilon_{0}]\right\\}\subset\Delta_{1}^{in},$
(55)

when viewed in the $(r_{1},y_{1},\epsilon_{1})$-coordinates, where
$\epsilon_{0}:=\delta\rho$, see also Lemma 3.1. Note that
$\gamma^{in}(c,\epsilon)\subset\Delta_{1}^{in}$ as long as we choose $\delta$
and $\rho$ sufficiently small. Thus we can use Proposition 5.1 and Lemma 5.2
in order to extend ${w_{\epsilon}^{u}}(p_{-})$ to $\Delta_{1}^{out}$.
Expanding in $\epsilon$ and using (51) together with the fact that
$\epsilon=\epsilon_{1}\rho$ for points on $\Delta_{1}^{in}$, we obtain

${Y_{u}(c,\epsilon)}=Y\left(1+\mathcal{O}(\epsilon^{N}),\rho^{-1}\epsilon,c\right)=c+\epsilon\int_{0}^{\delta}s^{-2}F_{1}(0,c,s,c){\mathrm{d}}s+\mathcal{O}(\epsilon^{2})$

as $\epsilon\to 0$, as required. ∎

We now derive an implicit equation for $Y_{s}(c,\epsilon)$.

######  Lemma 5.4.

For each $N\in\mathbb{N}$ there exists an $\epsilon_{0}=\epsilon_{0}(N)>0$
such that $Y_{s}(c,\epsilon)$ is $C^{N}$-smooth on
$c\in(1-\beta,1+\beta)\times[0,\epsilon_{0})$ and satisfies

$\frac{{Y_{s}(c,\epsilon)}^{2}-1}{2}-\epsilon\left[c\delta^{-1}f({Y_{s}(c,\epsilon)},\delta)-\int_{\delta}^{\infty}\left(cs^{-2}f_{1}(1,s)+\frac{s^{-4}}{2}{\mathrm{e}}^{-s^{-1}}\right){\mathrm{d}}s\right]+\mathcal{O}(\epsilon^{2})=0$
(56)

as $\epsilon\to 0$.

###### Proof.

We begin in $K_{2}$ coordinates $(\eta,\theta_{2},\epsilon)$, and look for an
expression for the intersection
${w_{\epsilon}^{s}}(P_{s})\cap\kappa_{12}(\Delta_{1}^{out})$, where
$P_{s}:=\\{(0,0,\epsilon):\epsilon\in[0,\epsilon_{0}]\\}$ denotes the line of
equilibria emanating from $p_{+}$. We use the fact that trajectories are
contained in constant level sets of the Hamiltonian function
$H(\theta_{2},\eta)$ defined in (23) when $\epsilon=0$, and consider now the
perturbed dynamics for $\epsilon\in[0,\epsilon_{0}]$. We have

$\frac{{\mathrm{d}}}{{\mathrm{d}}\theta_{2}}H=\epsilon\left(c\eta+\frac{\theta_{2}}{2}{\mathrm{e}}^{\theta_{2}}\right),$

using $H=H(\theta_{2},\eta(\theta_{2}))$. Along the stable manifold
$\eta(\theta_{2};c,\epsilon)=h^{s}(\theta_{2})+\mathcal{O}(\epsilon)$, where
$h^{s}(\theta_{2})$ is given by (24), regular perturbation theory implies

$\frac{{\mathrm{d}}}{{\mathrm{d}}\theta_{2}}H=\epsilon\left(ch^{s}(\theta_{2})+\frac{\theta_{2}}{2}{\mathrm{e}}^{\theta_{2}}\right)+\mathcal{O}(\epsilon^{2}),$

since we are working on compact subsets. Integrating this expression from
$\theta_{2}=0$ (where the saddle is) to $\theta_{2}=-1/\delta$ (where
$\kappa_{12}(\Delta_{1}^{out})$ is), we obtain

$H(-1/\delta,\eta^{out})-\frac{1}{2}=\epsilon\int_{0}^{-1/\delta}\left(ch^{s}(\theta_{2})+\frac{\theta_{2}^{2}}{2}{\mathrm{e}}^{\theta_{2}}\right){\mathrm{d}}\theta_{2}+\mathcal{O}(\epsilon^{2}),$
(57)

where $\eta^{out}$ denotes the $\eta$-coordinate of the intersection
${w_{\epsilon}^{s}}(P_{s})\cap\kappa_{12}(\Delta_{1}^{out})$ to be solved for.
Using (38), we have

$\displaystyle H(-1/\delta,\eta^{out})$
$\displaystyle=\frac{(y_{1}^{out})^{2}}{2}-\frac{\partial}{\partial\eta}H(-1/\delta,f_{1}(y_{1}^{out}))cr_{1}^{out}+\mathcal{O}((r_{1}^{out})^{2})$
$\displaystyle=\frac{(y_{1}^{out})^{2}}{2}-cr_{1}^{out}f_{1}(y_{1}^{out},\delta)+\mathcal{O}((r_{1}^{out})^{2})$

where $y_{1}^{out}={Y_{s}}(c,\epsilon)$ is given in terms of $\eta^{out}$ by
the right-hand side of (36) evaluated at $(r_{1}^{out},\eta^{out},\delta)$,
where $r_{1}^{out}=\delta^{-1}\epsilon$. Using (39) and (57), we obtain the
equation

$\frac{(y_{1}^{out})^{2}}{2}-\frac{1}{2}-\epsilon\left[c\delta^{-1}f(y_{1}^{out},\delta)-\int_{\delta}^{\infty}\left(cs^{-2}f_{1}(1,s)+\frac{s^{-4}}{2}{\mathrm{e}}^{-s^{-1}}\right){\mathrm{d}}s+\mathcal{O}(\epsilon)\right]=0,$

as required. ∎

In order to identify the zeros of $\mathcal{D}(c,\epsilon)$, it suffices to
consider solutions to the bifurcation equation

$\widetilde{B}(c,\epsilon):=B({Y_{u}}(c,\epsilon),c,\epsilon)=0,$ (58)

where the function $B(X,c,\epsilon)$ is defined using the implicit equation
(56), i.e.

$B(X,c,\epsilon):=\frac{X^{2}-1}{2}-\epsilon\left[c\delta^{-1}f(X,\delta)-\int_{\delta}^{\infty}\left(cs^{-2}f_{1}(1,s)+\frac{s^{-4}}{2}{\mathrm{e}}^{-s^{-1}}\right){\mathrm{d}}s\right]+\mathcal{O}(\epsilon^{2}).$
(59)

The following result establishes a solution to (58).

######  Lemma 5.5.

For each $N\in\mathbb{N}$ there exists an $\epsilon_{0}=\epsilon_{0}(N)>0$ and
a $C^{N}$-smooth function $\overline{c}:[0,\epsilon_{0})\to[1,\infty)$ such
that $\widetilde{B}(\overline{c}(\epsilon),\epsilon)=0$ for all
$\epsilon\in[0,\epsilon_{0})$. In particular,

${\overline{c}(0)={1},\quad\overline{c}^{\prime}(0)=\int_{-\infty}^{0}\left(1-h^{s}(x)\right){\mathrm{d}}x-1\approx
0.34405>0.}$

###### Proof.

Using (58) and (59), we obtain $\widetilde{B}(1,0)=0$ and
$\widetilde{B}_{c}^{\prime}(1,0)=1\neq 0$. Thus, the existence of the
$C^{N}$-smooth function $\epsilon\mapsto\overline{c}(\epsilon)$ satisfying
$\overline{c}(0)=1$ follows from the $C^{N}$-smoothness of $\widetilde{B}$ and
the implicit function theorem.

It remains to calculate $\overline{c}^{\prime}(0)$. Since
$\overline{c}^{\prime}(0)=-\widetilde{B}^{\prime}_{\epsilon}(1,0)/\widetilde{B}^{\prime}_{c}(1,0)=-\widetilde{B}^{\prime}_{\epsilon}(1,0)$
by implicit differentiation, it suffices to determine
$\widetilde{B}^{\prime}_{\epsilon}(1,0)$. From (58), using the expression for
$Y_{u}(c,\epsilon)$ in (54), we obtain

$\widetilde{B}_{\epsilon}^{\prime}(1,0)=\int_{0}^{\delta}s^{-2}F_{1}(0,1,s,1){\mathrm{d}}s-\delta^{-1}f(1,\delta)+\int_{\delta}^{\infty}\left(s^{-2}f_{1}(1,s)+\frac{s^{-4}}{2}{\mathrm{e}}^{-s^{-1}}\right){\mathrm{d}}s,$

where we also used the fact that $\overline{c}(0)=1$. In order to simplify the
right-hand side, we use the definitions of $f_{1}$ and $F_{1}$ given in (36)
and (41), together with
$\int_{0}^{\infty}\frac{s^{-4}}{2}{\mathrm{e}}^{-s^{-1}}{\mathrm{d}}s=1$. We
obtain

$\widetilde{B}_{\epsilon}^{\prime}(1,0)=1+\left[-\delta^{-1}f_{1}(1,\delta)+\frac{1}{2}\int_{0}^{\delta}\frac{s^{-4}{\mathrm{e}}^{-s^{-1}}}{f_{1}(1,s)}{\mathrm{d}}s+\int_{\delta}^{\infty}s^{-2}f_{1}(1,s){\mathrm{d}}s\right].$
(60)

Differentiating with respect to $\delta$ shows that (60) is independent of
$\delta$. Therefore,

$\begin{split}\widetilde{B}_{\epsilon}^{\prime}(1,0)&=\lim_{\delta\to
0}\widetilde{B}_{\epsilon}^{\prime}(1,0)\\\ &=1+\lim_{\delta\rightarrow
0}\left(\int_{\delta}^{\infty}s^{-2}f_{1}(1,s){\mathrm{d}}s-\delta^{-1}f_{1}(1,\delta)\right)\\\
&=1+\int_{0}^{\infty}s^{-2}\left(f_{1}(1,s)-f_{1}(1,0)\right){\mathrm{d}}s\\\
&=1+\int_{-\infty}^{0}\left(h^{s}(\theta_{2})-1\right){\mathrm{d}}\theta_{2},\end{split}$

as required. ∎

###  5.3 Completing the proof of Theorem 4.1

Lemma 5.5 shows the existence of a $C^{N}$-smooth minimal wave speed function
$\overline{c}(\epsilon)$ satisfying (26) and Assertion (ii) of Theorem 4.1. To
see that $\overline{c}$ is $C^{\infty}$, we first notice that $\overline{c}$
is defined on a domain $[0,\epsilon_{0})$ with $\epsilon_{0}>0$ small enough
(we can apply Lemma 5.5 with $N=1$). Clearly $\overline{c}$ is $C^{\infty}$
for any $\epsilon\in(0,\epsilon_{0})$, since the original ZFK system (6) is
$C^{\infty}$. We therefore only have to show that $\overline{c}$ is
$C^{\infty}$ at $\epsilon=0$. For this purpose, we use the uniqueness of
$\overline{c}(\epsilon)$ for $\epsilon\in(0,\epsilon_{0})$ and Lemma 5.5 to
conclude that $\overline{c}$ is $C^{N}$ on $[0,\epsilon_{0}(N))$ for some
$0<\epsilon_{0}(N)\ll 1$ for any $N\in\mathbb{N}$. The claim now follows.

Next, we notice that all heteroclinic connections with
$c\neq\overline{c}(\epsilon)$ are weak, due to the uniqueness of the strong
heteroclinic connection. In particular, for $c<\overline{c}(\epsilon)$ the
heteroclinic connections enter $\theta<0$. This shows Assertion (i) of Theorem
4.1. To prove Assertion (iii), we just have to show that the connections limit
to $\Gamma(c)$ in Hausdorff distance. For this purpose, we first notice that
for $c>\overline{c}(\epsilon)$ the stable manifold $w_{\epsilon}^{s}(p_{+})$
lies below the strong unstable manifold $w_{\epsilon}^{u}(p_{-})$. Moreover,
$w_{0}^{s}{(p_{+})}$ and $\Gamma^{1}(c)$ (as a trajectory of (10)) have
${q}:(r_{1},\eta,\epsilon_{1})=(0,1,0)$ as their $\alpha$-limit resp.
$\omega$-limit sets in the $K_{1}$-chart, as shown in Figure 7. It then
follows from Proposition 5.1 that the intersection of
$w_{\epsilon}^{s}(p_{+})$ with $\Delta^{in}=\\{(\theta,\eta):\theta=1-\rho\\}$
is $\mathcal{O}(\epsilon)$-close to $\Gamma^{2}(c)\cap\Delta^{in}$. Here we
abuse notation slightly by using $\Delta^{in}$ for the blown-down version of
(43). The result then follows from Fenichel theory and Proposition 3.2.

##  6 Summary and outlook

A primary aim of this work was to demonstrate the suitability of GSPT and
geometric blow-up for the study of dynamical systems with singular exponential
nonlinearities. We would like to emphasise and reiterate the following
methodological point, which is supported by our analysis of the ZFK equation
(1) as well as the analysis of electrical oscillator models in [18]:

Smooth dynamical systems with singular exponential nonlinearities can often be
formulated as singular perturbation problems with a piecewise-smooth singular
limit, after which the geometry and dynamics can be analysed using GSPT and
geometric blow-up.

This means that different limiting problems are to be expected on different
regions of the phase space which are separated by switching manifolds, along
which the system loses smoothness. Because the system is smooth for
$\epsilon>0$, one expects to be able to find a smooth connection between these
regions by inserting an inner ‘rescaling regime’, where another limiting
problem arises, after geometric blow-up; we refer again to [24]. In the case
of system (6), the singular exponential nonlinearity in the reaction term
$\omega(\theta,\epsilon^{-1})$ leads to a non-smooth singular limit in a
‘normalised system’ which is obtained from (6) by dividing the the right-hand
side by the positive expression in (28), since

$\frac{\omega(\theta,\epsilon^{-1})}{\epsilon^{-2}\left(1+\frac{1}{2}{\mathrm{e}}^{\epsilon^{-1}(1-\theta)}\right)}=\begin{cases}0,&\theta<1,\\\
\theta(1-\theta),&\theta>1,\end{cases}$

as $\epsilon\to 0$. This leads to different limiting problems on three
distinct regimes: $\theta<1$, $\theta=1+\mathcal{O}(\epsilon)$ and $\theta>1$.
In order to identify and characterise travelling wave solutions to the ZFK
equation (1), it sufficed to consider the dynamics in two regimes only,
namely, in the convective-diffusive zone wih $\theta<1$ bounded away from
$\Sigma=\\{\theta=1\\}$, and a reactive-diffusive zone with
$\theta=1+\mathcal{O}(\epsilon)$. After extending the phase space and applying
the blow-up transformation (30), we were able to smoothly connect the analyses
in each of these zones, which we performed separately in Section 3. The
connection problem in chart $K_{1}$ within the blow-up is smooth, even though
the singular limit in the original system (6) is not, showing that the blow-up
has ‘resolved’ the loss of smoothness.

Our main result is Theorem 4.1, which describes a $(c,\epsilon)$-family of
heteroclinic orbits in system (6) which correspond to travelling waves the
original ZFK equation (1). In essence, Theorem 4.1 is the rigorous and
geometrically informative counterpart of existing results on travelling waves
in (a number of closely related variants of) the ZFK equation, which have been
obtained formally using HAEA in e.g. [5, 6, 7, 42, 43]. We also moved beyond
the existing formal results by proving that the minimal wave speed function
$\overline{c}(\epsilon)$ is $C^{\infty}$ for all $\epsilon\in[0,\epsilon_{0})$
(Proposition 5.1 and the normal form in [9] was crucial for this), and
providing an asymptotic expansion for the slow manifold $S_{\epsilon}$ (see
again Proposition 3.2) which plays an important role in the construction of
the non-minimal waves for $c>\overline{c}(\epsilon)$. The latter feature is
not likely to be of significance for the ZFK problem in particular, given that
one expects the system to select the minimal wave speed
$\overline{c}(\epsilon)$ (recall Remark 4.2). Nevertheless, the proof of
Proposition 3.2 is of independent interest because it provides a rigorous and
dynamical systems based approach to the determination of asymptotic expansions
for flat slow manifolds.

This work paves the way for more complicated analyses, particularly because of
the emphasis on the development of transferable methods for study of dynamical
systems with exponential nonlinearities in general. We also left questions of
stabilitiy and selection completely open, see e.g. Remark 4.2. With regard to
problems in combustion theory, a natural continuation of this work would be to
consider the case of planar flames with non-unity Lewis number, in which case
one considers two species reaction-diffusion equations with different
diffusivity constants as in [7, eqn. (2.1.4)]. These and other related
problems remain for future work.

## References

  * [1] M. Avery and A. Scheel, Universal selection of pulled fronts, Communications of the American Mathematical Society, 2 (2022), pp. 172–231.
  * [2] C. Bonet-Revés and T. M-Seara, Regularization of sliding global bifurcations derived from the local fold singularity of Filippov systems, Discrete and Continuous Dynamical Systems- Series a, 36 (2016), pp. 3545–3601.
  * [3] E. Bossolini, M. Brøns, and K. U. Kristiansen, Singular limit analysis of a model for earthquake faulting, Nonlinearity, 30 (2017), p. 2805.
  * [4] M. Brøns, Relaxation oscillations and canards in a nonlinear model of discontinuous plastic deformation in metals at very low temperatures, in Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 461, The Royal Society, 2005, pp. 2289–2302.
  * [5] J. D. Buckmaster and G. S. S. Ludford, Theory of laminar flames, (No Title), (1982).
  * [6] W. B. Bush, Asymptotic analysis of laminar flame propagation: review and extension, International Journal of Engineering Science, 17 (1979), pp. 597–613.
  * [7] P. Clavin and G. Searby, Combustion waves and fronts in flows: flames, shocks, detonations, ablation fronts and explosion of stars, Cambridge University Press, 2016.
  * [8] T. de Carvalho, D. J. Tonon, et al., Generic bifurcations of planar Filippov systems via geometric singular perturbations, Bulletin of the Belgian Mathematical Society-Simon Stevin, 18 (2011), pp. 861–881.
  * [9] P. De Maesschalck and S. Schecter, The entry–exit function and geometric singular perturbation theory, Journal of Differential Equations, 260 (2016), pp. 6697–6715.
  * [10] F. Dumortier, J. Llibre, and J. C. Artés, Qualitative theory of planar differential systems, Springer, 2006.
  * [11] F. Dumortier, N. Popović, and T. J. Kaper, A geometric approach to bistable front propagation in scalar reaction–diffusion equations with cut-off, Physica D: Nonlinear Phenomena, 239 (2010), pp. 1984–1999.
  * [12] B. Erickson, B. Birnir, and D. Lavallée, A model for aperiodicity in earthquakes, Nonlinear processes in geophysics, 15 (2008), pp. 1–12.
  * [13] Y. Estrin and L. P. Kubin, Criterion for thermomechanical instability of low temperature plastic deformation, Scripta Metallurgica, 14 (1980), pp. 1359–1364.
  * [14] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, Journal of Differential Equations, 31 (1979), pp. 53–98.
  * [15] L. Glass and R. Edwards, Hybrid models of genetic networks: Mathematical challenges and biological relevance, Journal of theoretical biology, 458 (2018), pp. 111–118.
  * [16] D. Hester, The nonlinear theory of a class of transistor oscillators, IEEE Transactions on Circuit Theory, 15 (1968), pp. 111–117.
  * [17] L. Ironi, L. Panzeri, E. Plahte, and V. Simoncini, Dynamics of actively regulated gene networks, Physica D: Nonlinear Phenomena, 240 (2011), pp. 779–794.
  * [18] S. Jelbart, K. U. Kristiansen, P. Szmolyan, and M. Wechselberger, Singularly perturbed oscillators with exponential nonlinearities, Journal of Dynamics and Differential Equations, 34 (2022), pp. 1823–1875.
  * [19] S. Jelbart, K. U. Kristiansen, and M. Wechselberger, Singularly perturbed boundary-equilibrium bifurcations, Nonlinearity, 34 (2021), p. 7371.
  * [20] , Singularly perturbed boundary-focus bifurcations, Journal of Differential Equations, 296 (2021), pp. 412–492.
  * [21] C. K. Jones, Geometric singular perturbation theory, in Dynamical systems, vol. 1609 of Lecture Notes in Mathematics, Springer, 1995, pp. 44–118.
  * [22] K. U. Kristiansen, Blowup for flat slow manifolds, Nonlinearity, 30 (2017), p. 2138.
  * [23] , A new type of relaxation oscillation in a model with rate-and-state friction, arXiv preprint, (2019).
  * [24] , The regularized visible fold revisited, Journal of Nonlinear Science, (2020).
  * [25] K. U. Kristiansen and S. J. Hogan, On the use of blowup to study regularizations of singularities of piecewise smooth dynamical systems in $\mathbb{R}^{3}$, SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 382–422.
  * [26] , Regularizations of two-fold bifurcations in planar piecewise smooth systems using blowup, SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 1731–1786.
  * [27] , Resolution of the piecewise smooth visible–invisible two-fold singularity in $\mathbb{R}^{3}$ using regularization and blowup, Journal of Nonlinear Science, 29 (2019), pp. 723–787.
  * [28] K. U. Kristiansen and A. Sarantaris, On a tropicalization of planar polynomial ODEs with finitely many structurally stable phase portraits, arXiv preprint arXiv:2305.18002, (2023).
  * [29] K. U. Kristiansen and P. Szmolyan, Relaxation oscillations in substrate-depletion oscillators close to the nonsmooth limit, Nonlinearity, 34 (2021), p. 1030.
  * [30] P. Le Corbeiller, Two-stroke oscillators, IRE Transactions on Circuit Theory, 7 (1960), pp. 387–398.
  * [31] J. Llibre, P. R. Silva, and M. A. Teixeira, Regularization of discontinuous vector fields on $\mathbb{R}^{3}$ via singular perturbation, Journal of Dynamics and Differential Equations, 2 (2007), pp. 309–331.
  * [32] A. Machina, R. Edwards, and P. van den Driessche, Singular dynamics in gene network models, SIAM Journal on Applied Dynamical Systems, 12 (2013), pp. 95–125.
  * [33] Z. Miao and N. Popović, The effect of a cut-off on a model of invasion with dispersive variability, To appear in Proceedings of the Workshop on Topics in Multiple Time Scale Dynamics, Banff 2022, (2024).
  * [34] E. Plahte and S. Kjøglum, Analysis and generic properties of gene regulatory networks with graded response functions, Physica D: Nonlinear Phenomena, 201 (2005), pp. 150–176.
  * [35] N. Popović, Front speeds, cut-offs, and desingularization: A brief case, in Fluids and Waves: Recent Trends in Applied Analysis: Research Conference, May 11-13, 2006, the Universtiy of Memphis, Memphis, TN, vol. 440, American Mathematical Soc., 2007, p. 187.
  * [36] , A geometric analysis of front propagation in a family of degenerate reaction-diffusion equations with cutoff, Zeitschrift für angewandte Mathematik und Physik, 62 (2011), pp. 405–437.
  * [37] N. Popović, A geometric classification of traveling front propagation in the nagumo equation with cut-off, in Journal of Physics: Conference Series, vol. 268, IOP Publishing, 2011, p. 012023.
  * [38] N. Popović, A geometric analysis of front propagation in an integrable nagumo equation with a linear cut-off, Physica D: Nonlinear Phenomena, 241 (2012), pp. 1976–1984.
  * [39] S. Portisch, A novel approach to dimension reduction in enzyme kinetics, Master’s thesis, TU Wien, 2017.
  * [40] M. G. Velarde, Nonequilibrium cooperative phenomena in physics and related fields, vol. 116, Springer Science & Business Media, 2012.
  * [41] M. Wechselberger, Geometric singular perturbation theory beyond the standard form, in Frontiers in Applied Dynamical Systems: Reviews and Tutorials, Springer International Publishing, 2020.
  * [42] F. A. Williams, Combustion theory, CRC Press, 1985.
  * [43] Y. B. Zeldovich and D. Frank-Kamenetskii, Theory of thermal flame propagation (in russian), J Phys Chem, 12 (1938), pp. 100–105.