Weakly Nonlinear Analysis of Peanut-Shaped Deformations for Localized Spots of Singularly Perturbed Reaction-Diffusion Systems

Spatially localized patterns arise in a diverse range of applications including, the ferrocyanide-iodate-sulphite (FIS) reaction (cf. [14], [15]), the chloride-dioxide-malonic acid reaction (cf. [6]), certain electronic gas discharge systems [1], fluid-convection phenomena [10], and the emergence of plant root hair cells mediated by the plant hormone auxin (cf. [2]), among others. One qualitatively novel feature in many of these settings is the observation that spatially localized spot-type patterns can undergo a seemingly spontaneous self-replication process.

Many of these observed localized patterns, most notably those in chemical physics and biology, are modeled by nonlinear reaction-diffusion (RD) systems. In [25], where the two-component Gray-Scott RD model was used to qualitatively model the FIS reaction, full PDE simulations revealed a wide variety of highly complex spatio-temporal localized patterns including, self-replicating spot patterns, stripe patterns, and labyrinthian space-filling curves (see also [27], [19] and [20]). This numerical study showed convincingly that in the fully nonlinear regime a two-component RD system with seemingly very simple reaction kinetics can admit highly intricate solution behavior, which cannot be described by a conventional Turing stability analysis (cf. [31]) of some spatially uniform base state. For certain three-component RD systems in the limit of small diffusivity, Nishiura et. al. (cf. [23], [29]) showed from PDE simulations and a weakly nonlinear bifurcation analysis that a subcritical peanut-shaped instability of a localized radially symmetric spot plays a key role in understanding the dynamics of traveling spot solutions. These previous studies, partially motivated by the pioneering numerical study of [25], have provided the impetus for developing new theoretical approaches to analyze some of the novel dynamical behaviors and instabilities of localized patterns in RD systems in the “far-from-equilibrium” regime [21]. A survey of some novel phenomena and theoretical approaches associated with localized pattern formation problems are given in [36], [21] and [10]. The main goal of this paper is to use a weakly nonlinear analysis to study the onset of spot self-replication for certain two-component RD systems in the so-called “semi-strong” regime, characterized by a large diffusivity ratio between the solution components.

The derivation of amplitude, or normal form, equations using a multi-scale perturbation analysis is a standard approach for characterizing the weakly nonlinear development of small amplitude patterns near bifurcation points. It has been used with considerable success in physical applications, such as in hydrodynamic stability theory and materials science (cf. [5], [38]) and in biological and chemical modeling through RD systems defined in planar spatial domains and on the sphere (cf. [38], [17], [18], [3]). However, in certain applications, the effectiveness of normal form theory is limited owing to the existence of subcritical bifurcations (cf. [3]) or the need for an extreme fine-tuning of the model parameters in order to be within the range of validity of the theory (cf. [43]). In contrast to the relative ease in undertaking a weakly nonlinear theory for an RD system near a Turing bifurcation point of the linearization around a spatially uniform or patternless state, it is considerably more challenging to implement such a theory for spatially localized steady-state patterns. This is owing to the fact that the linearization of the RD system around a spatially localized spot solution leads to a singularly perturbed eigenvalue problem in which the underlying linearized operator is spatially heterogeneous. In addition, various terms in the multi-scale expansion that are needed to derive the amplitude equation involve solving rather complicated spatially inhomogeneous boundary value problems. In this direction, a weakly nonlinear analysis of temporal amplitude oscillations (breathing instabilities) of 1-D spike patterns was developed for a class of generalized Gierer-Meinhardt (GM) models in [37] and for the Gray-Scott and Schnakenberg models in [9]. A criterion for whether these oscillations, emerging from a Hopf bifurcation point of the linearization, are subcritical or supercritical was derived. A related weakly nonlinear analysis for competition instabilities of 1-D steady-state spike patterns for the GM and Schnakenberg models, resulting from a zero-eigenvalue crossing of the linearization, was developed in [16]. Finally, for a class of coupled bulk-surface RD systems, a weakly nonlinear analysis for Turing, Hopf, and codimension-two Turing-Hopf bifurcations of a patterned base-state was derived in [24].

The focus of this paper is to develop and implement a weakly nonlinear theory to analyze branching behavior associated with peanut-shaped deformations of a locally radially symmetric steady-state spot solution for certain singularly perturbed RD systems. Previous numerical simulations of the Schnakenberg and Brusselator RD systems in [13], [28] and [32] (see also [30]) have indicated that a non-radially symmetric peanut-shape deformation of the spot profile can, in certain cases, trigger a fully nonlinear spot self-replication event. The parameter threshold for the onset of this shape deformation linear instability has been calculated in [13] and [28] for the Schnakenberg and Brusselator models, respectively. We will extend this linear theory by using a multi-scale perturbation approach to derive a normal form amplitude equation characterizing the local branching behavior associated with peanut-shaped instabilities of the spot profile. From a numerical evaluation of the coefficients in this amplitude equation we will show that a peanut-shaped instability of the spot profile is always subcritical for both the Schnakenberg and Brusselator models. This theoretical result supports the numerical findings in [13], [28] and [32] that a peanut-shaped instability of a localized spot is the trigger for a fully nonlinear spot-splitting event, and it solves an open problem discussed in the survey article [39].

The dimensionless Schnakenberg model in the two-dimensional unit disk $\Omega=\{\mathbf{x}:|\mathbf{x}|\leq 1\}$ is formulated as

$$v_{t}=\varepsilon^{2}\Delta v-v+uv^{2},\quad\quad\tau u_{t}=D\Delta u+a-\varepsilon^{-2}uv^{2},\quad\mathbf{x}\in\Omega\,,$$

(1.1)

with $\partial_{n}v=\partial_{n}u=0$ on $\partial\Omega$. Here $\varepsilon\ll 1$, $D=\mathcal{O}(1)$, $\tau=\mathcal{O}(1)$, and the constant $a>0$ is called the feed-rate. For a spot centered at the origin of the disk, the contour plot in Fig. 1 of $v$ at different times, as computed numerically from (1.1), shows a spot self-replication event as the feed-rate $a$ is slowly ramped above the threshold value $a_{c}\approx 8.6$. At this threshold value of $a$ the spot profile becomes unstable to a peanut-shaped deformation (see §3 for the linear stability analysis).

Rigorous analytical results for the existence and linear stability of localized spot patterns for the Schnakenberg model (1.1) in the large $D$ regime $D=\mathcal{O}(\nu^{-1})$, where $\nu={-1/\log\varepsilon}$, are given in [41] and for the related Gray-Scott model in [40] (see [42] for a survey of such rigorous results). For the regime $D=\mathcal{O}(1)$, a hybrid analytical-numerical approach, which has the effect of summing all logarithmic terms in powers of $\nu$, was developed in [13] to construct quasi-equilibrium patterns, to analyze their linear stability properties, and to characterize slow spot dynamics. An extension of this hybrid methodology applied to other RD systems was given in [4], [28], [30], [32] and [2], and is surveyed in [39].

We remark that the mechanism underlying the self-replication of 1-D localized patterns is rather different than the more conventional symmetry-breaking mechanism that occurs in 2-D. In a one-dimensional domain, the self-replication behavior of spike patterns has been interpreted in terms of a nearly-coinciding hierarchical saddle-node global bifurcation structure of branches of multi-spike equilibria, together with the existence of a dimple-shaped eigenfunction of the linearization near the saddle-node point (see [22], [8], [7], [12], [35], [11], [19], [27] and the references therein).

The outline of this paper is as follows. For the Schnakenberg RD model (1.1), in §2 we use the method of matched asymptotic expansions to construct a steady-state, locally radially symmetric, spot solution centered at the origin of the unit disk. In §3 we perform a linear stability analysis for non-radially symmetric perturbations of this localized steady-state, and we numerically compute the threshold conditions for the onset of a peanut-shaped instability of a localized spot. Although much of this steady-state and linear stability theory has been described previously in [13], it provides the required background context for describing the new weakly nonlinear theory in §4. More specifically, in §4 we develop and implement a weakly nonlinear analysis to characterize the branching behavior associated with peanut-shaped instabilities of a localized spot. From a numerical evaluation of the coefficients in the resulting normal form amplitude equation we show that a peanut-shaped deformation of a localized spot is subcritical. By using the bifurcation software pde2path [34], the weakly nonlinear theory is validated in §4.1 by numerically computing an unstable non-radially symmetric steady-state spot solution branch that emerges from the peanut-shaped linear stability threshold of a locally radially symmetric spot solution. In §5 we perform a similar multi-scale asymptotic reduction to derive an amplitude equation characterizing the weakly nonlinear development of peanut-shaped deformations of a localized spot for the Brusselator RD model, originally introduced in [26]. From a numerical evaluation of the coefficients in this amplitude equation, which depend on a parameter in the Brusselator reaction-kinetics, it is shown that peanut-shaped linear instabilities are always subcritical. This theoretical result predicting subcriticality is again validated using pde2path [34]. In §6 we summarize a few key qualitative features of our hybrid analytical-numerical approach to derive the amplitude equation, and we discuss a few possible extensions of this work.

We use the method of matched asymptotic expansions to construct a steady-state single spot solution centered at $\mathbf{x}_{0}=\mathbf{0}$ in the unit disk. In the inner region near $\mathbf{x}=0$, we set

$$v=\sqrt{D}\,V(\mathbf{y})\,,\quad u={U(\mathbf{y})/\sqrt{D}}\,,\quad\mbox{where}\quad\mathbf{y}=\varepsilon^{-1}\mathbf{x}\,.$$

(2.2)

In the inner region, for $\mathbf{y}\in\mathbb{R}^{2}$, the steady-state problem is

$$\Delta_{\mathbf{y}}V-V+UV^{2}=0\,,\qquad\Delta_{\mathbf{y}}U-UV^{2}+\frac{a\varepsilon^{2}}{\sqrt{D}}=0\,.$$

(2.3)

We seek a radially symmetric solution in the form $V=V_{0}(\rho)+o(1)$ and $U=U_{0}(\rho)+o(1)$, where $\rho=|\mathbf{y}|$. Upon neglecting the $\mathcal{O}(\varepsilon^{2})$ terms, we obtain the core problem

$$\begin{split}&\Delta_{\rho}V_{0}-V_{0}+U_{0}V_{0}^{2}=0\,,\quad\Delta_{\rho}U_{0}-U_{0}V_{0}^{2}=0\,,\quad\mbox{where }\Delta_{\rho}\equiv\partial_{\rho\rho}+\rho^{-1}\partial_{\rho}\,,\\ &U_{0}^{\prime}(0)=V_{0}^{\prime}(0)=0\,;\quad V_{0}\to 0\,,\quad U_{0}\sim S\log\rho+\chi(S)+o(1)\,,\quad\mbox{as}\,\,\,\rho\to\infty\,,\end{split}$$

(2.4)

In particular, we must allow $U_{0}$ to have far-field logarithmic growth whose strength is characterized by the parameter $S>0$, which will be determined below (see (2.8)) in terms of the feed rate parameter $a$. The $\mathcal{O}(1)$ term in the far-field behavior depends on $S$, and is denoted by $\chi(S)$. It must be computed numerically from the BVP (2.4). A plot of the numerically-computed $\chi$ versus $S$ is shown in Fig. 2. By integrating the $U_{0}$ equation in (2.4), we obtain the identity that

$\displaystyle S=\int_{0}^{\infty}U_{0}V_{0}^{2}\,\rho\,\mathrm{d}\rho\,.$

(2.5)

In the limit $\varepsilon\to 0$, the term $\varepsilon^{-2}uv^{2}$ in the outer region can be represented, in the sense of distributions, as a Dirac source term using the correspondence rule

$$\varepsilon^{-2}uv^{2}\to 2\pi\sqrt{D}\left(\int_{0}^{\infty}U_{0}V_{0}^{2}\rho\,d\rho\right)\,\delta(\mathbf{x})=2\pi S\sqrt{D}\,\delta(\mathbf{x})\,,$$

(2.6)

where (2.5) was used. As a result, the outer problem for $u$ in (1.1) is

$$\Delta u=-\frac{a}{D}+\frac{2\pi S}{\sqrt{D}}\delta(\mathbf{x})\,,\quad\mathbf{x}\in\Omega\,;\qquad\partial_{n}u=0\,,\quad\mathbf{x}\in\partial\Omega\,.$$

(2.7)

We integrate (2.7) over the disk and use the Divergence theorem and $|\Omega|=\pi$, to get

$$S=\frac{a|\Omega|}{2\pi\sqrt{D}}=\frac{a}{2\sqrt{D}}\,.$$

(2.8)

To represent the solution to (2.7) we introduce the Neumann Green’s function $G(\mathbf{x};\mathbf{x}_{0})$ for the unit disk, which is defined uniquely by

$$\begin{split}&\Delta G=\frac{1}{|\Omega|}-\delta(\mathbf{x}-\mathbf{x}_{0})\,,\quad\mathbf{x}\in\Omega\,;\qquad\partial_{n}G=0\,,\quad\mathbf{x}\in\partial\Omega;\\ &\int_{\Omega}G\,d\mathbf{x}=0\,,\quad G\sim-\frac{1}{2\pi}\log|\mathbf{x}-\mathbf{x}_{0}|+R_{0}+o(1)\,,\quad\mbox{as}\quad\mathbf{x}\to\mathbf{x}_{0}\,,\end{split}$$

(2.9)

where $R_{0}$ is the regular part of the Green’s function. The solution to (2.7) is

$$u=-\frac{2\pi S}{\sqrt{D}}G(\mathbf{x};\mathbf{0})+\bar{u}\,,$$

(2.10)

where $\bar{u}$ is a constant to be determined below by asymptotic matching the inner and outer solutions. The Neumann Green’s function with singularity at the origin is

$$G(\mathbf{x};\mathbf{0})=-\frac{1}{2\pi}\log|\mathbf{x}|+\frac{|\mathbf{x}|^{2}}{4\pi}-\frac{3}{8\pi}\,.$$

(2.11)

Therefore, by using (2.11) in (2.10), the outer solution $u$ satisfies

$$u=\frac{S}{\sqrt{D}}\log|\mathbf{x}|+\frac{3S}{4\sqrt{D}}+\bar{u}+\mathcal{O}(|\mathbf{x}|^{2})\,,\quad\mbox{as}\quad\mathbf{x}\to\mathbf{0}\,.$$

(2.12)

By using the far-field behavior of the inner solution $U$ in (2.4), we obtain for $\rho\gg 1$ that

$$u=\frac{U}{\sqrt{D}}\sim\frac{1}{\sqrt{D}}\left[S\log|\mathbf{x}|+\frac{S}{\nu}+\chi(S)\right]\,,\quad\mbox{where}\quad\nu\equiv-\frac{1}{\log\varepsilon}\,.$$

(2.13)

From an asymptotic matching of (2.12) and (2.13), we identity $\bar{u}$ as

$$\bar{u}=\frac{1}{\sqrt{D}}\left(\chi(S)+\frac{S}{\nu}-\frac{3S}{4}\right)\,.$$

(2.14)

Upon substituting (2.14) and (2.11) into (2.10) we conclude that the outer solution is

$$u=\frac{1}{\sqrt{D}}\left(S\log|\mathbf{x}|-\frac{S|\mathbf{x}|^{2}}{2}+\chi(S)+\frac{S}{\nu}\right)\,,\qquad\mbox{where}\quad S=\frac{a}{2\sqrt{D}}\,.$$

(2.15)

In this section, we perform a linear stability analysis of the steady-state one-spot solution in the unit disk. For convenience, we will represent a column vector by the notation $(u_{1},u_{2})\mbox{ or }\begin{pmatrix}u_{1}\\ u_{2}\end{pmatrix}.$ For a steady-state spot centered at the origin, we will formulate the linearized stability problem in the quarter disk, defined by $\Omega_{+}=\{\mathbf{x}=(x,y)\,:\,\,|\mathbf{x}|<1,\,x\geq 0,\,y\geq 0\}$.

Let $v_{e},\,u_{e}$ be the steady-state spot solution centered at the origin. We introduce the perturbation

$$v=v_{e}+e^{\lambda t}\phi\,,\quad u=u_{e}+e^{\lambda t}\eta\,,$$

(3.17)

into (1.1) and linearize. This leads to the singularly perturbed eigenvalue problem

$$\varepsilon^{2}\Delta\phi-\phi+2u_{e}v_{e}\phi+v_{e}^{2}\eta=\lambda\phi\,,\qquad D\Delta\eta-\varepsilon^{-2}(2u_{e}v_{e}\phi+v_{e}^{2}\eta)=\tau\lambda\eta\,,$$

(3.18)

with $\partial_{n}\phi=\partial_{n}\eta=0$ on $\partial\Omega$.

In the inner region near $\mathbf{x}=\mathbf{0}$ we introduce

$$\begin{pmatrix}\phi\\ \eta\end{pmatrix}=\mathrm{Re}(e^{im\theta})\begin{pmatrix}\Phi(\rho)\\ N(\rho)/D\end{pmatrix}\,,\quad\mbox{where}\quad\rho=|\mathbf{y}|=\varepsilon|\mathbf{x}|\,,\quad\theta=\mbox{arg}(\mathbf{y})\,,$$

(3.19)

with $m=2,3,\ldots$. With $v_{e}\sim\sqrt{D}V_{0}$ and $u_{e}\sim{U_{0}/\sqrt{D}}$, we neglect the $\mathcal{O}(\varepsilon^{2})$ terms to obtain the eigenvalue problem

$$\mathcal{L}_{m}\begin{pmatrix}\Phi\\ N\end{pmatrix}+\begin{pmatrix}-1+2U_{0}V_{0}&V_{0}^{2}\\ -2U_{0}V_{0}&-V_{0}^{2}\end{pmatrix}\begin{pmatrix}\Phi\\ N\end{pmatrix}=\lambda\begin{pmatrix}1&0\\ 0&0\end{pmatrix}\begin{pmatrix}\Phi\\ N\end{pmatrix}\,,$$

(3.20)

where the operator $\mathcal{L}_{m}$ is defined by $\mathcal{L}_{m}\mathbf{\Phi}=\partial_{\rho\rho}\mathbf{\Phi}+\rho^{-1}\partial_{\rho}\mathbf{\Phi}-m^{2}\rho^{-2}\mathbf{\Phi}$. We seek eigenfunctions of (5.79) with $\Phi\to 0$ and $N\to 0$ as $\rho\to\infty$. An unstable eigenvalue of this spectral problem satisfying $\mbox{Re}(\lambda)>0$ corresponds to a non-radially symmetric spot-deformation instability.

For each angular mode $m=2,3,\ldots$, the eigenvalue $\lambda_{0}$ of (3.20) with the largest real part is a function of the source strength $S$. To determine $\lambda_{0}$ we discretize (3.20) as done in [13] to obtain a finite-dimensional generalized eigenvalue problem. We calculate $\lambda_{0}$ numerically from this discretized problem, with the results shown in the right panel of Fig. 3. In the left panel of Fig. 3 we show the quarter-disk geometry.

For the angular mode $m=2$, we find that $\mathrm{Re}(\lambda_{0})=0$ when $S=S_{c}\approx 4.3022$, which agrees with the result first obtained in [13]. At this critical value of $S$, we define

$$V_{c}(\rho)\equiv V_{0}(\rho\,;S_{c}),\quad U_{c}(\rho)\equiv U_{0}(\rho\,;S_{c})\,,\quad M_{c}\equiv\begin{pmatrix}-1+2U_{c}V_{c}&V_{c}^{2}\\ -2U_{c}V_{c}&-V_{c}^{2}\end{pmatrix}\,,$$

(3.21)

so that there exists a non-trivial solution, labeled by $\mathbf{\Phi}_{c}\equiv(\Phi_{c},N_{c})$, to

$$\mathcal{L}_{2}\mathbf{\Phi}_{c}+M_{c}\mathbf{\Phi}_{c}=\mathbf{0}\,.$$

(3.22)

For $m=2$, we have that $\Phi_{c}\to 0$ exponentially as $\rho\to\infty$ and $N_{c}=\mathcal{O}(\rho^{-2})$ as $\rho\to\infty$. As such, we impose $\partial_{\rho}N_{c}\sim-2{N_{c}/\rho}$ for $\rho\gg 1$. Since (3.22) is a linear homogeneous system, the solution is unique up to a multiplicative constant. We normalize the solution to (3.22) using the condition

$$\int_{0}^{\infty}\Phi_{c}^{2}\,\rho\,d\rho=1\,.$$

(3.23)

A plot of the numerically computed inner solution $V_{c}$ and $U_{c}$ is shown in Fig. 4.

Next, for $S=S_{c}$, it follows that there exists a nontrivial solution $\mathbf{\Phi_{c}}^{*}=(\Phi_{c}^{*},N_{c}^{*})$ to the adjoint problem

$$\mathcal{L}_{2}\mathbf{\Phi_{c}^{*}}+M_{c}^{T}\mathbf{\Phi_{c}^{*}}=\mathbf{0}\,,\quad\Phi_{c}^{*}\to 0\,,\quad\partial_{\rho}N_{c}^{*}\sim-\frac{2N_{c}^{*}}{\rho}\quad\mbox{as}\quad\rho\to\infty\,,$$

(3.24)

for which we impose the convenient normalization condition $\int_{0}^{\infty}(\Phi_{c}^{*})^{2}\rho\,d\rho=1$.

In Fig. 5 we plot the numerically computed nullvector $\Phi_{c}$ and $N_{c}$, satisfying (3.22), as well as the adjoint $\Phi_{c}^{*}$ and $N_{c}^{*}$, satisfying (3.24).

In this section we derive the amplitude equation associated with the peanut-splitting linear stability threshold for the Schnakenberg model. This amplitude equation will show that this spot shape-deformation instability is subcritical.

To do so, we first introduce a small perturbation around the linear stability threshold $S_{c}$ given by $S=S_{c}+\kappa\sigma^{2}$, where $\kappa=\pm 1$. In this way, the obtain the Taylor expansion $\chi(S)=\chi(S_{c})+\kappa\chi^{\prime}(S_{c})\sigma^{2}+\mathcal{O}(\sigma^{4})$. Then, we introduce a slow time scale $T=\sigma^{2}t$. As such, the inner problem in terms of $V={v/\sqrt{D}}$ and $U=\sqrt{D}u$ for $\mathbf{y}\in\mathbb{R}^{2}$ is

	$$\sigma^{2}V_{T}=\Delta_{\mathbf{y}}V-V+UV^{2}\,,\qquad\frac{\sigma^{2}\varepsilon^{2}\tau}{D}U_{T}=\Delta_{\mathbf{y}}U-UV^{2}+\frac{a\varepsilon^{2}}{\sqrt{D}}\,,$$		(4.34a)
for which we impose $V\to 0$ exponentially as $\rho\to\infty$, while
	$$U\sim\left(S_{c}+\kappa\sigma^{2}\right)\log\rho+\chi(S_{c})+\sigma^{2}\left[\kappa\chi^{\prime}(S_{c})+\mathcal{O}(1)\right]+\ldots\,,\quad\mbox{as}\quad\rho=|\mathbf{y}|\to\infty\,.$$		(4.34b)

In (4.34), we expand $V=V(\rho,\phi,T)$ and $U=U(\rho,\phi,T)$ as

$$V=V_{0}+\sigma V_{1}+\sigma^{2}V_{2}+\sigma^{3}V_{3}+\ldots\,,\quad U=U_{0}+\sigma U_{1}+\sigma^{2}U_{2}+\sigma^{3}U_{3}+\ldots\,,$$

(4.35)

where $V_{0},\,U_{0}$ is the radially symmetry core solution, satisfying (2.4). Furthermore, we assume that

$$\sigma^{3}\gg\mathcal{O}(\varepsilon^{2})\,,$$

(4.36)

so that the $\mathcal{O}(\varepsilon^{2})$ terms in (4.34a) are asymptotically smaller than terms of order $\mathcal{O}(\sigma^{k})$ for $k\leq 3$.

We then substitute (4.35) into (4.34) and collect powers of $\sigma$. From the $\mathcal{O}(1)$ terms, we obtain that $V_{0}$ and $U_{0}$ satisfy

	$\displaystyle\Delta_{\rho}$	$\displaystyle V_{0}-V_{0}+U_{0}V_{0}^{2}=0\,,\quad\Delta_{\rho}U_{0}-U_{0}V_{0}^{2}=0\,,$		(4.37a)
		$\displaystyle V_{0}\to 0,\quad U_{0}\sim S_{c}\log\rho+\mathcal{O}(1)\,,\quad\mbox{as}\quad\rho\to\infty\,.$		(4.37b)

From the far-field condition (4.37b), we can identify that $V_{0}$ and $U_{0}$ are the core solution with $S=S_{c}$. In other words, we have

$$V_{0}=V_{c}(\rho)\,,\quad U_{0}=U_{c}(\rho)\,.$$

(4.38)

From collecting $\mathcal{O}(\sigma)$ terms, and setting $V_{0}=V_{c}$ and $U_{0}=U_{c}$, we find that $\mathbf{V}_{1}=(V_{1},U_{1})$ satisfies

$$\Delta_{\mathbf{y}}\mathbf{V}_{1}+M_{c}\,\mathbf{V}_{1}=\mathbf{0}\,,\quad\mbox{where}\quad M_{c}=\begin{pmatrix}-1+2U_{c}V_{c}&V_{c}^{2}\\ -2U_{c}V_{c}&-V_{c}^{2}\end{pmatrix}\,.$$

(4.39)

We conclude that $\mathbf{V}_{1}$ is related to the eigenfunction solution to (3.22). We introduce the amplitude function $A=A(T)$, while writing $\mathbf{V}_{1}$ as

$$\mathbf{V}_{1}=A\cos(2\phi)\begin{pmatrix}\Phi_{c}\\ N_{c}\end{pmatrix}\,,$$

(4.40)

where $\Phi_{c}$ and $N_{c}$ satisfy (3.22) with normalization (3.23).

By collecting $\mathcal{O}(\sigma^{2})$ terms we readily obtain that $\mathbf{V}_{2}=(V_{2},U_{2})$ on $\mathbf{y}\in\mathbb{R}^{2}$ satisfies

	$$\Delta_{\mathbf{y}}\mathbf{V}_{2}+M_{c}\mathbf{V}_{2}=F_{2}\,\mathbf{q}\,,$$		(4.41a)
where we have defined $F_{2}$ and $\mathbf{q}$ by
	$$F_{2}\equiv 2V_{c}V_{1}U_{1}+U_{c}V_{1}^{2}\,,\qquad\mathbf{q}\equiv\begin{pmatrix}-1\\ 1\end{pmatrix}\,.$$		(4.41b)

By using (4.40) for $V_{1}$ and $U_{1}$, together with the identity $2\cos^{2}\phi=1+\cos(2\phi)$, we can write $F_{2}$ as

$$F_{2}=A^{2}F_{20}+A^{2}F_{20}\cos(4\phi)\,,\qquad F_{20}=\frac{1}{2}\left(U_{c}\Phi_{c}^{2}+2V_{c}\Phi_{c}N_{c}\right)\,.$$

(4.42)

This suggests a decomposition of the solution to (4.41a) in the form

$$\mathbf{V}_{2}=\mathbf{V}_{20}(\rho)+A^{2}\,\mathbf{V}_{24}(\rho)\cos(4\phi)\,,$$

(4.43)

where the problems for $\mathbf{V}_{20}$ and $\mathbf{V}_{24}$ are formulated below.

Firstly, we define $\mathbf{V}_{24}=(V_{24},U_{24})$ to be the radially symmetric solution to

	$$\mathcal{L}_{4}\mathbf{V}_{24}+M_{c}\mathbf{V}_{24}=F_{20}\,\mathbf{q}\,,$$		(4.44a)
where $\mathcal{L}_{m}\mathbf{V}_{24}=\partial_{\rho\rho}\mathbf{V}_{24}+\rho^{-1}\partial_{\rho}\mathbf{V}_{24}-m^{2}\rho^{-2}\mathbf{V}_{24}$, for which we can impose that
	$$V_{24}\to 0\,,\quad U_{24}=\mathcal{O}(\rho^{-4})\,\quad\longrightarrow\quad\partial_{\rho}\,U_{24}\sim-\frac{4}{\rho}\,U_{24}\,,\quad\mbox{as}\quad\rho\to\infty\,.$$		(4.44b)

Next, we define $\mathbf{V}_{20}=(V_{20},U_{20})$ to be the solution to

	$$\Delta_{\rho}\mathbf{V}_{20}+M_{c}\mathbf{V}_{20}=A^{2}F_{20}\,\mathbf{q}\,.$$		(4.45a)
We can impose $V_{20}\to 0$ exponentially as $\rho\to\infty$. As indicated in (4.34b), we have
	$$U_{2}\sim\kappa\log\rho+\mathcal{O}(1)\,,\quad\mbox{as}\quad\rho\to\infty\,.$$		(4.45b)

Since $U_{24}=\mathcal{O}(\rho^{-4})\ll 1$ as $\rho\to\infty$, we must have $U_{20}\sim\kappa\log\rho+\mathcal{O}(1)$.

Next, we decompose $\mathbf{V}_{20}$ by first observing that $\mathbf{W}_{2H}\equiv(\partial_{S}V_{c},\partial_{S}U_{c})$ is a radial solution to the homogeneous problem

$$\Delta_{\rho}\mathbf{W}_{2H}+M_{c}\mathbf{W}_{2H}=\mathbf{0}\,,\quad\mathbf{W}_{2H}\sim\left(0,\log\rho+\chi^{\prime}(S_{c})\right)\,,\quad\mbox{as}\quad\rho\to\infty\,.$$

(4.46)

This suggests that it is convenient to introduce the following decomposition to isolate the two sources of inhomogeneity in (4.45):

$$\mathbf{V}_{20}=\kappa\mathbf{W}_{2H}+A^{2}\hat{\mathbf{V}}_{20}\,,$$

(4.47)

where $\hat{\mathbf{V}}_{20}=(\hat{V}_{20},\hat{U}_{20})$ is taken to be the radial solution to

$$\Delta_{\rho}\hat{\mathbf{V}}_{20}+M_{c}\hat{\mathbf{V}}_{20}=F_{20}\,\mathbf{q}\,,\quad\hat{V}_{20}\to 0\,,\quad\partial_{\rho}\hat{U}_{20}\to 0\,,\quad\mbox{as}\quad\rho\to\infty\,.$$

(4.48)

In Appendix A we discuss in detail the derivation of the far-field condition for $\hat{U}_{20}$ imposed in (4.48). Moreover, since $\hat{U}_{20}\to U_{20\infty}\neq 0$ as $\rho\to\infty$, at the end of Appendix A we show how this fact can be accounted for in a simple modification of the outer solution given in (2.15).

In view of (4.47) and (4.44), the solution to (4.41a), as written in (4.43), is

$$\mathbf{V}_{2}=\kappa\mathbf{W}_{2H}+A^{2}\left[\hat{\mathbf{V}}_{20}+\mathbf{V}_{24}\cos(4\phi)\right]\,.$$

(4.49)

In the left and right panels of Fig. 6 we plot the numerically computed solution to (4.48) and (4.44), respectively.

The solvability condition, which yields the amplitude equation for $A$, arises from the $\mathcal{O}(\sigma^{3})$ problem. At this order, we find that $\mathbf{V}_{3}=(V_{3},U_{3})$ satisfies

	$$\Delta_{\mathbf{y}}\mathbf{V}_{3}+M_{c}\mathbf{V}_{3}=F_{3}\,\mathbf{q}+\partial_{T}V_{1}\,\mathbf{e}_{1}\,,$$		(4.50a)
where we have defined $F_{3}$ and $\mathbf{e}_{1}$ by
	$$F_{3}\equiv 2V_{c}V_{1}U_{2}+U_{1}V_{1}^{2}+2V_{c}U_{1}V_{2}+2U_{c}V_{1}V_{2}\,,\qquad\mathbf{e}_{1}\equiv\begin{pmatrix}1\\ 0\end{pmatrix}\,.$$		(4.50b)

Upon substituting (4.40) and (4.49) into $F_{3}$, we can write $F_{3}$ in (4.50b) in terms of a truncated Fourier cosine expansion as

	$$F_{3}=(\kappa g_{1}A+g_{2}A^{3})\cos(2\phi)+g_{3}A^{3}\cos(6\phi)\,,$$			(4.51a)
where $g_{1},\,g_{2}$ and $g_{3}$ are defined by
	$\displaystyle g_{1}$	$\displaystyle=2\Phi_{c}\partial_{S}(V_{c}U_{c})+N_{c}\partial_{S}(V_{c}^{2})\,,$		(4.51b)
	$\displaystyle g_{2}$	$\displaystyle=2V_{c}\Phi_{c}\hat{U}_{20}+V_{c}\Phi_{c}U_{24}+\frac{3}{4}\Phi_{c}^{2}N_{c}+(V_{c}N_{c}+U_{c}\Phi_{c})(2\hat{V}_{20}+V_{24})\,,$		(4.51c)
	$\displaystyle g_{3}$	$\displaystyle=\frac{1}{4}N_{c}\Phi_{c}^{2}+V_{c}\Phi_{c}U_{24}+(V_{c}N_{c}+U_{c}\Phi_{c})V_{24}\,.$		(4.51d)

In this way, the solution $\mathbf{V}_{3}=(V_{3},U_{3})$ to (4.50a) satisfies

$$\Delta\mathbf{V}_{3}+M_{c}\mathbf{V}_{3}=(\kappa g_{1}A+g_{2}A^{3})\cos(2\phi)\,\mathbf{q}+g_{3}A^{3}\cos(6\phi)\,\mathbf{q}+A^{\prime}\Phi_{c}\cos(2\phi)\,\mathbf{e}_{1}\,,$$

(4.52)

where $A^{\prime}\equiv{dA/dT}$. The right-hand side of this expression suggests that we decompose $\mathbf{V}_{3}$ as

	$$\mathbf{V}_{3}=\mathbf{W}_{2}(\rho)\cos(2\phi)+\mathbf{W}_{6}(\rho)\cos(6\phi)\,,$$			(4.53a)
so that from (4.52) we obtain that $\mathbf{W}_{2}$ and $\mathbf{W}_{6}$ are radial solutions to
	$\displaystyle\mathcal{L}_{2}\mathbf{W}_{2}+M_{c}\mathbf{W}_{2}$	$\displaystyle=(\kappa g_{1}A+g_{2}A^{3})\,\mathbf{q}+A^{\prime}\Phi_{c}\,\mathbf{e}_{1}\,,$		(4.53b)
	$\displaystyle\mathcal{L}_{6}\mathbf{W}_{6}+M_{c}\mathbf{W}_{6}$	$\displaystyle=g_{3}A^{3}\mathbf{q}\,.$		(4.53c)

We now impose a solvability condition for the solution to (4.53b). Recall from (3.24) that there is a non-trivial solution $\mathbf{\Phi}_{c}^{*}=(\Phi_{c}^{*},N_{c}^{*})$ to $\mathcal{L}_{2}\mathbf{\Phi}_{c}^{*}+M_{c}^{T}\mathbf{\Phi}_{c}^{*}=\mathbf{0}$.

As in the derivation of the eigenvalue expansion in (3.32), we have

$$\int_{0}^{\infty}\mathbf{\Phi}_{c}^{*}\cdot(\mathcal{L}_{2}\mathbf{W}_{2}+M_{c}\mathbf{W}_{2})\,\rho\,d\rho=0\,.$$

(4.54)

This yields that

$$\int_{0}^{\infty}\left(\mathbf{\Phi}_{c}^{*}\cdot\mathbf{q}\right)(\kappa g_{1}A+g_{2}A^{3})\,\rho\,\mathrm{d}\rho=-A^{\prime}\int_{0}^{\infty}\mathbf{\Phi}_{c}^{*}\cdot(\Phi_{c}\,\mathbf{e}_{1})\,\rho\,d\rho\,,$$

(4.55)