Snaking branches of planar BCC fronts in the 3D Brusselator

TL;DR

This paper investigates complex snaking bifurcation structures of planar fronts in the 3D Brusselator model, revealing localized patterns and transitions between different spatial structures using numerical continuation methods.

Contribution

It applies advanced numerical bifurcation analysis to the 3D Brusselator, exploring snaking branches and localized solutions between BCCs and tubes, with theoretical insights from Maxwell points.

Findings

Identification of snaking branches of planar fronts

Approximation of localized BCCs and embedded structures

Analysis of moving fronts between lamellas and tubes

Abstract

We present results of the application of the numerical continuation and bifurcation package pde2path to the 3D Brusselator model, focusing on snaking branches of planar fronts between body centered cubes (BCCs) and the spatial homogeneous solution, and on planar fronts between BCCs and tubes (also called square prisms). These solutions also yield approximations of localized BCCs, and of BCCs embedded in a background of tubes (or vice versa). Additionally, we compute some moving fronts between lamellas and tubes. To give some theoretical background, and to aid the numerics for the full system, we use the Maxwell points for the cubic amplitude system over the BCC lattice.

Tables

Table 1. Table 1: Landau coefficients and other data (see Remark 2.2 ) for ( 16 ), A = 2 𝐴 2 A=2 .

$R$	$\delta$	$q$	$\alpha$	$\beta$	$\gamma$	$\epsilon =B_c-B_f$	$B_M$	${\tilde{B}}_M$
1	2/3	0	-3	-6	-9	0	NA	NA
0.75	0.665	0.29	-1.2	-5.95	-7.15	0.01	3.054	3.21
0.52	0.532	0.515	1.88	-5.1	-3.2	0.05	2.266	NA
0.4	0.476	0.625	4.7	-3.78	0.93	0.18	1.82	NA

Equations

$${\partial }_t{u}_1$$

$${\partial }_t{u}_2$$

$${\partial }_{\nu }{u}_1={\partial }_{\nu }{u}_2=0\text{ on }\partial \Omega ,$$

$$A=2,D_1=1,D_2=(A/R{)}^2,$$

$$B<B_c=(1+R{)}^2,$$

$$\displaystyle\begin{array}[]{lll}\vec{k}^{(1)}=-\vec{k}^{(7)}=\frac{k_{c}}{\sqrt{2}}(1,1,0),&\vec{k}^{(2)}=-\vec{k}^{(8)}=\frac{k_{c}}{\sqrt{2}}(0,1,1),&\vec{k}^{(3)}=-\vec{k}^{(9)}=\frac{k_{c}}{\sqrt{2}}(1,0,1),\\ \vec{k}^{(4)}=-\vec{k}^{(10)}=\frac{k_{c}}{\sqrt{2}}(1,-1,0),&\vec{k}^{(5)}=-\vec{k}^{(11)}=\frac{k_{c}}{\sqrt{2}}(0,1,-1),&\vec{k}^{(6)}=-\vec{k}^{(12)}=\frac{k_{c}}{\sqrt{2}}(-1,0,1),\end{array}$$

$${\partial }_{t}w=Lw+N(w),$$

$$w(\vec{x},t)=\sum _{j=1}^6{A}_j(t)e^{\text{i}{\vec{k}}^{(j)}\cdot \vec{x}}\Phi +\text{c.c.}+\text{h.o.t},$$

$$\dot{A_i}=f_i(A_1,A_2,A_3,A_4,A_5,A_6),i=1,\dots ,6.$$

$$\begin{array}{rl}{f}_1=& \lambda A_1+q(A_2{\overline{A}}_6+A_3{A}_5)+c_{31}|A_1{|}^2{A}_1+c_{32}(|A_2{|}^2+|A_3{|}^2+|A_5{|}^2+|A_6{|}^2)A_1\\ & +c_{33}|A_4{|}^2{A}_1+c_{34}(A_2{A}_4{A}_5+A_3{\overline{A}}_4{\overline{A}}_6),\\ f_2=& \lambda A_2+q(A_1{A}_6+A_3{\overline{A}}_4)+c_{31}|A_2{|}^2{A}_2+c_{32}(|A_1{|}^2+|A_3{|}^2+|A_4{|}^2+|A_6{|}^2)A_2\\ & +c_{33}|A_5{|}^2{A}_2+c_{34}(A_1{\overline{A}}_4{\overline{A}}_5+A_3{A}_5{A}_6),\end{array}$$

$$\lambda =\delta (B-B_c)\text{ with }\delta =\frac{A^2}{(A^2-R^2)(R+1)},q=\frac{2A(1-R)}{A^2-R^2},c_{31}=\frac{8-38R-5R^2+8R^3}{9R(A^2-R^2)}.$$

$$\underset{R\to 1}{\lim }{c}_{32}/c_{31}=\underset{R\to 1}{\lim }{c}_{33}/c_{31}=\underset{R\to 1}{\lim }{c}_{34}/c_{31}=2$$

$$(A_1,A_2,A_3,A_4,A_5,A_6)=(A_1,A_2,A_2,A_1,A_2,A_2),$$

$$\displaystyle\begin{array}[]{rl}\dot{A}_{1}&=\lambda A_{1}+2qA_{2}^{2}+(\alpha A_{1}^{2}+2\beta A_{2}^{2})A_{1},\\ \dot{A}_{2}&=\lambda A_{2}+2qA_{1}A_{2}+(\beta A_{1}^{2}+\gamma A_{2}^{2})A_{2},\end{array}$$

$$\alpha =c_{31}+c_{33},\beta =2c_{32}+c_{34},\gamma =c_{31}+2c_{32}+c_{33}+c_{34}.$$

$$\displaystyle\begin{aligned} w_{{\rm BCC}}=2A^{\text{bcc}}\biggl{[}&\cos(\kappa(x+y))+\cos(\kappa(y+z)+\cos(\kappa(x+z))\\ &+\cos(\kappa(x-y))+\cos(\kappa(y-z))+\cos(\kappa(-x+z))\biggr{]}\Phi+\text{h.o.t.}\end{aligned},$$

$$B_f=B_c+\frac{q^2}{(\alpha +2\beta )\delta }.$$

$$w_{\mathrm{tube}}=2A^{\text{tube}}[\cos (\kappa (x+y))+\cos (\kappa (x-y))]\Phi +\text{h.o.t.},$$

$${\vec{A}}_{\pm }^{\text{bcc}}:=(A_{\pm }^{\text{bcc}},A_{\pm }^{\text{bcc}})\text{ and }{\vec{A}}_{\pm }^{\text{tube}}:=(A_{\pm }^{\text{tube}},0).$$

$$\displaystyle\begin{array}[]{rl}\dot{A}_{1}&=-d_{4}{\partial}_{z}^{4}A_{1}+\lambda A_{1}+aA_{2}^{2}+b(\frac{3}{2}A_{1}^{2}+6A_{2})A_{1},\\ \dot{A}_{2}&=d_{2}{\partial}_{z}^{2}A_{2}+\lambda A_{2}+aA_{1}A_{2}+b(3A_{1}^{2}+\frac{9}{2}A_{2}^{2})A_{2}.\end{array}$$

$$E_{\text{kin}}=d_4\left[{\partial }_z^3{A}_1{\partial }_z{A}_1-\frac{1}{2}({\partial }_z^2{A}_1{)}^2\right]+d_2({\partial }_z{A}_2{)}^2$$

$$F=\frac{1}{2}\lambda (A_1^2+2A_2^2)+2q{A}_1{A}_2^2+\frac{\alpha }{4}{A}_1^4+\beta A_1^2{A}_2^2+\frac{\gamma }{2}{A}_2^4$$

$$F({\vec{A}}^{\text{bcc}})=0$$

$$F({\vec{A}}^{\text{bcc}})=F({\vec{A}}^{\text{tube}})$$

$$\mathrm{Res}:=\Vert G(u){\Vert }_{\infty }\stackrel{!}{<}10^{-8},$$

$$w=4A^{\ast }[\cos (\kappa (x+y))+\cos (\kappa (x-y))+\cos (\kappa y+\tilde{\kappa }z)]\Phi +\text{h.o.t.}$$

$$\displaystyle u_{1}(x)=\left\{\begin{array}[]{ll}u_{{\rm BCC},1}&z>5\\ u_{{\rm tube},1}&z\leq 5\end{array}\right.,$$

$$w_{\mathrm{lam}}\sim \cos (\kappa z)\text{ with }\kappa =1.5k_c,$$

Full text

Keywords: Localized 3D patterns; Brusselator; snaking; numerical continuation and bifurcation

1 Introduction

Turing patterns [Tur52] are stationary solutions of Reaction–Diffusion PDE systems that bifurcate from a homogeneous steady state which becomes unstable with respect to perturbations with a critical wave number $k_{c}\neq 0$. If the bifurcation is subcritical and the bifurcating branch stabilizes in a fold, then this gives bistability between the homogeneous state and the larger amplitude patterns in the subcritical regime, and this often yields the existence of localized patches of periodic patterns. These localized patterns exist in extended regions in parameter space [Pom86], and come in ’snaking’ branches which move back and forth in parameter space. This mechanism is well studied in the one–dimensional and two–dimensional cases (1D and 2D, respectively), see, e.g., [BK06, BK07, BKL*+*09, ALB*+*10, KUW19] for mainly numerical results, and [CK09, DMCK11, KC13, dW19] for analysis using the Ginzburg-Landau formalism and beyond all order asymptotics.

1D patterns extended homogeneously into a second and third direction are also solutions over 2D and 3D domains and are then referred to as stripes and lamellas, respectively. Typical genuine 2D patterns are squares and hexagons, and extended homogeneously in a third direction these yield (square and hexagon, respectively) tubes, while the simplest genuine 3D periodic patterns are cubes (or balls). Numerically, Turing patterns in 3D have so far mostly been studied by direct numerical simulation (DNS, aka numerical time integration) [WBD97, HSO07]. Additionally, some experimental results (and associated DNS for the Lengyel Epstein model) are reported in [BJVE11]. See also [AGH*+*05, GH08] for further results and discussion.

In [UW14] we numerically study planar fronts between stripes and hexagons in the 2D Schnakenberg model, using the package pde2path [UWR14, Uec19c]. Such fronts (or heteroclinic connections) can be naturally extended to localized patterns (or homoclinic cycles) by gluing together fronts and backs (i.e., considering heteroclinic cycles). See also, e.g., [Wet16, Wet18] for various further results on localized 2D patterns in different reaction-diffusion systems, and the Swift–Hohenberg equation as another prototype pattern forming system. Moreover, [Wet16] also contains a brief outlook on branches of 3D patterns, and some results on 3D patterns including localized patterns are also explained in [Uec19a], with detailed explanations on the background and usage of pde2path.

In a similar fashion as in [UW14] in 2D, here we study 3D planar fronts between cubes and the homogeneous steady state, and between cubes and tubes, in the Brusselator model [PL68]

[TABLE]

The chemical concentrations $u_{1}=u_{1}(\vec{x},t)$ and $u_{2}=u_{2}(\vec{x},t)$, with spatial coordinate $\vec{x}=(x_{1},x_{2},x_{3})\in\Omega\subset{\mathbb{R}}^{3}$ and time $t\geq 0$, correspond to an activator and inhibitor, respectively, $D_{1}$ and $D_{2}$ are their diffusivities, $A$ and $B$ are parameters, and $\Delta$ is the Laplacian. We let $u=(u_{1},u_{2})$, and instead of the coordinates $(x_{1},x_{2},x_{3})$ of $\vec{x}$ we also write $(x,y,z)$. Moreover, (1) has to be complemented by suitable boundary conditions (BCs) on ${\partial}\Omega$, and we will restrict to homogeneous Neumann BCs, i.e.,

[TABLE]

where ${\partial}_{\nu}$ denotes the outer normal derivative. For the initial value problem we also need to prescribe initial conditions $u|_{t=0}=u_{0}$.

Homogeneous steady states of (1) are given by $u_{1}^{*}=A$ and $u_{2}^{*}=B/A$. We fix the parameters

[TABLE]

where $R$ is used as a convenient unfolding parameter, see below, and consider $B$ as the primary bifurcation parameter. The steady state $u^{*}=(u_{1}^{*},u_{2}^{*})$ is stable for

[TABLE]

where a Turing bifurcation occurs with critical wave number $k_{c}=\sqrt{R}$. Our focus will be on solution branches corresponding to a so called body centered cube (BCC) lattice. Close to bifurcation, these may be described by a system of equations for six amplitudes $A=(A_{1},\ldots,A_{6})\in{\mathbb{C}}^{6}$, see (10) in §2. This amplitude system has a variety of steady solutions $A^{*}$, which in the original system (1) correspond to, e.g., lamellas, tubes, and cubes, henceforth called BCCs. These solution branches of the amplitude equations have been classified and discussed in detail in [CK97, CK99], and the stability of the associated solutions of the original system close to onset has been studied in [CK01]. See also [GS02, Hoy06] for textbook expositions of the underlying and very important symmetry perspective. Additionally, the (ODE) amplitude system can be formally extended to a (1D PDE) modulation equation system by assuming a slow dependence of the amplitudes on one spatial coordinate. The steady modulation equations have a spatially conserved quantity (the potential energy) which thus defines Maxwell points for heteroclinics between different fixed points $A^{*}$. These results can then be used to identify parameter regimes for the search of snaking branches of steady fronts between BCCs and $u^{*}$, and between BCCs and tubes. For this, to keep the numerics inexpensive we choose small $\Omega$, e.g. boxes $\Omega=(-l,l)^{2}\times(-l_{z},l_{z})$, where $l=\pi/k_{c}$ $l_{z}=4l$. Additionally, we briefly illustrate that near to but outside the snaking region the dynamics of close by solutions show a stick–slip motion, and give examples of other moving fronts, for instance between lamellas and cubes.

Remark 1.1**.**

Our results are not specific to the Brusselator model (1), but can be expected for all 3D pattern forming systems with a subcritical Turing bifurcation, or, more generally, systems with a bistability of patterns and the homogeneous solutions, or a bistability of different patterns. Similar results are provided for the (quadratic–cubic) 3D Swift–Hohenberg equation in [Uec19a, §3]. There we also give detailed explanations on a number of issues that arise for numerical branch continuation and bifurcation in 3D pattern forming systems, including:

• •

The algorithm for branch switching at branch points of higher multiplicity, which naturally arise in 3D due to symmetries, see also [Uec19a].

• •

Tricks, including remarks on the choice of meshes, to avoid uncontrolled ’branch jumping’, which is a major issue in particular in 3D due to the multitude of different branches close to each other.

Here we focus on the Brusselator model (1) as a standard reaction–diffusion model. $\rfloor$

Acknowledgment. The work of DW was supported by the DFG under Grant No. 264671738.

2 The amplitude formalism

We briefly review the BCC amplitude equations for (1) close to the primary bifurcation from $u^{*}$, focusing on the bifurcating branches pertaining to Neumann BCs over cuboids. Amplitude equations for 3D pattern forming systems are derived and discussed in [CK97, CK99] for three lattices with cubic symmetry, namely the simple cubic, the face-centered cubic, and the body-centered cubic (BCC).

2.1 Derivation of the amplitude system

For the BCC lattice, the critical wave vectors are given by

[TABLE]

where $k_{c}=\sqrt{R}$ is the critical wave number. Setting $u=u^{*}+w$, $w=w(\vec{x},t)\in{\mathbb{R}}^{2}$, yields

[TABLE]

where $L$ is the linearization around $u^{*}$ and $N(w)$ denotes the nonlinear terms. We make the ansatz

[TABLE]

where $A_{j}\in{\mathbb{C}}$, and $\Phi=(\Phi_{1},\Phi_{2})\in{\mathbb{C}}^{2}$ is the critical eigenvector, independent of $j$ due to the rotational invariance of the Laplacian, and normalized to $\Phi_{1}=1$. In (9), c.c. means the complex conjugate of the preceding terms, and h.o.t. denotes higher order terms, which turn out to be nonlinear terms in the $A_{j}$. Plugging (9) into (8), sorting wrt. to the modes ${\rm e}^{{\rm i}\vec{k}^{(j)}\cdot\vec{x}}$, first solving for uncritical modes at, e.g., $e_{0}$, $e_{2k_{1}}$ and so on, we obtain the amplitude equations

[TABLE]

Their general form, dictated by symmetry [CK97], is

[TABLE]

and the remaining $f_{j},j=3,\ldots,6$ (we shall not need them explicitly) also follow from symmetry. The BCC lattice supports three-wave interactions, e.g., $k^{(1)}=k^{(2)}-k^{(6)}=k^{(3)}+k^{(5)}$, which explains the occurrence of the quadratic terms $q(A_{2}{\overline{A}}_{6}+A_{3}A_{5})$ in $f_{1}$. For the coefficients $\lambda$, $q$ and $c_{31}$ we have some analytic formulas which follow from, e.g., [VdWDB92], see also [CK99], namely

[TABLE]

Similar formulas can be derived for $c_{32},c_{33}$ and $c_{34}$, but we refrain from doing so, and instead will use numerical values computed by the pde2path tool ampsys [UW19], which is designed to do such computations with minimal user input. 333In [UW19] we apply our tool to a variety of models and wave vector lattices, and cases with known coefficients such as (12) are useful for checks of the implementation. Conversely, using ampsys on known cases is helpful to make sure that scalings of amplitudes are taken care of correctly. For instance, the formulas for $q$, $c_{31}$ in [CK99] are different from (12) because they are based on a normalization of the critical eigenvector $\Phi=(\Phi_{1},\Phi_{2})$ which is different from our normalization with $\Phi_{1}=1$.

Quadratic and cubic terms are considered to be of the same order for the derivation of (10), and this formally requires $q$ to be small, which means $R\approx 1$, cf. Remark 2.2. Moreover, as in [CK99],

[TABLE]

and $\lim_{R\to 1}c_{31}=-1$ for the choice $A=2$, which we fix in the numerics. However, we shall be interested in $1-R={\cal O}(1)$, and the deviations of $c_{32},c_{33}$ and $c_{34}$ from $2c_{31}$ turn out to be significant in this case.

The bifurcation diagrams for (10) with $B$ close to $B_{c}$, have been discussed in detail in [CK97, CK99]. Here we restrict to those branches that fulfill Neumann BC on cuboid domains of the form $\Omega=(-l_{x},l_{x})\times(-l_{y},l_{y})\times(-l_{z},l_{z})$ with $l_{x}=n_{1}l,l_{y}=n_{2}l,l_{z}=n_{3}l$, $n_{j}\in{\mathbb{N}}/2$ and $l=\sqrt{2}\pi/k_{c}$. By (7) and (9), this restricts (modulo phase-shifts, i.e., spatial translations by $nl$ for some $n\in{\mathbb{N}}/2$) the admissible solutions to

[TABLE]

where $A_{1},A_{2}\in{\mathbb{R}}$ fulfill

[TABLE]

with the effective coefficients

[TABLE]

In Table 1 we list these coefficients (together with further data explained below) for some chosen values of $R$, for which we shall also run numerics on the full system (1).

Remark 2.1**.**

We decrease $R$ rather far from $R=1$ to $R=0.4$, and, moreover, will use (16) for $B-B_{c}={\cal O}(1)$. Each operation alone makes the applicability of (16) quite questionable. However, some of the interesting results will occur in the strongly subcritical regime $R=0.4$, and applying (16) with care we get good predictions for these, while other effects cannot be captured. See, e.g., the discussion of Fig. 2 below, and [BMvS09] for a related general discussion about the use of amplitude equations for subcritical bifurcations. $\rfloor$

2.2 Steady solutions

The solution $A_{1}=A_{2}=A^{\text{bcc}}_{\pm}$ of (16) with $A^{\text{bcc}}_{\pm}=-\frac{q}{\alpha+2\beta}\pm\sqrt{\frac{q^{2}}{(\alpha+2\beta)^{2}}-\frac{\lambda}{\alpha+2\beta}}$ yields BCCs for (1)

in the form

[TABLE]

and naturally phase shifts in $A_{j}$ correspond to translations in $x,y,z$, which we shall not distinguish from (18). The BCC branch bifurcates supercritically if $q=0$ ($R=1$) and transcritically if $q\neq 0$ ($R<1$), and has a fold in

[TABLE]

The system (10) is equivariant under $q\mapsto-q$ and $A_{i}\mapsto-A_{i}$, $i=1,\ldots,6$, and this is naturally inherited by (16). Depending on the sign of $q$, one direction of the BCCs has maxima of $u_{1}$ in the centers of the balls, and we call these ’hot’ balls, while in the other direction we have ’cold’ balls, see Fig. 1(d,e). This classification is analogous to ’spots’ and ’gaps’ in the 2D case.

The other primary solution branch of (16) (i.e., of (10) and compatible with Neumann BCs) yields ’tubes’ (called ’squares’ in [CK97] and in much of the literature), i.e., $A_{1}=A^{\text{tube}}_{\pm}:=\pm\sqrt{-\frac{\lambda}{\alpha}},A_{2}=0$, and these bifurcate in pitchforks. The associated solutions $w$ of (8) are spatially homogeneous in $z$ direction, namely

[TABLE]

see Fig. 1(b). From $\alpha$ in Table 1 we readily see that the bifurcation of the tubes changes from super– to subcritical for $R$ decreasing from $0.75$ to $0.52$.

For simplicity, we shall also denote the vectors $\vec{A}=(A_{1},A_{2})$ similarly, i.e.,

[TABLE]

For $R\neq 1$, both of these branches (families of branches, via symmetries), $w_{{\rm BCC}}$ and $w_{{\rm tube}}$, are unstable close to bifurcation. However, on the level of the amplitude equations, the BCCs stabilize after the fold, and the tubes at an $O(1)$ distance from onset, while the stable BCC branch again destabilizes at $O(1)$ distance, and there is an unstable mixed mode branch connecting the BCCs and the tubes.

In Fig. 1(a) we illustrate the branching behavior for $R=0.75$ of three nontrivial branches of (16). Additional to the BCCs (black) and tubes (red) there is the mixed mode branch (orange) connecting tubes and BCCs between the points where these gain/loose stability, i.e., in their bistable range. In particular, in the amplitude equations we get two bistabilities: (a) bistability of $A=0$ and the $A^{\text{bcc}}_{+}$ part below onset, and (b) bistability of $A^{\text{bcc}}_{+}$ and $A^{\text{tube}}$ at $O(1)$ distance above onset. This does in general not mean that the associated solutions inherit these in the full system [CK01]. However, this turns out roughly to be the case over sufficiently small domains (which may be significantly extended in $z$). This motivates our main aim, i.e., to find snaking branches of localized BCCs, or more precisely, of fronts (a) between BCCs and $u=u^{*}$, and (b) between BCCs and tubes.

Remark 2.2**.**

(a) The reduced amplitude equations (16) have exactly the same structure as the amplitude equations for the three modes $(A_{1},A_{2},A_{3})$ in the 2D case on a planar hexagonal lattice when restricted to the subspace $A_{2}=A_{3}$, see, e.g., [UW14, §3.1]. The 3D tubes and BCCs thus correspond to 2D stripes and hexagons, respectively, and the stabilities within the amplitude system are also equivalent. However, the stability and bifurcation structures of the associated solutions in the original 2D vs 3D systems will in general be rather different.

(b) The 7th column of Table 1 indicates how the ’subcriticality’ $\varepsilon:=B_{c}-B_{f}=q^{2}/((\alpha+2\beta)\delta)$, cf. (19), of the BCC branch increases with decreasing $R$. In, e.g., [CK09, DMCK11] it is explained (for Swift–Hohenberg models) that the snaking width of branches of fronts connecting a subcritical pattern and [math] is exponentially small in this subcriticality $\varepsilon$, which for (1) means that we expect that $|B_{l}-B_{r}|\sim c_{1}\varepsilon^{-1}{\rm e}^{-c_{2}/\varepsilon}$ where $c_{1},c_{2}$ are constants and $B_{l}$ and $B_{r}$ denote the left and right ends of the snaking range, respectively. Relatedly, the steepness of the fronts scales as $\varepsilon$, and hence the required domain length as $1/\varepsilon$. See also [UW14] for numerical illustrations of this phenomenon, and [dW19] for further references and a transfer of the results from [CK09, DMCK11] to reaction diffusion systems. Thus, if we assume that snaking branches between BCCs and $u=u^{*}$ exist, then Table 1 also indicates that finding these should be more robust and less expensive at smaller $R$. The 8th column gives (for $R<1$) the approximate Maxwell point between BCCs and [math] (see §2.3), and the 9th column the one between BCCs and tubes (which for the used parameters only exists in the second row). $\rfloor$

2.3 Maxwell points in the amplitude system

The amplitude system (16) also already contains the information to derive a necessary condition for fronts between BCCs and zero, or BCCs and tubes to exist on the level of the amplitude equations (16). If we assume a slow $z$ dependence of the amplitudes $A_{1},A_{2}$, then we can formally derive an extension of (16) to

[TABLE]

The second order coefficient $d_{2}$ is determined as $d_{2}=-\frac{1}{2}(n^{(2)}\cdot(0,0,{\partial}_{k}))^{2}\mu_{1}(k)$ where $n^{(j)}=\vec{k}^{(j)}/\|\vec{k}^{(j)}\|_{2}$, cf., e.g., [Pis06, §4.6]. For the mode $A_{1}{\rm e}^{{\rm i}\vec{k}^{(1)}\cdot\vec{x}}$ we have $n^{(1)}\cdot(0,0,{\partial}_{k})=0$ and hence must expand the dispersion relation $\mu_{1}(k)$ to 4th order around $\vec{k}^{(1)}$, yielding $d_{4}=\frac{1}{4!}{\partial}_{k}^{4}\mu_{1}(k_{c})$.

The system (24) has the conserved quantity $E=E_{\text{kin}}+F$, i.e., $\frac{d}{dz}E(A(z))=0$, where

[TABLE]

can be considered as a kinetic energy, and

[TABLE]

as a potential energy. Thus, a necessary condition for the existence of steady front solutions of (24), connecting, e.g., $\vec{A}^{\text{bcc}}_{+}$ at $z=-\infty$ with $\vec{A}=(0,0)$ at $z=\infty$, or $\vec{A}^{\text{bcc}}_{+}$ at $z=-\infty$ with $\vec{A}^{\text{tube}}$ at $z=\infty$, is that the limit states (where ${\partial}_{z}=0$) have the same potential energy, i.e.,

[TABLE]

These equalities only hold at specific points, the Maxwell points. As already said, for (1) we use $B$ as a bifurcation parameter, for fixed $R=0.75,R=0.52$ and $R=0.4$ (and, for completeness, $R=1$), yielding the coefficients $\delta,q,\alpha,\beta,\gamma$ from Table 1. In Fig. 2 we plot, for these values, $F$ for the $\vec{A}^{\text{bcc}}$ and $\vec{A}^{\text{tube}}$ branches as a function of $B$, which defines $\lambda$ via $\lambda=\delta(B-B_{c})$. This illustrates four main things, and in the following section we use these (formal) results from the amplitude system (16) as a guide to search for the associated solutions of the original system (1):

• •

For $R=1$, the bifurcations of the BCCs and tubes are supercritical, and the $F$–plots show that no fronts between any of $\vec{A}^{\text{bcc}}_{+},\vec{A}^{\text{tube}}$ and $\vec{A}=0$ can exist. Thus we also do not expect such fronts in the full system, at least not near onset where we expect the amplitude equations to make good predictions.

• •

For $R=0.75$, there is a Maxwell point $B_{M}$ near $B=3.054$ for a front between $\vec{A}^{\text{bcc}}_{+}$ and $(0,0)$. However, the subcriticality is very weak and thus we should expect finding the associated fronts (if they exist) in the original system (1) to be a delicate and expensive task, cf. Remark 2.2(b). For $R=0.75$ there additionally exists a Maxwell point $\tilde{B}_{M}$ near $B=3.21$ for a front between $\vec{A}^{\text{bcc}}_{+}$ and $\vec{A}^{\text{tube}}$, and thus the possibility of steady fronts between $w_{{\rm BCC}}$ and $w_{{\rm tube}}$ in this parameter regime. This prediction, in particular the quantitative value for $\tilde{B}_{M}$, should only be considered as a hint as we are relatively far from onset.

• •

For decreasing $R$, the fold $B_{f}(R)$ and the Maxwell points $B_{M}(R)$ move farther away from $B_{c}(R)$ (cf. Table 1). Thus, if steady fronts between $w_{{\rm BCC}}$ and $w_{{\rm tube}}$ exist, then they should be easier to find at smaller $R$.

• •

On the other hand, at $R=0.52$ the pitchfork for $\vec{A}^{\text{tube}}$ is subcritical, and we should not expect the associated branch to make any reasonable predictions away from onset. In particular, while the sub/vs supercritical branching behavior of $w_{{\rm tube}}$ is correctly predicted, the $w_{{\rm tube}}$ branch in (1) has a fold rather close to onset, and this behavior can only be resolved by 5th order amplitude equations, which we do not consider here.

3 Results for the full system

To illustrate/corroborate some results from the amplitude formalism, and to find snaking branches for (1), we use pde2path [UWR14, Uec19c]. The approach is motivated by (and the results essentially similar to) the results on snaking branches for 1D and 2D problems in [UW14, Uec19a], but as already noted in the Introduction, the 3D case does present a number of significant numerical challenges. Besides the obvious issue of higher numerical costs due to more degrees of freedom (DoF) in the (discretized) 3D case, these challenges mainly include the branch switching at branch points of high multiplicity, and, in particular over non-small small domains, problems with undesired ’branch jumping’ due to many solution branches close to each other (more than in 1D and 2D). See [Uec19a, §3] for details on how we deal with these problems. Here we only remark that:

• •

The branch switching proceeds by (numerically) deriving and solving the pertinent algebraic bifurcation equations, which are essentially equivalent to the amplitude equations. No specific knowledge of the structure of the bifurcation problem is needed for this, but the user can (and should) use the symmetries to make a selection of branches to be continued.

• •

To have reliable and fast numerics we stick to rather small domains; in particular, for fronts we extend small $x$–$y$–squares in $z$–direction, i.e., choose long and slender bars. The typical number of DoF in the results below is on the order of $10^{5}$, and the residual tolerance is

[TABLE]

where $G(u)$ is the FEM discretization of the right hand side of (1). Typical runtimes for continuation of branches on an I7 laptop are, e.g., about 10-20 min for 50 points, including the stability computation. In detail, e.g., for the cubic domains in Fig. 3 we have 83.000 mesh points, which on a rectangular grid would correspond to about 44 points in each spatial direction. However, for symmetry reasons a staggered grid of ’symmetry type 1’ (see [Uec19a, §3.4]) is used. A zoom of a typical mesh for the long and slender bars used to compute fronts is given in Fig.3(h).

• •

To check accuray we essentially recomputed selected solutions on finer grids, typically with double the DoF. When first interpolating

a given solution to a solution $u^{\text{interp}}$ on the finer grid and then using a Newton–loop to obtain a solution $u^{\text{fine}}$ fulfilling (29) on that grid, then in all cases $\|u^{\text{interp}}-u^{\text{fine}}\|_{\infty}<0.005$. Moreover, a few continuation steps on the finer grid yield the same behaviour as on the original grid. Additionally, we used adaptive mesh–refinement on selected solutions (see Fig. 3 for an example), with $\|u-u^{{\rm refined}}\|_{\infty}<0.005$, and without visible changes in the solution structure. Thus we believe that our meshes are sufficiently fine that the numerical solutions are converged and reflect the true PDE behaviour.

3.1 Fronts between BCCs and $u^{*}$

First we seek fronts between BCCs and $u^{*}$. Following Remark 2.2(b), this should be easier for smaller $R$ than for $R$ close to 1. In Fig. 3 we start with $R=0.52$ which is roughly the largest value for which we find a snaking branch of a front between BCCs and $u^{*}$ on a domain as in Fig.3(f–h). We come back to $R=0.75$ in §3.2 where we consider the bistability range of BCCs and tubes.

Figure 3(a) shows the bifurcation diagram of BCCs and tubes on a small cube (8 times the minimal domain, i.e., $\Omega=(-l,l)^{3}$ with $l=\sqrt{2}\pi/k_{c}$), and (b)–(d) shows sample plots. The BCC branch qualitatively (and also quantitatively) agrees with the predictions from (16). The tubes (magenta) bifurcate subcritically (as predicted), but in the full system have a fold close to the bifurcation and hence can only be approximated by the amplitude system close to onset. In (e) we illustrate adaptive mesh refinement from a solution as in (b), for graphical reasons showing 1/8 of the computational domain and starting with a relatively coarse uniform mesh with $n_{p}=10351$ grid points (half the number of mesh points in each direction compared to the computations in (a)–(d)), and then adapting to $n_{p}=17258$. See [Uec19b] for the setup of 3D mesh adaptation in pde2path, based on the package trullekrul [Jen17], and [Uec19a, §3.4] for examples of the usage for localized 3D patterns in the Swift–Hohenberg equation.

In (f)–(h) we focus on the subcritical range on a long and slender bar $\Omega=(-l_{x},l_{x})\times(-l_{y},l_{y})\times(-l_{z},l_{z})$ with $l_{x}=l_{y}=\sqrt{2}\pi/(2k_{c})$ and $l_{z}=4\sqrt{2}\pi/k_{c}$, and with $1.8\times 10^{5}$ DoF. The snaking red branch bifurcates from the BCCs (black branch) close to onset, and corresponds to a front between BCCs and $u^{*}$. In the snaking region around $B=2.25$, which is reasonably close to the Maxwell point prediction $B_{M}=2.266$, it alternates between stable and unstable parts, and in each pair of folds an additional layer of BCCs is added. In (i) we illustrate the (uniform) meshing in (f)–(h) on the subdomain $(l_{x}/2,l_{x})\times(l_{y}/2,l_{y})\times(-l_{z},-2l_{z}/3)\subset\Omega$ corresponding to $\frac{1}{4}\frac{1}{4}\frac{1}{6}=\frac{1}{96}$ of the full domain at the bottom front right.

To increase the narrow snaking region in Fig. 3, we lower $R$ further to $R=0.4$ in Fig. 4. The branch of localized BCCs bifurcates from the BCC branch near onset as before, but is now to the right of the Maxwell point prediction $B_{M}=1.82$. However, during the snaking the localized BCCs significantly change their wave lengths in $z$ direction and terminate in a pitchfork bifurcation on a branch corresponding to cubes of the form (modulo a phase shift in $x,y$)

[TABLE]

with $\kappa=k_{c}/\sqrt{2}$ and $\tilde{\kappa}=9\kappa/8$. Such shifts to patterns with slightly different $|\vec{k}|$ (sideband patterns), are also known from 1D and 2D, cf., e.g., [UW14]. They are analyzed (aided by numerics) in detail for the 1D quadratic–cubic Swift–Hohenberg equation in [BBKM08] on a finite domain, including a relation to the Eckhaus instabilities of the periodic branches, and the phenomena explained there also occur here: The termination of the snake at the upper end where the pattern almost fills the domain depends critically on the domain size and the strength of the subcriticality, i.e., the snaking width. For stronger subcriticality, the local wave number $k_{\text{loc}}$ of the patterns in the snake varies more strongly, and then typically the snake bifurcating from the primary pattern ($|\vec{k}|=k_{c}$) tends to terminate on the branch of periodic patterns which is most subcritical, i.e., here the grey branch, which has the smallest $B$ value in its fold.

This also indicates that for small $R$ we should not expect the primary BCCs to be the “most stable” pattern. Instead, the stability range of “distorted” BCCs like in Fig. 4(d) may extend to significantly lower $B$ than that of the primary BCCs. In Fig. 4 (like in Fig. 3 and in similar figures below) we have a quasi 1D situation due to the small domain size in $x$ and $y$. Moreover, also the domain size in $z$ is not really large, and together with the Neumann BCs this restricts the allowed wave vectors to a still rather small set. Over larger domains, the variety of (stable) patterns and associated possible snaking away from onset rather quickly becomes excessive, and for instance the stability ranges of different distorted BCCs will deserve a dedicated study.

3.2 Fronts between BCCs and tubes

In §2.3 we explained that for $R=0.75$ the amplitude equations (16) also predict mixed mode branches (orange branch in Fig. 1(a)), which suggests the existence of fronts between BCCs and tubes, cf. Fig.2(b). However, these occur at ${\cal O}(1)$ distance from onset, and hence such predictions should be taken with caution. In Fig. 5(a) we show the BCCs (black), tubes (magenta) and mixed modes (orange) for (1) with $R=0.75$ over the cube $\Omega=(-l,l)^{3}$, $l=\sqrt{2}\pi/k_{c}$. This confirms the predictions from Fig. 1 over this small domain, and we may extend these periodic patterns over the boundaries to obtain the same patterns and branches over larger domains. However, it turns out that even the reliable continuation of the BCC branch to $B-B_{c}={\cal O}(1)$ over extended domains is a delicate task, and requires fine meshes and strict settings for the algorithm pmcont designed to mitigate undesired branch switching, see [Uec19a, §3.3]. See also the remarks at the end of §3.1.

Therefore we proceed differently to explore the range $B\approx\tilde{B}_{M}$ for fronts between BCCs and tubes over long and slender bars, aiming to find snaking branches of fronts between BCCs and tubes. Figure 5(b) shows an initial condition (IC) of the form

[TABLE]

composed of the primary BCCs above the interface at $z=5$, and tubes below, while $u_{2}$ is simply set to the homogeneous value $u_{2}=u_{2}^{*}=B/A$. The domain is $\Omega=(-l,l)^{2}\times(-l_{z},l_{z})$, where $l=\pi/(\sqrt{2}k_{c})$ and $l_{z}=12l$. Additionally, we choose rather carefully (see below) the value $B=3.3633$. Starting with this initial condition, direct numerical simulation (DNS) slowly decreases the residual ${\rm Res}(u(t))$ defined in (29), see (c1). However, this decrease is in general not monotonous, and the BCCs actually change their wave vector. Nevertheless, after the transient (at, e.g., $t=25$) we can run a Newton loop on the stationary problem, and converge to the solution $u_{f}$ illustrated in (d). This is a (stable, as it is essentially obtained from DNS) stationary front between BCCs on top and tubes at the bottom, but similarly to Fig. 4(d), the BCCs are clearly not the primary BCCs belonging to (7), but (rather strongly) distorted, i.e., of the form (30) with $\tilde{\kappa}=1.5\kappa$. Next we continue $u_{f}$ in $B$, and obtain the (narrow and short) snake shown in (c2). In one direction (brown part), the spots recede (sample plots (e,f)) as the parameter varies, and in the other direction (red part), the spots expand (sample plot (g)). In both directions, the branch eventually reconnects to the mixed mode branch between the tubes and the BCCs with $\tilde{\kappa}=1.5\kappa$.

The snake in Fig. 5 is rather narrow, and the starting point was obtained by a careful choice of $B$ for the DNS. In Fig. 6 we illustrate the “typical” behavior of DNS for ICs of the form (33), which also explains the idea how to find $B$ for Fig. 5. We use the same domain and IC as in Fig. 5. For $B=3.6$, in (a,b), the initial dynamics is very similar to that in Fig. 5, i.e., the solution evolves towards a (distorted) BCC-tubes fronts. However, once the solution is “near the snake” from Fig. 5, the BCC part continues to grow in time. If we were close enough to the snake of steady fronts, on a sufficiently large domain, then we would expect “stick–slip” motion. See, e.g., [BK06, §III.B] where this is analyzed semi–analytically for fronts between patterns and the trivial solution outside the pinning (snaking) region in the (1D) Swift–Hohenberg equation. The motion is slow when the moving front passes near a steady solution at a fold of the corresponding snake, and afterwards moves quickly to near the next fold. The transition time from one fold to the next can be formally derived to be ${\cal O}(\delta^{-1/2})$ with $\delta$ the distance from the snaking region, and the associated full formula shows excellent agreement with the numerics in [BK06]. See also [Llo19, Llo20] for numerical analysis of the depinning of fronts in the planar Swift–Hohenberg equation. Here, the domain is not quite long enough, and for clarity we chose $B$ not very close to the snake, and thus we do not really see the stick–slip effect, but just roughly periodic variations in the residual, and for a rather short transient.

If we use the solution from $t=200$ to start a Newton loop for the steady problem, then this gives convergence to the (distorted) BCC solution. Alternatively, continuing the DNS we also converge to this BCC after a very long transient. On the other hand, in (c,d) we choose $B=3.3$ to the left of the snake, and obtain convergence to the tubes, and the same happens (faster) at the Maxwell point prediction $\tilde{B}_{M}\approx 3.21$. If snaking branches containing stable steady fronts exist, then such simulations give a hint for the right parameters to find them, and that is how we found the $B$ value for Fig. 5 with some trial and error. In particular, the Maxwell point prediction $\tilde{B}_{M}\approx 3.21$ from Table 1 was rather “far off”. A certain deviation was expected a priori as we are at ${\cal O}(1)$ distance from criticality. Additionally, and a posteriori, we see that $\tilde{B}_{M}\approx 3.21$ was irrelevant as it is the Maxwell point prediction for fronts between BCCs and tubes, and not the $\tilde{\kappa}=1.5\kappa$ distorted BCCs and tubes obtained in the DNS. We can, e.g., a posteriori change the $z$ wave number of the BCCs in the IC to start closer to $u_{f}$, which then also allows to directly go to $u_{f}$ by a Newton loop for the steady problem. However, Fig. 5 illustrates that very good initial guesses are often not necessary, and instead rather poor initial guesses can be first improved by DNS.

In summary, the predictions from §2.3 are useful as they motivate the search for fronts and give hints for good parameter regimes. Of course, finding such fronts via DNS needs the existence of steady localized patterns of the desired form, and is easier and more robust if the (desired) snake is wide. On the other hand, in pattern forming systems such as (1) we may expect a (large) variety of (stable) steady patterns far from onset, and this increases the chances to converge towards some localized patterns. One more example is given in the next section, where at $R=1$ lamellas enter the game.

3.3 $R=1$: The comeback of the Lamellas

In Fig. 7 we illustrate some typical results for (1) at $R=1$. In (a) we show the BD over a short bar. The BCCs and tubes now both bifurcate in supercritical pitchforks, with the tubes stable. Additionally, we show the next two bifurcating branches. The orange branch consists of elongated BCCs, and the green branch are $\kappa=1.5k_{c}$ lamellas

[TABLE]

which we did not consider in §2 as they do not bifurcate at the primary bifurcation, but (on the given domain) at the third bifurcation point from the $u=u^{*}$ branch. These lamellas become stable at $B=B_{l}\approx 4.06$, on this domain, but similarly also on much longer domains.

Thus, we now have a bistable range between tubes and lamellas, and the lamellas turn out to play a crucial role in the DNS, as illustrated in (c,d). We set $B=4.2$, and again use an initial condition of type (33). Though there are no lamellas in the initial condition, the solution initially (till $t=100$, say), relaxes to a ’double–front’ from lamellas to tubes with a distinct ’cubes-like’ interface in between. This front then propagates downwards in a roughly periodic fashion (see the lower time series in (c)) but with essentially fixed shape.444The fixed shape appears to be another effect of locking due to the periodic pattern, as for double fronts between homogeneous states one would generically expect the middle state to expand or shrink. See, e.g., [CM99] for another striking example of such double-fronts, namely a ’roll belt’ ahead of hexagons invading the zero solution in a damped Kuramoto-Sivashinsky equation. A similar behavior occurs at other values of $B$ ($B>B_{l}$) and other initial conditions, i.e., for $R=1$ the lamellas always win on domains of the type considered here.

4 Discussion

We numerically studied patterns in the 3D Brusselator over boxes with Neumann BCs, specifically aiming at snaking branches of steady fronts between patterns, which can also be seen as approximations of localized patterns. The basic idea is as in 1D and 2D, namely to look for bifurcations from subcritical branches of patterns, or from mixed mode branches. However, the numerical challenges are significant. In 3D, pattern forming systems allow a much larger variety of steady patterns than in 1D or 2D. The problem is already quite complicated near onset, but on “nice domains” (e.g., small cuboids with Neumann BCs) the main branches can be found from (simplified and reduced) amplitude equations. Farther from onset, there typically is a multitude of patterns, in particular if the domain is not very small, and this makes (numerical) continuation and bifurcation analysis (much) harder than in 1D or 2D, essentially due to uncontrolled branch jumping in the continuation.

Therefore we focused on the simplest situations of small domains in the form of long but slender rods, with an underlying BCC lattice, and thus on specific localized patterns, namely localized BCCs, and fronts between BCCs and tubes (or localized BCCs embedded in a background of tubes or vice versa). Over larger domains we expect a huge variety of additional localized patterns, similar to but still extending the 2D examples in, e.g., [UW14, Wet18]. However, even for the minimal domains used, the search for localized patterns via continuation and bifurcation (as we did for the localized BCCs in Figs. 3 and 4) is rather delicate. Thus, to obtain starting points for the continuation of BCC-to-tubes fronts we found it more robust and efficient to use DNS, with the Maxwell point of the amplitude system as a guide for promising parameter regimes. Finally, we gave one example of a moving front between lamellas and tubes. It should be interesting to see whether any such localized patterns can be realized experimentally, as, e.g., the 3D Turing patterns presented in [BJVE11].