2. General setting and numerical solutions for (IBVP)

To begin this section, we briefly recall fundamental preliminary results concerning (IBVP), which can be obtained directly from application of the theory developed in [Reference Hamel and Ryzhik17]. We first observe, with key details established in [Reference Hamel and Ryzhik17] (Theorem 1.3), that (IBVP) has a unique, global (on $\overline{D}_{\infty }$ ) solution, which, for each $T\gt 0$ , depends continuously on intial data, throughout $\overline{D}_T$ . Moreover, $||u(\cdot, t)||_{\infty }$ is uniformly bounded on $[0,\infty )$ by a constant depending only upon $A$ and $D$ (although, for given $A\ge 0$ , we record from [Reference Hamel and Ryzhik17] (Remark 1.4) that this constant blows up as $D\to 0^+$ ). In addition, it is trivially established, via the strong maximum principle and comparison theorem, that, for any $t_0$ sufficiently large,

(10) \begin{align} 0\lt u(x,t)\le (t+t_0)^{-\frac{1}{2}}e^{(t+t_0)}e^{-x^2/4D(t+t_0)} \quad \forall \, (x,t)\in D_T, \end{align}

which provides useful information in later sections. However, we may observe immediately from (10) that should an identifiable wavefront location and structure develop in the solution to (IBVP) as $t\to \infty$ , say at location $|x|\sim S(t)$ (beyond which $u$ is exponentially small as $t\to \infty$ ), then,

(11) \begin{align} 0\lt S(t)\leq 2\sqrt{D}t-\frac{1}{2}\sqrt{D}\log t+ O(1), \end{align}

as $t\to \infty$ . We should remark here that Bouin, Henderson and Ryzhik [Reference Bouin, Henderson and Ryzhik6] establish an asymptotic estimate giving wavefront location at large- $t$ for the solution to the nonlocal Fisher-KPP equation, when the initial data is front-like and localised to the (wlog) left half-line, and in the case of the top hat kernel, this gives a form which satisfies the above inequality, with the coefficient of the $\log t$ term being $-\frac{3}{2}$ , as in the Bramson correction for the same evolution problem with the classical local Fisher-KPP equation. We, tentatively, anticipate that in the present situation, $S(t)$ will be similarly asymptotic to this form as $t\to \infty$ .

In the remainder of this section, we begin our principal study of the detailed qualitative and quantitative features of the solution to (IBVP), and in particular how these properties respond to decreasing the diffusion parameter $D$ . We will see that significant changes in structure and mechanism occur, particularly as $D$ decreases from moderately small ( $\sim 10^{-3}$ ) to extremely small ( $\sim 10^{-6}$ ). With this in mind, we develop a numerical scheme to approximate solutions to (IBVP) and use this to investigate the qualitative and quantitative structure of solutions to (IBVP), with particular attention to the structure of the solution as $t\to \infty$ . We note that similar evolutionary computations have been made and presented in [Reference Nadin, Perthame and Tang22], concerning a stability question relating to permanent form travelling wave structures connecting two unstable states of equation (1). Here our aim is to form a preliminary approach to elucidating the qualitative and quantitative properties of the evolution problem (IBVP). For the first set of numerical solutions that we present, we take $x_0=\frac{1}{2}$ and $g\;:\;{\mathbb{R}}\to{\mathbb{R}}$ as

(12) \begin{align} g(x) & = \begin{cases} (1-2x)^2(1+2x)^2, & |x|\leq \frac{1}{2} \\ 0, & |x|\gt \frac{1}{2}. \end{cases} \end{align}

We discretise $u$ on a uniform spatial grid of $N$ points, truncated to $0 \leq x \leq L$ , approximating the second derivative using central finite differences. The convolution term is evaluated using the trapezium rule, in other words assuming a linear variation of $u$ between grid points, dealing carefully with cases where the edge of the support of the top hat kernel lies between grid points. We also take into account the symmetry of the solution about $x=0$ and assume that $u=0$ for $x\gt L$ . In the simulations discussed below, we take $L = 10$ and $N = 1000$ . Timestepping is done using the midpoint method (second order accurate), with time step chosen adaptively so that the maximum change in $u$ at each step is below $10^{-2}$ .

Figure 1.

The numerical solution of (IBVP) for various values of $D$ with $A =0.01$ (solid black line), along with the minimum speed travelling wave (broken blue line).

Figure 1 shows the solution of (IBVP) when $t = 9/2 \sqrt{D}$ , which is just before the wavefront reaches the edge of the truncated domain, with $A = 0.01$ . For all values of $A$ that we investigated, indeed for all localised initial inputs of $u$ that we tried, the solution was qualitatively similar to those shown in Figure 1. In each case, a wavefront propagates in the positive $x$ -direction. Also shown is the corresponding minimum speed travelling wave solution (see Section 6 for details of the minimum speed travelling wave solution), calculated numerically using the same finite difference method to set up the discretised equations and ‘fsolve’ in Matlab to solve them. The equilibrium state $u=1$ is temporally unstable for $D \lt \Delta _1 \approx 0.00297$ whilst temporally stable for $D$ larger than this critical value (see Section 3 below). For $D \gt \Delta _1$ , the minimum speed travelling wave solution emerges, which leaves the equilibrium state $u=1$ in its wake. For $0\lt D \lt \Delta _1$ , however, a non-propagating stationary, spatially periodic state is left in the wake of the wavefront. For moderately small values of $D$ , as shown in Figure 1, the wavelength of this periodic state is close to $0.7$ , which is close to the most unstable wavelength in the linear stability theory for the equilibrium state $u=1$ , but as we will see, this only remains so at moderately small values of $D$ , with the wavelength selection mechanism and value changing significantly for much smaller values of $D$ (see Section 3, but also Figure 5 and subsequent discussion below). Movies of the numerical solutions illustrated in Figures 1 and 4 can be found here.

The temporal periodicity of the creation of the stationary spatially periodic state behind the wavefront is also reflected in a weak periodic variation in the position of the wavefront (defined to be the largest value of $x$ at which $u=\frac{1}{2}$ ), but with an average speed of approximately $2\sqrt{D}$ , the minimum wavespeed, which is clearly illustrated in Figure 2 (see Sections 3 and 6, regarding the specific notion of minimum wavespeed in the present context). For the solutions shown in Figure 2 with $D=0.003$ and $D=0.002$ , the formation of a stationary periodic state behind the wavefront does not cause oscillations in the position of the wavefront, whereas for $D = 0.001$ , the magnitude of the oscillation is large enough that there is a weak effect. Although choosing a smaller value of $u$ at which to define the position of the wavefront would eliminate these traces of the oscillations, as $D \to 0$ , the value of $u$ required to do so becomes exponentially small. It, therefore, seems reasonable to use $u = \frac{1}{2}$ as our qualitative definition of the position of the wavefront, and thereby retain an indication of the fundamentally oscillatory nature of the observable solution in our record of the progress of the travelling wave.

Figure 2.

The numerically calculated position of the wavefront for various values of $D$ . The broken line has slope $2 \sqrt{D}$ , the minimum wavespeed.

For smaller values of $D$ than those used in Figures 1 and 2, we find that the creation of the humps, which ultimately form the stationary and spatially periodic state to the rear of the wavefront, initiates, periodically in $t$ , just ahead of the wavefront and, as we shall see below, can be related to the dynamics of that part of the solution profile that becomes exponentially small with distance ahead of the wavefront. It is now this periodic mechanism that selects the final spatial wavelength of the stationary periodic state which forms at the rear of the wavefront. In order to accurately compute this part of the solution, we use a different numerical method, solving instead for $W = \log u$ , which satisfies the evolution equation

(13) \begin{equation} \frac{\partial W}{\partial t} = D\frac{\partial ^2 W}{\partial x^2} + D \left (\frac{\partial W}{\partial x}\right )^2 + 1 - \phi * e^W. \end{equation}

We also use an FFT to calculate the convolution term (and periodic boundary conditions with periodicity large enough that the effect on the travelling wave dynamics is negligible) and use five-point stencils for the derivatives for greater accuracy. In addition, since we cannot take the logarithm of an initial condition with compact support, we instead use a Gaussian of width $w\gt 0$ as the initial condition, so that

(14) \begin{equation} g(x) = e^{-x^2/w^2}. \end{equation}

Figure 3 shows a snapshot of the evolving solution for $D = 10^{-8}$ , $A=0.01$ and $w = 0.1$ . Although the creation of a spatially periodic steady state behind the wavefront is qualitatively similar to that shown for larger values of $D$ in Figure 1, we can now see how the behaviour ahead of the wavefront leads to the creation of the humps, which have unit weight, width of $O(D^{1/2})$ and height of $O(D^{-1/2})$ for $D \ll 1$ (see Section 4.3 and Figure 6). Far ahead of the wavefront, the solution evolves as

(15) \begin{equation} u(x,t) \sim \frac{w A}{\sqrt{w^2+4Dt}} e^t\exp \left ({-} \frac{x^2}{w^2+4Dt}\right )\,\,\mbox{for $x^2 \gg w^2+4Dt \geq O(1)$.} \end{equation}

Note that (15) complies with the upper bound given by (10) with $t_0 = w^2/4D$ .

Figure 3.

The numerical solution of (IBVP) for gaussian initial data with $A=0.01$ and width $w = 0.04$ , and $D = 10^{-8}$ , when $t = 6000$ . The upper panel shows $\log _{10} u$ . New spikes are initiated ahead of the wave at the point where $u$ is close to $10^{-700}$ , which can only be captured accurately by solving for $\log u$ instead of $u$ .

In a spatial interval of unit width, centred on the leading hump, the term $1-\phi *u$ is small. However, ahead of this interval, it abruptly changes to one because of the finite support of the kernel. This has the effect of abruptly turning on the term $e^t$ in (15), which manifests itself in the solution shown in Figure 3 as a rapid change in $\log _{10}u$ in the region a distance $\frac{1}{2}$ ahead of the rightmost hump (at $x \approx 3.65$ , with the hump at $x \approx 3.15$ ), where $\log _{10} u$ grows until a new hump is formed. The dynamics of this mechanism can most easily be observed in the animation that can be found here. This hump formation process, generated by the dynamics of the exponentially small part of the solution ahead of the wavefront, is similar to that studied in [Reference Billingham4] for a nonlocal reaction-diffusion equation with a different local reaction term. For $D\ll 1$ , this hump forms when the logarithm of the far field solution (15) becomes small, and we can approximate this location as the point where the logarithm is zero, namely $x = x_f(t)$ , with

(16) \begin{equation} x_f(t) = \sqrt{w^2+4Dt} \sqrt{t + \log (Aw) - \frac{1}{2} \log (w^2+4Dt)}. \end{equation}

Figure 4 shows the position of the wavefront along with $x_f(t)$ for various values of $D$ . As can be seen, $x_f(t)$ gives an excellent prediction of the position of the wavefront, and we observe that it conforms with the wavefront bound given in (11). We also note that this hump formation process naturally leads to a stationary, spatially periodic solution behind the wavefront with wavelength close to $\frac{1}{2}$ . Figure 5 shows how the wavelength varies with $D$ , consistent with this observation. Also shown is the most unstable wavelength according to the linearised theory of Section 3, which is close to the wavelength observed for moderately small values of $D$ , when the periodic states emerges behind the wavefront, but not for smaller values of $D$ , where the hump formation mechanism described above takes over and creates this state, and determines its wavelength to be approaching $\frac{1}{2}$ as $D\to 0^+$ . Finally, it is anticipated that additional details relating to this hump formation mechanism may be made available by a more detailed development of the large- $t$ asymptotic form (15) via the method of matched asymptotic coordinate expansions (see Leach and Needham [Reference Leach and Needham20]) or the approximation methods detailed in van Saarloos [Reference van Saarloos26]. However, for the present paper, we regard the above details as sufficient to elucidate this particular mechanism. Many of the additional distinct features described above, particularly when the diffusivity is very small, have not been reported in earlier literature, and their illumination is the focus of the following sections.

Figure 4.

The numerically calculated position of the wavefront for $D = 10^{-4}$ , $10^{-5}$ , $10^{-6}$ , $10^{-7}$ , $10^{-8}$ and $10^{-9}$ , with $w = 0.1$ and $A=0.01$ . This position is defined as the largest value of $x$ at which $u = \frac{1}{2}$ and, because the solution propagates through the formation of discrete spikes, is not a continuous function of time, $t$ , and takes the form of a step function. The broken line is the function $x_f(t)$ , defined in (16).

Figure 5.

The wavelength of the spatially periodic steady state left behind the wavefront, calculated numerically as a function of $D$ . The broken line shows the most unstable wavelength given by the linearised theory.

Figure 6.

The inverse of the height of the spikes behind the wavefront, calculated numerically as a function of $D$ . The broken line is 100/ $\sqrt{D}$ .

3. Equilibrium states, stability characteristics and linear evolution

The nonlocal PDE (5), with (9), has two equilibrium states. The unreacted state with

(17) \begin{align} u(x,t)=0\quad \forall \,{(x,t)}\in \overline{D}_{\infty }, \end{align}

and the fully reacted state

(18) \begin{align} u{(x,t)} =1\quad \forall \,{(x,t)}\in \overline{D}_{\infty }. \end{align}

We have seen in Section 2 that these equilibrium states play a key role in the large- $t$ structure of the solution to (IBVP). Of particular significance is the temporal stability of these equilibrium states. In this section, we examine the linearised temporal stability of each of these equilibrium states in detail and consider the consequences thereof in terms of spatio/temporal evolution. To this end, we formulate a linearised initial value problem. We write,

(19) \begin{align} u(x, t)=u_{e}+\delta{\overline{u}}(x, t),\quad (x, t) \in \overline{D}_{\infty }, \end{align}

with $\delta \ll 1$ and $u_e=1$ when considering the fully reacted state or $u_e=0$ when considering the unreacted state. On substituting from (19) into (5), and neglecting terms of $O(\delta ^2)$ as $\delta \to 0$ , we obtain a linear evolution equation for $\overline{u}$ , namely

(20) \begin{align} {{\overline{u}}_t}=D{\overline{u}}_{xx} +L({\overline{u}}),{\quad \text{on }{D}_{\infty }}, \end{align}

where

(21) \begin{align} L({\overline{u}})={\begin{cases}{\overline{u}};& u_e=0\\ -\int _{x-\frac{1}{2}}^{x+\frac{1}{2}}{\overline{u}}(y,t)dy;& u_e=1 \end{cases}}, \end{align}

after using (9). The linearised initial value problem is then composed of (20), with associated initial and far field conditions,

(22) \begin{align} &{\overline{u}} (x,0)=\overline{g}(x)\quad \forall \,x\in{\mathbb{R}}, \end{align}

(23) \begin{align} &{\overline{u}}{(x,t)}\to 0 \text{ as } |x|\to \infty \text{ uniformly on }\overline{D}_T\ (\text{each}\ T\gt 0). \end{align}

Here $\overline{g}\in C^1({\mathbb{R}}) \cap L^\infty ({\mathbb{R}})$ and non-negative, $\|{\overline{g}}\|_{\infty }=1$ whilst $\textrm{supp}(\overline{g})\subseteq [{-}x_0,x_0]$ $(x_0\gt 0)$ . This problem will be referred to as $\text{(LIVP)}_{0}$ when $u_e=0$ , and $\text{(LIVP)}_{1}$ when $u_e=1$ . We remark that aspects of the following linearised theory have been discussed in [Reference Banerjee, Kuznetsov, Udovenko and Volpert2] in relation to equilibrium state temporal stability. Here we develop this further in terms of detailing the full qualitative structure associated with the linearised evolution problems, which provides additional and significant information towards the analysis of (IBVP). We now consider $\text{(LIVP)}_{0}$ and $\text{(LIVP)}_{1}$ in turn.

3.1 Analysis of (LIVP) $_0$

We first seek elementary solutions to (20) and (21) in the form,

(24) \begin{align} {\overline{u}}{(x,t)}=e^{ikx-w t},\quad{(x,t)}\in{\overline{D}}_{\infty }, \end{align}

with $k\in{\mathbb{R}}$ and $w\in{\mathbb{C}}$ . On substitution from (24) in (20) and (21), we obtain the dispersion relation

(25) \begin{align} w=w_0(k)=Dk^2-1\quad \forall \,k\in{\mathbb{R}}, \end{align}

and so, as expected, $\text{(LIVP)}_{0}$ is nondispersive ( $w_0(k)\in{\mathbb{R}}$ $\forall$ $k\in{\mathbb{R}}$ ), and moreover, for each $D\gt 0$ ,

(26) \begin{align} w_0(k)\lt 0\quad \forall \,k\in \left (\frac{-1}{\sqrt{D}},\frac{1}{\sqrt{D}}\right ), \end{align}

with

(27) \begin{align} w_0(0)=-1=\inf _{k\in{\mathbb{R}}}w_0(k). \end{align}

We immediately conclude that the equilibrium state $u_e=0$ is temporally unstable, according to the linearised theory, at each $D\gt 0$ (this conclusion can readily extended to apply to the fully nonlinear and nonlocal PDE (5) with (9), via an application of the parabolic comparison theorem to the operator $N(w)\,:\!=w_t-Dw_{xx}-w$ on ${D}_T$ ; for brevity we omit the details). The Fourier integral theorem allows us to write down the solution to $\text{(LIVP)}_{0}$ , and this can then be estimated, via steepest descents, to obtain,

(28) \begin{align} {\overline{u}}{(x,t)} \sim \frac{\sqrt{\pi }}{\sqrt{D}t^{\frac{1}{2}}}\widehat{g} (0)e^te^{-x^2/4Dt}, \end{align}

as $t\to \infty$ uniformly for $x\in{\mathbb{R}}$ , with

\begin{align*} \widehat{g}(0)=\frac{1}{2\pi }\int _{-x_0}^{x_0}\overline{g}(s)ds \quad (\gt 0). \end{align*}

We observe from (28) that, as $t\to \infty$ , there are two symmetric ‘wavefronts’ where

(29) \begin{align} |x|\sim 2\sqrt{D}t\text{ as } t\to \infty, \end{align}

and behind the wavefronts $\overline{u}$ is growing exponentially in $t$ , whilst ahead of the wavefronts $\overline{u}$ is decaying exponentially in $t$ . This structure is consistent with (11) and the numerical solutions to (IBVP) in Section 2. It is also consistent with the rigorous result in [Reference Bouin, Henderson and Ryzhik6] for the corresponding evolution problem when the initial data is front-like and localised to the left half-line.

3.2 Analysis of (LIVP) $_1$

We seek elementary solutions to (20) and (21) in the form of (24), with again $k\in{\mathbb{R}}$ and $w\in{\mathbb{C}}$ . This now leads directly to the dispersion relation

(30) \begin{align} w=w_1(k)=Dk^2+\frac{2}{k}\sin \frac{1}{2}k\quad \forall \,k\in{\mathbb{R}} \end{align}

and we observe that $w_1(k)$ is an even function of $k$ . Moreover, $w_1(k)\in{\mathbb{R}}$ for all $k\in{\mathbb{R}}$ and so $\text{(LIVP)}_{1}$ is nondispersive. In further analysing (30), it is convenient to introduce the function $\Delta\; :\;{\mathbb{R}}^+\to{\mathbb{R}}$ such that

(31) \begin{align} \Delta (X) =-\frac{2}{X^3}\sin \frac{1}{2}X\quad \forall \,X\in{\mathbb{R}}^+. \end{align}

The zeros of $\Delta (X)$ are at

(32) \begin{align} X=2n\pi, \quad n\in \mathbb{N}, \end{align}

whilst the turning points are at

(33) \begin{align} X=\delta _n,\quad n\in \mathbb{N}, \end{align}

with $2n\pi \lt \delta _n\lt 2(n+1)\pi$ , and

(34) \begin{align} \delta _n\sim (2n+1)\pi \text{ as } n\to \infty . \end{align}

Furthermore, $\delta _n$ is a local maximum when $n$ is odd and a local minimum when $n$ is even. At each local maximum point, we write,

(35) \begin{align} \Delta _r =\Delta (\delta _{2r-1}),\quad r=1,2,\dots . \end{align}

We observe that,

(36) \begin{align} 0\lt \Delta _{r+1}\lt \Delta _r,\quad r=1,2,\dots \end{align}

and

(37) \begin{align} \Delta _r\sim \frac{2}{\pi ^3 (4r-1)^3}\text{ as }r\to \infty, \end{align}

whilst a straightforward numerical calculation gives ${\Delta }_1\approx 0.00297$ .

We can now readily interpret the dispersion relation (30). We first observe that for

(38) \begin{align} D\gt \Delta _1, \end{align}

then $w_1(k)\gt 0$ for all $k\in{\mathbb{R}}$ and so the equilibrium state $u_e=1$ is temporally asymptotically stable. However, for

(39) \begin{align} 0\lt D\lt \Delta _1, \end{align}

then

(40) \begin{align} \inf _{k\geq 0}w_1(k)=w_1(k_m)=k_m^2(D-\Delta (k_m))\lt 0, \end{align}

and so the equilibrium state $u_e=1$ is now temporally unstable. We note that $k=k_m$ is uniquely determined as the smallest positive root of the transcendental equation,

(41) \begin{align} w^{\prime}_1(k) \equiv 2k(D-\Delta (k)) - k^2\Delta '(k) \equiv \left(2Dk^3-2\sin \frac{1}{2}k +k\cos \frac{1}{2}k\right)k^{-2}=0, \end{align}

and it is straightforward to establish from this that,

(42) \begin{align} k_m\to \begin{cases} \delta _1\text{ as }&D\to \Delta _1^-\\ k_0\text{ as } &D\to 0^+ \end{cases}, \end{align}

where $k=k_0$ is the smallest positive root of the equation

(43) \begin{align} \tan \frac{1}{2}k=\frac{1}{2}k, \end{align}

so that $2\pi \lt k_0\lt 3\pi$ , and

(44) \begin{align} w_1(k_m)\to \frac{2}{k_0}\sin \frac{1}{2}k_0\text{ as }D\to 0^+. \end{align}

We note that in both cases,

\begin{align*} w_1(k)\sim Dk^2\text{ as }|k|\to \infty, \end{align*}

and so $\text{(LIVP)}_{1}$ is well-posed. Figure 5 shows the wavelength of the most unstable mode as a function of $D$ .

The solution to $\text{(LIVP)}_{1}$ is readily obtained as

(45) \begin{align} {\overline{u}}{(x,t)}=\int _{-\infty }^\infty \widehat{g}(k)e^{-w_1(k)t}e^{ikx}dk, \end{align}

for ${(x,t)}\in \overline{D}_{\infty }$ , with $\widehat{g}$ being the Fourier transform of $\overline{g}$ . We obtain from (45), via Laplace’s method, that,

(46) \begin{align} {\overline{u}}{(x,t)}\sim \frac{\sqrt{\pi }2^{3/2}|\widehat{g}(k_m)|}{(w^{\prime\prime}_1(k_m))^{\frac{1}{2}}t^{\frac{1}{2}}}e^{-w_1(k_m)t} \cos{(k_m x+\arg ( \widehat{g}(k_m))}, \end{align}

for $|x|=O(1)$ as $t\to \infty$ . It should be noted that further regions in the large- $t$ structure of $\overline{u}$ are required when $|x|=O(t)$ as $t\to \infty$ , which gives the transition into the spatially exponentially decaying far field for $|x|\gg O(t)$ . We observe, from (46), that when $0\lt D\lt \Delta _1$ , the solution to $\text{(LIVP)}_{1}$ evolves into an exponentially growing harmonic periodic state with spatial wave number $k_m$ , which depends upon $D$ , and temporal exponential growth rate $|w_1(k_m)|$ . This periodic state is stationary, and evolves when $|x|=O(1)$ as $t\to \infty$ , behind transition into the far field when $|x|=O(t)$ as $t\to \infty$ .

In relation to (IBVP), the above analyses indicate that when $t$ is large, two symmetric permanent form wavefronts propagate to left and right, with asymptotic propagation speeds of $\pm 2\sqrt{D}$ . However, in the region to the rear of the two wavefronts at $|x|\sim 2\sqrt{D}t$ , the large- $t$ structure to $\text{(LIVP)}_{1}$ indicates that, as $t\to \infty$ in (IBVP), and specifically when $|x|=O(1)$ as $t\to \infty$ , there are two possible steady spatial structures which emerge as a consequence of wavefront passage, depending upon $D$ . With $u\;:\;\overline{D}_{\infty }\to{\mathbb{R}}$ being the solution to (IBVP), these possibilities are

1. (P1) when $D\gt \Delta _1$ , then $u{(x,t)}\to 1$ as $t\to \infty$ , uniformly with $|x|=O(1)$

2. (P2) when $0\lt D\lt \Delta _1$ , then $u{(x,t)}\to P(x)$ as $t\to \infty$ , uniformly with $|x|=O(1)$ . Here $P\;:\;{\mathbb{R}}\to{\mathbb{R}}$ is a steady periodic solution to equation (5) with (9), which is positive, has $1\in \textrm{Im} (P)$ and has wavelength which (i) is close to $2\pi /k_m(D)$ when $D$ is close to $\Delta _1$ , and is selected well behind the emerging wavefronts by a maximum linear growth rate mechanism on the now weakly linearly unstable equilibrium state $u=1$ and (ii) moves away from this value as $D$ decreases, decreasing and drifting towards the minimum wavelength of $1/2$ (see Figure 5 and Section 4) as $D\to 0^+$ , with this wavelength now selected via the hump formation mechanism, which operates to the immediate rear of the wavefronts (see Section 2).

We recall that both (P1) and (P2) are supported by the numerical solutions to (IBVP) presented in Section 2 in conjunction with the detailed linearised theory of this section.

As a consequence of the theory in this section, and in particular to further enable the development of theory to support and extend the conjecture (P2), the next natural key step is to investigate the existence and structural nature of the positive periodic steady solutions to the nonlocal equation (5) with (9), which are conjectured to emerge in the large- $t$ development of the solution to (IBVP). The starting point for this study begins by investigating the emergence of periodic steady solutions via steady state bifurcations from the equilibrium solution $u_e=1$ . Thereafter, the detailed structure of these fully nonlinear and nonlocal periodic solutions is developed in detail for significant asymptotic limits. In particular, a fully developed theory is obtained in the small $D$ limit, when the nonlocal length scale is much larger than the diffusion length scale, leading to spatially periodic states consisting of separated, periodically distributed, localised humps, characterised by nonlocal effects which are regulated by weak diffusion.

4. Positive periodic steady states

We begin by considering the steady state form of the nonlocal equation (5) with (9), and in particular, we seek to identify steady state bifurcations from the equilibrium state $u_e=1$ , which give rise to periodic steady states. For the general form of the nonlocal Fisher-KPP equation, under quite general conditions on the kernel, generic considerations concerning local bifurcations to small amplitude bounded steady states on the real line, as the associated equilibrium state loses its temporal stability, have been developed in [Reference Faye and Holzer11] and a detailed consideration of the spectral theory on the real line of the linear operators thereby encountered has been made by Volpert and Vougalter [Reference Volpert, Vougalter, Lewis, Maini and Petrovsky29]. By restricting attention to bifurcations to periodic steady states, we effectively restore compactness to the operators and can then develop a standard weakly nonlinear bifurcation theory. Indeed, under the periodic restriction, these operator theoretic results conform with the basic weakly nonlinear theory outlined below for the current specific situation, which we then develop in detail onto the global bifurcation branch. Any steady state to the nonlocal equation (5) with (9), in the present context, is a function $F\;:\;{\mathbb{R}}\to{\mathbb{R}}$ , with $F\in C^{2}({\mathbb{R}})\cap L^\infty ({\mathbb{R}})$ , and such that,

(47) \begin{align} DF_{xx}+F\left (1-\int _{x-\frac{1}{2}}^{x+\frac{1}{2}}F(y)dy\right )=0,\quad x\in{\mathbb{R}}. \end{align}

We restrict attention to considering the existence of positive periodic steady states, which oscillate about the equilibrium state $u_e=1$ . Specifically, we will fix $D\gt 0$ , and consider the bifurcation to periodic steady states from $u_e=1$ , with fundamental wavelength $\lambda$ , as the bifurcation parameter. As a preliminary, we first observe, via bootstrapping in (47), that when $F\in C^2({\mathbb{R}}) \cap L^\infty ({\mathbb{R}})$ is a steady state, then in fact, $F\in C^\infty ({\mathbb{R}})$ . Moreover (see for example [Reference Hale16]), this can be improved to

(48) \begin{align} F\in C^\omega ({\mathbb{R}}) \cap L^\infty ({\mathbb{R}}). \end{align}

Now, let $F=F_p(x,\lambda, D)$ be a positive, periodic, steady state at diffusivity $D$ and with fundamental wavelength $\lambda \gt 0$ . We define

(49) \begin{align} \alpha \equiv \max _{x\in [0,\lambda )}F_p(x,\lambda, D)-\min _{x\in [0,\lambda )}F_p(x,\lambda, D). \end{align}

the peak-to-trough magnitude of this periodic steady state. In general, we anticipate that $\alpha =\alpha (\lambda, D)$ . For fixed $D\gt 0$ , we now suppose that a bifurcation to periodic steady states occurs from $u_e=1$ as $\lambda$ passes through $\lambda =\lambda _b$ $(\gt 0)$ . To determine the possible values of $\lambda _b$ , we consider those values $\lambda =\lambda _b$ when the linearised form of (47) has solution

(50) \begin{align} F(x)=1+\alpha \cos \left(\frac{2\pi }{\lambda _b}x+\phi \right),\quad x\in{\mathbb{R}}, \end{align}

as $\alpha \to 0^+$ , with $\phi$ being a constant phase. On substitution from (50) into (47) and allowing $\alpha \to 0^+$ , a non-trivial solution requires that $\lambda _b$ satisfies the transcendental equation, with fixed $D\gt 0$ ,

(51) \begin{equation} D =-\frac{\lambda _b^3}{4\pi ^3}\sin \frac{\pi }{\lambda _b} =\Delta \left (\frac{2\pi }{\lambda _b}\right ), \end{equation}

with $\Delta (\cdot )$ as introduced in (31). We can now interpret the root structure of (51), using (32)–(38). In particular, for each $r=1,2,\dots$ , then, with

(52) \begin{align} \Delta _{r+1}\lt D\lt \Delta _r, \end{align}

equation (51) has exactly $2r$ positive roots, which we label as $\lambda _i^\pm (D)$ for $i=1,2,\dots, r$ , with

(53) \begin{align} \frac{1}{2i}\lt \lambda _i^-(D)\lt \frac{2\pi }{\delta _{2i-1}}\lt \lambda _i^+(D)\lt \frac{1}{(2i-1)}, \end{align}

for each $i=1,\dots, r$ . We note that, for fixed $k\in \mathbb{N}$ , then $\lambda ^\pm _k(D)$ exist and are continuous for $D\in [0,\Delta _k)$ . Moreover, the following properties are readily established:

(54) \begin{align} & \lambda _k^\pm (D)\to \frac{2\pi }{\delta _{2k-1}}\text{ as } D\to \Delta _k^-, \end{align}

(55) \begin{align} & \lambda _k^{+({-})}(D)\text{ is monotone decreasing (increasing) with } D\in [0,\Delta _k), \end{align}

(56) \begin{align} &\lambda _k^+(D)\to \frac{1}{(2k-1)}^-\text{ as }D\to 0^+, \end{align}

(57) \begin{align} &\lambda _k^-(D)\to \frac{1}{2k}^+\text{ as }D\to 0^+. \end{align}

Now, fix $r\in \mathbb{N}$ , and fix $\Delta _{r+1}\lt D\lt \Delta _r$ . At this fixed $D$ , there are thus $2r$ local bifurcation points to periodic steady states, namely at critical wavelengths

(58) \begin{align} \lambda =\lambda _k^\pm (D),\quad k=1,2,\dots, r. \end{align}

Further lengthy, but straightforward calculations (standard weakly nonlinear asymptotics as the bifurcation point is approached) establish (which can be made rigorous in a standard way within the classical local bifurcaton theory, as developed in, for example, Kuznetsov [Reference Kuznetsov19]) that each bifurcation is a steady state pitchfork bifurcation from $u_e=1$ , being subcritical at $\lambda =\lambda _k^+(D)$ and supercritical at $\lambda =\lambda _k^-(D)$ $(k=1,2,\dots, r)$ . Additional investigation (numerical continuation of the weakly nonlinear bifurcation branches) reveals that, for each $k$ , the two bifurcated branches at $\lambda _k^\pm (D)$ are connected in the $(\lambda, \alpha )$ bifurcation diagram. We can represent this curve as $\alpha =\alpha (\lambda ;D)$ for $\lambda _k^-(D)\leq \lambda \leq \lambda _k^+(D)$ (for each $k=1,2,\dots, r$ ). Close to the bifurcation points $\lambda =\lambda _k^\pm (D)$ a very standard weakly nonlinear local bifurcation theory is readily developed (details are omitted for brevity as they follow a now standard routine) to establish that, as should be expected for a nondegenerate bifurcation,

(59) \begin{align} \alpha (\lambda, D)=C_k^\pm (D)\left |{\lambda -\lambda _k^\pm (D)}\right |^{\frac{1}{2}}+O(\left |{\lambda -\lambda _k^\pm (D)}\right |), \end{align}

as $\lambda \to (\lambda _k^+(D))^-$ or as $\lambda \to (\lambda _k^-(D))^+$ , respectively. Here, $C_k^\pm (D)$ is a positive constant. The bifurcated periodic steady states have the form

(60) \begin{align} u=F_p(x,\lambda, D)=1+\alpha (\lambda, D)\cos \left(\frac{2\pi }{\lambda }x+\phi \right)+O(\alpha ^2(\lambda, D)), \end{align}

as $\lambda \to \lambda _k^\pm (D)$ respectively, with $\alpha (\lambda, D)$ as in (59), and $\phi \in [0,2\pi )$ being an arbitrary phase (representing translational invariance of (47) in $x$ ). In addition, numerical investigation confirms that, with $D$ fixed, then $\alpha (\lambda, D)$ has a single stationary point for $\lambda \in [\lambda _k^-(D),\lambda _k^+(D)]$ , which is a maximum point.

Figure 7.

The $(\lambda, \alpha )$ bifurcation diagram with $D = 3\times 10^{-6}$ , $10^{-5}$ , $3\times 10^{-5}$ , $10^{-4}$ , $3\times 10^{-4}$ , $10^{-3}$ , $2 \times 10^{-3}$ . The amplitude, $\alpha$ , increases as $D$ decreases.

The numerically calculated $(\lambda, \alpha )$ bifurcation diagram is given in Figure 7 for various values of $D$ . It is now straightforward to construct the bifurcation locus in the $(\lambda, D)$ plane. This is given directly by the intersection of the curve

\begin{align*} D=-\frac{\lambda ^3}{4\pi ^3}\sin \frac{\pi }{\lambda }, \end{align*}

with the positive quadrant of the $(\lambda, D)$ plane. This intersection is in the form of a countably infinite sequence of ‘tongue like’ curves based on the $\lambda -$ axis, and with monotone decreasing ‘height’. The $i^{th}$ ‘tongue’ has base points at $\lambda =1/(2i-1)$ and $\lambda =1/2i$ , with ‘height’ $D=\Delta _i$ . We label the set of interior points of the $i^{th}$ ‘tongue’ as $\Omega _i$ , and its boundary as $\partial \Omega _i=\overline{\Omega _i}\setminus \Omega _i$ . We write

(61) \begin{align} \Omega =\bigcup _{i=1}^\infty \Omega _i,\quad \partial \Omega =\bigcup _{i=1}^\infty \partial \Omega _i. \end{align}

At each point ${(\lambda, D)}\in \Omega$ , there is a unique (up to translation in $x$ ) and nontrivial periodic steady state, with fundamental wavelength $\lambda$ and amplitude $\alpha (\lambda, D)$ . Since the equation (47) is invariant under the transformation $x\mapsto -x$ , this allows for the existence of a translation in $x$ , so that at each ${(\lambda, D)}\in \Omega$ we may select a representative periodic steady state $u=F_p(x,\lambda, D)$ which is an even function of $x\in{\mathbb{R}}$ , so that

(62) \begin{align} &F_p({-}x,\lambda, D)=F_p(x,\lambda, D)\quad \forall x\in{\mathbb{R}}, \end{align}

(63) \begin{align} &F^{\prime}_p(0,\lambda, D)=0, \end{align}

and

(64) \begin{align} F^{\prime}_p\left (\pm \frac{1}{2}\lambda, \lambda, D\right )=0. \end{align}

A plot of $\Omega$ in the $(\lambda, D)$ plane is shown in Figure 8. The numerically calculated maximum value of $u=F_p(x,\lambda, D)$ , which we denote by $u_{max}$ , is shown in Figure 9 for the first four tongues of $\Omega$ . This is a more suitable measure of amplitude than $\alpha$ for the purposes of this semi-logarithmic plot, since $\alpha$ is zero at the edges of the region where the nonlinear solutions exist.

Figure 8.

The region $\Omega$ lies below these curves or tongues.

Figure 9.

Semilogarithmic plots of $u_{max}$ in the first four tongues of $\Omega$ .

We now examine the structure of the periodic steady states with ${(\lambda, D)}\in \overline{\Omega }$ . Firstly, on the bifurcation locus, with ${(\lambda, D)}\in \partial \Omega \cap ({\mathbb{R}}^+)^2$ we have

(65) \begin{align} F_p(x,\lambda, D)=1\quad \forall x\in{\mathbb{R}}, \end{align}

whilst it follows from (59) and (60) that

(66) \begin{align} F_p(x,\lambda, D)=1+O(|\lambda -\lambda _0|^{\frac{1}{2}}+|D-D_0|^{\frac{1}{2}})\text{ as }(\lambda, D)\to (\lambda _0,D_0), \end{align}

uniformly for $x\in{\mathbb{R}}$ , with $(\lambda _0,D_0)\in \partial \Omega \cap ({\mathbb{R}}^+)^2$ and ${(\lambda, D)}\in \Omega$ . In terms of regularity, it follows from (48) and (47) (see for example, [Reference Hale16]) that

(67) \begin{align} F_p\in C^{\omega, 0,0}({\mathbb{R}}\times (\overline{\Omega }\cap ({\mathbb{R}}^+)^2)). \end{align}

We can now establish the following bounds:

1. (B1) For any ${(\lambda, D)}\in \Omega$ , then,

   \begin{equation*}\inf _{x\in {\mathbb {R}}}F_p(x,\lambda, D)\gt 0.\end{equation*}

Proof. Set $(\lambda ^*,D^*)\in \Omega$ , then $(\lambda ^*,D^*)\in \Omega _i$ for some $i\in \mathbb{N}$ . The point $(\lambda ^*,D^*)$ can then be connected, in $\Omega _i$ , to a point $(\lambda _0,D_0)\in \partial \Omega _i\cap ({\mathbb{R}}^+)^2$ via a straight line segment $L$ . Now, suppose that

(68) \begin{align} \inf _{x\in{\mathbb{R}}}F_p(x,\lambda ^*,D^*)\leq 0. \end{align}

We recall from (65) that

(69) \begin{align} \inf _{x\in{\mathbb{R}}}F_p(x,\lambda _0,D_0)=1. \end{align}

Thus, from (67), it follows that there exists $(\lambda _1,D_1)\in L\cap \Omega _i$ , such that

(70) \begin{align} \inf _{x\in{\mathbb{R}}} F_p(x,\lambda _1,D_1)=0. \end{align}

As a consequence (since $F_p$ is periodic in $x\in{\mathbb{R}}$ , together with (67)) there exists $x_1\in{\mathbb{R}}$ such that

(71) \begin{align} F_p(x_1,\lambda _1,D_1)=0. \end{align}

It then follows, via (70), (71), (67) and (47), that

(72) \begin{align} F_p^{(n)}(x_1,\lambda _1,D_1)=0\quad \forall \,n\in \mathbb{N}. \end{align}

It is then a consequence of (72) with (67) that,

\begin{align*} F_p(x,\lambda _1,D_1)=0\quad \forall \,x\in{\mathbb{R}}. \end{align*}

However, $(\lambda _1,D_1)\in \Omega _i$ , and we have established earlier in this section that as such $F_p(x,\lambda _1,D_1)$ must be a non-trivial periodic function of $x\in{\mathbb{R}}$ , with fundamental wavelength $\lambda _1$ , and so cannot be the trivial function. Thus we arrive at a contradiction. The result follows.

We next have

1. (B2) For any ${(\lambda, D)}\in \Omega$ , then,

   \begin{equation*}\sup _{x\in {\mathbb {R}}}F_p(x,\lambda, D)\gt 1\end{equation*}

Proof. Set ${(\lambda, D)}\in \Omega$ , then ${(\lambda, D)}\in \Omega _i$ for some $i\in \mathbb{N}$ . Suppose that

(73) \begin{align} \sup _{x\in{\mathbb{R}}}F_p(x,\lambda, D)\leq 1. \end{align}

Then since $F_p(x,\lambda, D)$ is a non-trivial periodic function in $x\in{\mathbb{R}}$ , with fundamental wavelength $\lambda$ , it follows that there exists $x^*\in{\mathbb{R}}$ such that

(74) \begin{align} 0\lt F_p(x^*,\lambda, D)=\inf _{x\in{\mathbb{R}}}F_p(x,\lambda, D)\lt 1, \end{align}

with the left-hand inequality following from (B1). From (73) and (74), it follows that

(75) \begin{align} \int _{x^*-\frac{1}{2}}^{x^*+\frac{1}{2}}F_p(y,\lambda, D) dy\lt 1. \end{align}

Also, from (47), we have

(76) \begin{align} F^{\prime\prime}_p(x^*,\lambda, D)=-F_p(x^*,\lambda, D)\left (1-\int _{x^*-\frac{1}{2}}^{x^*+\frac{1}{2}}F_p(y,\lambda, D) dy\right )\lt 0, \end{align}

via (74) and (75). However, via (74) and (67), $x=x^*$ is a local minimum point for $F_p(x,\lambda, D)$ , and so

\begin{equation*}F^{\prime\prime}_p(x^*,\lambda, D)\geq 0,\end{equation*}

which contradicts (76). Thus,

\begin{equation*}\sup _{x\in {\mathbb {R}}}F_p(x,\lambda, D)\gt 1,\end{equation*}

as required.

Correspondingly, we now have

1. (B3) For any ${(\lambda, D)}\in \Omega$ , then,

   \begin{equation*}\inf _{x\in {\mathbb {R}}}F_p(x,\lambda, D)\lt 1.\end{equation*}

To complete this section, we next consider in detail the structure of the periodic steady states in the $(\lambda, D)$ plane as $D\to 0^+$ . We begin by examining this limit for the periodic steady states in the principal tongue $\Omega _1$ and then give the corresponding results for $\Omega _i$ , $i=2,3,4,\dots$ .

4.1 Asymptotic structure as $D\to 0^+$ with $(\lambda, D)\in \Omega _1$

We consider the asymptotic structure of $F_p(x,\lambda, D)$ with $(\lambda, D)\in \Omega _1$ as $D\to 0^+$ . This is achieved via a detailed application of the method of matched asymptotic expansions. Due to the evenness of $F_p$ in $x$ , we need only consider this structure for $x\in \left [0,\frac{1}{2}\lambda \right ]$ , with $\lambda \in \left (\frac{1}{2},1\right )$ fixed as $D\to 0^+$ . We recall that, $F_p(x,\lambda, D)\gt 0$ for all $x\in \left [0,\frac{1}{2}\lambda \right ]$ , whilst evenness requires

(77) \begin{align} F_{p}^{\prime }(0, \lambda, D)=F_{p}^{\prime }\left (\frac{1}{2} \lambda, \lambda, D\right )=0. \end{align}

An examination of (47), together with a consideration of numerical solutions, dictates that, for $\lambda \in \left (\frac{1}{2},1\right )$ fixed, the asymptotic structure of $F_p(x,\lambda, D)$ , as $D\to 0^+$ , develops into a three region structure as follows:

Region I (support region)

\begin{align*} x\in [0,S-O(D^\gamma )), \,\, F_p=O(1)^+\text{ as }D\to 0^+, \end{align*}

with $0\lt S\lt \frac{1}{2}\lambda$ and $\gamma \gt 0$ to be determined. In general, $S$ will depend upon both $\lambda$ and $D$ , with $S=O(1)^+$ as $D\to 0^+$ , and $\lambda \in \left (\frac{1}{2},1\right )$ . With this in mind, we expand $S(\lambda, D)$ as

(78) \begin{align} S(\lambda, D)=\bar{S} (\lambda ) +D^\gamma S_1(\lambda )+ o(D^\gamma ),\text{ as } D\to 0^+, \end{align}

with $\lambda \in \left (\frac{1}{2},1\right )$ , which allows for weak displacements in the location of region II below.

Region II (transition region)

\begin{align*} x\in \left (S-O(D^\gamma ),S+O(D^\gamma )\right ),\,\, F_p=O(D^\delta )^+\text{ as }D\to 0^+, \end{align*}

with $\delta \gt 0$ to be determined.

Region III (exponential region)

\begin{align*} x\in \Big(S+O(D^\gamma ),\frac{1}{2}\lambda \Big],\,\, F_p=O(E(D))^+\text{ as } D\to 0^+, \end{align*}

with $E(D)$ indicating terms exponentially small in $D$ as $D\to 0^+$ .

The above structure requires, for the change in structure across region II, that the leading order interval length of region III (which is the leading order separation of consecutive support regions) must have the same length as the half span of the nonlocal term, which is $\frac{1}{2}$ , and so,

\begin{align*} \bar{S}+\frac{1}{2}=\lambda -\bar{S}, \end{align*}

which gives,

(79) \begin{align} \bar S =\frac{1}{2}\left (\lambda -\frac{1}{2}\right )=\bar S(\lambda ), \end{align}

and we note from this that

\begin{align*} 0\lt \bar S(\lambda )\lt \frac{1}{4}, \end{align*}

for $\lambda \in \left (\frac{1}{2},1\right )$ , with,

(80) \begin{align} \bar{S}(\lambda )-\frac{1}{2}\lt -\bar{S}(\lambda )\quad \forall \, \lambda \in \left (\frac{1}{2},1\right ). \end{align}

With (78)–(80), then equation (47) requires that

(81) \begin{align} \int _{x-\frac{1}{2}}^{x+\frac{1}{2}} F_p(y,\lambda, D)dy=\int _{-\frac{1}{2}}^{\frac{1}{2}} F_p(y,\lambda, D)dy+ O(E(D))\text{ as } D\to 0^+, \end{align}

for $x\in \left [0,\bar{S}(\lambda )-O(D^\gamma )\right )$ , whilst

(82) \begin{align} \int _{-\frac{1}{2}}^{\frac{1}{2}} F_p(y,\lambda, D)dy=1+\alpha _1 D+\beta _1 D^r +o(D^r), \end{align}

as $D\to 0^+$ , with the constants $r\gt 1$ , $\alpha _1$ and $\beta _1$ to be determined.

We begin in region I and expand in the form

(83) \begin{align} F_p(x,\lambda, D)=F_0(x,\lambda )+D^m F_1(x,\lambda )+o(D^m)\text{ as } D\to 0^+, \end{align}

with $x\in [0,S(\lambda )-O(D^\gamma ))$ and $m\gt 0$ to be determined. We substitute into (47), to obtain, at leading order as $D\to 0^+$ ,

(84) \begin{align} &F^{\prime\prime}_0-\alpha _1 F_0=0,\quad 0\lt x\lt \bar S(\lambda ), \end{align}

(85) \begin{align} &F^{\prime}_0(0,\lambda )=F_0(\bar S(\lambda ),\lambda )=0, \end{align}

(86) \begin{align} &F_0(x,\lambda )\gt 0,\quad 0\leq x\lt \bar S(\lambda ), \end{align}

with the second condition in (85) required by asymptotic matching (via Van Dyke’s principle, [Reference Van Dyke25]) between region I and region II. Finally, (82) requires

(87) \begin{align} \int _0^{\bar S(\lambda )}F_0(y,\lambda ) dy=\frac{1}{2}, \end{align}

using the evenness of $F_p(x,\lambda, D)$ , and its asymptotic form in regions II and III. Now, the problem (84)–(86) is a classical Sturm–Liouville eigenvalue problem, with eigenvalue $-\alpha _1$ , whilst condition (86) requires that $-\alpha _1$ is the lowest eigenvalue. This may be treated directly, to obtain,

(88) \begin{align} \alpha _1=-\frac{\pi ^2}{\left (\lambda -\frac{1}{2}\right )^2}, \end{align}

with

(89) \begin{align} F_0(x,\lambda )=B \cos \frac{\pi x}{\left (\lambda -\frac{1}{2}\right )},\quad 0\leq x\leq \bar S(\lambda ), \end{align}

and $B\gt 0$ an arbitrary constant. It remains to apply condition (87), which gives

(90) \begin{align} B=\frac{\pi }{(2\lambda -1)}, \end{align}

and so,

(91) \begin{align} F_0(x,\lambda )=\frac{\pi }{(2\lambda -1)}\cos \frac{\pi x}{\left (\lambda -\frac{1}{2}\right )},\quad 0\leq x\leq \bar S(\lambda ), \end{align}

and $\bar{S}(\lambda )=\frac{1}{2}\left (\lambda -\frac{1}{2}\right )$ . Via (81), (82) and (88), we also have

(92) \begin{align} \int _{x-\frac{1}{2}}^{x+\frac{1}{2}} F_p(y,\lambda, D)dy=2\int _0^{\frac{1}{2}}F_p(y,\lambda, D)dy=1-\frac{\pi ^2}{\left (\lambda -\frac{1}{2}\right )^2}D+\beta _1 D^r+o(D^r), \end{align}

as $D\to 0^+$ , with $x\in [0,S(\lambda )-O(D^\gamma ))$ . This completes the leading order structure in region I. We now move on to consider region II, the transition region. In region II we introduce the scaled coordinate $X=O(1)$ as $D\to 0^+$ , given by,

(93) \begin{align} x=S(\lambda, D)+D^\gamma X, \end{align}

It then follows from the expansion in region I, on moving into region II, via (91), (83) with (87), that $F_p=O(D^\gamma )$ as $D\to 0^+$ in region II. Thus, we expand in the form

(94) \begin{align} F_p(X,\lambda, D)=D^\gamma \bar F_0(X,\lambda )+o(D^\gamma )\text{ as } D\to 0^+, \end{align}

with $X=O(1)$ . We remark here that the term at $O(D^\gamma )$ in $S(\lambda, D)$ could be removed by a constant shift in the coordinate $X$ . However, it is convenient for later matching purposes to leave this shift in $S$ ; this also leads to a parameter free equation at leading order, with all parameters shifted into the matching condition. We next consider the form of the nonlocal term in (47) in region II, which can be written as,

(95) \begin{align} \notag &\int _{x-\frac{1}{2}}^{x+\frac{1}{2}} F_p(y,\lambda, D)dy =\int _{S(\lambda, D)-\frac{1}{2}+D^\gamma X}^{S(\lambda, D)+\frac{1}{2}+D^\gamma X} F_p(y,\lambda, D)dy \nonumber \\ &=\int _{S(\lambda, D)-\frac{1}{2}+D^\gamma X}^{S(\lambda, D)-\frac{1}{2}} F_p(y,\lambda, D)dy+\int _{S(\lambda, D)-\frac{1}{2}}^{S(\lambda, D)+\frac{1}{2}-O(1)} F_p(y,\lambda, D)dy+\int _{S(\lambda, D)+\frac{1}{2}-O(1)}^{S(\lambda, D)+\frac{1}{2}+D^\gamma X} F_p(y,\lambda, D)dy. \end{align}

Now, using the periodicity and evenness of $F_p$ , and its structure in each of regions I-III, we have firstly that

(96) \begin{align} \int _{S(\lambda, D)-\frac{1}{2}+D^\gamma X}^{S(\lambda, D)-\frac{1}{2}} F_p(y,\lambda, D)dy=O(E(D)D^\gamma ), \end{align}

as $D\to 0^+$ with $X=O(1)$ , whilst,

(97) \begin{align} \int _{S(\lambda, D)-\frac{1}{2}}^{S(\lambda, D)+\frac{1}{2}-O(1)} F_p(y,\lambda, D)dy=1-\frac{\pi ^2}{\left (\lambda -\frac{1}{2}\right )^2}D+\beta _1 D^r+o(D^r), \end{align}

as $D\to 0^+$ . The final term in (95) may now be written as

(98) \begin{align} &\int _{S(\lambda, D)+\frac{1}{2}-O(1)}^{S(\lambda, D)+\frac{1}{2}+D^\gamma X} F_p(y,\lambda, D)dy = D^\gamma \int _{Y=-O(D^{-\gamma })}^{Y=X} F_p\left (S(\lambda, D)+\frac{1}{2}+D^\gamma Y,\lambda, D\right )dY \notag \\&= D^\gamma \int _{Y=-O(D^{-\gamma })}^{Y=X} F_p(\lambda -S(\lambda, D)+D^\gamma Y,\lambda, D)dY \notag \\&= D^\gamma \int _{Y=-O(D^{-\gamma })}^{Y=X} F_p({-}S(\lambda, D)+D^\gamma Y,\lambda, D)dY\notag \\ &= D^\gamma \int _{Y=-O(D^{-\gamma })}^{Y=X} F_p(S(\lambda, D)-D^\gamma Y,\lambda, D)dY \sim D^{2\gamma }\int _{Y=-\infty }^{Y=X} \bar F_0({-}Y,\lambda )dY \notag \\&=D^{2\gamma }\int _{z=-X}^{z=\infty } \bar F_0(z,\lambda )dz, \end{align}

as $D\to 0^+$ with $X=O(1)$ . Using (96)–(98) in (95), we finally arrive at

(99) \begin{equation} \int _{S(\lambda, D)-\frac {1}{2}+D^\gamma X}^{S(\lambda, D)+\frac {1}{2}+D^\gamma X} F_p(y,\lambda, D)dy =1-\frac{\pi ^2}{\left (\lambda -\frac{1}{2}\right )^2}D+\beta _1 D^r+D^{2\gamma }\int _{z=-X}^{z=\infty } \bar F_0(z,\lambda )dz+o(D^{2\gamma },D^r,E(D)), \end{equation}

as $D\to 0^+$ with $X=O(1)$ . We now substitute from (93), (94) and (99) into equation (47). To obtain a non-trivial balance at leading order requires us to choose

(100) \begin{align} \gamma =\frac{1}{4}. \end{align}

At leading order, equation (47) then becomes,

(101) \begin{align} \Psi _{\bar X \bar X}-\Psi \int _{w=-\bar X}^{w=\infty }\Psi (w)dw=0,\quad \bar X\in{\mathbb{R}}, \end{align}

which is both nonlinear and nonlocal. In (101), we have introduced the simple scalings

(102) \begin{align} \bar F_0=(2\lambda -1)^{-\frac{3}{2}}\Psi, \,\,\,X =(2\lambda -1)^{\frac{1}{2}}\bar X. \end{align}

Equation (101) is now completed by matching conditions between region II and region I (as $\bar{X}\to -\infty$ ) and region III (as $\bar{X}\to \infty$ ). The matching process is straightforward (following Van Dyke’s asymptotic matching principle, [Reference Van Dyke25]), and leads to the two boundary conditions,

(103) \begin{align} &\Psi (\bar X)\to 0\text{ as } \bar X\to \infty, \end{align}

(104) \begin{align} &\Psi (\bar X)\sim -2\pi ^2(\bar X+l)\text{ as } \bar X\to -\infty . \end{align}

Here

(105) \begin{align} S_1(\lambda )=(2\lambda -1)^{\frac{1}{2}} l. \end{align}

It is also required that

(106) \begin{align} \Psi (\bar X)\gt 0\quad \forall \,\bar X\in{\mathbb{R}}. \end{align}

We observe that the problem (101), (103), (104) and (106) is a nonlinear, nonlocal eigenvalue problem, with real eigenvalue $l$ . It is also free of parameters. A numerical investigation of this problem establishes that a solution exists if and only if $l=l^*$ , and the solution is unique. Here

(107) \begin{align} l^*\approx -\frac{3.493}{2\pi ^2}, \end{align}

and a numerically determined graph of $\Psi (\bar{X})$ against $\bar{X}$ , calculated using a simple finite difference method and the trapezium rule, is shown in Figure 10. It is also straightforward to establish that

(108) \begin{align} \Psi (\bar X)=-2\pi ^2(\bar X+l^*)+\Psi _{-\infty }\bar X^{-2}e^{-\frac{1}{2}\pi \bar{X}^2}(1+o(1)), \end{align}

as $\bar{X}\to -\infty$ , whilst,

(109) \begin{align} \Psi (\bar X)=\Psi _{\infty } e^{-\frac{1}{2}\pi \bar{X}^2}(1+o(1)), \end{align}

as $\bar{X} \to +\infty$ . Here, $\Psi _{\infty }$ ( $\gt 0$ ) and $\Psi _{-\infty }$ are parameter free, globally determined constants. In principle, these can be determined from the numerical solution shown in Figure 10, but in practice, since these corrections are exponentially small, we are not able to resolve them with sufficient accuracy.

Figure 10.

The numerically calculated solution of (101) subject to (102). The broken line is $\Psi = -2 \pi ^2 (\bar{X}+l^*)$ .

We now return to region I. To allow matching with region II, we now require

(110) \begin{align} m=\frac{1}{4}, \end{align}

whilst the most structured balance in equation (47) dictates that

(111) \begin{align} r=m+1=\frac{5}{4}. \end{align}

The problem for $F_1$ is then,

(112) \begin{align} F^{\prime\prime}_{1}+\frac{\pi ^{2}}{\left (\lambda -\frac{1}{2}\right )^{2}} F_{1}&=\beta _1 F_{0}(x, \lambda )\notag \\ &=\frac{\beta _1\pi }{(2\lambda -1)}\cos \frac{\pi x}{\left (\lambda -\frac{1}{2}\right )},\quad 0\lt x\lt \frac{1}{2}\left (\lambda -\frac{1}{2}\right ), \end{align}

(113) \begin{align} &F^{\prime}_1 (0,\lambda )=F_1\left (\frac{1}{2}\left (\lambda -\frac{1}{2}\right ),0\right )=0, \end{align}

(114) \begin{align} &\int _0^{\frac{1}{2}\left (\lambda -\frac{1}{2}\right )} F_1(y,\lambda )dy=0. \end{align}

The second of conditions (113) arises from matching with region II (with both expansions taken up to $O(D^{\frac{1}{4}})$ in regions I and II). A solution to (112) and (113) exists if and only if

(115) \begin{align} \beta _1=0, \end{align}

after which

(116) \begin{align} F_1(x,\lambda )=C\cos \frac{\pi x}{\left (\lambda -\frac{1}{2}\right )}, \end{align}

with $C$ an arbitrary real constant. However, substitution from (116) into (114), then requires $C=0$ , and so,

(117) \begin{align} F_1(x,\lambda )=0\quad \forall \,x\in \left [0,\frac{1}{2}\left (\lambda -\frac{1}{2}\right )\right ]. \end{align}

Thus, in region I we have

(118) \begin{align} F_p(x,\lambda, D)=\frac{\pi }{(2\lambda -1)}\cos \frac{\pi x}{\left (\lambda -\frac{1}{2}\right )}+o(D^{\frac{1}{4}}), \end{align}

as $D\to 0^+$ with $x\in \left [0,\frac{1}{2}\left (\lambda -\frac{1}{2}\right )\right )$ .

We finally move into the exponential region, region III. In this region, $F_p$ is exponentially small in $D$ as $D\to 0^+$ and develops as a WKB expansion. For brevity we omit details, but obtain

(119) \begin{align} F_p(x,\lambda, D)=A^*(\lambda )e^{-\Phi _0(\lambda )D^{-\frac{1}{2}}}\cosh \left(\frac{\Phi (x)}{D^\frac{1}{2}}\right)(1+o(1)), \end{align}

as $D\to 0^+$ , with $x\in \left (\frac{1}{2}(\lambda -\frac{1}{2})+O(D^{\frac{1}{4}}),\frac{1}{2}\lambda \right ]$ . Here,

(120) \begin{align} \Phi (x)=\frac{1}{\sqrt 2}\int _x^{\frac{1}{2}\lambda }\left (1-\sin \left(\frac{\pi w}{\left (\lambda -\frac{1}{2}\right )}\right)\right )^{\frac{1}{2}}dw\quad \forall \,x\in \left [\frac{1}{2}\left (\lambda -\frac{1}{2}\right ),\frac{1}{2}\lambda \right ), \end{align}

and

(121) \begin{align} \Phi _0(\lambda )=\Phi \left (\frac{1}{2}\left (\lambda -\frac{1}{2}\right )\right ). \end{align}

The constant $A^*(\lambda )$ is independent of $D$ and is related to the constant $\Psi _{\infty }$ and $\lambda$ , via matching expansion (119) as $x\to \frac{1}{2}\left (\lambda -\frac{1}{2}\right )^-$ with expansion (94) as $X\to \infty$ , which gives

(122) \begin{align} A^*(\lambda )=\frac{2\Psi _{\infty }}{(2\lambda -1)^{\frac{3}{2}}}, \end{align}

with details omitted for brevity. This completes the asymptotic structure as $D\to 0^+$ , with fixed $\lambda \in \left (\frac{1}{2},1\right )$ .

We make the following key observations. First, up to terms exponentially small in D, as $D\to 0^+$ ,

(123) \begin{align} \textrm{supp}(F_p(x,\lambda, D))=[{-}S(\lambda, D),S(\lambda, D)], \end{align}

with

(124) \begin{align} S(\lambda, D)=\frac{1}{2}\left (\lambda -\frac{1}{2}\right )-\frac{3.493}{\sqrt 2\pi ^2}\left (\lambda -\frac{1}{2}\right )^{\frac{1}{2}}D^{\frac{1}{4}}+o(D^{\frac{1}{4}})\,\text{ as } D\to 0^+, \end{align}

whilst,

(125) \begin{align} \alpha (\lambda, D)=\frac{\pi }{(2\lambda -1)}+o(D^{\frac{1}{4}})\,\text{ as } D\to 0^+, \end{align}

and occurs at $x=0$ . In addition,

(126) \begin{align} \int _{-S(\lambda, D)}^{+S(\lambda, D)} F_p(y,\lambda, D)dy=1-\frac{\pi ^2}{\left (\lambda -\frac{1}{2}\right )^2}D+o(D^{\frac{5}{4}}), \end{align}

as $D\to 0^+$ . For comparison, we present in Figure 11, a numerically determined representation of $F_p$ , with $\lambda =0.75$ and various values of $D$ . This is in excellent agreement with the asymptotic form developed above.

Figure 11.

The solution $F_p$ , calculated numerically for $\lambda = \frac{3}{4}$ and $D = 10^{-3}$ , $10^{-4}$ , $10^{-5}$ , $10^{-6}$ , $10^{-7}$ and $10^{-8}$ . The broken line shows the leading order outer solution given by (118).

Finally, we observe from (123)to (126) in particular, that the asymptotic structure developed above, for $F_p(x,\lambda, D)$ as $D\to 0^+$ , with fixed $\lambda \in \left (\frac{1}{2},1\right )$ , becomes nonuniform as $\lambda \to \frac{1}{2}$ , and, specifically, when $\lambda =\frac{1}{2}+O(D^{\frac{1}{2}})$ as $D\to 0^+$ . For the present, we continue by considering the asymptotic structure of $F_p(x,\lambda, D)$ as $D\to 0^+$ in each of the remaining tongues $\Omega _i$ , $i=2,3,\dots$ .

4.2 Asymptotic structure as $D\to 0^+$ with $(\lambda, D)\in \Omega _i$ , $i=2,3,\dots$

We now consider the asymptotic form of $F_p(x,\lambda, D)$ with ${(\lambda, D)}\in \Omega _i$ ( $i\geq 2$ ) as $D\to 0^+$ . We, again, need only consider the structure for $x\in \left [0,\frac{1}{2}\lambda \right ]$ , with now $\lambda \in \left (\frac{1}{2i},\frac{1}{2i-1}\right )$ fixed, as $D\to 0^+$ . The structure is again as in the case of $\Omega _1$ , and for brevity, we restrict attention to the key support region, labelled as region I. The analysis follows exactly that of the previous subsection, and we therefore present only the salient features. In this case, we obtain,

(127) \begin{align} &\bar S (\lambda )= \frac{1}{2}i\left (\lambda -\frac{1}{2i}\right ), \end{align}

(128) \begin{align} & F_0(x,\lambda )=\frac{\pi }{(2i-1)(2i\lambda -1)}\cos \frac{\pi x}{i\left (\lambda -\frac{1}{2i}\right )}\quad \forall \,x\in [0,\bar S(\lambda )], \end{align}

(129) \begin{align} & (2i-1)\int _{-S(\lambda, D)}^{S(\lambda, D)}F_p(y,\lambda, D)dy=1-\frac{\pi ^2}{i^2\left (\lambda -\frac{1}{2i}\right )^2}D+o(D), \end{align}

as $D\to 0^+$ , whilst,

(130) \begin{align} \alpha (\lambda, D)=\frac{\pi }{(2i-1)(2i\lambda -1)}+o(1)\text{ as } D\to 0^+, \end{align}

with $\lambda \in \left (\frac{1}{2i},\frac{1}{2i-1}\right )$ fixed. We observe that, as before, this structure becomes nonuniform when $\lambda =\frac{1}{2i}+o(1)$ as $D\to 0^+$ .

4.3 Asymptotic structure as $D\to 0^+$ with $\lambda =\frac{1}{2} +O(D^{\frac{1}{2}})$

We return to $(\lambda, D)\in \Omega _1$ , and now develop the structure to $F_p(x,\lambda, D)$ with $\lambda =\frac{1}{2} +O(D^{\frac{1}{2}})$ and $x\in \left [0,\frac{1}{2}\lambda \right ]$ as $D\to 0^+$ . We write

(131) \begin{align} \lambda =\frac{1}{2}+D^{\frac{1}{2}}\bar \lambda, \end{align}

with $\bar{\lambda }=O(1)^+$ as $D\to 0^+$ . An examination of (123)–(125), with (131), then requires,

(132) \begin{align} x=O(D^{\frac{1}{2}})^+ \text{ and } F_p=O(D^{-\frac{1}{2}})^+, \end{align}

as $D\to 0^+$ in the support region, with

(133) \begin{align} F_p=O(E(D)), \end{align}

when

(134) \begin{align} x\in (O(D^{\frac{1}{2}}),\frac{1}{4}+\frac{1}{2}\bar{\lambda } D^{\frac{1}{2}}], \end{align}

in the exponential region. Following (132), in the support region we introduce

(135) \begin{align} \hat{X}=\frac{x}{D^{\frac{1}{2}}}=O(1)^+, \end{align}

as $D\to 0^+$ , and expand in the form

(136) \begin{align} F_p(\hat{X},{\bar{\lambda }},D)=v(\hat{X},{\bar{\lambda }})D^{-\frac{1}{2}}+o(D^{-\frac{1}{2}}), \end{align}

as $D\to 0^+$ with $\hat{X}=O(1)^+$ . After a careful consideration of the nonlocal term, substitution from (131), (135) and (136) into (47) gives, at leading order, the nonlocal, nonlinear equation,

(137) \begin{align} v_{\hat{X} \hat{X}} +v\left (1-2\int _{-\infty }^\infty v(s,{\bar{\lambda }})ds +\int _{\hat{X}-{\bar{\lambda }}}^{\hat{X}+{\bar{\lambda }}}v(s,{\bar{\lambda }})ds \right )=0;\quad \hat{X}\gt 0, \end{align}

which must be solved subject to the conditions,

(138) \begin{align} &v_{\hat{X}}(0,{\bar{\lambda }})=0, \end{align}

(139) \begin{align} &v(\hat{X},{\bar{\lambda }})\to 0\text{ as } \hat{X}\to \infty, \end{align}

(140) \begin{align} &v(\hat{X},{\bar{\lambda }})\gt 0 \quad \forall \, x\geq 0, \end{align}

(141) \begin{align} & v({-}\hat{X},{\bar{\lambda }})=v(\hat{X},{\bar{\lambda }}) \quad \forall \hat{X}\geq 0. \end{align}

Here, the boundary condition (139) accommodates asymptotic matching to the exponential region.

Figure 12.

The numerical solution (bold curves) of (137) to (141) for $\bar{\lambda } = 0.5$ , $1$ , $10$ , $50$ and $100$ . The broken line is the asymptotic solution for $\bar{\lambda }\ll 1$ for ${\bar{\lambda }} = 0.5$ and $1$ , given by (156), whilst the dash-dotted line is the asymptotic solution for ${\bar{\lambda }} \gg 1$ for ${\bar{\lambda }} = 50$ and $100$ , which comes from rescaling (118).

Figure 13.

The maximum value, $v(0,\bar{\lambda })$ of the numerically calculated solution for (137) to (141). The broken lines show the predicted asymptotic behaviour as $\bar{\lambda } \to 0$ and $\bar{\lambda } \to \infty$ .

The problem (137)–(141) has been considered numerically, which has demonstrated that a unique solution exists for each ${\bar{\lambda }}\gt 0$ . We observe that this solution has

(142) \begin{align} v(\hat{X},{\bar{\lambda }})\sim v_{\infty } ({\bar{\lambda }}) e^{-\sigma _{\infty }({\bar{\lambda }})\hat{X}}\text{ as } \hat{X} \to \infty, \end{align}

with

(143) \begin{align} \sigma _{\infty }({\bar{\lambda }})=4\int _0^\infty v(s,{\bar{\lambda }})ds-1\gt 0, \end{align}

and $v_{\infty }({\bar{\lambda }})\gt 0$ a globally dependent constant. Figure 12 shows the numerical solution of (137) to (141) for various values of $\bar{\lambda }$ . The solution becomes larger and narrower as $\bar{\lambda }$ increases, but then becomes smaller and fatter as $\bar{\lambda }$ increases past about five (see also Figure 13). Also shown is a comparison with the asymptotic solutions for $\bar{\lambda }$ both large and small, with which there is excellent agreement. Figure 14 shows the numerically calculated periodic solution of the full problem for $D = 10^{-3}$ and $\lambda = \frac{1}{2} + 10 \sqrt{D}$ , along with the corresponding asymptotic solution for $\bar{\lambda } = 10$ and the solution (118). The agreement is excellent, and we can also see that the $\lambda = O(1)$ asymptotic solution, (118), is not a good approximation to the solution if $\lambda - \frac{1}{2}$ is of $O(\sqrt{D})$ , as expected. The maximum of $v$ is at $\hat{X}=0$ and a graph of $v(0,{\bar{\lambda }})$ against $\bar{\lambda }$ is shown in Figure 13, which shows that there is a single maximum at ${\bar{\lambda }} \approx 5.8$ , with, as required,

(144) \begin{align} v(0,{\bar{\lambda }})&\sim \frac{\pi }{2{\bar{\lambda }}}\text{ as }{\bar{\lambda }}\to \infty, \nonumber \\v(0,{\bar{\lambda }})&\to 0^+\text{ as }{\bar{\lambda }}\to 0^+ \end{align}

We note that

(145) \begin{align} \alpha ({\bar{\lambda }},D)=D^{-\frac{1}{2}}v(0,{\bar{\lambda }})+o(D^{-\frac{1}{2}}), \end{align}

as $D\to 0^+$ . The amplitude given by $\alpha ({\bar{\lambda }},D)$ is indistinguishable from the amplitude given by the numerical solution of the full problem, shown in Figure 7, for $D$ less than about $10^{-4}$ .

Figure 14.

The numerically calculated periodic solution of the full problem for $D = 10^{-3}$ and $\lambda = \frac{1}{2} + 10 \sqrt{D}$ , along with the corresponding asymptotic solution for $\bar{\lambda } = 10$ (broken line) and the solution (118), (dash-dotted line).

To complete this subsection, we examine the problem (137)–(141) in the limit ${\bar{\lambda }}\to 0^+$ . A balancing of terms in (137) leads us to write,

(146) \begin{align} \tilde{X}={\bar{\lambda }}\hat{X}=O(1)^+\text{ as }{\bar{\lambda }}\to 0^+, \end{align}

and expand in the form

(147) \begin{align} v(\tilde{X},{\bar{\lambda }})={\bar{\lambda }} \tilde{v}(\tilde{X})+o({\bar{\lambda }}) \text{ as }{\bar{\lambda }}\to 0^+. \end{align}

Using (146) and (147), the nonlocal term in (137) becomes

(148) \begin{align} \int _{\hat{X}-{\bar{\lambda }}}^{\hat{X}+{\bar{\lambda }}}v(s,{\bar{\lambda }})ds =\int _{\hat{X}-{\bar{\lambda }}}^{\hat{X}+{\bar{\lambda }}}({\bar{\lambda }} \tilde{v}({\bar{\lambda }} s)+o({\bar{\lambda }}))ds =\int _{{\bar{\lambda }}\hat{X}-{\bar{\lambda }}^2}^{{\bar{\lambda }}\hat{X}+{\bar{\lambda }}^2}(\tilde{v}(w)+o(1))dw =2{\bar{\lambda }}^2\tilde{v}(\tilde{X})+o({\bar{\lambda }}^2), \end{align}

as ${\bar{\lambda }}\to 0^+$ . The most structured balance at leading order then requires,

(149) \begin{align} \int _{-\infty }^\infty v(\tilde{X},{\bar{\lambda }})d\tilde{X}=\frac{1}{2}{\bar{\lambda }}+\tilde{I}{\bar{\lambda }}^3+o({\bar{\lambda }}^3), \end{align}

as ${\bar{\lambda }}\to 0^+$ with $\tilde{I}$ a constant to be determined. On substitution from (146)–(149) into (137)–(141), we arrive at the leading order problem,

(150) \begin{align} &\tilde{v}_{\tilde{X}\tilde{X}}-2\tilde{v}(\tilde{I}-\tilde{v})=0,\, \tilde{X}\gt 0, \end{align}

(151) \begin{align} &\tilde{v}_{\tilde{X}}(0)=0, \end{align}

(152) \begin{align} &\tilde{v}(\tilde{X})\to 0\text{ as } \tilde{X}\to \infty, \end{align}

(153) \begin{align} &\tilde{v}(\tilde{X})\gt 0\quad \tilde{X}\geq 0, \end{align}

(154) \begin{align} &\int _0^\infty \tilde{v}(w)dw=\frac{1}{4}. \end{align}

This problem is a nonlinear eigenvalue problem, with eigenvalue $\tilde{I}\in{\mathbb{R}}$ , which can be solved directly. There is a single eigenvalue given by

(155) \begin{align} \tilde{I}=\frac{1}{72}, \end{align}

with the associated eigenfunction uniquely determined as

(156) \begin{align} \tilde{v}(\tilde{X})=\frac{1}{48}\text{sech}^2\left (\frac{1}{12}\tilde{X}\right )\quad \forall \,\tilde{X}\geq 0. \end{align}

Thus, (147) and (156) give

(157) \begin{align} v(0,{\bar{\lambda }})\sim \frac{1}{48}{\bar{\lambda }}\,\text{ as }{\bar{\lambda }}\to 0^+, \end{align}

which is shown in Figure 13, and (154) and (155) lead to

(158) \begin{align} \int _{-\infty }^\infty v(\hat{X},{\bar{\lambda }})d\hat{X}=\frac{1}{2}+\frac{1}{72}{\bar{\lambda }}^2 +o({\bar{\lambda }}^2), \end{align}

as ${\bar{\lambda }}\to 0^+$ . Equations (156)–(158) are in excellent agreement with the numerical solution to problem (137)–(141), when $\bar{\lambda }$ is small

A final remark in this section relates to the temporal stability of the periodic steady states identified with ${(\lambda, D)}\in \Omega$ . In this respect, we have performed a limited stability analysis of $u=F_p(x,\lambda, D)$ with ${(\lambda, D)}\in \Omega$ when $D$ is small by considering the exponentially small part of the solution. Without presenting details, this analysis reveals that each such periodic steady state is locally, temporally, asymptotically stable to perturbations in the region between the $O(1)$ parts of the solution, where we would expect any instability to occur (extended intervals where $u$ remains uniformly very small, when long enough, can destabilise as they inherit the linear instability characteristics of the zero equilibrium). In relation to (IBVP), we observe that the results of this section support the possibility (P2), for the case when $0\lt D\lt \Delta _1$ .

We readily verify that the wavelength selected by (IBVP) has, when $D$ is close to $\Delta _1$ ,

(159) \begin{align} \lambda \approx \lambda _m(D)=\frac{2\pi }{k_m(D)}\in \Omega _1, \end{align}

and so we have, in (P2),

(160) \begin{align} P(x) = F_p(x,\lambda _m(D),D),\, x\in \mathbb{R}. \end{align}

However, this wavelength decreases as $D$ decreases, approaching the limiting value of $\frac{1}{2}$ as $D\to 0^+$ , and so now in (P2),

(161) \begin{align} P(x) = F_p(x,\lambda (D),D),\, x\in \mathbb{R}, \end{align}

where now $\lambda (D)\to \frac{1}{2}^+$ as $D\to 0^+$ , which therefore has the structure given in detail in Section 4.3. In particular, the periodic solution selected therefore is of the nature of a spike, with base width of $O(D^{\frac{1}{2}})$ and height of $O(D^{-\frac{1}{2}})$ . This is in full accord with the numerical solution of Section 2. A final point to note is that the detailed analysis of the existence and structure of the family of steady periodic solutions in the first tongue $\Omega _1$ on the $(\lambda, D)$ plane given in this section not only provides substance to the conjecture (P2) in the present context, but also plays a crucial role in studying the corresponding Cauchy problem on a closed finite spatial interval (with either Dirichlet or Neumann endpoint conditions), which we present in the second of this series of papers, [Reference Needham and Billingham23].

We now move on to consider bifurcations to periodic travelling waves from the equilibrium solution $u_e=1$ . This enables us to investigate in more depth the structure of the solution to (IBVP) when $t$ is large, now in the region of the propagating wavefronts identified in Section 3. In particular, the next section first rules out the possibility of nondecaying, permanent form periodic waves travelling with, and to the rear of, the principal wavefronts.

5. Positive periodic travelling waves

We seek to identify the bifurcation to spatially periodic travelling wave solutions from the equilibrium state $u_e=1$ . At fixed $D\gt 0$ , a periodic travelling wave solution, with propagation speed $v=v_b\neq 0$ will bifurcate from $u_e=1$ , at wavelength $\lambda _b\gt 0$ if and only if, $(\lambda _b,v_b)$ satisfies the complex equation,

(162) \begin{align} 4 \pi ^{2} D-2 i \pi v_{b} \lambda _{b}+\frac{\lambda _{b}^{3}}{\pi } \sin \frac{\pi }{\lambda _{b}} = 0 \end{align}

which is readily obtained from equation (5) with (9), when written in the travelling wave coordinate $z=x-v t$ , and linearised about $u_e=1$ . It is immediate that (162) has no solutions which have $v_b\neq 0$ and $\lambda _b\gt 0$ . We can therefore conclude that there are no bifurcations to periodic travelling waves from the equilibrium state $u_e=1$ .

We are now able to use this result in moving on to consider transitional permanent form travelling waves, which accommodate the transition from the equilibrium state $u_e=0$ ahead of the wavefront to the equilibrium state $u_e=1$ or to a spatially periodic permanent form travelling wave, to the rear of the wavefront.

6. Transitional travelling waves

A transitional travelling wave (or wavefront), in the present context, is defined as: a travelling wave of permanent form, propagating with constant propagation speed $v\gt 0$ . The permanent wave form is nonnegative, and ahead of the wavefront the wave form approaches the equilibrium state $u_e=0$ , whilst to the rear of the wavefront, the wave form either approaches the equilibrium state $u_e=1$ or approaches a non-trivial, positive periodic travelling wave (with the same propagation speed $v\gt 0$ ). We anticipate that such structures will play a key role in the large- $t$ evolution of the solution to (IBVP), when $|x|=O(t)$ as $t\to \infty$ . This motivates the present consideration of transitional travelling waves. Following Section 5, we can immediately eliminate the existence of transitional travelling waves of the second type described above, and we therefore concentrate on the existence of transitional travelling waves of the first type, which we will henceforth refer to as (TPTW) solutions to equation (5) with (9). For a different treatment of permanent form transitional travelling waves to (1), that lead up to the concept of generalised travelling waves, the reader is referred to [Reference Apreutesei, Bessonov, Volpert and Vougalter1] and [Reference Banerjee, Kuznetsov, Udovenko and Volpert2] together with the references therein. In addition, in this context, numerical computations of transitional travelling waves have been considered in [Reference Nadin, Perthame and Tang22], relating particularly to such structures connecting two unstable steady states.

We introduce the travelling coordinate $z=x-vt$ with $v\gt 0$ being the constant propagation speed. We represent a (TPTW) solution with propagation speed $v\gt 0$ as $u=u_T(z,v)$ , with $u_T(\cdot, v)\in C^2({\mathbb{R}})$ satisfying,

(163) \begin{align} & D u^{\prime\prime}_T+v u^{\prime}_T+u_T\left (1-\int _{z-\frac{1}{2}}^{z+\frac{1}{2}}u_T(s,v)ds\right )=0,\quad z\in{\mathbb{R}}, \end{align}

(164) \begin{align} &u_T(z,v)\geq 0 \quad \forall \, z\in{\mathbb{R}}, \end{align}

(165) \begin{align} &u_{T}(z, v) \rightarrow \left \{\begin{array}{lll} 1 & \text{ as } & z \rightarrow -\infty \\ 0 & \text{ as } & z \rightarrow +\infty . \end{array}\right . \end{align}

The problem (163)–(165) will be henceforth refered as (BVP). We begin with some qualitative results.

1. (T1) Let $u_T(\cdot, v)\;:\;{\mathbb{R}}\to{\mathbb{R}}$ be a (TPTW), then $u_T(\cdot, v)\in C^\omega ({\mathbb{R}})$ .

1. (T2) Let $u_T(\cdot, v)\;:\;{\mathbb{R}}\to{\mathbb{R}}$ be a (TPTW), then $u_T(z,v)\gt 0$ for all $z\in{\mathbb{R}}$ .

1. (T3) The existence of a (TPTW) requires $v\geq 2\sqrt D$ .

Proof. Let $u_T(\cdot, v)\;:\;{\mathbb{R}}\to{\mathbb{R}}$ be a (TPTW) with propagation speed $v\gt 0$ . Via (163)–(165) and the linearisation theorem (see for example [Reference Hale16]), we must have

(166) \begin{align} &D u_{T}^{\prime \prime }+v u_{T}^{\prime }+u_{T}=0, \quad z\gg 1, \end{align}

(167) \begin{align} &u_T(z,v)\to 0\,\text{ as } z\to \infty, \end{align}

(168) \begin{align} &u_T(z,v)\gt 0 \,\text{for} \, z\gg 1, \end{align}

with the strict inequality in (168) following from (T2). As a consequence of (166) and (167), there exists constants $A_{\infty }$ and $B_{\infty }$ , not both zero, such that

(169) \begin{align} u_T(z,v)\sim A_{\infty } e^{\lambda _+(v)z}+B_{\infty } e^{\lambda _-(v)z}, \end{align}

with

(170) \begin{align} \lambda _{\pm }(v)=-\frac{1}{2 D}\left (v \pm \left (v^{2}-4 D\right )^{1 / 2}\right ), \end{align}

and the obvious modification when $v = 2\sqrt{D}.$ It then follows from (169) and (170) with (168) that the propagation speed $v\geq 2 \sqrt{D}$ , as required.

In the remainder of this section, we suppose that $u_T\;:\;{\mathbb{R}}\to{\mathbb{R}}$ is a (TPTW), which has the minimum possible propagation speed

(171) \begin{align} v=2\sqrt{D}. \end{align}

which we refer to as $u_T(z)$ . Our intention is now to examine in detail the form of $u_T(z)$ as $z\to \ -\infty$ , using the linearisation theorem (see, for example, [Reference Hale16]). We remark here that Fang and Zhao [Reference Fang and Zhao10] establish that there exists a critical value of $D$ such that $u_T(z)$ is monotone if and only if $D$ exceeds this critical value. We extend this here by obtaining details of this critical value and establishing the existence of a second, and lower, critical value of $D$ , for which the behaviour of $u_T(z)$ as $z\to -\infty$ , although not globally monotone, is ultimately monotone exponential decay for $D$ above this second critical value, but becomes harmonically oscillating exponential decay for $D$ below this critical value. First we write,

(172) \begin{align} u_T(z)=1+{\overline{u}}(z), \end{align}

after which ${\overline{u}} (z)$ must satisfy

(173) \begin{align} &D{\overline{u}}''+2\sqrt{D}{\overline{u}}'-\int _{z-\frac{1}{2}}^{z+\frac{1}{2}}{\overline{u}}(s)ds=0,\quad ({-}z)\gg 1, \end{align}

(174) \begin{align} &{\overline{u}}(z)\to 0\, \text{ as } z\to -\infty . \end{align}

The form of ${\overline{u}}(z)$ as $z\to -\infty$ is determined by examining elementary solutions in the form

(175) \begin{align} {\overline{u}} (z)= e^{\sigma z}, \end{align}

with $\sigma \in{\mathbb{C}}$ to be determined, and $\mathrm{Re}(\sigma )\gt 0$ , to satisfy condition (174) (both the real and imaginary parts of (175) provide real-valued solutions to the real linear equation (173)). On substitution from (175) into (173), we find that $\sigma \in{\mathbb{C}}$ allows (175) to solve (173) and (174) if and only if

(176) \begin{align} D\sigma ^2 + 2 \sqrt{D} \sigma -\frac{2}{\sigma } \sinh{\frac{1}{2}\sigma }=0, \text{ with } \mathrm{Re}(\sigma )\gt 0. \end{align}

For each fixed $D\gt 0$ , the transcendental equation (176) has a countably infinite number of roots in the complex plane, which we label as

(177) \begin{align} \sigma _n (D)\in{\mathbb{C}} \text{ for each } n\in \mathbb{Z}\setminus \{0\}, \end{align}

and, in addition

(178) \begin{align} \sigma _n(\cdot )\in PC^1(\overline{{\mathbb{R}}}^+)\cap C(\overline{{\mathbb{R}}}^+). \end{align}

with,

(179) \begin{align} &\mathrm{Re}(\sigma _{\pm (2r-1)}(D))\gt 0,\, r\in \mathbb{N}, \nonumber \\&\mathrm{Re}(\sigma _{\pm 2r}(D))\lt 0,\, r\in \mathbb{N}. \end{align}

We next observe that

(180) \begin{align} \sigma _n(D) = 2n\pi i -({-}1)^n 8 n^2 \pi ^2 \sqrt{D}+ O(D)\,\text{ as } D\to 0^+ \end{align}

with $n\in \mathbb{Z}\setminus \{0\}$ . Also,

(181) \begin{align} &\sigma _{\pm (2r-1)}(D)=2 \log \frac{1}{2}D \pm 4 (r-1)\pi i +o(1), \nonumber \\&\sigma _{\pm 2r}(D)=-2 \log \frac{1}{2}D \pm 4 (r-1)\pi i +o(1), \end{align}

as $D\to \infty$ , with $r\in \mathbb{N}\setminus \{1\}$ . The full locus of $\sigma _n(D)$ is sketched in Figure 15 for various values of $n$ . We observe from (176) that, in fact,

(182) \begin{align} \sigma _{-n}(D)=\bar{\sigma }_n(D), \end{align}

for $n\in \mathbb{Z}\setminus \{\pm 1,\pm 2\}$ .

Figure 15.

The locus of the eigenvalues $\sigma _n(D)$ in the complex plane. Note that the imaginary part has been scaled by a factor of $2 \pi$ . These start from $\sigma _n = 2 \pi n i$ when $D=0$ (shown as circles) and move through the complex plane as $D$ increases. For $n = \pm 1$ , $\pm 2$ , $\sigma _n$ reaches the real axis at the points marked with a square, with $\sigma _1$ , $\sigma _2 \to 0$ as $D \to \infty$ , and $\sigma _{-1}\to \infty$ , $\sigma _{-2} \to - \infty$ as $D \to \infty$ . Both $\sigma _{-1}$ and $\sigma _{-2}$ are shown as broken lines.

We now consider the cases $n= \pm 1,\pm 2$ . We concentrate on $\sigma _{\pm 1}(D)$ . For $0\lt D\lt D_+$ , we have that $\sigma _{-1}(D)=\bar \sigma _{1}(D)$ and $\mathrm{Im}(\sigma _1(D))\gt 0$ . At $D=D_+$ , $\sigma _{-1}(D_+)=\sigma _1(D_+)=\sigma _+\in{\mathbb{R}}^+$ . For $D\gt D_+$ , $0\lt \sigma _1(D)\lt \sigma _+\lt \sigma _{-1}(D)$ , with

(183) \begin{align} &\sigma _1(D)\sim \frac{(\sqrt{2}-1)}{\sqrt{D}}, \nonumber \\&\sigma _{-1}(D)\sim 2\log \frac{1}{2}D, \end{align}

as $D\to \infty$ , as can be seen in Figure 15. In the above, $(\sigma _+,D_+)\in ({\mathbb{R}}^+)^2$ is the solution to the transcendental equations

(184) \begin{align} & D\sigma ^3 +2 \sqrt{D}\sigma ^2 -2 \sinh \frac{1}{2}\sigma =0, \nonumber \\& 3D \sigma ^2 + 4\sqrt{D} \sigma - \cosh \frac{1}{2} \sigma =0, \end{align}

which gives, via numerical approximation,

(185) \begin{align} (\sigma _+,D_+)\approx (4.437,\ 2.824 \times 10^{-2}). \end{align}

We are now able to interpret these results in relation to a (TPTW) with minimum propagation speed $v=2\sqrt D$ . In particular, translation invariance in $z$ can and henceforth will be assumed to be, fixed so that

(186) \begin{align} u_T(z)\sim z e^{-\frac{1}{\sqrt{D}}z}\,\text{ as } z\to \infty, \end{align}

via (169) and (170), whilst, for $0\lt D\lt D_+$ , there are global constants $A_{\infty }\neq 0$ and $\phi _{\infty }\in [0,2\pi )$ such that

(187) \begin{align} u_T(z)\sim 1-A_{\infty } e^{a(D)z }\cos (b(D)z+\phi _{\infty })\,\text{ as } z\to -\infty, \end{align}

with

(188) \begin{align} a(D)=\mathrm{Re}(\sigma _1(D)), \quad b(D)=\mathrm{Im}(\sigma _1(D)), \end{align}

which are both positive. In particular,

(189) \begin{align} & a(D) \sim \begin{cases} 8 \pi ^{2} \sqrt{D} &\text{ as } D \rightarrow 0^{+}, \\ \sigma _{+}-O\left (\left (D_{+}-D\right )\right ) &\text{ as } D \rightarrow D_{+}^{-} \end{cases} \end{align}

(190) \begin{align} & b(D) \sim \begin{cases} 2 \pi +O(D) &\text{ as } D \rightarrow 0^{+}, \\ O((D_{+}-D)^{\frac{1}{2}}) &\text{ as } D \rightarrow D_+^-. \end{cases} \end{align}

However, for $D\gt D_+$ , there is a global constant $A_{\infty }'\neq 0$ such that

(191) \begin{align} u_T(z)\sim 1-A_{\infty }' e^{\sigma _1(D)z}\,\text{ as } z\to -\infty, \end{align}

recalling that $0\lt \sigma _1(D)\lt \sigma _+$ and decreasing for $D\in (D_{\infty },\infty )$ , with

(192) \begin{align} \sigma _1(D)\sim \begin{cases} \sigma _+-O((D-D_+)^{\frac{1}{2}})&\text{ as } D\to D^+_+ \\ \frac{(\sqrt 2 -1)}{\sqrt D} & \text{ as } D\to \infty . \end{cases} \end{align}

Principally, we see that a (TPTW) with minimum speed is monotone as $z\to -\infty$ when $D\in (D_+,\infty )$ , but has decaying oscillations as $z\to -\infty$ when $D\in (0,D_+)$ . Numerical solutions of (IBVP) with $v=2\sqrt{D}$ and $D\gt 0$ are shown in Figure 1 and are consistent with the existence of a (unique by fixing translation invariance with $u_T(0)=\frac{1}{2}$ ) (TPTW) at minimum propagation speed $v = 2\sqrt{D}$ , and that this is the travelling wave generated in the initial value problem (5) to (7) for $D\gt \Delta _1$ . For $D\lt \Delta _1$ , as discussed in Section 3, for moderate values of $D$ , the minimum speed (TPTW) provides a good approximation to the solution of (IBVP) at the wavefront, which moves at speed $2\sqrt{D}$ , ahead of a stationary periodic state.