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

HAL
open science

# Fast traveling waves for an epidemic model

Emmanuel Dormy, Andrew M Soward

► To cite this version:

Emmanuel Dormy, Andrew M Soward. Fast traveling waves for an epidemic model. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2025, 481 (2306), 10.1098/rspa.2024.0562. hal-05402189

HAL Id: hal-05402189

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

Submitted on 6 Dec 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
PROCEEDINGS A

royalsocietypublishing.org/journal/rspa

Research

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

Check for
updates

Article submitted to journal

**Subject Areas:**

mathematical modelling

**Keywords:**

epidemic model, asymptotic

expansion, traveling wave solutions

**Author for correspondence:**

Emmanuel Dormy

e-mail: emmanuel.dormy@ens.fr

Andrew Soward

**e-mail:**

andrew.soward@newcastle.ac.uk

# Fast traveling waves for an epidemic model

E. Dormy¹ and A. M. Soward²

¹Département de Mathématiques et Applications,
UMR-8553, École Normale Supérieure, CNRS, PSL
University, 75005 Paris, FRANCE

² School of Mathematics and Statistics, Newcastle
University, Newcastle upon Tyne NE1 7RU, UK

Asymptotic solutions are investigated for the traveling wave consisting of infectives $I(x - ct)$ propagating at speed $c$ into a region of uninfected susceptibles $S = S_+$, on the basis that $S_+$ is large. In the moving frame, three domains are identified. In the narrow leading frontal region, the infectives terminate relative abruptly. Conditions ahead (increasing $x$) of the front control the speed $c$ of the front advance. In the trailing region (decreasing $x$), the number of infectives decay relatively slowly. Our asymptotic development focuses on the dependence of $I$ on $S$ in the central region. There, the apparently simple problem is complicated by the presence of both algebraic and logarithmic dependencies. Still, we can construct an asymptotic expansion to a high order of accuracy that embeds the trailing region solution. A proper solution in the frontal region is numerical, but here the central region solution works well too.

We also investigated numerically the evolution from an initial state to a traveling wave. Following the decay of transients, the speed adopted by the wave is fast, though the slowest of those admissible. The asymptotic solutions are compared with the numerical solutions and display excellent agreement.
# 1. Introduction

We investigate a model based on an extension of the classical Kermack and McKendrick [1], or SIR, model for epidemics. In this so-called 'compartment model', a population N is divided into three compartments: those susceptible of catching a disease, those currently infected and those removed (i.e. either recovered or dead, but not susceptible to catching the disease). The time evolution of infectives I(t) and of susceptibles S(t) is then governed by the two coupled ODEs

$$ \frac{dI}{dt} = \beta IS - \alpha I, \qquad \frac{dS}{dt} = -\beta IS, \qquad (1.1a,b) $$

where $\beta$ controls the infections through interactions between susceptibles and infectives and $\alpha$ the recovery rate (or, in some cases, the death rate). The third compartment $R(t)$ is readily recovered as $N - S - I$, where from (1.1a, b)

$$ \frac{d}{dt}(N - R) = \frac{d}{dt}(I + S) = -\alpha I = \frac{\alpha}{\beta} \frac{d}{dt} \ln S, \qquad (1.1c) $$

which on integration shows that

$$ I + S - (\alpha/\beta) \ln S = \text{const.}, \qquad (1.1d) $$

where the constant is fixed by the initial values $S(0)$ and $I(0)$. Regarding $I$ as a function of $S$, the time dependence follows on integration of (1.1b). The integrable system (1.1) was initially validated against data from the plague death in the island of Bombay in 1905-1906. Since then, it has proved very useful in higher dimensional extensions that concern periods of epidemics and pandemics (including the recent Covid-19 pandemic, e.g. [2,3]). A recent investigation involving infectious periods is provided by Fowler and Hollingsworth [4] and see references therein.

## (a) Spatial dependence

We consider an extension of the above model, which allows for a spatial dimension $x$, introduced by Källén [5] to study the spatial propagation of an epidemic. On denoting the line density of infectives by $u(x, t)$ and susceptibles by $v(x, t)$, their evolution satisfies the coupled equations

$$ \frac{\partial u}{\partial t} = d \frac{\partial^2 u}{\partial x^2} + \beta uv - \alpha u, \qquad \frac{\partial v}{\partial t} = -\beta uv, \qquad (1.2a,b) $$

where the parameter $d$ measures spatial diffusion of infectives and $\alpha, \beta, d$ are all constants. Here, we have adopted the notation of Hosono & Ilyas [6] eq. (1.1), but note that a similar reduction of equations in the style (1.1) to Källén [5] eq. (1.1) is given by Källén et al [7] eqs.(3), (4). We solve (1.2) as an initial value problem for $t > 0$ with initial data $u(x, 0) = u_0(x)$ (a localised density, or patch, of infectives, i.e. $u_0(x) \to 0$, as $|x| \to \infty$) and $v(x, 0) = v_+^-$ (a constant, independent of $x$).

With the line densities $u$ and $v$ both having the dimensions of 1/length, the nonlinear balance of $\beta uv$ and $\alpha u$ suggests new dimensionless variables $f(> 0)$, $g(> 0)$ and $\phi$, defined by

$$ [u, v] = (\alpha/\beta) [f, g], \qquad g = e^{\phi}, \qquad \text{equivalently} \qquad \phi = \ln g. \qquad (1.3a-c) $$

We note that that $\beta/\alpha$ has the dimension. of length, while two further useful lengths are

$$ \ell = \sqrt{d/\alpha} \qquad \text{and} \qquad \mathcal{L} = \sqrt{d/(\beta v_+)}, \qquad (1.4a,b) $$

suggested by the balance, in (1.2a), of $d\partial^2 u/\partial x^2$ with $\alpha u$ and $-\beta v_+ u$ respectively. With the change of variables (1.3) and the introduction of $\ell$ (1.4a), the governing equations (1.2) reduce
to

$$ \frac{1}{\alpha} \frac{\partial f}{\partial t} = \ell^2 \frac{\partial^2 f}{\partial x^2} + f(e^{\phi} - 1), \qquad \frac{1}{\alpha} \frac{\partial \phi}{\partial t} = -f. \qquad (1.5a,b) $$

Together (1.5a,b) imply

$$ \frac{1}{\alpha} \frac{\partial}{\partial t} (f + e^{\phi} - \phi) = \ell^2 \frac{\partial^2 f}{\partial x^2}, \qquad (1.5c) $$

a construction suggested by (1.1c).

For the special x-independent case $d=0$ ($\ell=0$) considered by Kermack and McKendrick [1], (1.5c) has the integral $f = F(\phi)$, where

$$ F(\phi) = \mathcal{D} + \phi - e^{\phi}, \qquad (1.6a) $$

essentially (1.1d), in which $\mathcal{D}$ is the constant of integration, fixed by the initial conditions, $f = f_0$, $\phi = \phi_0 = \ln g_0$ (say) at $t = 0$. The actual time dependence is given implicitly by integrating (1.5b)

$$ \alpha t = - \int_{\phi_0}^{\phi} \frac{d\phi}{F}. \qquad (1.6b) $$

For our more general x-dependent case, we adopt the initial conditions

$$ (\beta/\alpha)u_0(x) = f(x, 0) = f_0(x), \qquad (\beta/\alpha)v_+ = g(x, 0) = g_+ = \text{const.} \qquad (1.7a,b) $$

with $f_0(x) \to 0$, as $|x| \to \infty$. Källén [5] investigated traveling waves, moving at speed $c$. These consist of a localised patch of infectives $u$ propagating into a region of uninfected susceptibles $v_+$. He proved the existence of a lower bound

$$ \bar{c} = \begin{cases} 2\sqrt{d\beta v_+(1-\gamma)}, & \text{where} \\ 2\alpha\ell\sqrt{g_+-1}, & \end{cases} \qquad \frac{1}{\gamma} = g_+ = \frac{\beta v_+}{\alpha} = \left(\frac{\ell}{\mathcal{L}}\right)^2, \qquad (1.8a,b) $$

on the speed $c(>\bar{c})$ of the patch, albeit expressed here in the notation of [6], which on demanding that $\bar{c}$ be real only exists when $\gamma < 1$. Note that $\gamma^{-1} = R_0$ is known as the 'basic reproduction number' or 'basic reproduction rate' (e.g. [8] p.322), which can be easily estimated at the beginning of a pandemic [9]. Hosono & Ilyas [6] showed that introducing a diffusive term into (1.2b) does not modify this minimal wave speed (i.e. the diffusion of the infectives $u$ alone controls spatial propagation of the epidemic). We emphasise that both [5,6] identify $\gamma = \alpha/(\beta v_+)$ (1.8b) (Källén uses $r$, see [5] Theorem B, p. 582) as the key parameter, which determines when traveling waves of permanent form (of concern to us) are possible, namely

$$ \gamma < 1 \qquad \text{equivalently} \qquad \ell > \mathcal{L}, \quad g_+ > 1, \quad \phi_+ > 0, \qquad (1.9) $$

to avoid inadmissible complex $\bar{c} \propto \sqrt{g_+-1}$. The approach in [6] was extended to more general systems by Ai & Huang [10].

Of course, there are many studies of traveling wave systems similar to ours. Background material may be found in Murray's books [8,11]. We mention the investigations [12-15] of the Fisher-Kolmogorov-Petrowskii-Piscounov (FKPP) model [16,17], as well as the detailed studies of [18-20].

## (b) Organisation

Our article is organised as follows.

In §2, we formulate the traveling wave (speed $c$) problem and introduce the dimensionless parameter, $\epsilon = (\alpha\ell/c)^2$ (2.1b), upon which we again rescale. In §2(a), we summarise known features of the solutions and, in §2(b), apply simple phase plane arguments, similar in spirit to §2 of Källén [5], to establish the existence of the solution as well as its nature. In §2(c), we emphasise how the conditions at the leading edge of the wave control the realised wave speed, by appealing to the Dee & Langer [21] (see also van Saarloos [22]) condition, which concerns the propagation of finite amplitude disturbances into undisturbed unstable regions (see, e.g., [23]).
In §§3 and 4, we consider a large number of susceptibles $v_+ \gg \alpha/\beta$, equivalently $\ell \gg \mathcal{L}$ by (1.8b). For that, the lower bound $\bar{c}$ on $c$, measured in units of $\alpha\ell$ by $\bar{c}^{-1/2} = \bar{c}/\alpha\ell \gg 1$ (see (2.25a), (3.1d,e)), is fast (see (1.8a)). So, as the solutions of initial value problems in §5 confirm, the emergent traveling wave is generally the slowest $c = \bar{c}$ of the admissible fast waves speeds $c (\ge \bar{c})$. Rather than limit attention to $\bar{c}$, we consider neighbouring values of $c$ and adopt $\epsilon = (\alpha\ell/c)^2 = O(\bar{c})$ (2.1b) as our small expansion parameter. Our asymptotic problem based on $\epsilon \ll 1$ is complicated because of both algebraic and logarithmic dependence on $v$, hinted at by (1.1c,d).

In §3, we find that the propagating patch of infectives, $u(x - ct)$, exhibits three asymptotic domains. In the tail, the solution, described in §3(a), is simple; at any particular point $x$, spatial diffusion is negligible and the remaining balance $\partial u/\partial t = -\alpha u$ leads to exponential decay of infectives on the time scale $T_- = \alpha^{-1} (3.8e)$. The tail is embedded in the central (or intermediate) region solution considered in §3(b). No analytic solution appears to exist in the frontal region of width $O(\mathcal{L})$ (see (1.4b)) beyond which the susceptibles remain undisturbed, $v = v_+$. Nevertheless, we attempt a power series expansion in §3(c), also embedded in the central region solution.

In §4, various approximate solutions are described. Encouraged by our §3 findings that the tail and front are embedded in our intermediate solution of §3(b), we extend that analysis to higher orders in §4(a). We find that the ensuing series solution works well almost everywhere except close to where the susceptibles are undisturbed ahead of the wave. The implications of this weakness are discussed briefly in §4(b). To rectify the defect we consider the heuristic solutions (4.16a,b) in §4(c), which are surprisingly good everywhere.

In §5, we solve an initial value problem numerically. The results confirm our §4(a) analytic predictions, as highlighted by our comparisons in various figures. Furthermore, the numerical results emphasise that the intermediate and frontal domains occupy the same spatio-temporal range, and in that sense there are only two domains: the tail and the front. In §6, we conclude with a few brief remarks.

## 2. The traveling wave problem

We consider traveling wave solutions dependent on $x - ct$ alone, that move with speed

$$c(> 0), \qquad (2.1a)$$

which, w.l.o.g., is in the positive $x$-direction. We will focus attention on fast waves, $c \gg \alpha\ell$, for which we introduce the small dimensionless parameter

$$\epsilon = \epsilon(c) = (\alpha\ell/c)^2 = \alpha d/c^2 \ll 1 \qquad (2.1b)$$

and use it to rescale (1.3a), (1.6a):

$$(f, g) = \epsilon^{-1}(f, g), \qquad F = \epsilon^{-1}F. \qquad (2.1c,d)$$

It is unfortunate that our non-dimensionalisation involves $c$, which is part of the solution. Nevertheless, our numerical solutions of initial value problems reported in §5 reveal that the traveling wave that emerges, as $t \to \infty$, has the slowest permissible speed $c = \bar{c}$ (1.8a), defined by the problem's data. For our analytic development it is not necessary to be that precise, rather we simply assume that $\epsilon = O(\bar{c})$, where, here and henceforth, the overline is used to identify values at $c = \bar{c}$, e.g., $\bar{c} = \epsilon(\bar{c})$. On use of (2.1b,c), the relations (1.3a,b) become

$$[u, v] = (\alpha/\beta)\epsilon^{-1}(f, g) = (c^2/\beta d)[f(\xi), g(\xi)], \qquad g = \epsilon e^{\phi(\xi)} \qquad (2.2a,b)$$

with $\xi$ defined by

$$x - ct = \ell\epsilon^{1/2}\xi = (d/c)\xi. \qquad (2.2c)$$

Likewise the useful formula (1.6a) may be written in the style of (2.1d) as

$$F = \epsilon F(\phi) = \epsilon(D + \phi - e^{\phi}), \qquad \text{in which} \quad D = \epsilon D. \qquad (2.3a,b)$$
Differentiation of (2.3a) gives

$$F' = \epsilon F' = \epsilon(1 - e^{\phi}) \ge 0, \quad \text{when} \quad \phi \le 0, \qquad (2.4a,b)$$

where here and henceforth the ' denotes differentiation with respect to $\phi$. Furthermore, since $F'' = -\epsilon e^{\phi} < 0$ everywhere, F achieves its maximum

$$F(0) = D - \epsilon = \epsilon(D - 1) \quad \text{at} \quad \phi = 0, \quad \text{where} \quad g = \epsilon \qquad (2.5a-c)$$

and $F' = F'(0) = 0$.

In terms of our new variables (2.1) and (2.2), equations (1.5a,b) (equivalently (1.2a,b)) become

$$\frac{d^2 f}{d\xi^2} + \frac{df}{d\xi} = f F' = \frac{dF}{d\xi}, \quad \frac{d\phi}{d\xi} = f. \qquad (2.6a,b)$$

Integration of (2.6a) yields

$$\frac{df}{d\xi} + f = F \qquad \text{equivalently} \qquad f f' + f = F, \qquad (2.7a,b)$$

where $D = \epsilon D$ in the definition (2.3) of F is determined by the boundary conditions (2.8a) below.

As $\xi \to \pm\infty$, we demand that $f \to 0$, i.e., $f_{\pm} = 0$, where here and henceforth we denote such terminal values by the subscripts $\pm$. Furthermore, by (2.7b), the vanishing of $f$ is mimicked by $F$, i.e., $F_{\pm} = 0$. In turn, use of (2.3a) determines

$$e^{\phi_{-}} - \phi_{-} = D = e^{\phi_{+}} - \phi_{+}. \qquad (2.8a)$$

This has distinct real solutions $\phi_{\pm}$ with the property,

$$\phi_{-} < 0 < \phi_{+}, \quad \text{provided} \quad D > 1. \qquad (2.8b,c)$$

Furthermore, we may integrate (2.6b) to obtain $\xi(\phi)$ in the form

$$\xi = \int_{0}^{\phi} \frac{d\phi}{f(\phi)} = \int_{\epsilon}^{g} \frac{dg}{gf(\phi(g))} \qquad (\xi' > 0), \qquad (2.9)$$

where, w.l.o.g., we have chosen the $\xi$-origin $\xi = 0$ to be located, where $\phi = 0$ (equivalently $g = \epsilon$) inside the range (2.8b). It is important to appreciate that during the transient evolution of the wave from the initial state (1.7), the boundary condition (1.7b) on $g$, as $x \to \infty$, remains unaltered, keeping the value $g_{+} = \epsilon^{-1}g_{-}$. In contrast, as $\xi \to -\infty$, the initial value is forgotten and the new value $g_{-} = \epsilon^{-1}g_{+}$, determined by (2.8a), is attained behind the wave but ahead of any transients.

We conclude with two integral constraints. For the total number

$$U = \int_{-\infty}^{\infty} u \, dx = \frac{c}{\beta} U, \qquad (2.10a)$$

infected at any instant, we may use (2.6b) and (2.8a) to determine

$$U = \int_{-\infty}^{\infty} f \, d\xi = \phi_{+} - \phi_{-} \equiv \llbracket \phi \rrbracket = \llbracket g \rrbracket, \qquad (2.10b)$$

where here and henceforth, $\llbracket \cdots \rrbracket$ denotes the jump between the terminal values, as $\xi \to \pm\infty$. In a similar spirit, integration of (2.7b) with respect to $\phi$ yields the additional identity

$$\begin{aligned} U_1 &= \int_{-\infty}^{\infty} f^2 d\xi = \int_{\phi_{-}}^{\phi_{+}} F d\phi \\ &= \epsilon(D\llbracket\phi\rrbracket + \frac{1}{2}\llbracket\phi^2\rrbracket - \llbracket g \rrbracket) = \frac{1}{2}\epsilon U(g_{+} + g_{-} - 2) = -\frac{1}{2}U(F_{+}' + F_{-}'). \end{aligned} \qquad (2.11)$$

Generally, we use $\phi$, rather than $\xi$, as our independent variable, envisaging $\xi = \xi(\phi)$ as in (2.9).
(a) The solution structure

The function $f(\phi)$ is defined on the range $\phi_- < \phi < \phi_+$ (2.8b) and is maximum at the turning point $\phi = \phi_{\max}$, where $f' = 0$; there, values will, henceforth, be identified by the subscript $\max$. We differentiate (2.7b) to obtain

$$f f'' + (f')^2 + f' = F' \equiv \epsilon(1 - e^{\phi}). \qquad (2.12a)$$

Regarded as a quadratic equation for $f'$, (2.12a) has the solution

$$f' = \frac{1}{2} \left[ -1 + \sqrt{1 + 4(F' - f f'')} \right], \qquad (2.12b)$$

where the sign of the root taken is fixed by the requirement that $f' > 0$ at the end-point $\phi = \phi_-$, at which $f = 0$. The reality of $f'$, defined by (2.12b), implies

$$f f'' \le \frac{1}{4} + F', \quad \text{also} \quad f_{\max} f'_{\max} = F'_{\max}. \qquad (2.13a,b)$$

The condition $f''_{\max} < 0$ for a maximum (recall $f > 0$) requires $F'_{\max} < 0$, which in turn implies $\phi_{\max} > 0$ ($\xi_{\max} > 0$). Significantly, elsewhere, we have

$$f \ge F \quad \text{with} \quad f' \le 0 \quad \text{on} \quad \phi \ge \phi_{\max} \quad (\xi \ge \xi_{\max}). \qquad (2.13c)$$

Obviously $f$ is bounded above by $f_{\max} < F(0) = D$. The above features are evident in figures 1–4 below.

Close to the endpoints $\phi_{\pm}$, where $f_{\pm} = 0$, the system (2.6) becomes linear with solutions $f \propto \exp(k_{\pm}\xi)$. Accordingly, by (2.12b), $f'_{\pm} = (f^{-1}df/d\xi)_{\pm} = k_{\pm}$ take the values

$$k_{\pm} = \frac{1}{2} \left[ -1 + \sqrt{1 + 4F'_{\pm}} \right], \qquad (2.14a)$$

or, in its more primitive form derived from (2.12a) again with $f = f_{\pm} = 0$,

$$k_{\pm}^2 + k_{\pm} = F'_{\pm} = \epsilon(1 - e^{\phi_{\pm}}) \quad \Rightarrow \quad k_{\pm} < F'_{\pm}, \qquad (2.14b,c)$$

since $k_{\pm}^2 \ge 0$, consistent with the requirements (2.13c). By (2.2c), $k_{\pm}$ has the dimensional form

$$k_{\pm} = (c/d)k_{\pm}, \quad \text{equivalently} \quad \ell k_{\pm} = \epsilon^{-1/2}k_{\pm}. \qquad (2.14d)$$

Noting that $\pm k_{\pm}$ denote the dimensional decay rates as $\xi \to \pm\infty$, we may define the e-folding lengths/times at fixed $t/x$ by

$$L_{\pm} = \mp 1/k_{\pm} = \mp \epsilon^{1/2} \ell/k_{\pm}, \quad T_{\pm} = L_{\pm}/c, \qquad (2.14e,f)$$

respectively, with $L_-$ for the tail, and $L_+$ for the front.

(b) Phase plane solution

Writing (2.7b), in the form

$$\frac{df}{d\phi} = \mathcal{R}(\phi, f) \equiv \frac{F - f}{f} \qquad (F = F(\phi)) \qquad (2.15)$$

its solution is readily understood in the $\phi-f$ phase plane. The key boundary condition is the location $(\phi_-, 0)$ of the left-hand end point, determined by (2.8a). On increasing $\phi$ from $\phi_-$, the required solution emerges tangent to the line $f = k_-(\phi - \phi_-)$. It continues to increase remaining below the F-curve, achieving its maximum on intersecting the F-curve at $\phi = \phi_{\max} > 0$. On further increase of $\phi$, the solution $f$ now remains above the F-curve. To establish that the solution reaches the required end point $\phi_+$, we argue in the following paragraph that, in addition to being bounded below by the F-curve, the solution is bounded above by the line $f = k_+(\phi - \phi_+)$. Trapped in that wedge-like region the only point at which the solution can terminate is $(\phi_+, 0)$. The argument is a variant to that given by Källén in the proof of their Theorem B on p.853 in §2 of [5], in which their $h(g)$ after rescaling is our $F(\ln g)$.
The essential idea is to show that, if the solution escaped from the wedge across the straight line $f = k_+(\phi - \phi_+)$, the solution slope $|f'|$ must be smaller in magnitude than $|k_+|$. We now show that cannot happen by considering

$$ \mathcal{R}(\phi, k_+(\phi - \phi_+)) - k_+ = \left[ \frac{F}{k_+(\phi - \phi_+)} - 1 \right] - k_+ = \frac{F - F'_+(\phi - \phi_+)}{k_+(\phi - \phi_+)} \quad (2.16) $$

since $1 + k_+ = F'_+/k_+$ (see (2.14b)). Further, since $F'' = -\epsilon e^{\phi} < 0$, F is a concave function

$$ (F - F_+) - F'_+(\phi - \phi_+) < 0. \quad (2.17a) $$

So as $F_+ = 0$, (2.16) and (2.17a) establish

$$ \mathcal{R}(\phi, k_+(\phi - \phi_+)) < k_+ \quad \text{on} \quad \phi < \phi_+. \quad (2.17b) $$

This shows that on increasing $\phi$ the solution f-curve could cross the line $f = k_+(\phi - \phi_+)$ into (but not out of) the wedge $k_+(\phi - \phi_+) > f(\phi) > F(\phi)$ ($\phi_{max} < \phi < \phi_+$).

### (c) The front $\phi = \phi_+$

At $\phi = \phi_+$, in addition to $k_+$ (see (2.14a)), there is a second solution $k_+^\dagger = \frac{1}{2}[-1 - \sqrt{1 + 4F'_+}]$. Indeed, in the $\xi$-f plane, there are two linear solutions $f \propto \exp(k_+\xi)$ and $f \propto \exp(k_+^\dagger\xi)$ as $\xi \to \infty$. However, as $-k_+^\dagger > -k_+$, the $k_+^\dagger$ solution decays faster than the $k_+$ solution leaving only the $k_+$ solution. In the $\phi$-f plane, the repeated root

$$ k_+^\dagger = k_+ = \bar{k}_+ \equiv -\frac{1}{2} \quad (2.18a) $$

simply identifies the degenerate member of the family of solutions with tangents $f = k_+(\phi - \phi_+)$ at $\phi = \phi_+$. According to (2.14a), this degeneracy occurs when

$$ 4F'_+ = 4\bar{F}'_+ \equiv -1, \quad (2.18b) $$

$$ \text{equivalently} \quad \epsilon = \bar{\epsilon}, \quad \text{where} \quad \bar{\epsilon} = \frac{1}{4}(e^{\phi_+} - 1)^{-1} \quad (2.18c,d) $$

on use of (2.4a). The two real solutions $k_+$ and $k_+^\dagger$ are only possible when

$$ 1 > -4F'_+ = F'_+/\bar{F}'_+ = \epsilon/\bar{\epsilon} = (\bar{c}/c)^2 \quad \Rightarrow \quad \begin{cases} \epsilon < \bar{\epsilon}, \\ c > \bar{c}, \end{cases} \quad (2.19a,b) $$

where we have noted that $\epsilon c^2 = (\alpha \ell)^2$ (see (2.1b)). The result shows that the limiting case $c = \bar{c}$ identifies the slowest possible traveling wave. Furthermore, (2.3b), (2.8a) and (2.18d) imply

$$ \bar{D} = \bar{\epsilon}D = \frac{e^{\phi_+} - \phi_+}{4(e^{\phi_+} - 1)} = \frac{1}{4} - \frac{\ln g_+ - 1}{4(g_+ - 1)} \quad (2.20a) $$

$$ \text{with} \quad \bar{D} \uparrow \frac{1}{4} \quad \text{as} \quad g_+ \to \infty. \quad (2.20b) $$

To set the above results in context, we consider small amplitude disturbances $f \propto \exp(kx - pt)$ ahead of the front, for which we can linearise (1.5a) on the basis that $1 - e^{\phi} \approx 1 - e^{\phi_+} = F'_+$ (see (2.4a)). Then, $p = p(k)$ satisfies the dispersion relation

$$ p = \alpha(-\ell^2 k^2 + F'_+), \quad (2.21a) $$

where, in conventional wave notation, $k/i$ and $p/i$ are respectively the complex wave number and complex frequency. From (2.21a), we deduce that the phase and group velocities are

$$ c_p = \frac{p/i}{k/i} = \frac{p}{k}, \quad c_g = \frac{\partial(p/i)}{\partial(k/i)} = \frac{\partial p}{\partial k} = -2\alpha\ell^2 k. \quad (2.21b,c) $$

At the head of the front, where $k = k_+ = \epsilon^{-1/2}\ell^{-1}k_+$, $F' = F'_+ = \epsilon^{-1}F'_+$, (2.21b,c) determine

$$ c_p = c(-k_+ + F'_+/k_+), \quad c_g = -2ck_+, \quad (2.22a,b) $$
$$c_g - c_p = c(-k_+ - F'_+/k_+). \qquad (2.22c)$$

Moreover, the result $c^2F'_+ = -\frac{1}{4}\bar{c}^2$ (see (2.19a)) may be substituted into (2.14a) to obtain

$$\frac{1}{2}c_g = -ck_+ = \frac{1}{2}[c - \sqrt{c^2 - \bar{c}^2}], \quad cF'_+/k_+ = \frac{1}{2}[c + \sqrt{c^2 - \bar{c}^2}], \qquad (2.23a,b)$$

where we have noted that $c_g = -2ck_+$ (2.22b). Substitution of these results into the right-hand sides of (2.22a,c) determines

$$c_p = c, \qquad c_g - c_p = -\sqrt{c^2 - \bar{c}^2} \qquad (2.24a,b)$$

$$\text{implying that} \quad 0 < c_g < c, \quad \text{when} \quad c > \bar{c}. \qquad (2.24c,d)$$

It follows the slowest wave speed $c = \bar{c}$, meets the Dee and Langer [21] criterion $c_g = c_p = \bar{c}$ for identifying realised front speed, while (2.23a) evaluated at $c = \bar{c}$ recovers the repeated root values $k_+ = k_+^1 = -\frac{1}{2}$ (2.18a). For other front speeds, the condition $c_g < c$ (see (2.24c,d)) indicates that ‘wave energy’ cannot be transported ahead of the wave to maintain its form, and the front speed relaxes to the slowest available speed $\bar{c}$. Indeed the numerical solution of initial value problems confirm that this Dee and Langer limit $c = \bar{c}$ is normally reached after all transients have decayed.

To summarise, (2.1b) and (2.18d) determine

$$\bar{c}^{1/2} \bar{c} = \alpha \ell, \qquad \bar{c} = \frac{1}{4}(g_+ - 1)^{-1}, \qquad (2.25a,b)$$

which together recover Källén's result (1.8a) for the slowest wave speed $\bar{c}$, specifically Theorem B of [5], p. 852, where $c_0$ corresponds to our $\bar{c}$, while Hosono & Ilyas in their Main Theorem [6], p. 936, used $c^*$. In addition, using (2.18a), the e-folding length and time (2.14e,f) at the front are

$$\bar{L}_+ = -\frac{1}{k_+} = -\frac{\bar{c}^{1/2} \ell}{k_+} = 2\bar{c}^{1/2} \ell = \frac{2d}{\bar{c}}, \qquad \bar{T}_+ = \frac{2d}{\bar{c}^2} = \frac{\bar{L}_+^2}{2d} \qquad (2.26a,b)$$

$$\text{with} \quad \bar{L}_+ \downarrow \mathcal{L} = \sqrt{d/(\beta v_+)} \quad \text{as} \quad g_+ \to \infty. \qquad (2.26c)$$

### 3. The large $g_+$ limit ($\ell \gg \mathcal{L}$). (I) Various regions

Henceforth we limit our attention to a large initial value $v_+$ of $v$ in the sense that

$$v_+ \gg \alpha/\beta \quad \text{implying} \quad \ell \gg \mathcal{L} \quad \text{and} \quad g_+ \gg 1 \qquad (3.1a-c)$$

by (1.8b). Then (2.25a,b) show, under the approximation $g_+ - 1 \approx g_+$, that

$$\bar{c} \ll 1 \quad \text{and} \quad \frac{\bar{c}}{2\alpha\ell} = \frac{1}{2\bar{c}^{1/2}} \approx \sqrt{g_+} = \frac{\ell}{\mathcal{L}} \gg 1, \qquad (3.1d,e)$$

which means that the slowest wave, speed $\bar{c}$, is fast when measured in units of $\alpha\ell$. Though those slowest waves are our primary interest, it is not necessary to be that precise in our asymptotic development. Rather, we will simply assume that $c$ is close $\bar{c}$ with the consequence that

$$\epsilon = O(\bar{c}) \ll 1 \quad \text{and} \quad \epsilon \mathcal{D} = \mathcal{D} = O(\bar{D}) = O(1), \qquad (3.2a,b)$$

remembering the equivalent requirements $\epsilon < \bar{c}$ and $c > \bar{c}$ (2.19b).

With $g_+ = e^{\phi_+}$ and $\phi_+ = \ln g_+$ prescribed by the solution to the related initial value problem, the relation (2.8a) determines

$$\mathcal{D}/\epsilon \equiv \mathcal{D} = e^{\phi_+} - \phi_+ \quad \Rightarrow \quad 0 < \phi_+ \approx \ln(\mathcal{D}/\epsilon) = O(\ln(\epsilon^{-1})). \qquad (3.3a,b)$$

Likewise, with $\mathcal{D}$ given by (3.3a), the value of $\phi_-$ is the second solution of (2.8a), namely

$$\mathcal{D}/\epsilon = \mathcal{D} = e^{\phi_{-}} - \phi_{-} \quad \Rightarrow \quad 0 > \phi_{-} \approx -\mathcal{D}/\epsilon = O(\epsilon^{-1}) \qquad (3.3c,d)$$

(see also the remarks following (2.9)). In view of the large magnitudes $|\phi_\pm|$ of the endpoint locations (3.3b,d), we seek asymptotic solutions for $f(\phi)$ on the three asymptotic regimes $-\phi = O(\epsilon^{-1})$, $\phi = O(1)$ and $\phi = O(\ln(\epsilon^{-1}))$, respectively, in §§3(a-c) below.
(a) The tail: $-\phi = O(\epsilon^{-1}) \gg 1$, $-\xi = O(\epsilon^{-1})$

On the range $-\phi = O(\epsilon^{-1}) \gg 1$, we neglect all exponentially small terms involving $e^{\phi}$, as in (3.3d), so that (2.3a), (2.4a) are approximated by

$$F = D + \epsilon\phi, \qquad F' = \epsilon. \tag{3.4a,b}$$

Hence, the endpoint, at which $F_- = 0$, is simply $\phi_- = -D/\epsilon$ (3.3d), and there (2.14a) determines

$$f'_{-} = k_{-} = \frac{1}{2}[-1 + \sqrt{1 + 4\epsilon}] \tag{3.5a}$$

$$= \epsilon - \epsilon^2 + 2\epsilon^3 + \dots \quad \text{for} \quad \epsilon \ll 1. \tag{3.5b}$$

Thus, $F'_{-} - f'_{-} \approx \epsilon^2$, with the consequence that $F'_{-} > f'_{-} > 0$ as claimed by (2.14c).

For the simplified $F = D + \epsilon\phi$ (3.4a), the equation $f(f' + 1) = F$ (2.7b) has the linear solution

$$f = f^{[0]} \equiv a + e\phi = k_{-}(\phi - \phi_{-}^{[0]}). \tag{3.6a}$$

It intersects the $\phi$-axis at

$$\phi_{-}^{[0]} = -a/e = -D \approx \phi_{-} \tag{3.6b}$$

with only an exponentially small error, and where

$$e(e+1) = \epsilon, \qquad e = k_{-} = O(\epsilon), \qquad a/D = e/\epsilon = k_{-}/\epsilon. \tag{3.6c-e}$$

Despite the small slope $k_- = \epsilon + O(\epsilon^2)$ of $f^{[0]}$ (3.6a), it grows linearly over the long $\phi = O(\epsilon^{-1})$ range, becoming $O(1)$ when $\phi = O(1)$.

The spatial dependence (2.9) based on $f^{[0]}$ is

$$a(\xi - \xi_{-}) = \int_{0}^{\phi} \frac{a \, d\phi}{a + e\phi} = \frac{a}{e} \ln \left( 1 + \frac{e}{a} \phi \right), \tag{3.7}$$

where the constant $\xi_- = O(\epsilon)$ is fixed by matching with the intermediate solution (3.12b,c) below. The inverse of (3.7) is

$$\phi = -(a/e)[1 - e^{e(\xi - \xi_{-})}] \quad \text{or} \quad \phi - \phi_{-}^{[0]} = (D/\epsilon)e^{k_{-}(\xi - \xi_{-})} \tag{3.8a,b}$$

on use of (3.6d,e) (or solve (2.7b) directly, noting $F_- = 0$ and using (3.6c,d)). Then (3.6a) becomes

$$f^{[0]} = ae^{e(\xi - \xi_{-})} = ae^{k_{-}(\xi - \xi_{-})}, \tag{3.8c}$$

in accord with the decay rate $k_-$ identified above (2.14a). According to (2.14e, f) with $k_- \approx \epsilon$ (see (3.6d,e)), the e-folding length and time at fixed $t$ and $x$ respectively are

$$L_{-} = \ell/\epsilon^{1/2} = c/\alpha, \qquad T_{-} = 1/\alpha, \tag{3.8d,e}$$

namely the free decay time at a fixed position when diffusion is neglected.

Finally $g = e^{\phi}$ becomes

$$g = g_{-} \exp[(D/\epsilon)e^{k_{-}(\xi - \xi_{-})}] = g_{-} \exp[e^{-1}f^{[0]}], \qquad g_{-} = e^{-D/\epsilon}, \tag{3.9a,b}$$

showing that the rise of $g$ is far more abrupt than that of $f^{[0]}$, as confirmed by the numerical results portrayed in figure 5b below. The disconcerting double exponential is natural, as $g$ collapses on the short frontal length scale $L_+ = O(\epsilon L_-)$, since $L_+/L_- = k_-/k_+$ (see (2.14e)) and noting that $k_+ = O(1)$, $k_- = O(\epsilon)$, i.e., $L_+ \ll L_-$.

(b) The central ‘intermediate’ region: $\phi = O(1)$, $\xi = O(1)$

When $-\phi = O(\epsilon^{-1})$, the term $g = e^{\phi}$ in (2.3a) is exponentially small and negligible. However, as $-\phi$ decreases, $g$ increases becoming $O(1)$, when $\phi = O(1)$, and must be retained. It is easy to show,
by direct substitution into (2.7b), that the solution correct to $O(\epsilon)$ is

$$f = f^{[1]} = f^{[0]} - \epsilon\kappa g = a + e\phi - \kappa g \quad \text{with} \quad \kappa = (a+1)^{-1}. \quad (3.10a,b)$$

A formal derivation is provided in §4(a) below. We compare, in passing, $F = D + \epsilon\phi - g$ (2.3a) with $f^{[1]}$, for which the structure is the same but the value of the coefficients differ (see (3.6c,e)) leading to the difference

$$f^{[1]} - F = -e(a + e\phi) + a\kappa g. \quad (3.11)$$

Evidently the straight line part $f^{[0]} = a + e\phi$ of $f^{[1]}$ lies very slightly below the corresponding $D + \epsilon\phi$ of $F$, but meeting at $\phi = \phi^{[0]}$. On increasing $\phi$, the exponential difference $a\kappa g$ causes the curves to cross when $\phi = O(1)$, after which $f^{[0]}$ lies above $F$ by an $O(1)$ amount. However, note that $f^{[1]}$ vanishes at $\phi^{[1]}$ (say) which differ from $\phi_{\pm}$ (see (4.11) below). Of course, these end-points lie outside the domain of validity of the 'intermediate' approximations.

The spatial dependence (2.9) based on $f^{[1]}$ (3.10a) yields

$$\xi = \int_{0}^{\phi} \frac{d\phi}{a + e\phi - \kappa g} = \int_{0}^{\phi} \left( 1 - \frac{\epsilon\kappa e^{\phi}}{a + e\phi} \right)^{-1} \frac{d\phi}{a + e\phi}. \quad (3.12a)$$

Taking the binomial expansion of the first inverse in the second representation and then integrating yields

$$7a\xi = \frac{a}{e} \ln \left( 1 + \frac{e\phi}{a} \right) + \frac{\epsilon\kappa}{a} (e^{\phi} - 1) + O(\epsilon^2) \quad \text{for } g = e^{\phi} = O(1). \quad (3.12b)$$

Its first ln-term recovers large negative $\phi$ solution (3.7), which on matching determines

$$\xi_{-} = -\kappa\epsilon/a^2. \quad (3.12c)$$

For our $\phi = O(1)$ range, (3.12b) reduces, on expanding the ln-term, to

$$a\xi = \phi - \frac{1}{2}(e/a)\phi^2 + (\kappa\epsilon/a)(e^{\phi} - 1) + O(\epsilon^2) \quad (\phi = O(1)). \quad (3.13a)$$

Its inverse

$$\phi = a\xi + \frac{1}{2}ea\xi^2 - (\kappa\epsilon/a)(e^{a\xi} - 1) \quad (\xi = O(1)) \quad (3.13b)$$

clearly identifies the transition in space from $\phi$ increasing algebraically, when $\xi < 0$, to a subtracted exponential correction of increasing magnitude, when $\xi > 0$.

As $\phi$ increases in size, but before entering the frontal domain $\phi = O(\ln \epsilon^{-1})$ (§3(c)), i.e.,

$$1 \ll e^{\phi} \ll \epsilon^{-1} \quad \text{equivalently} \quad \epsilon \ll g \ll 1, \quad (3.14)$$

we may neglect $e\phi$ in the denominator of (3.12a), which becomes approximately

$$a\xi = \int_{0}^{\phi} \frac{d\phi}{1 - (\epsilon\kappa/a)e^{\phi}} = \phi - \ln[1 - (\epsilon\kappa/a)e^{\phi}] + a\xi_{-} + O((\epsilon e^{\phi})^2). \quad (3.15a)$$

On the basis that $\epsilon e^{\phi} \ll 1$ (3.14), equation (3.15a) has the expansion

$$a\xi = \phi + (\epsilon\kappa/a)(e^{\phi} - 1) + O((\epsilon e^{\phi})^2), \quad (3.15b)$$

which matches correctly with (3.13a) up to an error $O(\epsilon e^{\phi^2})$, small compared to the terms $O(\epsilon e^{\phi})$ retained in (3.15).
On writing (3.15a) in the form $-a(\xi - \xi_{-}) = \ln[e^{-\phi} - (\epsilon\kappa/a)]$, the inverse is

$$ \frac{g}{\epsilon} = g = \frac{e^{a(\xi - \xi_{-})}}{1 + (\epsilon\kappa/a)e^{a(\xi - \xi_{-})}} + O((\epsilon e^{\phi})^2) \quad \text{for} \quad a\xi \ll O(\ln(\epsilon^{-1})), \quad (3.16a) $$

as required by $e^{\phi} \ll \epsilon^{-1}$ (3.14) (recall that $\phi \approx a\xi$ (3.15b)). The result determines

$$ f \approx a - \kappa g = \frac{a}{1 + (\epsilon\kappa/a)e^{a(\xi - \xi_{-})}} + O((\epsilon e^{\phi})^2). \quad (3.16b) $$

Though (3.16a,b) are only valid for $(\epsilon\kappa/a)e^{a\xi} \ll 1$, the formulae suggest exponential decay to the terminal values

$$ f \downarrow 0 \qquad g \uparrow a/\kappa \qquad \text{as} \qquad \xi \to \infty. \quad (3.17) $$

Of course, (3.17) takes the solution beyond its domain of validity identified in (3.16a). The true terminal value determined by $F \approx D - g = 0$ is $g_{+} \approx a$ (see (3.22a) below).

### (c) The front region: $\phi = O(\ln(\epsilon^{-1}))$, $\xi = O(1)$

The right-hand end point $\phi = \phi_{+} > 0$, at which $f = 0$, occurs when

$$ \phi_{+} = \ln(g_{+}/\epsilon) = \ln g_{+}, \quad (3.18a) $$

setting the scale

$$ \phi = O(\ln(\epsilon^{-1})) \gg 1, \qquad g = O(1) \quad (3.18b,c) $$

of the large $\phi$-range. We briefly consider solutions based on the assumptions

$$ e = \epsilon + O(\epsilon^2), \qquad a = D + O(\epsilon), \qquad \epsilon\phi = O(\epsilon \ln(\epsilon^{-1})) \ll 1 \quad (3.19a-c) $$

(see (3.6c,e) and (3.18a)). The smallness of $\epsilon\phi$, guaranteed by (3.19c), means that it may be neglected in $F = D + \epsilon\phi - g$ (2.3a) leaving $F = a - g$, correct to leading order. Then, under the assumption that $f$ is a function of $g$, equation (2.15) reduces to

$$ \frac{df}{dg} = \frac{f'}{g'} = \frac{a - g - f}{gf}. \quad (3.20) $$

The required solution of (3.20) in the neighborhood of the singular point $f = a, g = 0$ is

$$ f - a = -\kappa g, \qquad \text{where} \qquad \frac{df}{dg} = -\kappa \left( = -\frac{1}{a+1} \right). \quad (3.21a,b) $$

On increasing $g$, a numerical solution may be found, which necessarily terminates at the singular point $f = 0, g = a$, in the neighborhood of which (3.20) is solved by

$$ f = k(g - a)/a \qquad \text{where} \qquad k = a \frac{df}{dg} = \frac{1}{2}[-1 + \sqrt{1 - 4a}]. \quad (3.22a,b) $$

Valid solutions require real $k$, and so

$$ 0 < a \approx D < \frac{1}{4} \quad (3.22c) $$

c.f. the slowest wave speed condition $\bar{D} \approx \frac{1}{4}$ (2.20b). A power series solution

$$ f = a - \kappa g + \sum_{n=2}^{\infty} a_n (\kappa g)^n \quad (3.23a) $$

may be attempted, in which the leading term $a - \kappa g$ agrees with $f^{[1]} = a + e\phi - \kappa g$ (3.10a) on omitting $\epsilon\phi$. That omission leads to the approximate $\xi$-solution discussed between (3.14) and
(3.17) of the previous section §3(b). The other coefficients $a_n$ are determined in §4(a) below from

$$a_n = a_{00}^{[n]} \qquad (3.23b)$$

as given by (A 3)–(A 7) for $n = 2, \dots, 5$, in which the other coefficients $a_{mr}^{[n]}$ are irrelevant here whenever either $m$ or $r$ are non-zero. As the series (3.23a) is embedded in the composite expansion considered in the next §4(a), it is premature to discuss convergence here.

Finally, we recall that, for the critical case, the e-folding distance at the front given by (2.26a) is $\bar{L}_+ = 2\bar{\epsilon}^{1/2} \ell = 2d/\bar{c}$. This $\xi = O(1)$ length-scale is the same as that for the central ‘intermediate’ region. Indeed, from a spatial point of view, the front encompasses the two domains. For that reason $\bar{L}_+$ gives only a qualitative estimate of the front width. So instead, we estimate the distance from where $f$ takes its maximum value $f_{\max}$ to where $f$ essentially vanishes. To that end, we consider the slope $-df/d\xi$ at its maximum $-(df/d\xi)_{\max}$ (Henceforth, values at the corresponding $\xi = \xi_{\max}$ will be identified by $f_{\max}$) with the inflection point property $(d^2f/d\xi^2)_{\max} = 0$. On the basis that $g_{\max}$ like $g_+$ is large compared to unity, (2.6a) determines

$$- \left( \frac{1}{f} \frac{df}{d\xi} \right)_{\max} = \bar{\epsilon} (g_{\max} - 1) \approx \epsilon g_{\max} \qquad (3.24a)$$

at $\xi_{\max}$, Then on approximating the front by the straight line

$$f - f_{\max} = \left( \frac{df}{d\xi} \right)_{\max} (\xi - \xi_{\max}), \qquad (3.24b)$$

we estimate the width of the front from the maximum $f = f_{\max}$ to $f = 0$ as

$$\Delta\xi = - \frac{\aleph}{\bar{\epsilon}g_{+}}, \quad \text{where} \quad \aleph = \frac{f_{\max}g_{+}}{f_{\max}g_{\max}}. \qquad (3.25a,b)$$

Measured in dimensional units of distance $x - ct$ (2.2c), $\Delta\xi$ becomes

$$\Delta(x - ct) = \bar{\epsilon}^{1/2} \ell \Delta\xi = \frac{\aleph \bar{\epsilon}^{1/2} \ell}{\bar{\epsilon} g_{+}} \approx 4 \aleph \bar{\epsilon}^{1/2} \ell = 2 \aleph \bar{L}_{+} = 2 \aleph \mathcal{L}, \qquad (3.25c)$$

where we have noted (2.25), (2.26) as they apply in our large $g_+$ limit. The plausible estimates $f_{\max}/f_{\max} \sim g_+/g_{\max} \sim 2$ give $\aleph = 4$, and whence

$$\Delta(x - ct) \sim 8\mathcal{L}, \qquad (3.25d)$$

somewhat larger than $\mathcal{L}$.

## 4. The large $g_+$ limit ($\ell \gg \mathcal{L}$). (II) Approximate solutions

In §4(a), we construct a composite solution $f(\phi)$ on the entire $\phi$-range based on the asymptotic theory of §3, while in §4(b) we consider briefly the corresponding end point values of $\xi$. In §4(c), we propose a simple heuristic solution which captures the essence of the numerical solution in the large $g_+$ limit.

### (a) Composite solution valid for all $\phi$

In §3(b), we considered the main domain $\phi = O(1)$ and found the lowest order solution $f = f^{[1]}(\phi) = a + e\phi - \kappa g$ (3.10a) correct to $O(\epsilon)$. The leading order terms for the $\xi(\phi)$ expansions were also found but, being cumbersome, we investigate them no further. Instead, we limit attention to $f(\phi)$ and leave $\xi(\phi)$ to be determined at each order by the integral (4.12a) in §4(b) below.
We extend (3.10a) for $f^{[1]}$ to higher orders by considering the series

$$f = a + e\phi + \mathcal{F}(\phi), \quad \mathcal{F} = \sum_{n=1}^{\infty} (\kappa g)^n a_n(\phi). \qquad (4.1a,b)$$

On neglecting the term $e\phi$ in (4.1a) and assuming the coefficients $a_n$ to be constants, the representation (4.1a,b) reduces to the large $g$ approximation (3.23a). This embedding of the front solution of §3(c), and for that matter the tail solution of §3(a), explain our description of (4.1a,b) as a composite solution. Differentiation of (4.1) gives

$$f' = e + \mathcal{B}(\phi), \quad \mathcal{B} = \sum_{n=1}^{\infty} (\kappa g)^n b_n(\phi), \quad b_n = a'_n + n a_n. \qquad (4.1c-e)$$

Substitution of (4.1a,c) into $f(f' + 1) = D + \epsilon\phi - g$ (see (2.3a), (2.7b)) and noting (3.6c,e) yields

$$(a + e\phi)\mathcal{B} + (1 + e)\mathcal{F} + C = 0, \qquad (4.2a)$$

where

$$\mathcal{C} = \mathcal{F}\mathcal{B} + g = \sum_{n=1}^{\infty} (\kappa g)^n c_n(\phi), \quad c_n = \begin{cases} \kappa^{-1} & (n=1), \\ \sum_{s=1}^{n-1} a_{n-s}(\phi) b_s(\phi) & (n \ge 2). \end{cases} \qquad (4.2b,c)$$

Further substitution of the power series (4.1b,d), noting (4.1e), and equating the coefficients of various powers of $g$ yields

$$(a + e\phi)(a'_n + n a_n) + (1 + e)a_n = -c_n \quad (n \ge 1). \qquad (4.3)$$

We solve (4.3) subject to the boundary conditions that $a_n \to 0$ ($n \ge 1$), as $\phi \to -\infty$ (i.e., solutions, that grow exponentially, are rejected). To begin, we considered the lowest order case $n=1$ and expressed the solution of the linear first order ode (4.3) with $c_1 = \kappa^{-1}$ for $a_1(\phi)$ as an integral, which we evaluated by the method of steepest descent to obtain an asymptotic power series expansion. The approach may be applied to the successive cases $n > 1$ for $a_n(\phi)$ but becomes increasingly cumbersome. In all cases, the asymptotic power series takes the form

$$\begin{bmatrix} a_n(\phi) \\ b_n(\phi) \\ c_n(\phi) \end{bmatrix} = \sum_{m=0}^{\infty} e^m \sum_{r=0}^{m} \begin{bmatrix} a_{mr}^{[n]} \\ b_{mr}^{[n]} \\ c_{mr}^{[n]} \end{bmatrix} \phi^r \quad (n \ge 1), \qquad (4.4a)$$

with the curious termination of the r-sum at $r=m$, traceable to the product $e\phi$ in (4.3), and where

$$a_{mr}^{[n]} = 0 \quad \text{whenever} \quad r > m, \qquad (4.4b)$$

$$b_{mr}^{[n]} = n a_{mr}^{[n]} + (r+1)a_{m(r+1)}^{[n]}, \qquad (4.4c)$$

$$c_{mr}^{[1]} = \begin{cases} \kappa^{-1} & (m=r=0). \\ 0 & \text{otherwise.} \end{cases} \qquad (4.4d)$$

Furthermore, in Appendix A(a), we show that, for $n \ge 2$,

$$c_{mr}^{[n]} = \sum_{s=1}^{n-1} \sum_{q=0}^{r} \sum_{p=q}^{q+(m-r)} a_{(m-p)(r-q)}^{[n-s]} \left[ s a_{pq}^{[s]} + (q+1) a_{p(q+1)}^{[s]} \right] \qquad (4.5a)$$

may be expressed more concisely as

$$c_{mr}^{[n]} = n \Upsilon_{mr}^{[n]} + (r+1) \Upsilon_{m(r+1)}^{[n]}, \qquad (4.5b)$$
where

$$ \Upsilon_{mr}^{[n]} = \frac{1}{2} \sum_{s=1}^{n-1} \sum_{q=0}^{r} \sum_{p=q}^{q+(m-r)} a_{(m-p)(r-q)}^{[n-s]} a_{pq}^{[s]}. \quad (4.5c) $$

When $\phi = O(1)$, both g and e are $O(\epsilon)$ (see (2.2b) and (3.6d)), suggesting the reorganised expansion $f = f^{[N]} + O(\epsilon^{N+1})$ of (4.1a,b) with (4.4a), namely

$$ f^{[N]} = a + e\phi + \mathcal{F}^{[N]}(\phi), \quad \mathcal{F}^{[N]} = \begin{cases} 0 & (N=0), \\ \sum_{n=1}^{N} e^n f_n(\phi) & (N \ge 1) \end{cases} \quad (4.6a,b) $$

with the value of the integer N taken to obtain our chosen order of accuracy, and where

$$ e^n f_n = \sum_{m=0}^{n-1} (\kappa g)^{n-m} e^m \left( \sum_{r=0}^{m} a_{mr}^{[n-m]} \phi^r \right) \quad (4.6c) $$

$$ = (\kappa g)^n a_{00}^{[n]} + \dots + e^{n-1} \kappa g (a_{(n-1)0}^{[1]} + \dots + a_{(n-1)(n-1)}^{[1]} \phi^{n-1}) \quad (n \ge 2). \quad (4.6d) $$

To illustrate the nature of the expansion, the first three terms are

$$ \mathbf{e}f_1 = \kappa g a_{00}^{[1]} \quad (a_{00}^{[1]} = -1 \quad \text{to agree with (3.10a)}), \quad (4.7a) $$

$$ \mathbf{e}^2 f_2 = (\kappa g)^2 a_{00}^{[2]} + \mathbf{e} \kappa g (a_{10}^{[1]} + a_{11}^{[1]} \phi), \quad (4.7b) $$

$$ \mathbf{e}^3 f_3 = (\kappa g)^3 a_{00}^{[3]} + \mathbf{e}(\kappa g)^2 (a_{10}^{[2]} + a_{11}^{[2]} \phi) + \mathbf{e}^2 \kappa g (a_{20}^{[1]} + a_{21}^{[1]} \phi + a_{22}^{[1]} \phi^2), \quad (4.7c) $$

consistent with (4.6d), from which we may construct

$$ \mathbf{f}^{[3]} = \mathbf{a} + e\phi + \kappa g [a_{00}^{[1]} + \mathbf{e}(a_{10}^{[1]} + a_{11}^{[1]}\phi) + \mathbf{e}^2(a_{20}^{[1]} + a_{21}^{[1]}\phi + a_{22}^{[1]}\phi^2)] \\ + (\kappa g)^2 [a_{00}^{[2]} + \mathbf{e}(a_{10}^{[2]} + a_{11}^{[2]}\phi)] + (\kappa g)^3 a_{00}^{[3]}. \quad (4.8) $$

Equating the coefficient of $(\kappa g^n)e^m\phi^r$ in equation (4.3) gives

$$ (na + 1)a_{mr}^{[n]} + a(r+1)a_{m(r+1)}^{[n]} + na_{(m-1)(r-1)}^{[n]} + (r+1)a_{(m-1)r}^{[n]} = -c_{mr}^{[n]}. \quad (4.9) $$

Taking all but the first term on the left-hand side to the right-hand side gives

$$ a_{mr}^{[n]} = -\kappa_n [\mathbf{a}(r+1)a_{m(r+1)}^{[n]} + na_{(m-1)(r-1)}^{[n]} + (r+1)a_{(m-1)r}^{[n]} + c_{mr}^{[n]}], \quad (4.10a) $$

in which

$$ \kappa_n = (na + 1)^{-1}, \quad (4.10b) $$

such that $\kappa_1 = \kappa$, and where

$$ c_{mr}^{[n]} = \begin{cases} c_{mr}^{[1]}, & \text{given by (4.4d),} \quad (n=1), \\ n\Upsilon_{mr}^{[n]} + (r+1)\Upsilon_{m(r+1)}^{[n]}, & \text{given by (4.5b),} \quad (n \ge 2). \end{cases} \quad (4.10c) $$

The coefficients in each expression for $e^n f_n$ ($n=1, 2, 3, \dots$) (4.6c,d) may be determined in sequence, beginning with $a_{00}^{[1]} = -\kappa c_{00}^{[1]} = -1$ (see (4.4d), (4.7a) and (A 3)) for $n=1$. For each subsequent $n$, the needed coefficients $a_{mr}^{[n-m]}$ in (4.6d), begin with $m=r=n-1$, namely $a_{(n-1)(n-1)}^{[1]} = -(-\kappa)^{n-1}$ on the right and finish with $m=r=0$, namely $a_{00}^{[n]} = -\kappa_n n \Upsilon_{00}^n$ on the left. The sequential style of the coefficients is illustrated by the explicit formulae for the first few $e^n f_n$ ($n=1, \dots, 5$) given in Appendix A(b).
(b) Convergence and end-point values $\phi_{\pm}^{[N]}$

The series solution (4.1) constructed in §4(a) is based on the assumption that $\phi = O(1)$ for $\epsilon \ll 1$. Nevertheless, even when $\phi \gg 1$, the power series (4.6) gives results that compare remarkably well with the numerical results reported in §5 below almost everywhere, except close to $\phi = \phi_+$.

It should not be overlooked that the end point values,

$$ \phi = \phi_{\pm}^{[N]}, \quad \text{at which} \quad f^{[N]}(\phi) = 0, \qquad (4.11a,b) $$

differ from the true values $\phi_{\pm}$, at which $f(\phi) = 0$. This has the consequence that

$$ \epsilon^{-1} g_{\pm}^{[N]} = g_{\pm}^{[N]} = e^{\phi_{\pm}^{[N]}} \qquad (4.11c) $$

also unfortunately differs from the true values $\epsilon^{-1}g_{\pm} = g_{\pm}$; the matter of $g_{+}^{[N]}$ is explored in §5(b) (see (5.9) and figures 5, 6).

Furthermore, the wave coordinate $(c/d)(x - ct) = \xi = \int_{0}^{\phi} f^{-1} d\phi$ (see (2.2c) and (2.9)) becomes

$$ \xi^{[N]}(\phi) \equiv \int_{0}^{\phi} \frac{d\phi}{f^{[N]}(\phi)}. \qquad (4.12a) $$

In application to the truncated series (4.6a), we need the solution

$$ \phi = \phi^{[N]}(\xi), \quad \text{which solves} \quad \xi^{[N]}(\phi) = \xi, \qquad (4.12b,c) $$

to determine $f^{[N]}$, as a function of $\xi$, in the form

$$ f^{[N]} = f^{[N]}(\phi^{[N]}(\xi)). \qquad (4.12d) $$

(c) A heuristic solution

We consider (2.7b) expressed in the iterated form

$$ f = \frac{F}{1 + f'} \quad \text{where} \quad f' = \frac{1}{2} \left[ -1 + \sqrt{1 + 4(F' - f f'')} \right] \qquad (4.13a,b) $$

(2.12b). In the small $\epsilon$ limit our numerical results reported in the next §5 indicate that $f$ remains close to $F$ suggesting that a reasonable first approximation of $f$ is

$$ f \approx F_0 \equiv F. \qquad (4.14) $$

Moreover, approximations based on small $f'$, in the denominator of (4.13a), gain support from our slowest wave result $\overline{F}' = -\frac{1}{4}$ at $\phi = \phi_+$, where we expect $|f'|$ to take its largest value. Unfortunately (2.14a,c) and (2.18b) indicate slightly larger values

$$ \frac{1}{4} = -\overline{F}'_+ < -F'_+ < -f'_+ < -\overline{f}'_+ = \frac{1}{2}. \qquad (4.15) $$

In this small $f'$ spirit, setting $f' = F'$ in (4.13a) provides an improved approximation

$$ f \approx F_1 \equiv \frac{F}{1 + F'}. \qquad (4.16a) $$

Furthermore, on setting $f f''$ on the right-hand side of (4.13b) equal to $FF''$, we obtain

$$ f \approx F_2 \equiv 2F / \left[ 1 + \sqrt{1 + 4(F' - F F'')} \right], \qquad (4.16b) $$

which has the additional merit that $F'_2 = k_+$ at $\phi = \phi_+$ as demanded by (2.14a) of a true solution.
It should not be overlooked that, for $\phi = O(1)$, the lowest order approximation of $f$ is $f = D$. So correct to lowest order (2.12a) reduces to

$$Df'' + f' = F' = \epsilon(1 - e^{\phi}) \quad (4.17a)$$

with the unique (bounded as $\phi \to -\infty$) solution

$$f' = \epsilon \left[ 1 - \frac{e^{\phi}}{D + 1} \right] \quad (4.17b)$$

in agreement at leading order with the asymptotic result $f^{[1]}' = e - \epsilon\kappa e^{\phi}$ (differentiate (3.10a)). From this point of view, even (4.16a), which employs $f' = F' = \epsilon(1 - e^{\phi})$ (4.17a), is a poor approximation to (4.17b), as $\phi$ increases through positive values. However, for small $\epsilon$, the value of $D$ is bounded above by $\bar{D} = \frac{1}{4}$ and the approximation $f' = F'$ made in (4.16a) is not too bad as our numerical results reported in the next §5 show.

## 5. Numerical solutions

We solved the equations (1.2) numerically subject to the initial conditions

$$v_0(x) = v_+ = (\alpha/\beta)g_+ = 1 \quad (g_+ = 1/\gamma), \quad (5.1a)$$

$$u_0(x) = \begin{cases} 0.3(1-x^2) & \text{if } |x| < 1, \\ 0 & \text{otherwise,} \end{cases} \quad (5.1b)$$

adopted by Hosono and Ilyas [6]. We considered the parameter values

$$\left. \begin{array}{l} d = 1, \\ \alpha = 0.1 \end{array} \right\} \implies \ell = \sqrt{\frac{d}{\alpha}} = \sqrt{10}, \quad (5.2a,b)$$

and generated two data sets by the choices

$$\beta = \begin{cases} 1/3, \\ 1 \end{cases} \implies \frac{\alpha}{\beta} = \begin{cases} 0.3 & \text{case I,} \\ 0.1 & \text{case II,} \end{cases} \quad (5.2c,d)$$

identified by the labels I and II. They determine

$$g_+ = \frac{\beta v_+}{\alpha} = \begin{cases} 10/3, \\ 10 \end{cases} \implies \mathcal{D} = g_+ - \ln g_+ \approx \begin{cases} 2.1294 & \text{I,} \\ 7.6971 & \text{II.} \end{cases} \quad (5.3a,b)$$

### (a) The traveling wave parameter values

In this section, we summarise the derived parameter values needed for comparison of our numerical results with the analytic traveling waves

$$[u, v] = \frac{\alpha}{\beta} [f, g] = \frac{\alpha}{\beta} \left[ \frac{1}{\epsilon} \frac{d\phi}{d\xi}, e^{\phi} \right] = \frac{\alpha}{\beta\epsilon} [f(\xi), g(\xi)], \quad (5.4)$$

where $\xi = (c/d)(x - ct)$, and $\epsilon^{1/2}c = \sqrt{\alpha d} = \alpha\ell$ relates $\epsilon$ to $c$ (see (2.1), (2.2)). Such waves evolve, as $t \to \infty$, from the solutions of the governing equations (1.2) with initial data (5.1) and parameter values (5.2), (5.3), after the transients have decayed. The realised values

$$c \approx \begin{cases} 0.96 & \text{I,} \\ 1.8974 & \text{II} \end{cases} \quad (5.5)$$

compare respectively very well with the slowest wave speeds

$$\bar{c} = 2\alpha\ell\sqrt{g_+ - 1} = \begin{cases} \sqrt{14/15} \approx 0.9661 & \text{I,} \\ \sqrt{18/5} \approx 1.8974 & \text{II,} \end{cases} \quad (5.6a)$$
leading to the corresponding

$$ \bar{\epsilon} = \left( \frac{\alpha l}{\bar{c}} \right)^2 = \frac{1}{4(g_+ - 1)} = \begin{cases} 3/28 \approx 0.1071 & \text{I,} \\ 1/36 \approx 0.0278 & \text{II.} \end{cases} \qquad (5.6b) $$

This agreeable situation, which is in accord with the Dee and Langer criterion explained in §2(c), allows us (up to graph plotting accuracy) to compare the numerical results with the analytic results for $[u, v] = (\alpha/\beta\bar{\epsilon})[f(\xi), g(\xi)]$ (5.4), for which (5.2d), (5.3a) and (5.6b) determine

$$ \frac{\alpha}{\beta \bar{\epsilon}} = \begin{cases} 2.8, \\ 3.6, \end{cases} \qquad g_+ = \frac{\bar{\epsilon} \beta v_+}{\alpha} = \begin{cases} 5/14 & \text{I,} \\ 5/18 & \text{II.} \end{cases} \qquad (5.7a,b) $$

Though $f = d\phi/d\xi$ in (5.4) is expressed as a function of $\xi$, it may be inverted to determine $\xi$ as a function of $\phi$ (see (2.9)). This relationship is needed to link the numerical results to those from our asymptotic development of §4 and appendix A, which builds on the following parameters:

$$ e = \frac{2\bar{\epsilon}}{1 + \sqrt{1 + 4\bar{\epsilon}}} \approx \begin{cases} 0.0976 & \text{I,} \\ 0.0270 & \text{II,} \end{cases} \qquad (5.8a) $$

which solves $e(e + 1) = \bar{\epsilon}$ (3.6c), and

$$ a = De \approx \begin{cases} 0.2079, \\ 0.2082, \end{cases} \qquad \kappa \approx \begin{cases} 0.8279 & \text{I,} \\ 0.8277 & \text{II,} \end{cases} \qquad (5.8b,c) $$

where we have used $a = e\bar{D}/\bar{\epsilon} = De$ (see (3.6e)), and $\kappa = 1/(a+1)$ (3.10b). The numerical values (5.8b) are compatible with the properties $a \uparrow \bar{D} \uparrow \frac{1}{4}$, as $\epsilon \to 0$ (see (2.20b)).

## (b) Comparison of the numerical and analytical results

Hosono and Ilyas [6] portray results for the case (a) in their figure 3, which except for axis scalings resembles our figure 5(a), case I. However, we see from our results portrayed for various $\beta$ in figure 6(b) that their choice $\gamma = g_+^{-1} = 0.3$ (see (5.3a), case I) is only just within the range of applicability of our asymptotic theory. Accordingly, for case II, we reduced $\gamma$ to 0.1.

In figures 1(a,b), we plot the numerical solution $f = f(\phi)$ together with their asymptotic approximations $f = f^{[N]} = f^{[0]} + \mathcal{F}^{[N]}$ (see (4.6a)). Taken correct to order $N$, the series $\mathcal{F}^{[N]}(4.6b)$ ($N \ge 1$) has coefficients $a_{mr}^{[n]}$ ($0 \le r \le m, n \le N$), introduced in (4.4a), which are given by (A 3)–(A 8) in appendix A(b) for $1 \le N \le 5$. Each $f^{[N]}$ intersects the the $\phi$-axis at two points $\phi = \phi_{\pm}^{[N]}$ (4.11a,b). The lowest order approximation $f^{[0]} = a + e\phi$ is appropriate to the large negative $\phi$-range, $-\phi = O(\epsilon^{-1}) \gg 1$, discussed in §3(a), and gives $\phi_{-}^{[0]} = -a/e = -\mathcal{D}$ (see (3.6b,e)). Generally, the neglect of exponentially small terms is justified with all $f^{[N]}$-curves terminating close to $\phi_{-}$. However, figure 1(a), shows $\phi_{-}^{[0]}$ differing significantly from $\phi_{-}$, supporting the view expressed in the previous paragraph that the value $\bar{\epsilon} = 3/28$ (5.6b), for case I $\beta = 1/3$, is not small enough for the neglect of exponentially small terms, required by asymptotic theory.

When $\phi = O(1)$ discussed in §3(b), the first (large) exponential correction to the straight line $f^{[0]}$ is $ef_1 = f^{[1]} - f^{[0]} = \mathcal{F}^{[1]} = -\epsilon\kappa e^{\phi}$ (see (3.10), also (4.7a) and the definition above it) leading to the curve $f^{[1]}$, which gives remarkably good agreement with the numerical solution. On adding successively $e^2f_2$ (4.7b), $e^3f_3$ (4.7c), ... to obtain sequentially $f^{[2]}$, $f^{[3]}$ (4.8), ... better agreement with the numerics is obtained everywhere at each higher order, as summarised compactly by figures 6.

The evident weakness of the asymptotic solutions for $\phi \gg 1$ in figure 1 (even more pronounced in the plots versus g in figure 3) prompted our alternative heuristic approximations $F_0 = F$ and $F_1, F_2$ (4.16a,b) described in §4(c). Plots of all three are illustrated in figures 2 together with f for comparison. The figures show that F itself is a reasonable approximation differing little from f, whose maximum at $\phi = \phi_{max}$ is close to $\phi = 0$, where F is maximum. The relative behaviour
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Plots of $f$ versus $\phi$. The full numerical solution (black) is compared to the reference curve $F$ (2.3a) (dashed red), both with the same end-points at $f = F = 0$. The others display the asymptotic solution $f = f^{[N]}$ (see (4.6)) with error $O(\epsilon^{N+1})$ at various truncation levels $N$. The lowest order linear solution $f^{[0]}$ is the highest straight line (dashed black). The higher order solutions modified exponentials are successively $f = f^{[1]}$ (magenta), $f^{[3]}$ (blue), and $f^{[5]}$ (red); (a) case I, $\beta = 1/3$, (b) case II, $\beta = 1$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** As in figure 1, with the full numerical solution $f$ (black) and the reference curve $F$ (dashed red), which is also the first heuristic approximation $F_0$. The higher order approximations (4.16a,b) are $F_1$ (red) and $F_2$ (green).</figcaption></figure>

of $f$ and $F$ were discussed at length in §§2(a),(b) but note, particularly their crossing at $\phi_{max}$ (see (2.13c)). Agreement with $f$, however, is significantly improved by $F_1$. The more complicated approximation $F_2$ has the merit, commented on below (4.16b), that by construction its slope at the end point $\phi_+$ is the same as that of $f$. The $F_1$, $F_2$ and $f$-curve plots are virtually indistinguishable.

As $g = \epsilon g (= \epsilon e^{\phi})$ is the primitive variable, we provide in figures 3 and 4 plots of $f$ versus $g$, rather than $\phi$. This alteration has other merits, as fore and aft of $\phi = 0$, namely $g = \bar{\epsilon}$, the $\phi$ and $g$ scales are very different. Specifically, relative to the $\phi$-plots of figures 1 and 2, the $g$-plots of figures 3 and 4 compress (expand) the horizontal scale for $\phi < 0$ ($> 0$). This permits us to examine left-hand ($\phi < 0$) behaviour in greater detail on figures 1 and 2, particularly near $\phi_-$, and right-hand ($g > \bar{\epsilon}$) behaviour on figures 3 and 4, particularly near $g_+$, where the end points $g = g_+^{[N]}$, at which $f^{[N]} = 0$, (see (4.11)), clearly differ.

The asymmetries fore and aft visible on figures 1, 2 and 3, 4, are consolidated in figure 5 with plots of $u$ and $v$ versus the space coordinate $x - ct = (d/c)\xi$ in the frame moving with the wave. Essentially, as the wave speed $\bar{c}$ increases, the $u$ profile develops an ever increasing asymmetry with the formation of a front at its leading edge, illustrated well by the small $\alpha/\beta = 0.1$, case II, in figure 5(b). The feature has the following interpretation: As the wave of infectives $u$ travels to the
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Same as in figure 1, but plots are versus $g = \bar{e}e^{\phi}$ instead of $\phi$. The end point g-values for the numerics and F, given by (5.7b) are (a) case I, $g_+ = 5/14$, (b) case II, $g_+ = 5/18$, while the corresponding $g_+^{[N]}$ (see (4.11c)) exceed these values.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Same as in figure 2, but plots are versus $g = \bar{e}e^{\phi}$ instead of $\phi$.</figcaption></figure>

right into the hoard of susceptibles $v_+$ ($\xi > 0$), almost all are rapidly infected until $u$ is maximised at $\phi = \phi_{max}$ close to $\xi = 0$, where $v/v_+$ is very small. Thereafter, the decay of infectives $u$ is slow.

Concerning the asymptotic solutions, we recall that our $\xi$-origin $\xi = 0$ is located at $\phi = 0$ (close to $\phi_{max}$), where $g = e^{\phi} = 1$. There, $v = \alpha/\beta = 0.3$, for case I, and $v = 0.1$, for case II, portrayed in figures 5(a,b) respectively. Since $\phi^{[N]} = 0$ at $\phi = 0$ and $g^{[N]} = e^{\phi^{[N]}}$ (see (4.11) and (4.12)), all asymptotically truncated solutions $v^{[N]} = (\alpha/\beta)g^{[N]}$ intersect the $v$-curve at $\xi = 0$. This is confirmed by an extremely close inspection of figures 5, which also shows that $v^{[N]} \ge v$ for $\xi \ge 0$ for the values $N = 3$ and 5 selected.

We assess the failure of $v_+^{[N]}$ to drop to the true $v_+$ by displaying the error ratio

$$ \varepsilon^{[N]} = (v_+ - v_+^{[N]})/v_+ \qquad (v_+ = 1 \text{ see (5.1a)}) \qquad (5.9a) $$

in figures 6. The log-log plot figure 6(a) clearly shows systematic but slow convergence on increasing the order $N$ at fixed $\beta = \frac{1}{3}$ (case I) and $\beta = 1$ (case II). The almost linearity of both sets of points is approximated well by the power law fit

$$ \varepsilon^{[N]} \propto N^{-\sqrt{2}}, \qquad (5.9b) $$

with the intriguing index $-\sqrt{2}$. Figure 6(b) shows how the accuracy of the solution varies little on increasing $\beta$ but the unreliability of the results suggested below $\beta = \frac{1}{3}$ supports our view that $\beta = \frac{1}{3}$ is as close to the lower limit of the reliability of our large $v_+$ asymptotics.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Plots of *u* and *v* versus *x* − *ct* with the full numerics (black). Our asymptotic description of *u* and *v* derived from *f* = $f^{[3]}$ (blue) and *f* = $f^{[5]}$ (red) are also shown. The various *u* plots are virtually indistinguishable. The same is true for the *v* curves, except near the right of the figures, where the asymptotic solutions fail to meet *v* = 1 boundary condition attained by the full numerics; (a) case I, $\beta$ = 1/3, (b) case II, $\beta$ = 1.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The error $\mathcal{E}$ (5.9a) of the asymptotic value of $v^{[N]}$ relative to the numerical value $v$, as $x \to \infty$ for orders $N=1$ (magenta), $N=2$ (green), $N=3$ (blue), $N=4$ (cyan), $N=5$ (red). (a) The $\mathcal{E}$ vs. $N$ points (both log scales) for our two cases I, $\beta=1/3$ ○'s, II, $\beta=1$ ☐'s. The dashed line is the power law fit (5.9b). (b) Curves of $\mathcal{E}$ vs. $\beta$ at various orders $N$, with the values at $\beta=1/3$ (case I) and $\beta=1$ (case II) identified by ○'s and ☐'s, as in (a).</figcaption></figure>

**6. Concluding remarks**

Our main conclusion is that the asymptotic series solution $f^{[N]}$ (4.6), for an appropriate choice of $N$, leads to a remarkably accurate approximation of $u$ as a function of $v$, provided that $v_+ \gtrsim \alpha/\beta$, i.e., $\mathcal{L} \lesssim \ell$. The largest $N(=5)$ considered shows good convergence for small $\gamma=0.1$, i.e. large basic reproduction number $R_0$ (see figures 6, $\beta=1$). Our form (4.8) for $f^{[3]}$ deserves comparison with the simpler expansion, equ. (14) of Canosa [14], for the FKPP model. El-Hachem et al. [15] provide phase plane plots in their figure 11 (cf. our figure 3) based on the formula (14) of [14] and remark on the accuracy of that series at large $c$.

Having obtained $u=u(v)$, the complete solution requires their dependence on the traveling wave coordinate $x-ct=\ell e^{1/2}\xi$ (2.2c). That is determined implicitly by the integral (2.9) which gives $\xi=\xi(v)$. The integration is awkward and inversion to obtain $v=v(\xi)$ and $u=u(\xi)$ worse. Some progress was possible at low order, particularly $N=1$ (see (3.12)–(3.16)). Generally, and to obtain the results for $N=3, 5$ portrayed on figure 5, we integrated (2.9) numerically. Of course, numerical integration of the governing equations (1.2), treated as an initial value problem, determines (after the transients have decayed) $u=u(\xi)$ and $v=v(\xi)$. All our analytic results
displayed in figures 1–6 adopt the choice of $c = \bar{c}$, predicted by the Dee and Langer condition, whose relevance to the numerical solutions is confirmed by comparing (5.6a) with (5.5).

Though $\phi$ has three ranges $-\phi \gg 1$ §3(a), $\phi = O(1)$ §3(b) and $\phi = O(\ln(\epsilon^{-1}))$ §3(c), it is remarkable that, for the spatio-temporal dependence on $\xi$, there are only two, $-\xi = O(\epsilon^{-1})$ §3(a) and $\xi = O(1)$ §3(b,c), such that the central and frontal regions blend into one. The situation is highlighted in figure 5(b), case II, where the collapse of succptibles $v$ identifies the frontal region, which terminates, where the infectives are maximised. Based on $\epsilon \ll 1$, we estimate the short front width ahead of that maximum to be in the order of $16\bar{\epsilon}^{1/2}\ell$ (see (3.25c,d), $\aleph = 4$), namely $8\mathcal{L} = 8\sqrt{d/(\beta v_+)}$ (see (1.4b) and (2.26a,c)). Behind the maximum, the infectives decay relatively slowly on the long length $L_- \approx \ell/\epsilon^{1/2}$ (3.8d), determined by the fast wave speed $\bar{c}$ and the exponential decay $\propto e^{-\alpha t}$ at fixed locations on the time scale $T_- = 1/\alpha$ (3.8e). For further clarification, we recall that the length $\ell$ is determined by the role of diffusion in the governing equations (see (1.4b)), while $\bar{\epsilon}^{1/2} = \alpha\ell/\bar{c} \approx \frac{1}{2}/\sqrt{g_+} = \frac{1}{2}\ell/\mathcal{L}$ (see (3.1e)) is a measure of the nonlinearity introduced via the boundary condition $v = v_+ = (\alpha/\beta)g_+$ (see (1.8)), with the consequence that $\bar{c} \propto \sqrt{v_+}$. Thus, the front shortens proportional to $1/\sqrt{v_+}$, while the tail lengthens $\propto \sqrt{v_+}$. Unlike the tail, the front exhibits a fully advective-diffusive balance.

Finally we remark, that, in real situations, the space in which epidemics evolve is fully 2-dimensional. If the infection originated in the vicinity of a point, we might expect the infection to spread with a roughly circular front. For that our 1-dimensional analysis would apply at a sufficiently long time. when the front is far from the focus of the waves and so locally straight.

**Acknowledgements.** This work was motivated by the Covid-19 pandemic and started during the lockdown periods in the UK and in France, in 2020, when both authors were in isolation. This research would not have been possible without the use of videoconferencing. Both authors are grateful to Raymond E. Goldstein and Steve M. Tobias for useful discussions in the course of this work.

# APPENDICES

## A. The composite expansion

$$ (a) \text{ Derivation of (4.5b): } c_{mr}^{[n]} = n\Upsilon_{mr}^{[n]} + (r+1)\Upsilon_{m(r+1)}^{[n]} \quad (n \ge 2) $$

We write $c_{mr}^{[n]}$ defined by (4.5a) in the alternative form

$$ c_{mr}^{[n]} = n\alpha_{mr}^{[n]} + (r+1)\beta_{m(r+1)}^{[n]} \quad (n \ge 2), \quad (A 1a) $$

where

$$ \alpha_{mr}^{[n]} = \frac{1}{n} \sum_{s=1}^{n-1} \sum_{q=0}^{r} \sum_{p=q}^{q+(m-r)} a_{(m-p)(r-q)}^{[n-s]} sa_{pq}^{[s]}, \quad (A 1b) $$

$$ \beta_{m(r+1)}^{[n]} = \frac{1}{r+1} \sum_{s=1}^{n-1} \sum_{q=0}^{r} \sum_{p=q+1}^{q+(m-r)} a_{(m-p)(r-q)}^{[n-s]} (q+1) a_{p(q+1)}^{[s]}, \quad (A 1c) $$

which, following the change of variables $r+1 \mapsto r$, $q+1 \mapsto q$, becomes

$$ \beta_{mr}^{[n]} = \frac{1}{r} \sum_{s=1}^{n-1} \sum_{q=0}^{r} \sum_{p=q}^{q+(m-r)} a_{(m-p)(r-q)}^{[n-s]} qa_{pq}^{[s]}. \quad (A 1d) $$

Here, we have left the vanishing $q=0$ term in the middle $q$-sum so that the summation ranges agree with those for $\alpha_{mr}^{[n]}$ (A 1b).

For the $\alpha_{mr}^{[n]}$ (A 1b), $\beta_{mr}^{[n]}$ (A 1d) and $\Upsilon_{mr}^{[n]}$ (4.5c), we interchange the summation variables $s \mapsto n-s$, $p \mapsto m-p$, $q \mapsto r-q$, and add the new sums to their respective originals to obtain

$$ \alpha_{mr}^{[n]} = \Upsilon_{mr}^{[n]} = \beta_{mr}^{[n]}, \quad (A 2) $$

which on substitution into (A 1a) yields (4.5b), as required.
(b) Sequential derivation of $a_{mr}^{[n]}$ (4.10a)

We collect together the coefficients $a_{mr}^{[n]}$ (see (4.10a,b)) in explicit form, as needed to define $e^n f_n$ at each order $e^n$ ($n = 1, \dots, 5$). Their evaluation may be executed consecutively, following the sequence listed below, where the right-hand sides of each successive formula is determined by the previously calculated coefficients.

ef₁ coefficient:

$$a_{00}^{[1]} = -1. \qquad (A3)$$

e²f₂ coefficients:

$$a_{11}^{[1]} = -\kappa a_{00}^{[1]} = \kappa, \qquad (A4a)$$

$$a_{10}^{[1]} = -\kappa [aa_{11}^{[1]} + a_{00}^{[1]}] = \kappa^2, \qquad (A4b)$$

$$a_{00}^{[2]} = -\kappa_2 2\Upsilon_{00}^{[2]} = -\kappa_2, \qquad \Upsilon_{00}^{[2]} = \frac{1}{2}(a_{00}^{[1]})^2. \qquad (A4c)$$

e³f₃ coefficients:

$$a_{22}^{[1]} = -\kappa a_{11}^{[1]} = -\kappa^2, \qquad (A5a)$$

$$a_{21}^{[1]} = -\kappa [2aa_{22}^{[1]} + a_{10}^{[1]} + 2a_{11}^{[1]}], \qquad (A5b)$$

$$a_{20}^{[1]} = -\kappa [aa_{21}^{[1]} + a_{10}^{[1]}], \qquad (A5c)$$

$$a_{11}^{[2]} = -\kappa_2 [2a_{00}^{[2]} + 2\Upsilon_{11}^{[2]}], \qquad \Upsilon_{11}^{[2]} = a_{00}^{[1]}a_{11}^{[1]}, \qquad (A5d)$$

$$a_{10}^{[2]} = -\kappa_2 [aa_{11}^{[2]} + a_{00}^{[2]} + 2\Upsilon_{10}^{[2]} + \Upsilon_{11}^{[2]}], \qquad \Upsilon_{10}^{[2]} = a_{00}^{[1]}a_{10}^{[1]}, \qquad (A5e)$$

$$a_{00}^{[3]} = -\kappa_3 3\Upsilon_{00}^{[3]}, \qquad \Upsilon_{00}^{[3]} = a_{00}^{[1]}a_{00}^{[2]}. \qquad (A5f)$$

e⁴f₄ coefficients:

$$a_{33}^{[1]} = -\kappa a_{22}^{[1]} = \kappa^3, \qquad (A6a)$$

$$a_{32}^{[1]} = -\kappa [3aa_{33}^{[1]} + (a_{21}^{[1]} + 3a_{22}^{[1]})], \qquad (A6b)$$

$$a_{31}^{[1]} = -\kappa [2aa_{32}^{[1]} + (a_{20}^{[1]} + 2a_{21}^{[1]})], \qquad (A6c)$$

$$a_{30}^{[1]} = -\kappa [aa_{31}^{[1]} + a_{20}^{[1]}], \qquad (A6d)$$

$$a_{22}^{[2]} = -\kappa_2 [2a_{11}^{[2]} + 2\Upsilon_{22}^{[2]}], \qquad \Upsilon_{22}^{[2]} = a_{00}^{[1]}a_{22}^{[1]} + \frac{1}{2}(a_{11}^{[1]})^2, \qquad (A6e)$$

$$a_{21}^{[2]} = -\kappa_2 [2aa_{22}^{[2]} + 2(a_{10}^{[2]} + a_{11}^{[2]}) + 2\Upsilon_{21}^{[2]} + 2\Upsilon_{22}^{[2]}], \qquad \Upsilon_{21}^{[2]} = a_{00}^{[1]}a_{21}^{[1]} + a_{11}^{[1]}a_{10}^{[1]}, \qquad (A6f)$$

$$a_{20}^{[2]} = -\kappa_2 [aa_{21}^{[2]} + a_{10}^{[2]} + 2\Upsilon_{20}^{[2]} + \Upsilon_{21}^{[2]}], \qquad \Upsilon_{20}^{[2]} = a_{00}^{[1]}a_{20}^{[1]} + \frac{1}{2}(a_{10}^{[1]})^2, \qquad (A6g)$$

$$a_{11}^{[3]} = -\kappa_3 [3a_{00}^{[3]} + 3\Upsilon_{11}^{[3]}], \qquad \Upsilon_{11}^{[3]} = a_{00}^{[1]}a_{11}^{[2]} + a_{11}^{[1]}a_{00}^{[2]}, \qquad (A6h)$$

$$a_{10}^{[3]} = -\kappa_3 [aa_{11}^{[3]} + a_{00}^{[3]} + 3\Upsilon_{10}^{[3]} + \Upsilon_{11}^{[3]}], \qquad \Upsilon_{10}^{[3]} = a_{00}^{[1]}a_{10}^{[2]} + a_{00}^{[2]}a_{10}^{[1]}, \qquad (A6i)$$

$$a_{00}^{[4]} = -\kappa_4 4\Upsilon_{00}^{[4]}, \qquad \Upsilon_{00}^{[4]} = a_{00}^{[1]}a_{00}^{[3]} + \frac{1}{2}(a_{00}^{[2]})^2. \qquad (A6j)$$
e⁵f₅ coefficients:

$$a_{44}^{[1]} = -\kappa a_{33}^{[1]} = -\kappa^4, \qquad (A 7a)$$

$$a_{43}^{[1]} = -\kappa [4aa_{44}^{[1]} + a_{32}^{[1]} + 4a_{33}^{[1]}], \qquad (A 7b)$$

$$a_{42}^{[1]} = -\kappa [3aa_{43}^{[1]} + a_{31}^{[1]} + 3a_{32}^{[1]}], \qquad (A 7c)$$

$$a_{41}^{[1]} = -\kappa [2aa_{42}^{[1]} + a_{30}^{[1]} + 2a_{31}^{[1]}], \qquad (A 7d)$$

$$a_{40}^{[1]} = -\kappa [aa_{41}^{[1]} + a_{30}^{[1]}], \qquad (A 7e)$$

$$a_{33}^{[2]} = -\kappa_2 [2a_{22}^{[2]} + 2\gamma_{33}^{[2]}], \qquad (A 7f)$$

$$a_{32}^{[2]} = -\kappa_2 [3aa_{33}^{[2]} + 2a_{21}^{[2]} + 3a_{22}^{[2]} + 2\gamma_{32}^{[2]} + 3\gamma_{33}^{[2]}], \qquad (A 7g)$$

$$a_{31}^{[2]} = -\kappa_2 [2aa_{32}^{[2]} + 2a_{20}^{[2]} + 2a_{21}^{[2]} + 2\gamma_{31}^{[2]} + 2\gamma_{32}^{[2]}], \qquad (A 7h)$$

$$a_{30}^{[2]} = -\kappa_2 [aa_{31}^{[2]} + a_{20}^{[2]} + 2\gamma_{30}^{[2]} + \gamma_{31}^{[2]}], \qquad (A 7i)$$

$$a_{22}^{[3]} = -\kappa_3 [3a_{11}^{[3]} + 3\gamma_{22}^{[3]}], \qquad (A 7j)$$

$$a_{21}^{[3]} = -\kappa_3 [2aa_{22}^{[3]} + 3a_{10}^{[3]} + 2a_{11}^{[3]} + 3\gamma_{21}^{[3]} + 2\gamma_{22}^{[3]}], \qquad (A 7k)$$

$$a_{20}^{[3]} = -\kappa_3 [aa_{21}^{[3]} + a_{10}^{[3]} + 3\gamma_{20}^{[3]} + \gamma_{21}^{[3]}], \qquad (A 7l)$$

$$a_{11}^{[4]} = -\kappa_4 [4a_{00}^{[4]} + 4\gamma_{11}^{[4]}], \qquad (A 7m)$$

$$a_{10}^{[4]} = -\kappa_4 [aa_{11}^{[4]} + a_{00}^{[4]} + 4\gamma_{10}^{[4]} + \gamma_{11}^{[4]}], \qquad (A 7n)$$

$$a_{00}^{[5]} = -\kappa_5 [5\gamma_{00}^{[5]}], \qquad (A 7o)$$

where

$$\gamma_{33}^{[2]} = a_{00}^{[1]}a_{33}^{[1]} + a_{11}^{[1]}a_{22}^{[1]}, \qquad (A 8a)$$

$$\gamma_{32}^{[2]} = a_{00}^{[1]}a_{32}^{[1]} + a_{11}^{[1]}a_{21}^{[1]} + a_{22}^{[1]}a_{10}^{[1]}, \qquad (A 8b)$$

$$\gamma_{31}^{[2]} = a_{00}^{[1]}a_{31}^{[1]} + a_{10}^{[1]}a_{21}^{[1]} + a_{20}^{[1]}a_{11}^{[1]}, \qquad (A 8c)$$

$$\gamma_{30}^{[2]} = a_{00}^{[1]}a_{30}^{[1]} + a_{10}^{[1]}a_{20}^{[1]}, \qquad (A 8d)$$

$$\gamma_{22}^{[3]} = (a_{00}^{[1]}a_{22}^{[2]} + a_{00}^{[2]}a_{22}^{[1]}) + a_{11}^{[1]}a_{11}^{[2]}, \qquad (A 8e)$$

$$\gamma_{21}^{[3]} = (a_{00}^{[1]}a_{21}^{[2]} + a_{00}^{[2]}a_{21}^{[1]}) + (a_{10}^{[1]}a_{11}^{[2]} + a_{10}^{[2]}a_{11}^{[1]}), \qquad (A 8f)$$

$$\gamma_{20}^{[3]} = (a_{00}^{[1]}a_{20}^{[2]} + a_{00}^{[2]}a_{20}^{[1]}) + a_{10}^{[1]}a_{10}^{[2]}, \qquad (A 8g)$$

$$\gamma_{11}^{[4]} = a_{00}^{[1]}a_{11}^{[3]} + a_{00}^{[2]}a_{11}^{[2]} + a_{00}^{[3]}a_{11}^{[1]}, \qquad (A 8h)$$

$$\gamma_{10}^{[4]} = a_{00}^{[1]}a_{10}^{[3]} + a_{00}^{[2]}a_{10}^{[2]} + a_{00}^{[3]}a_{10}^{[1]}, \qquad (A 8i)$$

$$\gamma_{00}^{[5]} = a_{00}^{[1]}a_{00}^{[4]} + a_{00}^{[2]}a_{00}^{[3]}. \qquad (A 8j)$$
References

1. Kermack, W.O., McKendrick, A.G. 1927 A contribution to the mathematical theory of epidemics. Proc. Roy. Soc A 115, 700-721. See https://doi.org/10.1098/rspa.1927.0118.

2. Ram, V., Schaposnik, L.P. 2021 A modified age-structured SIR model for COVID-19 type viruses. Sci. Rep. 11, 15194. See https://doi.org/10.1038/s41598-021-94609-3.

3. Odagaki, T. 2023 New compartment model for COVID-19. Sci. Rep. 13, 5409. See https://doi.org/10.1038/s41598-023-32159-6.

4. Fowler, A.C., Hollingsworth, T.D. 2015 Simple approximations for epidemics with exponential and fixed infectious periods. Bull. Math. Biol. 77(8), 1539-1555. See https://doi.org/10.1007/s11538-015-0095-3.

5. Källén, A. 1984 Thresholds and travelling waves in an epidemic model for rabies. Nonlinear Anal.-Theor. 8(8), 851-856. See https://doi.org/10.1016/0362-546X(84)90107-X.

6. Hosono, Y., Ilyas, B. 1995 Travelling waves for a simple diffusive epidemic model. Math. Mod. Meth. Appl. S. 5(7), 935-966. See https://doi.org/10.1142/S0218202595000504.

7. Källén, A., Arcuri, P.A., Murray, J.D. 1985 A simple model for the spatial spread and control of rabies. J. Theor. Biol. 116 3, 377-93. See https://doi.org/10.1016/S0022-5193(85)80276-9.

8. Murray, J.D. 2002 Mathematical Biology: I. An Introduction Interdisciplinary Applied Mathematics, Volume 17, Eds. S.S. Antman, J.E. Marsden, L. Sirovich, S. Wiggins. ISBN 0-387-95223-3, Springer-Verlag, New York, Berlin, Heidelberg.

9. Locatelli, I., Trächsel, B., Rousson, V. 2021 Estimating the basic reproduction number for COVID-19 in Western Europe. PLoS ONE 16(3): e0248731. See https://doi.org/10.1371/journal.pone.0248731.

10. Ai, S., Huang, W. 2005 Travelling waves for a reaction-diffusion system in population dynamics and epidemiology. P. Roy. Soc. Edinb. A135, 663-675. See https://doi.org/10.1017/S0308210500004054.

11. Murray, J.D. 2003 Mathematical Biology: II. Spatial Models and Biomedical Applications Interdisciplinary Applied Mathematics, Volume 17, Eds. S.S. Antman, J.E. Marsden, L. Sirovich, S. Wiggins. ISBN 0-387-95228-4, Springer-Verlag, New York, Berlin, Heidelberg.

12. Popović, N., Kaper, T.J. 2006 Rigorous asymptotic expansions for critical wave speeds in a family of scalar reaction-diffusion equations. J. Dyn, Differ. Equ 18(1), 103-139. See https://doi.org/10.1007/s10884-005-9002-1.

13. Kreten, F. 2022 Traveling waves of an FKPP-type model for self-organized growth. J. Math. Biol. 84, 42. See https://doi.org/10.1007/s00285-022-01753-z.

14. Canosa, J. 1973 On a nonlinear diffusion equation describing population growth. IBM J. Res. Dev. 17 307-313. See https://doi.org/10.1147/rd.174.0307.

15. El-Hachem, M., McCue, S.W., Simpson, M.J. 2022 Non-vanishing sharp-fronted travelling wave solutions of the Fisher-Kolmogorov model. Math. Med. Biol. 39(3), 226-250. See https://doi.org/10.1093/imammb/dqac004.

16. Fisher, R.A. 1937 The wave of advance of advantageous genes. Ann. Eugenics 7, 353-369. See https://doi.org/10.1111/j.1469-1809.1937.tb02153.x.

17. Kolmogorov, A.N., Petrovskii, I.G., Piscounov, N.S. 1937 Étude de l'équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique. Moscow University, Bull. Math. 1, 1-25.

18. Simitev, R.D., Biktashev, V.N. 2011 Asymptotics of conduction velocity restitution in models of electrical excitation in the heart. Bull. Math. Biol. 73(1), 72-115. See https://doi.org/10.1007/s11538-010-9523-6.

19. Bezekci, B., Idris, I., Simitev, R.D., Biktashev, V.N. 2015 Semianalytical approach to criteria for ignition of excitation waves. Phys. Rev. E 92, 042917. See https://doi.org/10.1103/PhysRevE.92.042917.

20. Ramaswamy, H., Oberai, A.A., Yortsos, Y.C. 2021 A comprehensive spatial-temporal infection model. Chem. Eng. Sci. 233, 116347. See https://doi.org/10.1016/j.ces.2020.

21. Dee, G., Langer, J.S. 1983 Propagating pattern selection. Phys. Rev. Lett. 50(6), 383-386. See https://doi.org/10.1103/PhysRevLett.50.383.

22. van Saarloos, W. 2003 Front propagation into unstable states. Phys. Rep. 386(2-6), 29-222. See https://doi.org/10.1016/j.physrep.2003.08.001.

23. Satnoianu, R.A., Merkin, J.H., Scott, S.K. 2000 Pattern formation in a differential-flow reactor model. Chem. Eng. Sci. 55 461-469. See https://doi.org/10.1016/S0009-2509(99)00340-1.