2.1. The attracting and repelling QSSs

The PDE (2.1) has an attracting QSS for all $\mu < -\delta $ , where $\delta>0$ is small and $\mathcal {O}(1)$ with respect to $\varepsilon $ , and solutions approach it at an exponential rate. Similarly, it has a repelling QSS for all $\mu> \delta $ , from which solutions diverge at an exponential rate. These QSS are given by

(2.2) $$ \begin{align} A_{\mathrm{QSS}} (x,y) &= - \sqrt{\varepsilon} \frac{I_a(x,y)}{\mu+i\omega_0} \nonumber\\ &\quad+ \varepsilon^{{3}/{2}} \bigg( \frac{I_a (x,y)+d(\mu+i\omega_0)\Delta I_a(x,y)}{(\mu+i \omega_0)^3} -\frac{(1+i\alpha)I_a^3(x,y)}{(\mu+i\omega_0)^2(\mu^2+\omega_0^2)} \bigg) \nonumber \\ &\quad+ \mathcal{O}(\varepsilon^{{5}/{2}}). \end{align} $$

Here, the $\mathcal {O}(\varepsilon ^{{5}/{2}})$ terms depend on $(x,y)$ and $\mu $ .

2.2. Solution of the linearized equation

In this section, we consider $\mu _0<0$ and solve the linearized equation,

$$ \begin{align*} \begin{aligned} A_t &= (\mu + i \omega_0) A + \sqrt{\varepsilon} I_a (x,y) + \varepsilon d \Delta A,\\ \mu_t &= \varepsilon. \end{aligned} \end{align*} $$

This equation may be simplified using an integrating factor. Let $B(x,y,\mu )=A(x,y,\mu )e^{{-(\mu + i \omega _0)^2}/{2\varepsilon }}$ . The equation for B is an inhomogeneous heat equation,

(2.3) $$ \begin{align} \sqrt{\varepsilon} B_\mu = \sqrt{\varepsilon}d \Delta B + I_a(x,y) e^{{-(\mu+i\omega_0)^2}/{2\varepsilon}}. \end{align} $$

The solution of (2.3) with initial data $B(x,y,\mu _0)=A_0(x,y)e^{{-(\mu _0+i\omega _0)^2}/{2\varepsilon }}$ is obtained by superimposing the homogeneous solution $B_h(x,y,\mu )$ and the particular solution $B_p(x,y,\mu )$ . For $\mu> \mu _0$ , we find

$$ \begin{align*} B_h(x,y,\mu) = \frac{e^{{-(\mu_0+i\omega_0)^2}/{2\varepsilon}}}{4\pi d (\mu-\mu_0)} \int_{\mathbf{R}^2} \textrm{exp} \bigg[ \frac{-(x-x')^2 - (y-y')^2}{4d(\mu - \mu_0)} \bigg] A_0(x',y')\,dx'\,dy'. \end{align*} $$

Here, we have used the fundamental solution of the heat equation in two dimensions,

$$ \begin{align*} \Phi(x,y,\mu)=\frac{1}{4\pi d \mu}e^{({-(x^2 + y^2)})/{4d\mu}}. \end{align*} $$

Also, for $\mu> \mu _0$ ,

$$ \begin{align*} \begin{split} B_p(x,y,\mu) &= \frac{1}{\sqrt{\varepsilon}} \int_{\mu_0}^\mu g(x,y,\mu-\tilde{\mu}) e^{{-(\tilde{\mu} + i \omega_0)^2}/{2\varepsilon}}\,d\tilde{\mu}, \\ g(x,y,\mu) &= \frac{1}{4d\pi\mu} \int_{\mathbf{R}^2} \mathrm{exp}\bigg[ \frac{-(x-x')^2 - (y-y')^2}{4d\mu}\bigg] I_a(x',y')\,dx'\,dy'. \end{split} \end{align*} $$

Transforming back to the original dependent variable using the integrating factor, we find that, for $\mu> \mu _0$ , the homogeneous solution is

(2.4) $$ \begin{align} A_h(x,y,\mu) = \frac{e^{(\mu^2 - \mu_0^2 + 2i \omega_0(\mu-\mu_0))/2\varepsilon}}{4\pi d (\mu-\mu_0)} \int_{\mathbf{R}^2} e^{({-(x-x')^2 - (y-y')^2})/({4d(\mu - \mu_0)})} A_0(x',y')\,dx'\,dy'. \end{align} $$

Furthermore, for $\mu> \mu _0$ , the particular solution is

$$ \begin{align*} \begin{aligned} A_p(x,y,\mu) &= \frac{e^{(\mu+i\omega_0)^2/2\varepsilon}}{\sqrt{\varepsilon}} \int_{\mu_0}^\mu g(x,y,\mu-\tilde{\mu}) e^{{-(\tilde{\mu} + i \omega_0)^2}/{2\varepsilon}}\,d\tilde{\mu}, \nonumber \\ g(x,y,\mu) &= \frac{1}{4\pi d\mu} \int_{\mathbf{R}^2} \mathrm{exp}\bigg[ \frac{-(x-x')^2 - (y-y')^2}{4d\mu}\bigg] I_a(x',y')\,dx'\,dy'. \end{aligned} \end{align*} $$

This completes the derivation of the solution of the linearized equation. We calculate $A_h$ for several different types of initial data $A_0(x,y)$ in Section 3.1, and we calculate $A_p$ for several different types of source terms in Section 3.2. Also, we observe here that it will be useful to distinguish between initial data $A_0(x,y)$ given for $\mu _0\le -\omega _0$ and for $-\omega _0 < \mu _0 < 0$ .

2.3. Solutions stay near the repelling QSS at least until $\mu =\omega _0$

In this section, we consider solutions with $\mu _0 \le -\omega _0$ . We show that not only do the solutions of (2.1) with $\mu _0\le -\omega _0$ remain close to the attracting QSS until the time of the instantaneous Hopf bifurcation at $\mu =0$ , but after the parameter crosses the instantaneous Hopf bifurcation they remain close to the repelling QSS as well, at least until the time $\mu =+\omega _0$ at all points $(x,y)$ for the functions $I_a(x,y)$ we consider.

We consider complex values of $\mu $ in a horizontal strip with mid-line on the real axis and of sufficient height (at least $3\omega _0$ ). This enables the use of classical methods of stationary phase and steepest descents (see, for example, [Reference Bender and Orszag6, Reference Kevorkian and Cole26, Reference Murray34]) to calculate $A_p(x,y,\mu )$ . We find, for $\delta \le \mu \le \omega _0$ ,

(2.5) $$ \begin{align} A_p(x,y,\mu)&= - \sqrt{\varepsilon} \frac{I_a(x,y)}{\mu + i \omega_0}+ \varepsilon^{{3}/{2}} \bigg( \frac{I_a(x,y)+d(\mu+i\omega_0)\Delta I_a(x,y)} {(\mu+i \omega_0)^3} \bigg) \nonumber \\ &\quad + \mathcal{O}\bigg(\frac{\varepsilon^{{5}/{2}}}{(\mu+i\omega_0)^5}\bigg)\nonumber \\ &\quad + ( \sqrt{2\pi} g(x,y,\mu+i\omega_0) +\mathcal{O}(\sqrt{\varepsilon}) ) e^{(\mu+i\omega_0)^2/2\varepsilon}. \end{align} $$

The calculation is presented in Appendix A.

The first and second terms are precisely the leading-order terms in the expansion of the repelling QSS for the linear CGL (compare with (2.2), where the QSSs are given for the cubic CGL equation). The third (remainder) term contains the higher-order terms in the asymptotic expansion of the repelling QSS, and continued integration by parts will yield them. The fourth term is exponentially small for $\mu \in [\delta ,\omega _0 -K\varepsilon ^r)$ , for some $K>0$ and any $0<r<1$ . It is a classic Stokes-type term. This term is not in the expansion of the repelling QSS (on $\mu>\delta $ ) to all orders. Rather, it is beyond all orders, $\mathcal {O}(e^{-{\omega _0^2}/{2\varepsilon }})$ , arising naturally from tracking solutions on (and near) the attracting QSS along a contour over the saddle point in the complex $\mu $ plane and into the regime of $\mathrm {Re} (\mu )>0$ . It is a measure of the exponentially small distance between the attracting and repelling QSS at $\mu =0$ .

Overall, formula (2.5) shows that, at all points $(x,y)$ , the solutions of the linear CGL equation with Gevrey regular data $A_0(x,y)$ given at $\mu _0 \le -\omega _0$ remain near the repelling QSS at least until $\mu = \omega _0$ to leading order.

Next, one also needs to include the nonlinear terms from (2.1). This was done for DHB in the one-dimensional CGL PDE in Section 6 of [Reference Goh, Kaper and Vo17]. There, we first used the same type of integrating factor (as used in (2.3) above) to rewrite the nonlinear PDE for $B(x,\mu )$ . Then, we split the solution into two parts: $B(x,\mu ) = B_p(x,\mu ) + b(x,\mu )$ , where we recall that $B_p$ is the particular solution of the linearized equation. The PDE for the remainder term $b(x,\mu )$ was converted into an integral equation. We showed formally that there is a mild solution of that integral equation, using an iterative method, and that the magnitude of $b(x,\mu )$ remains small at least until $\mu =\omega _0$ . Hence, the solution of the full nonlinear PDE remains near the repelling QSS at least until $\omega _0$ .

As the estimates of the mild formulation of the one-dimensional nonlinear equation only used $L^\infty \rightarrow L^\infty $ estimates for the heat semi-group, we expect that a similar formal analysis will hold for two spatial dimensions, as well. That is, based on decomposing $B(x,y,\mu ) = B_p(x,y,\mu ) + b(x,y,\mu )$ , we expect that $b(x,y,\mu )$ remains small at least until $\mu $ reaches $\omega _0$ . Fundamentally, the nonlinear terms remain exponentially small at least as long as the linear terms do.

3.1. The spatio-temporal memory surface

By writing the homogeneous solution $A_h(x,y,\mu )$ as a single exponential function, we define the spatio-temporal memory surface to be the set of points $(x,y,\mu _{\mathrm {mem}}(x,y))$ at which the real part of the argument of the exponential vanishes. In this section, we examine three different types of initial data: constant, Gaussian and periodic, in order to study how the spatio-temporal memory surface depends on the functional form of $A_0(x,y)$ .

For constant initial data $A_0(x,y)=1$ , one finds from (2.4) that

(3.1) $$ \begin{align} A_h(x,y,\mu)= e^{(\mu^2-\mu_0^2 + 2i\omega_0(\mu-\mu_0))/{2\varepsilon}}. \end{align} $$

It is independent of $(x,y)$ . Hence, $A_h(x,y)$ is exponentially small for all $\mu \in (\mu _0,-\mu _0)$ . At $\mu =-\mu _0$ , the real part of the exponential vanishes, which implies that the memory surface is a horizontal plane in the $(x,y,\mu )$ space: that is,

(3.2) $$ \begin{align} \mu_{\mathrm{mem}}(x,y)=-\mu_0 \quad \textrm{for all } (x,y), \end{align} $$

where we recall that $\mu _0<0$ . Then, for $\mu> - \mu _0$ , $A_h$ becomes exponentially large.

For Gaussian initial data $A_0(x,y)=e^{{-(x^2+y^2)}/{4\sigma }}$ , formula (2.4) shows that the homogeneous solution is given by

(3.3) $$ \begin{align} A_h(x,y,\mu) = \bigg[\frac{\sigma}{\sigma+d(\mu-\mu_0)}\bigg] e^{(\mu^2-\mu_0^2 +2i\omega_0(\mu-\mu_0))/2\varepsilon} e^{{-(x^2+y^2)}/({4(\sigma+ d(\mu-\mu_0))})}. \end{align} $$

At each $(x,y)$ , the real part of the argument of the exponential vanishes for

(3.4) $$ \begin{align} \frac{\mu^2 -\mu_0^2}{2\varepsilon} &+ \ln(\sigma) - \frac{1}{2}\ln ((\sigma+d_R(\mu-\mu_0))^2 + d_I^2(\mu-\mu_0)^2) \nonumber \\ &- \frac{(x^2+y^2)(\sigma + d_R(\mu-\mu_0))}{4[(\sigma+d_R(\mu-\mu_0))^2 + d_I^2 (\mu-\mu_0)^2]} = 0, \end{align} $$

where we recall that $d=d_R + i d_I$ in (2.1). Hence, to leading order asymptotically, we find that the spatio-temporal memory surface is parabolic in x and y,

(3.5) $$ \begin{align} \mu_{\mathrm{mem}}(x,y) &= \vert \mu_0 \vert - \frac{\varepsilon}{\vert \mu_0\vert} \bigg( \ln(\sigma) - \frac{1}{2}\ln ( (\sigma+2d_R\vert\mu_0\vert)^2 + 4d_I^2\vert\mu_0\vert^2) \bigg) \nonumber \\ &\quad + \frac{\varepsilon}{4\vert \mu_0 \vert} \bigg( \frac{(x^2+y^2)(\sigma + 2d_R\vert\mu_0\vert)}{(\sigma+2d_R\vert\mu_0\vert)^2 + 4d_I^2 \vert\mu_0\vert^2} \bigg) + \mathcal{O}(\varepsilon^2). \end{align} $$

For periodic initial data $A_0(x,y)=\cos (\pi(x-y)/L) \cos (\pi(x+y)/L)$ , the homogeneous solution for $\mu> \mu _0$ is

(3.6) $$ \begin{align} A_h(x,y,\mu) = e^{(\mu^2-\mu_0^2 +2i\omega_0(\mu-\mu_0))/{2\varepsilon}} e^{-{4\pi^2 d (\mu-\mu_0)}/{L^2}} \cos \frac{\pi}{L}(x-y)\cos\frac{\pi}{L}(x+y). \end{align} $$

Hence, at any point $(x,y)$ , the real part of the argument vanishes for

(3.7) $$ \begin{align} \frac{1}{2\varepsilon}(\mu^2 - \mu_0^2) - \frac{4\pi^2 d (\mu - \mu_0)}{L^2} + \ln \bigg\lvert\! \cos\frac{\pi}{L}(x-y)\cos\frac{\pi}{L}(x+y)\bigg\rvert = 0. \end{align} $$

Asymptotically to leading order, one finds logarithmic dependence on $A_0$ ,

(3.8) $$ \begin{align} \mu_{\mathrm{mem}}(x,y) = \vert \mu_0 \vert + \frac{8 \varepsilon \pi^2 d}{L^2} -\frac{\varepsilon}{\vert \mu_0\vert} \ln \bigg\lvert\! \cos\frac{\pi}{L}(x-y)\cos\frac{\pi}{L}(x+y)\bigg\rvert. \end{align} $$

These theoretical results for $\mu _{\mathrm {mem}}(x,y)$ with constant, Gaussian and periodic initial data are compared with simulations of (2.1) in Section 4.

3.2. The spatio-temporal buffer surface

In this section, we focus on the particular solution $A_p(x,y,\mu )$ with special emphasis on the fourth (last) term in the general formula (2.5). We consider $\delta \le \mu \le \omega _0$ and label this term G, where

(3.9) $$ \begin{align} G(x,y,\mu) = ( \sqrt{2\pi} g(x,y,\mu+i\omega_0) +\mathcal{O}(\sqrt{\varepsilon})) e^{(\mu+i\omega_0)^2/{2\varepsilon}}. \end{align} $$

The function G is the component of the particular solution $A_p(x,y,\mu )$ that measures the deviation from the repelling QSS (where we recall that the first three terms in (2.5) represent the QSS). It is generated by passage over the saddle point at $\mu =-i\omega _0$ in the complex $\mu $ -plane, as shown in Appendix A. It is present in all solutions that start from initial data at any $\mu _0<0$ , including solutions on the attracting QSS.

The spatio-temporal buffer surface is defined as the surface along which

(3.10) $$ \begin{align} \lvert G(x,y,\mu) \rvert = 1, \end{align} $$

to leading order. Here, we derive a general formula for the spatio-temporal buffer surface for smooth, bounded sources $I_a(x,y)$ , and asymptotic formulas for it with constant, Gaussian, periodic and stripe sources. We see that, in general, G ceases to be exponentially small and becomes $\mathcal {O}(1)$ in a spatially dependent manner. This implies that all solutions with initial data given at $\mu _0<0$ , including those on the attracting QSS, must leave an $\mathcal {O}(1)$ neighbourhood of the repelling QSS when $\mu $ reaches the buffer surface, irrespective of how large $|\mu _0|$ is; that is, irrespective of how far in the distant past the solutions approached the attracting QSS.

From (3.9), we find that the spatio-temporal buffer surface is given implicitly by

$$ \begin{align*} \mathrm{Re} \bigg[\!\!\ln ( \sqrt{2\pi} g(x,y,\mu+i\omega_0)) + \frac{1}{2\varepsilon} (\mu+i\omega_0)^2 \bigg] = 0 \end{align*} $$

for general smooth and bounded source terms $I_a(x,y)$ . We label the solution $\mu _{\mathrm {buf}}(x,y)$ . Its graph is the spatio-temporal buffer surface. To leading order,

(3.11) $$ \begin{align} \mu_{\mathrm{buf}}^2(x,y) = \omega_0^2 - \varepsilon \ln(2\pi) - 2\varepsilon \ln\lvert g(x,y,\omega_0+i \omega_0)\rvert. \end{align} $$

At each point $(x,y)$ , the buffer surface $\mu _{\mathrm {buf}}(x,y)$ determines the time at which the deviation of $A_p(x,y,\mu )$ from the repelling QSS ceases to be exponentially small.

We begin with constant source terms, $I_C(x,y)=c$ , and set $c=1$ without loss of generality. Substituting this into the integral for g in (2.5), we find

(3.12) $$ \begin{align} g(x,y,\mu-\tilde{\mu}) = 1\quad \textrm{for all } (x,y) \quad \textrm{and}\quad \mu>\tilde{\mu} \ge \mu_0. \end{align} $$

Hence, by (3.11), the buffer surface is to leading order

(3.13) $$ \begin{align} \mu_{\mathrm{buf}}(x,y) = \omega_0 - \frac{\varepsilon}{2\omega_0} \ln(2\pi)\quad \textrm{for all } (x,y). \end{align} $$

It is independent of $(x,y)$ . Hence, at all points in the domain, the particular solution ceases to be exponentially small uniformly at this time, $\mu _{\mathrm {buf}}$ .

Next, we consider Gaussian source terms, $I_G(x,y)=e^{{-(x^2+y^2)}/{4\sigma }}$ with $\sigma>0$ . Using (2.5), we find

$$ \begin{align*} g(x,y,\mu-\tilde{\mu}) = \bigg( \frac{\sigma}{\sigma+d(\mu-\tilde{\mu})} \bigg) \,\mathrm{exp} \bigg[ \frac{-(x^2+y^2)}{4(\sigma+ d(\mu-\tilde{\mu}))}\bigg]. \end{align*} $$

Hence, (3.11) implies that, to leading order, the buffer surface is

(3.14) $$ \begin{align} \mu_{\mathrm{buf}}(x,y) &= \omega_0 - \frac{\varepsilon}{2\omega_0} \ln \bigg(\frac{2\pi\sigma^2}{(\sigma+\omega_0(d_R-d_I))^2 + \omega_0^2 (d_R+d_I)^2}\bigg) \nonumber \\ &\quad + \frac{\varepsilon}{4\omega_0} \bigg[\frac{(\sigma+\omega_0(d_R-d_I))}{(\sigma+\omega_0(d_R-d_I))^2 + \omega_0^2 (d_R+d_I)^2}\bigg] (x^2+y^2). \end{align} $$

For periodic source terms, $I_P(x,y)=\cos (\pi(x-y)/L) \cos (\pi(x+y)/L)$ , (2.5) yields

$$ \begin{align*} g(x,y,\mu-\tilde{\mu}) = e^{-4\pi^2 d(\mu -\tilde{\mu})/L^2} \cos\frac{\pi}{L}(x-y)\cos\frac{\pi}{L}(x+y). \end{align*} $$

Hence, (3.11) implies that, to leading order, the buffer surface is

(3.15) $$ \begin{align} \mu_{\mathrm{buf}}(x,y) &= \omega_0 -\frac{\varepsilon}{2\omega_0}\ln(2\pi) +\frac{4\varepsilon\pi^2(d_R - d_I)}{L^2} \nonumber \\ &\quad -\frac{\varepsilon}{\omega_0} \ln\bigg\lvert\!\cos\frac{\pi}{L}(x-y) \cos\frac{\pi}{L}(x+y)\bigg\rvert. \end{align} $$

Finally, we let H denote the Heaviside step function ( $1$ on $x>0$ and $0$ on $x<0$ ) and study stripe source terms

$$ \begin{align*} I_S(x,y)=\sum_k [ H(x-x_k)-H(x-(x_k+h)) ]. \end{align*} $$

Here, h is the width of the stripe, $x_k$ denotes the left edge of the kth stripe and $x_{k+1}-x_k =\Delta> h$ for each k. We use (2.5) to derive

(3.16) $$ \begin{align} g(x,y,\mu-\tilde{\mu}) = \frac{1}{2} \sum_{k} \bigg( \operatorname{erf} \bigg( \frac{x_k+h-x}{\sqrt{4d(\mu-\tilde{\mu})}} \bigg) - \operatorname{erf} \bigg( \frac{x_k-x}{\sqrt{4d(\mu-\tilde{\mu})}} \bigg) \bigg), \end{align} $$

where $\operatorname {erf}z = {2}/{\sqrt {\pi }} \int _0^z e^{-t^2}\,dt$ is the error function. Hence, (3.11) implies that, to leading order, the buffer surface is given by

(3.17) $$ \begin{align} \mu_{\mathrm{buf}}(x,y) = \omega_0 - \frac{\varepsilon}{2\omega_0} \ln(2\pi) - \frac{\varepsilon}{\omega_0} \ln \lvert g(x,y,\omega_0(1+i)) \rvert, \end{align} $$

where g is given by (3.16). For these source terms, the spatio-temporal buffer surface is compared with the numerical solutions of the PDE (2.1) in Section 4.

4.1. DHB with constant source term

In the first representative simulation, we study DHB in (2.1) with constant source term $I_C(x,y)=1$ and Gaussian initial data

(4.1) $$ \begin{align} A_0(x,y) = c_1+c_2 e^{-(x^2+y^2)/4\sigma},\quad c_1,\quad c_2> 0, \end{align} $$

given at $\mu _0=-0.3$ , which we note is in $(-\omega _0,0)$ . The dynamics are shown in Figure 1, and a quantitative comparison between the escape surface for this solution and the memory surface calculated using the analysis of the previous section is shown in Figure 2. (Note that the solution is radially symmetric, because the initial data was chosen to be radially symmetric for simplicity in this first simulation; see also Figure 1(a).)

Figure 1

Delayed onset of oscillations in (2.1) with constant source term $I_C(x,y)=1$ and Gaussian initial data (4.1) given at $\mu _0 = -0.3$ . (a) Snapshot of $\operatorname {Re}u$ at $\mu = -\mu _0 = 0.3$ showing that the solution is radially symmetric. (b) Space-time evolution along $y=0$ . The red line indicates the instantaneous Hopf bifurcation. The temporal evolution of $\operatorname {Re}u$ (black curve) is compared with the QSS (blue curve) at the centre of the Gaussian at $(x,y)=(0,0)$ (c), at a radial distance of three spatial units from the origin at $(x,y) = ({3}/{\sqrt {2}},{3}/{\sqrt {2}})$ (d), and at a radial distance of 10 spatial units from the origin at $(x,y) = ({10}/{\sqrt {2}},{10}/{\sqrt {2}})$ (e). In each case, the numerical solution stays close to the QSS well past the instantaneous Hopf bifurcation at $\mu = 0$ , after which there is a hard onset to large-amplitude oscillations.

Figure 2

(a) In the three-dimensional $(x,y,\mu )$ space, the surface $\mu _{\mathrm {esc}}(x,y)$ is where the hard onset of the oscillations occurs. It has been obtained from direct numerical simulation of the PDE (2.1) with constant source term $I_C(x,y)=1$ and Gaussian initial data $A_0(x,y)=c_1+c_2 e^{{-(x^2+ y^2)}/{4\sigma }}$ with $(c_1,c_2,\sigma ) = (0.5,0.5,2.5)$ , given at $\mu _0=-0.3$ . (b) The memory surface $\mu _{\mathrm {mem}}(x,y)$ is given by (3.4). It gives the leading order asymptotics of the escape time for all $(x,y)$ . (c) The difference $|\mu _{\mathrm {esc}} (x,y) - \mu _{\mathrm {mem}}(x,y)|$ is shown in the projection onto the $(x,y)$ plane. Here, the parameters are $d=1$ , $\omega _0=0.5$ , $\alpha =0.2$ and $\varepsilon =0.01$ .

As expected, the solution stays near the repelling QSS at least until $\mu $ reaches $-\mu _0$ (see Figure 1(a) and (b)). Then, the delayed, post-Hopf, temporal oscillations are initiated at the origin of the domain, near $\mu =-\mu _0$ (see Figure 1(c)). As $\mu $ continues to slowly increase, the oscillations occur on successively larger disks about the origin (see Figure 1(d) and (e) for two other points). Then, once $\mu $ reaches approximately $0.327$ (which depends on the chosen threshold), there is a fairly rapid transition to large-scale oscillations, and these occur on the entire domain.

For this first simulation, we also show the escape surface in Figure 2(a). Below the escape surface, that is for all $\mathcal {O}(1)$ values $\mu < \mu _{\mathrm {esc}}(x,y)$ , the solution is near the repelling QSS. Then, at each point $(x,y)$ , the oscillations of amplitude $\mathcal {O}(\sqrt {\varepsilon })$ set in when $\mu $ reaches $\mu _{\mathrm {esc}}(x,y)$ . Furthermore, at each point $(x,y)$ , the amplitude of the oscillations becomes large ( $\mathcal {O}(1)$ ) as soon as $\mu $ is slightly beyond $\mu _{\mathrm {esc}}(x,y)$ , which confirms for all points $(x,y)$ what is shown for just three points in Figure 1.

For comparison, we show the memory surface $\mu _{\mathrm {mem}} (x,y)$ in Figure 2(b). The memory surface was computed as follows. First, we combined the homogeneous solutions (3.1) and (3.3) to determine the homogeneous solution, $A_h$ , with the initial data (4.1), $A_0(x,y)=c_1+c_2 e^{{-(x^2+ y^2)}/{4\sigma }}$ . Then, by enforcing the condition $| A_h | = 1$ , we obtained the leading order asymptotic relationship

(4.2) $$ \begin{align} \mu^2 - \mu_0^2 + 2\varepsilon \ln \bigg| c_1 + c_2 \frac{\sigma}{\sigma+d(\mu-\mu_0)} \exp \bigg( \frac{-(x^2+y^2)}{4(\sigma+d(\mu-\mu_0))} \bigg) \bigg| = 0 \end{align} $$

for the memory surface corresponding to this $A_0$ .

The difference $|\mu _{\mathrm {esc}} (x,y) - \mu _{\mathrm {mem}}(x,y)|$ is shown in Figure 2(c). At the centre of the Gaussian, the difference is small (with magnitude of approximately $10^{-5}$ , that is, $\mathcal {O}(\varepsilon ^{5/2})$ ). Then, in an annular region about the origin (red and dark red), the difference is slightly larger, due to nonlinear effects. Finally, both surfaces exhibit a fairly rapid transition into the regime (blue region in Figure 2(c), red in (a) and (b)) where they are essentially constant, since the Gaussian is tiny. Here, $\mu _{\mathrm {mem}}(x,y)=\sqrt {\mu _0^2 - 2 \varepsilon \ln c_1} \approx 0.322$ to leading order (as obtained from (4.2) for large $x^2 + y^2$ ), which agrees well with the value observed numerically in (2.1). Similar results were observed for solutions and the escape and memory surfaces with other values of $-\omega _0< \mu _0<0$ (data not shown).

Moreover, for the CGL with this (constant) source term, the buffer surface lies above the memory surface for all $(x,y)$ in three dimensions, since ${\mu _{\mathrm {buf}}(x,y) = \sqrt {\omega _0^2 - \varepsilon \ln (2\pi )}\!\approx \!0.481}$ for all $(x,y)$ . Hence, for all $(x,y)$ , the minimum in $\mu _{\mathrm {esc}}(x,y) \sim \mathrm {\min } \{ \mu _{\mathrm {mem}} (x,y), \mu _{\mathrm {buf}}(x,y) \}$ is given entirely by the memory surface in this simulation.

The second representative simulation is also with $I_C(x,y)=1$ . However, now the initial data given at $-\omega _0 < \mu _0 < 0$ is periodic,

$$ \begin{align*} A_0(x,y) = p_1 + p_2 \cos \frac{\pi}{L}(x-y) \cos \frac{\pi}{L}(x+y), \end{align*} $$

with $p_1>p_2>0$ , so that $A_0(x,y)$ is strictly positive everywhere. The results are shown in Figure 3.

Figure 3

(a) In the three-dimensional $(x,y,\mu )$ space, the surface $\mu _{\mathrm {esc}}(x,y)$ is where the hard onset of the oscillations occurs. It has been obtained from direct simulation of (2.1) with constant source term $I_C(x,y)=1$ and periodic initial data $A_0(x,y)=p_1+p_2 \cos (\pi(x-y)/L) \cos (\pi(x+y)/L)$ at $\mu _0=-0.3$ with $(p_1,p_2,L)=(1,0.5,25)$ . (b) The memory surface $\mu _{\mathrm {mem}}(x,y)$ given by (3.7). (c) The surface $|\mu _{\mathrm {num}} (x,y) - \mu _{\mathrm {mem}}(x,y)|$ shown in the projection onto the $(x,y)$ plane. The parameters are $d=1, \omega _0=0.5$ , $\alpha = 0.2$ and $\varepsilon =0.01$ .

The escape surface computed from the numerical simulations (Figure 3(a)) shows that the oscillations first occur at time $\mu \approx 0.281$ at the points $(x,y)=(j_1 L, j_2 L)$ , where $j_1,j_2=-1,0,1$ (dark blue), that is, at the maxima of $A_0(x,y)$ , where $\cos (\pi(x-y)/L) \cos (\pi(x+y)/L)=1$ . About each of those points, the oscillations set in on successively larger disks as $\mu $ slowly increases, until those disks collide at $\mu \approx 0.306$ (near the transition from yellow to orange). Then, as $\mu $ increases further, the escape surface consists of the four inverted paraboloid segments (orange and red). The peaks occur at time $\mu \approx 0.315$ at the points $(x,y)=( {j_1 L}/{2}, {j_2 L}/{2})$ , where $j_1,j_2=-1,1$ (dark red), that is, at the minima of $A_0(x,y)$ , where $\cos (\pi(x-y)/L) \cos (\pi(x+y)/L) = -1$ .

For comparison, the memory surface $\mu _{\mathrm {mem}}(x,y)$ is defined by

(4.3) $$ \begin{align} \mu^2-\mu_0^2 + 2\varepsilon \ln \bigg| p_1+p_2 \exp \bigg(-\frac{4\pi^2 d(\mu-\mu_0)}{L^2} \bigg) \cos \frac{\pi}{L}(x-y) \cos \frac{\pi}{L}(x+y) \bigg| = 0, \end{align} $$

and shown in Figure 3(b). It has global minima (dark blue) when

$$ \begin{align*}\mu_{\mathrm{mem}}=\sqrt{\lvert \mu_0 \rvert^2 - 2 \varepsilon \ln ( 1 + 0.5 e^{({-8\pi^2 d |\mu_0|})/{L^2}} ) } \approx 0.2866 \end{align*} $$

to leading order, with $\mu _0=-0.3$ . These occur at the points where $A_0(x,y)$ has its maxima. Also, the memory surface has global maxima (dark red) when $\mu \approx 0.3211$ , at the points where $A_0(x,y)$ has its minima. In between, it has the same conical shape qualitatively as the escape surface.

The difference between $\mu _{\mathrm {esc}}(x,y)$ and $\mu _{\mathrm {mem}}(x,y)$ is shown in Figure 3(c). The nonlinear terms cause $\mu _{\mathrm {esc}}$ to grow more steeply from the local minima, compared with the memory surface, so that the disks collide at a slightly later time than for the memory surface. Then, after the disks collide, the inverted paraboloids are slightly wider in $\mu _{\mathrm {esc}}(x,y)$ due to the nonlinear terms.

Finally, for this second representative simulation, we report that the buffer surface also lies above the memory surface for all $(x,y)$ in three dimensions. Indeed, with constant source $I_C(x,y)=1$ , $\mu _{\mathrm {buf}}(x,y) = \sqrt {\omega _0^2 - \varepsilon \ln (2\pi )}\approx 0.481$ . Hence, here also the minimum is given by the memory surface for all $(x,y)$ .

4.2. DHB with Gaussian source term

Gaussian source terms can be used to model spatially localized, radially symmetric inputs, such as a circular spot of visible light that shines on the reactor in a chemical pattern-forming experiment. A representative example of (2.1) with Gaussian source term is shown in Figure 4. Here, $I_G(x,y) = \exp ( -({x^2+y^2})/{4\sigma } )$ with $\sigma = 1$ , and the initial data at $\mu _0 = -0.75$ (which is chosen so that $\mu _0 < - \omega _0$ ) is

$$ \begin{align*} A_0(x,y) = \cos \frac{\pi}{L}(x-y)\cos \frac{\pi}{L} (x+y)\quad\text{with } L=25. \end{align*} $$

Figure 4

(a) In the three dimensional $(x,y,\mu )$ space, the surface $\mu _{\mathrm {esc}}(x,y)$ is where the hard onset of the oscillations occurs. It has been obtained from direct simulation of (2.1) with Gaussian source term $I_G(x,y)=\exp ( -({x^2+y^2})/{4\sigma } )$ and periodic initial data $A_0(x,y)=\cos (\pi(x-y)/L) \cos (\pi(x+y)/L)$ at $\mu _0=-0.75$ with $L=25$ . (b) The predicted escape surface is given by the minimum of $\mu _{\mathrm {buf}}(x,y)$ (parabolic part) and $\mu _{\mathrm {mem}}(x,y)$ (periodic part). (c) The surface $|\mu _{\mathrm {esc}} (x,y) - \min \{ \mu _{\mathrm {mem}}, \mu _{\mathrm {buf}} \} |$ shown in the projection onto the $(x,y)$ plane. The parameters are $d=1, \omega _0=0.5$ , $\alpha = 0.2$ and $\varepsilon =0.01$ .

From direct numerical simulations of the PDE, we find that the domain of the escape surface (Figure 4(a)) can be split into two distinct parts: the annular region $R_G = \{(x,y) \mid x^2+y^2 \leq 11.5^2 \}$ and its complement $R_P = \mathbb {R}^2 \backslash R_G$ . On $R_G$ , the earliest numerically detected escape occurs at the origin for $\mu \approx 0.4905$ . The escape surface is radially symmetric and increases on concentric rings, which thus creates a rotationally symmetric paraboloid. In contrast, for $(x,y)\in R_P$ , the escape surface is no longer radially symmetric, and is instead a periodic tile pattern, reflecting the initial data. The minima occur at $\mu \approx 0.738$ (yellow/green regions) and the maxima at $\mu \approx 0.812$ (red lines between yellow/green regions).

For this source term and initial condition, the buffer surface is given by (3.14) and the memory surface is given by (3.7). In Figure 4(b), we plot $\mathrm {\min } \{ \mu _{\mathrm {buf}}(x,y), \mu _{\mathrm {mem}}(x,y) \}$ . On $R_G$ , the minimum is given by the buffer surface, whereas it is given by the memory surface on $R_P$ . The buffer surface predicts that the earliest onset occurs at the origin at $\mu \approx 0.4908$ . The values of $\mu _{\mathrm {buf}}$ increase in a radially symmetric fashion until ${\mu \approx -\mu _0 = 0.75}$ on the boundary of $R_G$ . Then, for $(x,y)$ on $R_P$ , the memory surface predicts the onset. (Compare Figure 4(a) and (b).)

The difference $\lvert \mu _{\mathrm {esc}}(x,y)- \mathrm {\min }\{ \mu _{\mathrm {buf}}(x,y), \mu _{\mathrm {mem}}(x,y) \} \rvert $ is small throughout, as shown in Figure 4(c). The difference is especially small (of the order of the neglected terms in the asymptotic expansion) near the onset at the origin, which is at the minimum of the buffer surface, and also near the local minima of the memory surface.

4.3. DHB with stripe source term

Stripe source terms are also of fundamental interest. For example, chemical reactions irradiated with constant intensity light filtered through stripe masks can produce complex patterns (see [Reference Berenstein, Dolnik, Zhabotinsky and Epstein7, Reference Gentili and Micheau16] and references within).

Here, we use a simple model for a strictly positive stripe source term,

(4.4) $$ \begin{align} I_S(x,y)=1+\frac{1}{2}\sum_{k=-2}^2 [ H(x-x_k)-H(x-x_k-h) ], \end{align} $$

where h is the stripe width, $x_k$ denotes the left edge of the kth stripe and $x_{k+1}-x_k =\Delta> h$ . We set $x_0= -1.25$ , $h=2.5$ and $\Delta =10$ .

For the PDE (2.1), the surfaces $\mu _{\mathrm {esc}}(x,y)$ (calculated from direct numerical simulations) and $\mu _{\mathrm {buf}}(x,y)$ (calculated from (3.12) and (3.16)) are shown in Figure 5(a) and (b), along with the difference in Figure 5(c). From the plot of $\mu _{\mathrm {esc}}(x,y)$ , we see that the hard onset of oscillations occurs first at the points where $I_S(x,y)$ has its maximum, that is, inside the stripes. See the green stripes in the plot of $\mu _{\mathrm {esc}}(x,y)$ in Figure 5(a). The buffer curve $\mu _{\mathrm {buf}}(x,y)$ also has local minima inside the stripes (blue stripes in Figure 5(b)), and these stripes have approximately the same width as those of the escape surface. The difference in the green and blue stripes is of $\mathcal {O}(\varepsilon \sqrt {\varepsilon })$ .

In the complementary regions, the PDE simulations show that the duration of the DHB is longer. The escape surface lies just below $\omega _0$ , by about $2\varepsilon $ (see the yellow stripes in Figure 5(a)). In the complementary regions, the buffer surface transitions more gradually to its local maxima, which also occur at the local minima of the source term. Here, the difference is of $\mathcal {O}(\varepsilon \sqrt {\varepsilon })$ .

Figure 5

(a) The surface $\mu _{\mathrm {esc}}(x,y)$ , shown here in the projection onto the $(x,y)$ plane, is where the hard onset of oscillations occurs. It has been obtained from numerical simulation of (2.1) with stripe source term $I_S(x,y)$ (4.4) and constant initial data $A_0(x,y)=1$ given at $\mu _0=-1$ . (b) The buffer surface $\mu _{\mathrm {buf}}(x,y)$ defined by (3.11), with g given by the linear combination of (3.12) and (3.16). (c) The difference $|\mu _{\mathrm {esc}} (x,y) - \mu _{\mathrm {buf}}(x,y)|$ . The parameters are $d = 1, \omega _0=0.5$ , $\alpha = 0.2$ and $\varepsilon =0.01$ .

The narrow red strips between the green and yellow bands in Figure 5(a) correspond to the regions where the Heaviside function has a discontinuity. This discontinuity is smoothed out in the predicted buffer surface Figure 5(b), so that $\mu _{\mathrm {buf}}(x,y)$ is smooth and continuous across the whole domain. This difference is highlighted in Figure 5(c). There, the error is small and of $\mathcal {O}(\varepsilon \sqrt {\varepsilon })$ throughout the blue regions, and is of $\mathcal {O}(\varepsilon )$ in the thin strips where the Heaviside function has jump discontinuities.

We add that, in this simulation with the stripe source term (4.4) and the constant initial data given at $\mu _0=-1$ , the memory surface (3.2) lies well above the buffer surface (3.13) for all $(x,y)$ . Hence, the escape time is determined exclusively by the buffer surface, as reported.

4.4. The competition between $\mu _{\mathbf{mem}}(x,y)$ and $\mu _{\mathbf {buf}}(x,y)$

In the above simulations, we have shown that at each point $(x,y)$ , $\mu _{\mathrm {esc}}(x,y) \sim {\min }(\mu _{\mathrm {mem}}(x,y),\mu _{\mathrm {buf}}(x,y))$ to leading order. Therefore, in effect, there is a competition at every point $(x,y)$ between $A_h(x,y,\mu )$ and $A_p(x,y,\mu )$ in which the first one that ceases being exponentially small determines the maximal duration of the delay in the Hopf bifurcation to leading order.

The memory surface $\mu _{\mathrm {mem}}(x,y)$ depends to leading order on the initial time $\mu _0$ . Then, at $\mathcal {O}(\varepsilon )$ , it also depends on the logarithm of $\lvert A_h(x,y;\mu ) \rvert $ (recall, for example, equations (3.2), (3.5) and (3.8)).

The space-time buffer surface $\mu _{\mathrm {buf}}(x,y)$ (recall (3.11) in Section 3.2) is determined to leading order by $\omega _0$ , which is the frequency at the instantaneous Hopf bifurcation, and then at $\mathcal {O}(\varepsilon )$ by the logarithm of the Stokes term $g(x,y,\omega _0 (1+i))$ in (3.11). Indeed, as highlighted by the analysis in the complex $\mu $ plane in Appendix A, the Stokes line through the saddle point at $\mu =-i\omega _0$ crosses the negative $\mu $ -axis at $\mu =-\omega _0$ (see Figure A.1). Hence, the initial time $\mu _0$ is to the left of that here, and the contour used to find $A_p(x,y,\mu )$ for values of $\mu>0$ goes through that saddle. The saddle point method (applied in Appendix A) shows that the function $G(x,y,\mu )$ given by (3.9) arises owing to passage over the saddle at $-i\omega _0$ . For each $\mu>0$ , $G(x,y,\mu )$ measures the deviation of the solution $A_p$ from the repelling QSS. On $\mu> 0$ , it stays exponentially small at least until $\mu $ reaches $+\omega _0$ , where the other Stokes line through the saddle reaches the positive $\mu $ -axis. Then, as $\mu $ slowly increases beyond $\omega _0$ , $G(x,y,\mu )$ ceases to be exponentially small—and grows to become $\mathcal {O}(1)$ —in a spatially dependent manner. This is why the buffer surface is defined by (3.11), that is, when and where $\lvert G \rvert = 1$ .

Therefore, at $\mathcal {O}(1)$ , $\mu _0$ and $\omega _0$ determine the competition, and at $\mathcal {O}(\varepsilon )$ it is determined by the logarithmic terms. Moreover, for any finite value of $\varepsilon $ , there can be a changeover between the two terms for which one wins the competition, since, in general, the logarithmic terms in both expressions can grow due to their spatial dependence (see, for example, Figure 4).

5.1. Hopf bifurcation of the QSS

The Brusselator PDE (5.1) has a QSS,

$$ \begin{align*} \begin{aligned} u_{\mathrm{QSS}}(x,y) &= a(x,y) + \varepsilon \Delta \bigg( d_u a(x,y) + d_v \frac{b}{a(x,y)} \bigg) + \mathcal{O}(\varepsilon^2), \\ v_{\mathrm{QSS}}(x,y) &= \frac{b}{a(x,y)} + \frac{\varepsilon}{a^2(x,y)} \Delta \bigg((1-b) \frac{d_v b}{a(x,y)}-b d_u a(x,y) \bigg) + \mathcal{O}(\varepsilon^2). \end{aligned} \end{align*} $$

We perform a linear stability analysis about the QSS by setting $u=u_{\mathrm {QSS}} + \hat {u}$ and ${v = v_{\mathrm {QSS}} + \hat {v}}$ . The linearized system for $\hat {u}$ and $\hat {v}$ is

$$ \begin{align*} \left[\!\!\begin{array}{c} \hat{u}_t \\ \hat{v}_t \end{array}\!\!\right] = M \left[\!\!\begin{array}{c} \hat{u} \\ \hat{v} \end{array}\!\!\right] + \varepsilon \left[\!\!\begin{array}{cc} d_u & 0 \\ 0 & d_v \end{array}\!\!\right] \left[\!\! \begin{array}{c} \Delta \hat{u} \\ \Delta \hat{v} \end{array}\!\!\right], \end{align*} $$

where

$$ \begin{align*} M = \left[\!\!\begin{array}{cc} 2 u_{\mathrm{QSS}} v_{\mathrm{QSS}} - (1+b) & u_{\mathrm{QSS}}^2 \\ b - 2 u_{\mathrm{QSS}} v_{\mathrm{QSS}} & - u_{\mathrm{QSS}}^2 \end{array}\!\!\right]. \end{align*} $$

Hence, by setting $\mathrm {Tr}(M)=0$ , we find that the PDE (5.1) exhibits a spatially dependent Hopf bifurcation at $b = b_H(x,y) $ , where

(5.2) $$ \begin{align} b_H (x,y) = \frac{1 + a^2 + 2 \varepsilon a d_u \Delta a + \mathcal{O}(\varepsilon^2)} {1 + 2\varepsilon d_v ((1-a^2)/a) \Delta ({1}/{a}) + \mathcal{O}(\varepsilon^2) }, \end{align} $$

and we have written a for $a(x,y)$ . At each point $(x,y) \in \mathbf {R}^2$ , the QSS is linearly stable for $b < b_H(x,y)$ and linearly unstable for $b>b_H(x,y)$ . Moreover, the frequency, $\omega _H$ , at $b_H$ is $\omega _H = a(x,y)+\mathcal {O}(\varepsilon )$ .

5.2. Spatio-temporal dependence of the DHB in the Brusselator model

In this section, we present the results of numerical simulations showing the spatio-temporal dependence of the DHB in the Brusselator model (5.1). We take the initial values $b(t_0)< b_H(x,y)$ for all $(x,y)$ , so that the system slowly passes through the Hopf bifurcation $b_H(x,y)$ as b slowly increases (recall (5.2)). We work with a stripe source term,

(5.3) $$ \begin{align} a_S(x,y) = 1+ \varepsilon \alpha \sum_k [ H(x-x_k)-H(x-(x_k+h))], \end{align} $$

where $\alpha>0$ , $x_k$ denotes the left edge of the kth stripe, $x_{k+1}-x_k=\Delta $ , $h<\Delta $ is the width of the stripe and the stripes are sufficiently well separated.

In the simulations, the solutions rapidly approach the attracting QSS for $b < b_H(x,y)$ and, after the instantaneous Hopf point is crossed, they stay near the repelling QSS for long times, of $\mathcal {O}(1/\varepsilon )$ in the fast time t. This is the delay in the Hopf bifurcation, and the solutions leave a neighbourhood of the QSS in a spatially dependent manner. Moreover, they do so by making a large rapid jump, so that there is a hard onset of large-amplitude oscillations.

For the representative simulation shown in Figure 6, we set $\alpha =5$ , $h=0.1$ , $x_0=-0.05$ and $\Delta =0.3$ . The initial data is spatially periodic,

(5.4) $$ \begin{align} u_0(x,y)=a(x,y) + 0.5 \cos(\pi^2 x y),\quad v_0(x,t)= \frac{b_0}{a(x,y)} + \varepsilon. \end{align} $$

The initial value of the parameter is $b_0=b(t_0)=1.5$ (Figure 6(a)), below the Hopf point $b_H=2$ here. As b slowly increases, the solution rapidly approaches the attracting QSS, and the component of the solution given by the spatially periodic initial data decreases. By the time $b=1.7$ , the remnant of the initial data is faint (Figure 6(b)). Then, when b reaches $b_H=2$ , there is no visible trace of the initial data (Figure 6(c)).

Figure 6

DHB in the two-dimensional Brusselator with stripe source (5.3) and spatially periodic initial data (5.4) given at $b(t_0)=1.5$ . (a)–(i) The solution $u(x,y)$ rapidly converges to the source-dependent QSS, stays close to the QSS well beyond the instantaneous Hopf bifurcation at $b=2$ and exhibits a hard transition to oscillations. Note the different scales on the vertical colour bars in panels (a)–(i). (j) Space-time evolution along $y=0$ . The parameters are $\varepsilon = 0.01, d_u = 1 \times 10^{-3}$ , $d_v = 5 \times 10^{-4}$ , $\alpha =5$ , $h=0.1$ , $x_0=-0.05$ and $\Delta =x_{k+1}-x_k=0.4$ .

As b continues to increase, the solution remains near the repelling QSS (Figure 6(d), where $b=2.2$ ). It is not until $b=2.4$ that small amplitude temporal (post-Hopf) oscillations set in (Figure 6(e)). Also, the solution has regained a significant component given by the initial data. The oscillations rapidly gain large amplitude as b increases further (Figure 6(f)–(h), where $b=2.6$ , $b=2.8$ , $b=3.0$ ; note the changes in the vertical scales). Hence, the observed onset is a hard onset. Also, significant interactions occur between the components coming from the initial data and the stripe source.

Further increases in b result in spatio-temporal patterns, such as spiral waves (Figure 6(i)). These arise from the interaction of the oscillations (which occur beyond the memory surface), the stripe source and the diffusion.

A cross-section through $y=0$ is shown in Figure 6(j). The space-time plot shows that the solution stays close to the QSS (blue state) well past the instantaneous Hopf bifurcation (red line) and the large-amplitude oscillations do not set in until $b \approx 2.5$ . Cross-sections taken at other values of y are qualitatively the same, with some slight shifts in x.

The hard onset of oscillations occurs along the escape surface, $b_{\mathrm {esc}}(x,y)$ , which at each point $(x,y)$ is given by the value of b at which

$$ \begin{align*}|(u(x,y),v(x,y))-(u_{\mathrm{QSS}}(x,y),v_{\mathrm{QSS}}(x,y))| = \delta_{\mathrm{th}}.\end{align*} $$

This is the analogue of $\mu _{\mathrm {esc}}(x,y)$ . Plots of $b_{\mathrm {esc}}(x,y)$ are shown for two different initial values: $b(t_0)=0.1$ in Figure 7(a) and $b(t_0)=1.5$ in Figure 7(b). Here, the parameters are the same as in Figure 6. These escape surfaces reveal that the spatio-temporal buffer surface (defined by the particular solution of the equation linearized about the QSS) and the spatio-temporal memory surface (defined by the homogeneous solution) play similar roles here for (5.1). For initial values $b(t_0)$ far from $b_H=2$ , as in Figure 7(a), the escape surface is determined entirely by the stripe source and, hence, by the particular solution. The initial data has no influence on $b_{\mathrm {esc}}(x,y)$ , because the escape surface is almost invariant. Indeed, this was verified across four different choices of $b(t_0)$ . Then, for initial values of $b(t_0)$ close to $b_H$ (as in (b)), the escape surface is located below the buffer surface at all points. The homogeneous solution determines the exit time for all $(x,y)$ . The earliest escape times (dark blue regions) are where the initial data has its smallest amplitude (in the hyperbolic bands with the minima of the cosine term) and the light blue hyperbolic bands occur where the cosine term is maximal. In the regions of the stripes, the initial data is larger and the hard onset occurs later, although still before the buffer surface.

Figure 7

Contours of the escape surface, $b_{\mathrm {esc}}(x,y)$ , where the hard onset of oscillations occurs. The parameters are as in Figure 6 with different values of b at the initial time $t_0$ : $b(t_0)=0.1$ (a) and $b(t_0)=1.5$ (b). The escape surface in (a) depends only on the stripe source and has no memory of the initial data. Similar escape surfaces to (a) are observed for $b(t_0)=0, 0.2$ and $0.3$ . The escape surface in (b) depends on the initial profile.

Overall, we find that these primary features of DHBs in the two-dimensional Brusselator—the long, spatially dependent DHB; the hard, spatially dependent onset of large-amplitude oscillations; and the properties of the escape surface—are the same as those observed for DHBs in the CGL PDE in two dimensions (Section 4), as well as for DHBs in PDEs in one dimension [Reference Avitabile, Desroches, Veltz and Wechselberger3, Reference Goh, Kaper and Vo17, Reference Kaper and Vo25], in analytic ODEs [Reference Baer, Erneux and Rinzel4, Reference Hayes, Kaper, Szmolyan and Wechselberger20, Reference Neishtadt35–Reference Neishtadt and Treschev39, Reference Shishkova44, Reference Su46, Reference Su, Jones and Khibnik47] and in the applications cited in the Introduction.