1. Introduction

Patterns, which represent a kind of organized yet heterogeneous macroscopic structure in time or space or space-time, are widely investigated using the reaction diffusion equations [24, 8, 3]. In the natural world, many animals have fascinating color patterns on their skins, exemplified by the coloration of zebrafish and tigers [20]. Such patterns are one of the most obvious traits of animals and serve a variety of functions. For example, patterns have been successfully used in camouflage, warning, social aggregation, mate choice, adaptive radiation, and other strategies [23, 14, 18]. Given patterns' prominence and ecological functions, zebrafish patterns often are determined by natural selection and are of particular interest to biologists [16]. These patterns elicited long-standing interest from developmental and cell biologists as well: their accessibility to observation and manipulation has made them a classic and enduring system for studying basic genetic and cellular mechanisms [16, 22, 19].

In recent years, some research results have been brought to bear on zebrafish patterns. In [4, 5], the authors established three different models to investigate the molecular mechanisms of zebrafish patterns, and discussed one-dimensional patterns. Schnackenberg [17] studied the following reaction-diffusion zebrafish model,

$\left\{\begin{array}{lcr}\frac{\partial u}{\partial t}=\nabla^2u+\gamma(a-u+u^2v)\\ \frac{\partial v}{\partial t}=d\nabla^2v+\gamma(b-u^2v), \end{array}\right.$

(1)

where $u$ and $v$ denote the concentration of activator and substrate, respectively; $a$ and $b$ are the basic production rate of two substances; $\gamma$ is the reaction rate. In this work, Schnackenberg studied a two-component chemical reaction system that had to involve at least three reactions for exhibiting limit cycle behavior. Under this condition, possible candidates for chemical limit systems were selected by postulating their steady state to be an instable focus. [12] also gave an account of various biological pattern formation phenomena, including the zebrafish patterns. The patterns are formed in the early developmental stages. In particular, the number of stripes did not change in an animal's life time, even when the body size increased considerably. However, the phenomenon has not always existed in nature. The authors [10] proved that when fish are growing, their skin patterns change, and this change in their patterns could be explained by a simple reaction-diffusion system. Based on the molecular mechanisms, it was suggested that leopard gene production is a component of putative reaction-diffusion system [1, 9], which also showed that mutations in the zebrafish gene, leopard, changed the pattern from stripes to spots. All the pattern variations of leopard mutations could be generated in a simulation by changing a parameter value that corresponds to the reaction kinetics in a putative reaction-diffusion system. [11, 26] discussed the discovery that zebrafish patterns appeared independent of the internal tissues or the body structures, and these patterns were robust against perturbation. In [21], Hiroto et al. studied Turing instability, and they controlled the variation of one parameter to discuss the directionality of the tripes of fish in terms of the anisotropy of diffusion. Two zebra fish of the same species have two different spots on the skin, their hybrid offspring may be give rise to spot-stripe mixtures on the skin [13], then the authors investigated the intermediate patterns and obtained the numerical result, which explained the hybridization of zebrafish in nature very well.

As is well known, the above existing zebrafish spatio-temporal models in ecological networks only concentrate on self-diffusion. However, zebrafish gene products can freely diffuse within a cell, apart from the random motion of individuals, i.e. self-diffusion, one production of genes tends to diffuse in the direction of the other. More precisely, the movement of a gene at any particular location is influenced by the gradient of the concentration of the other at that location. In biomathematics, such a scenario can be well described by reaction-diffusion systems with cross-diffusion terms [27], but it has not been applied to study zebrafish patterns. Therefore, considering the effect of cross-diffusions between zebrafishes will provide a new insight into pattern formation.

In light of the above discussions, based on system (1), the zebrafish pattern model can be formulated by the following equation with cross-diffusion terms:

$\left\{\begin{array}{lcr} \frac{\partial u(t,x,y)}{\partial t}=d_{11}\nabla^2 u(t,x,y)+d_{12}\nabla^2 v(t,x,y)+R(a-u^2(t,x,y)+u^2(t,x,y)v(t,x,y))\\ \frac{\partial v(t,x,y)}{\partial t}=d_{21}\nabla^2 u(t,x,y)+d_{22}\nabla^2 v(t,x,y)+R(b-u^2(t,x,y)v(t,x,y)) \end{array}\right.$

(2)

for $t>0$, $(x,y) \in \Omega=[0,L]\times[0,L]$ with homogeneous Neumann boundary conditions

$\frac{\partial u}{\partial \nu}(t,x,y)=\frac{\partial v}{\partial \nu}(t,x,y)=0,\;\;t\geq0, (x,y) \in \partial \Omega,$

and initial conditions

$\left\{ \begin{array}{l} u (0,x,y) = u_{0} \ge 0\\ v (0,x,y) = v_{0} \ge 0\\ \end{array} \right.$

(3)

where $u$ and $v$ denote the concentration of activator and substrate, respectively; $a$ and $b$ are the basic production rate of two substances; $R$ is the reaction rate. $\Omega$ is a bounded domain in $R^n$ with smooth boundary $\partial \Omega$, $\nu$ is the outward normal vector. $\nabla$ denotes the Laplacian operator in $R^{n}$, namely, $\nabla^2=\frac{\partial^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}$. $d_{11},d_{22}>0$ are the self-diffusion coefficients, $d_{12},d_{21}$ are the cross-diffusion coefficients, which may take positive or zero values. A positive value shows that one species of molecules tends to diffuse in the direction of the other concentration species. In this paper, our main contributions are summarized as follows.

(1) Diffusion is a very important factor to selection of patterns, especially, cross-diffusions. In this work, to better describe the zebrafish patterns, we propose a new viewpoint to research zebrafish patterns, i.e. studying the spatial dynamics of model (2) with cross-diffusions.

(2) Through bifurcation theory analysis of the equilibrium points for the proposed model, we obtain the sufficient conditions of Hopf bifurcation and Turing instability, respectively. To determine the selection of Turing patterns, we deduce the amplitude equations based on multiple-scale analysis.

(3) It is easy to gain five categories of Turing patterns by analyzing amplitude equations. The associated characteristic equation is established based on random initial perturbation. Then through stability analysis, we obtain the corresponding stable range of these patterns.

(4) Numerical results show that the five categories of patterns are decided by the parameter of model (2) or the coefficients of cross-diffusions. It is important to note that cross-diffusions also can change the selection of patterns, which is one reason the coloration of zebrafish varies.

---

2. The model and the analysis

In this section, based on Hopf bifurcation and Turing bifurcation theory, we will discuss the dynamical behavior of model (2). Before our discussion, according to thermodynamics theroy one should consider the following fact

$d_{11}d_{22}-d_{12}d_{21}\geq 0,$

which implies that all eigenvalues of the diffusion matrix $D=\left( {\begin{array}{*{20}{c}} d_{11} & d_{12} \\ d_{21} & d_{22} \\ \end{array}} \right)$ are positive definite.

For simplicity, one can rewrite model (2) as the following form

$\left\{\begin{array}{lcr} \frac{\partial u}{\partial t}=d_{11}\nabla^2 u+d_{12}\nabla^2v+f(u,v)\\ \frac{\partial v}{\partial t}=d_{21}\nabla^2 u+d_{22}\nabla^2 v+g(u,v),\\ \end{array}\right.$

(4)

where $f(u,v)=R(a-u+u^2v)$, $g(u,v)=R(b-u^2v)$.

The steady state of this system is $E_{0}$, $E_{0}=(u^{\ast},v^{\ast})=(a+b,\frac{b}{(a+b)^2})$, $a+b>0$ and $b>0$ because of practical significance. Let us briefly recall here the results of the linear stability analysis around $E_{0}$. The Jacobian matrix corresponding to $E_{0}$ is

$J= \left( {\begin{array}{*{20}{c}} a_{11} & a_{12} \\ a_{21} & a_{22} \\ \end{array}} \right),$

(5)

where $a_{11}=\frac{\partial f}{\partial u}|_{(u^{\ast},v^{\ast})}=-R\frac{a-b}{a+b}$, $a_{12}=\frac{\partial f}{\partial v}|_{(u^{\ast},v^{\ast})}=R(a+b)^2$, $a_{21}=\frac{\partial g}{\partial u}|_{(u^{\ast},v^{\ast})}=-R\frac{2b}{a+b}$, $a_{22}=\frac{\partial g}{\partial v}|_{(u^{\ast},v^{\ast})}=-R(a+b)^2$.

Then we get the following linear equation

$\left\{\begin{array}{lcr} \frac{\partial u}{\partial t}=d_{11}\nabla^2 u+d_{12}\nabla^2v+a_{11}u+a_{12}v\\ \frac{\partial v}{\partial t}=d_{21}\nabla^2 u+d_{22}\nabla^2 v+a_{21}u+a_{22}v.\\ \end{array}\right.$

(6)

Expanding the variables in the Fourier space

$\left(\begin{array}{ccc} u\\ v\\ \end{array}\right)=\sum\limits_{k=0}^{\infty} \left(\begin{array}{ccc} c_{k}^{1}\\ c_{k}^{2}\\ \end{array}\right)\exp(\lambda_{k}t+ikr),$

(7)

where $\lambda_{k}$ is the growth rate of perturbations in time $t$, $i$ is the imaginary unit and $i^2=-1$. $r=(x,y)$ is the spatial vector in two dimensions. We substitute (7) into (6), then obtaining the characteristic equation:

$\lambda^2-tr_{k}\lambda+\Delta_{k}=0.$

(8)

According to (8), it is easy to show that the eigenvalues $\lambda_{k}$ as follows

$\lambda_{k}=\frac{tr_{k}\pm\sqrt{tr_{k}^2-4\Delta_{k}}}{2},$

(9)

where

$tr_{k}=a_{11}+a_{22}-(d_{11}+d_{22})k^2,$

(10)

$\Delta_{k}=(d_{11}d_{22}-d_{12}d_{21})k^4+(a_{12}d_{21}+a_{21}d_{12}-a_{11}d_{22}-a_{22}d_{11})k^2+a_{11}a_{22}-a_{12}a_{21}.$

(11)

Hopf bifurcation occurs when $Im(\lambda_{k})\neq 0$, $Re(\lambda_{k})=0$, at $k=0$, i.e. $a_{11}+a_{22}=0$. So the critical value of the Hopf bifurcation parameter $b_{H}$ satisfies the following equation

$b^3+3ab^2+(3a^2-1)b+a^3+a=0.$

(12)

The unbalance changes of phases, corresponding to Turing branches, are the transitions of model from the uniform state to the oscillatory state. After the process, the formed patterns are called Turing patterns. From above discussion, we can obtain the necessary conditions for causing Turing instability. For some $k$, we have

$\left\{\begin{array}{lcr} tr_{0}=-\frac{a-b}{a+b}-(a+b)^2<0\\ \Delta_{0}=(a+b)^2>0\\ \Delta_{k}=(d_{11}d_{22}-d_{12}d_{21})k^4+R((a+b)^2d_{21}-\frac{2b}{a+b}d_{12}+\frac{a-b}{a+b}d_{22}\\ +(a+b)^2d_{11})k^2+R^2(a+b)^2<0. \end{array}\right.$

(13)

(13) indicates that system (2) is unstable for some perturbations to the wave number. So getting $\Delta_{k}=0$ at the critical value, namely, Turing bifurcation occurs when $Im(\lambda_{k})=0$, $Re(\lambda_{k})=0$, at $k=k_{T}\neq 0$. Therefore, the critical value of the Turing bifurcation parameter $b_{T}$ satisfies the following equation

$(a+b)^3(d_{11}+d_{21})-2bd_{12}+(a-b)d_{22}+2\sqrt{d}(a+b)^2=0,$

(14)

where $d=d_{11}d_{22}-d_{12}d_{21}$.

When Turing patterns come into being, the wave number $k_{T}$ satisfies

$k_{T}^2=-\frac{(a+b)^2d_{21}-\frac{2b}{a+b}d_{12}+\frac{a-b}{a+b}d_{22}+(a+b)^2d_{11}}{2d}R.$

(15)

From the above analysis, it can be obtained that the simplified conditions of Turing instability for model (2)

$\left\{\begin{array}{lcr} a-b+(a+b)^3>0\\ (a+b)^3(d_{11}+d_{21})-2bd_{12}+(a-b)d_{22}+2\sqrt{d}(a+b)^2<0.\\ \end{array}\right.$

(16)

Remark 1. Eq. (16) shows that diffusion can damage stability. Moreover, if the basic production rate of the activator is greater than the basic production rate of the substrate, namely, $a\geq b$, then the cross-diffusion coefficient $d_{12}$ is the key factor to damage stability and form Turing bifurcation. While, if the basic production rate of the activator is less than the basic production rate of the substrate, namely, $a< b$, then the self-diffusion coefficient $d_{22}$ together with the cross-diffusion coefficient $d_{12}$ are both the key factors to damage stability and form Turing bifurcation.

At the Turing bifurcation threshold, the spatial symmetry of system (2) is broken, and the patterns are stationary in time and oscillatory in space with the corresponding wavelength $\lambda_{T}=\frac{2\pi}{k_{T}}$.

According to Hopf and Turing bifurcation conditions, we could get the Hopf bifurcation region and Turing instability region. They are shown in the Figure 1.

Now let us observe the real parts of eigenvalues when $b$ is decreasing. Figure 2 shows the change of the real parts of eigenvalues.

---

3. Turing pattern of zebrafish

---

3.1. Pattern selection

In this section, we perform extensive numerical simulations of the spatially extended model (2) in two-dimensional spaces, and the qualitative results are shown here. All our numerical simulations employ the non-zero initial and zero-flux boundary conditions with a system size of $100\times100$, The space step is $\triangle x=1$, $\triangle y=1$, and the time step is $\triangle t=0.01$.

Take $a=0.14$, $R=1$, $d_{11}=1$, $d_{12}=2$, $d_{21}=2$, $d_{22}=20$. By a simple calculation, we find $g_{1}=2.4494$, $g_{2}=14.6185$, $h=1.7567$, $\mu_{1}=-0.0243$, $\mu_{2}=0$, $\mu_{3}=0.051$, $\mu_{4}=0.4067$ in Appendix A and B. Obviously the parameter values of $g_{1}$, $g_{2}$, and $h$ have the following relations: $g_{2}>g_{1}>|h|>0$, otherwise it is necessary to include some other terms up to the fourth order or higher.

Initially, the entire system is placed in the stationary state $(u^{\ast},v^{\ast})$. We run the numerical simulations until they reach a stationary state or until they show a behavior that does not seem to change its characteristics any longer. In the numerical simulations, different types of dynamics are observed. It is also found that the distribution of $u$ and $v$ are always the same type. Consequently, our analysis of pattern formation can be restricted to one distribution. In this section, the distribution of $u$ is showed. In the following, it will show the Turing patterns for the parameters $(a,b)$ located in the Turing space. It is reasonable to take random small initial perturbation of the stationary state $(u^{\ast},v^{\ast})$.

Example 1. Let $b=1.4$ and other parameters unchanged. By calculating, we obtain $\mu=0.0267$ satisfying $\mu_{2}<\mu<\mu_{3}$. According to Conclusion (Ⅲ), there is only $H_{\pi}$ hexagon patterns. These facts are illustrated by the numerical simulations in Figure 3, which shows that the $H_{\pi}$ hexagon patterns prevail over the whole domain, and the dynamics of system (2) does not undergo any further changes. The numerical simulation is compatible with the theoretical analysis. This pattern can be seen in nature, such as Figure 4. From Figure 4, we find that if the basic production rate of the substrate is 1.4 in a zebrafish, then the spot pattern of the zebrafish is unchanged and it is always a hexagon pattern.

Example 2. Take $b=1.3$, and other parameters are fixed. It is easy to gain $\mu=0.0962$ satisfying $\mu_{3}<\mu<\mu_{4}$. According to Conclusion (Ⅱ), (Ⅲ), $H_{\pi}$ hexagon patterns and stripe patterns appear. Figure 5 explains the fact that the stationary stripe patterns and $H_{\pi}$ emerge at the same time. This phenomenon is called pinning effect [7]. Combine with Figure 3, it can be obtained that when $\mu>\mu_{3}$, $H_{\pi}$ hexagon patterns begin to break up into stripe patterns gradually. Moreover, this pattern can be searched in the skin of zebrafish (Figure 6). From Figure 6, we find that if the basic production rate of the substrate is 1.3 in a zebrafish, then the spot patterns of the zebrafish will be the hexagon patterns and the stripe patterns.

Example 3. Choose $b=1.2$, and keep other parameters as above. By the calculation, the parameter satisfies $\mu_{3}<\mu=0.1657<\mu_{4}$. According to Conclusion (Ⅱ), stripe patterns emerge. The fact is illustrated by the numerical simulations in Figure 7. The numerical simulation is compatible with the theoretical analysis. As above analysis, the $H_{\pi}$ hexagon patterns vanish. In other words, there are only the stripe patterns. In nature, we can search this pattern on the body of zebra fish (Figure 8). From Figure 8, we find that if the basic production rate of the substrate is 1.2 in a zebrafish, then the spot pattern of the zebrafish is always a stripe pattern.

Example 4. Let $b=0.96$, and other parameters as Example 1. By a simple calculation, we obtain $\mu=0.3326$, which satisfies $\mu_{3}<\mu<\mu_{4}$. According to Conclusion (Ⅱ), (Ⅲ), $H_{0}$ hexagon patterns and stripe patterns emerge. The fact is illustrated by the numerical simulations in Figure 9. The numerical simulation corresponds to the theoretical analysis. System (2) is bistable, that is to say, the two kinds of patterns exist at the same time. When $\mu$ is close to the critical value $\mu_{4}$, the stripe patterns begin to break up and system (2) transfers to $H_{0}$ hexagon patterns from the stripe patterns gradually. Actually, the beautiful pattern exists in nature, such as Figure 10. From Figure 10, we find that if the basic production rate of the substrate is 0.96 in a zebrafish, then the spot pattern of the zebrafish becomes the hexagon pattern and the stripe pattern.

Example 5. Take $b=0.8$, and fix other parameters. By calculating, the parameter satisfies $\mu_{4}<\mu=0.4438$. Figure 11 shows that there are only $H_{0}$ hexagon patterns. The numerical simulation can not correspond to the theoretical analysis. This phenomenon can not be explained by the amplitude equations. This is the reentry of the hexagon patterns and can be explained as follows: when the system gets away from the Turing critical bifurcation line, some of the primary slave modes turn into active modes. We can not adiabatically eliminate them when deducing the amplitude equations. On the contrary, we should add them into the amplitude equations. In addition, this pattern can be seen in nature, such as Figure 12. From Figure 12, we find that if the basic production rate of the substrate is 0.8 in a zebrafish, then the spot pattern of the zebrafish is only the hexagon pattern.

---

3.2. Impact of cross-diffusions on pattern selection

Example 6. To observe the impact of different cross-diffusions on the selection of patterns, we consider system (2) with $a=0.14$, $b=1.4$, $R=1$, $d_{11}=1$, $d_{22}=20$. In order to compare with system (2) having cross-diffusion, it is necessary to take system (2) without cross-diffusions into account, namely, $d_{12}=d_{21}=0$. It is easy to obtain $g_{1}=20.5941$, $g_{2}=42.9383$, $h=0.3092$, $\mu_{1}=-0.00022447$, $\mu_{3}=0.0039$, $\mu_{4}=0.0161$, $\mu=0.0157$. Obviously, $\mu_{3}<\mu<\mu_{4}$. From Figure 13, we can observe that the stationary stripe patterns and $H_{0}$ emerges.

Example 7. Take $d_{12}=d_{21}=1$, and other parameters show as above. By a simple calculation, it is easy to show that $g_{1}=9.0002$, $g_{2}=24.764$, $h=1.2476$, $\mu_{1}=-0.0066$, $\mu_{3}=0.0564$, $\mu_{4}=0.2679$, $\mu=0.000785$, $\mu_{2}<\mu<\mu_{3}$. According to Conclusion (Ⅲ), there are only $H_{\pi}$ hexagon patterns in Figure 14, which is different from Figure 13. The numerical simulation corresponds to the theoretical analysis. The result illustrates that cross-diffusions can change the selection of patterns. This pattern can be seen in nature, such as Figure 15. From Figure 15, we find that the cross-influence on the activator and the substrate in a zebrafish will result in the single hexagon pattern.

Example 8. Choose $d_{12}=2$, $d_{21}=1$, and other parameters do not change. By a simple calculation, $g_{1}=5.7966$, $g_{2}=21.3726$, $h=1.6727$, $\mu_{1}=-0.0144$, $\mu_{3}=0.0669$, $\mu_{4}=0.3802$, $\mu=0.11$, $\mu_{3}<\mu<\mu_{4}$. According to Conclusion (Ⅱ), there are only stripe patterns. The patterns in Figure 15 is obviously different with Figure 13. The numerical simulation is compatible with the theoretical analysis. The result also implies the effect of cross-diffusions on the selection of patterns. This pattern can be seen in nature, such as Figure 17. From Figure 17, we find that the cross-influence on the activator and the substrate in a zebrafish can also result in the single stripe pattern.

---

4. Conclusions and discussions

On the one hand. The amplitude equations for the active modes are established, which determines the stability of amplitudes towards uniform and inhomogeneous perturbations. It presents all five categories of Turing patterns by analyzing amplitude equations, which indicates that the model dynamics exhibits complex pattern. In addition, system (2) perfectly simulates the pattern on the body of zebrafish. More specifically, from Figure 3 to Figure 12, we can find: in the range $\mu_{2}<\mu<\mu_{3}$, $H_{\pi}$ hexagon patterns emerge; in the range $\mu_{3}<\mu<\mu_{4}$, $H_{\pi}$-hexagon-stripe mixtures$\rightarrow$ stripes$\rightarrow$ $H_{0}$-hexagon-stripe mixtures can be observed; in the range $\mu_{4}<\mu$, $H_{0}$ hexagon patterns appear. At the same time, it is worth noting that $\mu$ is not close to the critical points $\mu_{2}$ and $\mu_{3}$, however, the numerical results cannot correspond perfectly to our theoretical analysis. This means that our theoretical analysis is appropriate just for the adjacent domains of the critical points $\mu_{2}$ and $\mu_{3}$. In a word, we obtain that the sequence "$H_{\pi}$ hexagons $\rightarrow$ $H_{\pi}$-hexagon-stripe mixtures $\rightarrow$ stripes $\rightarrow$ $H_{0}$-hexagon-stripe mixtures $\rightarrow$ $H_{0}$ hexagons" emerges by increasing the values of $b$. That implies parameter can lead to the change of pattern.

On the other hand, we also discuss the effect of cross-diffusions on pattern selection. From Figure 13, Figure 14, and Figure 16, it is obvious that $H_{0}$ hexagon patterns and stripe patterns coexist when system (2) have not cross-diffusions. However, there are only stripe patterns or hexagon patterns when the cross-diffusion terms emerge and change. In other words, cross-diffusions can change the selection of patterns. The methods and results in the present paper may enrich the research of the pattern formation in the zebrafish model, or may be useful for other reaction-diffusion systems. For example, the pattern research in reaction-diffusion ecological system is common and the theory is mature. But, with the best of our knowledge, a few researchers can combine the theoretical patterns with the realistic patterns. That is to say, they can not find the realistic patterns to support their theoretical analysis. This paper makes up these lacks and we have give some real zebrafish to support our results.

---

Acknowledgment

The work is partially supported by National Natural Science Foundation of China under Grant 11571170 and 61174155. The work is also sponsored by Funding of Jiangsu Innovation Program for Graduate Education \break KYZZ$15\_0091$, the Fundamental Research Funds for the Central Universities.

---

Appendix A: Amplitude equations

Based on the above discussion, we still cannot determine the selection and competition of Turing patterns. In the following, we will discuss the pattern selection close to the onset $b=b_{T}$ (Although we have not given the exact expression of $b_{T}$, we can also discuss the amplitude equations near $b=b_{T}$ so long as $b_{T}$ being as an entirety. Moreover, in numerical simulations we can calculate $b_{T}$ accurately.), the eigenvalues associated with the critical modes are close to zero and they are slowly varying modes, but the off-critical modes relax quickly, so only perturbations with $k$ around $k_{T}$ need to be considered. Consequently, all dynamics can be reduced to the dynamics of active slow modes [2], so the stability and selection of different patterns close to the onset can be studied by investigating the amplitude equations of the system. Driving the coefficients of amplitude equations has two methods, one is symmetrical analysis and the other is standard multiple-scale analysis. Symmetrical analysis is very concise, however, it has many limits as well, especially on describing a partial differential equation, which has specific relationship between coefficients and parameters of system. As is well known, the standard multiple-scale analysis is one of the best way to derive amplitude equations [15, 6]. The method is based on the theory near the instability threshold. The basic state is unstable only with respect to perturbations with wave numbers close to the critical value $k_{T}$.

Let $\bar{u}=u-u^{\ast}$, $\bar{v}=v-v^{\ast}$ and drop the bars for simplicity of notations. System (2) can be rewritten as the following form

$\left\{\begin{array}{lcr} \frac{\partial u}{\partial t}=d_{11}\nabla^2 u+d_{12}\nabla^2v+a_{11}u+a_{12}v+R(v^{\ast}u^2+2u^{\ast}uv+u^2v)\\ \frac{\partial v}{\partial t}=d_{21}\nabla^2 u+d_{22}\nabla^2 v+a_{21}u+a_{22}v-R(v^{\ast}u^2+2u^{\ast}uv+u^2v),\\ \end{array}\right.$

(17)

transforming the positive constant steady state $E_{0}$ into the zero equilibrium. In the following, we take Taylor series expansion and omit the terms $o(x^3)$. The solutions of system (17) can be expanded as

$\left(\begin{array}{ccc} u\\ v\\ \end{array}\right)=\sum\limits_{i=1}^{3} \left(\begin{array}{ccc} A_{i}^{u}\\ A_{i}^{v}\\ \end{array}\right)e^{ik_{i}r}+c.c,$

(18)

where $c.c$ denotes the complex conjugate.

Setting

$c= \left(\begin{array}{ccc} u\\ v\\ \end{array}\right) , N= \left(\begin{array}{ccc} N_{1}\\ N_{2}\\ \end{array}\right).$

Model (17) can be converted to the following system

$\frac{\partial c}{\partial t}=Lc+N,$

(19)

where

$L= \left( {\begin{array}{*{20}{c}} a_{11}+d_{11}\nabla^2 & a_{12}+d_{12}\nabla^2 \\ a_{21}+d_{21}\nabla^2 & a_{22}+d_{22}\nabla^2 \\ \end{array}} \right),$

(20)

$N= \left(\begin{array}{ccc}R(v^{\ast}u^2+2u^{\ast}uv+u^2v)\\ -R(v^{\ast}u^2+2u^{\ast}uv+u^2v)\\ \end{array}\right).$

(21)

We only analyse the behavior of the controlled parameter which is close to the onset $b=b_{T}$ for this system. Then expanding $b$ in the following term with this method

$b_{T}-b=\varepsilon^2b_{2}+O(\varepsilon^3),$

(22)

where $\varepsilon$ is a small parameter.

By expanding the variable $c$ and the nonlinear term $N$ according to the small parameter, we obtain the following results:

$c= \left(\begin{array}{ccc} u\\ v\\ \end{array}\right)= \varepsilon\left(\begin{array}{ccc} u_{1}\\ v_{1}\\ \end{array}\right)+ \varepsilon^2\left(\begin{array}{ccc} u_{2}\\ v_{2}\\ \end{array}\right)+\cdots,$

(23)

$N=\varepsilon^2h^2+\varepsilon^3h^3+O(h^4),$

(24)

where $h^2$, $h^3$ correspond to the second and third order of $\varepsilon$ in the expansion of the nonlinear operator respectively. $L$ can be expanded as follows

$L=L_{T}+(b_{T}-b)M,$

(25)

where

$L_{T}= \left( {\begin{array}{*{20}{c}} a_{11}^{\ast}+d_{11}\nabla^2 & a_{12}^{\ast}+d_{12}\nabla^2 \\ a_{21}^{\ast}+d_{21}\nabla^2 & a_{22}^{\ast}+d_{22}\nabla^2 \\ \end{array}} \right),$

(26)

$M= \left( {\begin{array}{*{20}{c}} b_{11}& b_{12}\\ b_{21}& b_{22}\\ \end{array}} \right),$

(27)

where $a_{11}^{\ast}=a_{11}|_{b=b_{T}}, a_{12}^{\ast}=a_{12}|_{b=b_{T}}, a_{21}^{\ast}=a_{21}|_{b=b_{T}}, a_{22}^{\ast}=a_{22}|_{b=b_{T}}$, $b_{11}=$$-R\frac{2a}{(a+b_{T})^2}, b_{12}=-2R(a+b_{T}), b_{21}=R\frac{2a}{(a+b_{T})^2}, b_{22}=2R(a+b_{T})$.

The core of the standard multiple-scale analysis is separating the dynamical behavior of system (19) according to different time scale or spatial scale. It is needed to separate the time scale for model (17). Each time scale $T_{i}$ can be considered as an independent variable. $T_{i}$ corresponds to the dynamical behaviors of variable whose scale are $\varepsilon^{-i}$. Therefore, the derivative with respect to time converts to the following term

$\frac{\partial}{\partial t}=\varepsilon^2\frac{\partial}{\partial T_{2}}+O(\varepsilon^3).$

(28)

We substitute (20) $\sim$ (27) into (19) and expand (19) according to different orders of $\varepsilon$, obtaining three equations

$\varepsilon$:

$L_{T} \left(\begin{array}{ccc} u_{1}\\ v_{1}\\ \end{array}\right)=0.$

(29)

$\varepsilon^2$:

$L_{T} \left(\begin{array}{ccc} u_{2}\\ v_{2}\\ \end{array}\right)= \left(\begin{array}{ccc} -R(v^{\ast}u_{1}^2+2u^{\ast}u_{1}v_{1})\\ R(v^{\ast}u_{1}^2+2u^{\ast}u_{1}v_{1})\\ \end{array}\right).$

(30)

$\varepsilon^3$:

$L_{T} \left(\begin{array}{ccc} u_{3}\\ v_{3}\\ \end{array}\right)=\frac{\partial}{\partial T_{2}} \left(\begin{array}{ccc}u_{1}\\ v_{1}\\ \end{array}\right)-b_{2}M \left(\begin{array}{ccc}u_{1}\\ v_{1} \end{array}\right)- \left(\begin{array}{ccc}R(2v^{\ast}u_{1}u_{2}+2u^{\ast}u_{1}v_{2}+2u^{\ast}u_{2}v_{1}+u_{1}^2v_{1})\\ -R(2v^{\ast}u_{1}u_{2}+2u^{\ast}u_{1}v_{2}+2u^{\ast}u_{2}v_{1}+u_{1}^2v_{1})\\ \end{array}\right).$

(31)

For the first order of $\varepsilon$, since $L_{T}$ is the linear operator of system (2) close to the onset, $(u_{1},v_{1})^{T}$ is the linear combination of the eigenvectors that corresponds to the eigenvalue zero. Therefore, solving the first order of $\varepsilon$, we will obtain

$\left(\begin{array}{ccc} u_{1}\\ v_{1}\\ \end{array}\right)=\left(\begin{array}{ccc} S_{1}\\ 1\\ \end{array}\right)(W_{1}e^{ik_{1}r}+W_{2}e^{ik_{2}r}+W_{3}e^{ik_{3}r})+c.c,$

(32)

where $S_{1}=-\frac{a_{12}^{\ast}-d_{12}k_{T}^2}{a_{11}^{\ast}-d_{11}k_{T}^2}$, $|k_{j}|=k_{T},(j=1,2,3)$. $W_{j}$ is the amplitude of the mode $e^{ik_{j}r}$ when the system is under the first-order perturbation. Its form is determined by the perturbation term of the high order.

For the second order of $\varepsilon$, we have

$L_{T} \left(\begin{array}{ccc} u_{2}\\ v_{2}\\ \end{array}\right)= \left(\begin{array}{ccc} -R(v^{\ast}u_{1}^2+2u^{\ast}u_{1}v_{1})\\ R(v^{\ast}u_{1}^2+2u^{\ast}u_{1}v_{1})\\ \end{array}\right)=\left(\begin{array}{ccc}F_{u}\\ F_{v} \end{array}\right).$

(33)

According to the Fredholm condition, to guarantee the existence of the nontrivial solution of the equation, the vector function of the right-hand side of (31) must be orthogonal with the zero eigenvectors of operator $L_{T}^{+}$. $L_{T}^{+}$ is the adjoint and transposed operator of $L_{T}$. In the system, the zero eigenvectors of operator $L_{T}^{+}$ are

$\left(\begin{array}{ccc}1\\ S_{2}\\ \end{array}\right)e^{-ik_{j}r}+c.c (j=1,2,3),$

(34)

where $S_{2}=-\frac{a_{11}^{\ast}-d_{11}k_{T}^2}{a_{21}^{\ast}-d_{21}k_{T}^2}$.

The orthogonality condition is

$(1,S_{2}) \left(\begin{array}{ccc}F_{u}^{i}\\ F_{v}^{i}\\ \end{array}\right)=0,$

(35)

where $F_{u}^{i}$, $F_{v}^{i}$, respectively, represent the coefficients corresponding to $e^{ik_{i}r}$ in $F_{u}$, $F_{v}$.

We might as well take $e^{ik_{1}r}$ as an example

$(1,S_{2}) \left(\begin{array}{ccc}-R(2v^{\ast}S_{1}^2\bar{W_{2}}\bar{W_{3}}+4u^{\ast}S_{1}\bar{W_{2}}\bar{W_{3}})\\ R(2v^{\ast}S_{1}^2\bar{W_{2}}\bar{W_{3}}+4u^{\ast}S_{1}\bar{W_{2}}\bar{W_{3}})\\ \end{array}\right)=0,$

(36)

Solving the above equation and obtaining $\bar{W}_{2}\bar{W}_{3}=0$, in a similar way, $\bar{W}_{3}\bar{W}_{1}=0$, $\bar{W}_{1}\bar{W}_{2}=0$.

We suppose

$\begin{array}{l} \left( {\begin{array}{*{20}{c}} {{u_2}}\\ {{v_2}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{X_0}}\\ {{Y_0}} \end{array}} \right) + \sum\limits_{i = 1}^3 {\left( {\begin{array}{*{20}{c}} {{X_i}}\\ {{Y_i}} \end{array}} \right)} {e^{i{k_i}r}} + \sum\limits_{i = 1}^3 {\left( {\begin{array}{*{20}{c}} {{X_{ii}}}\\ {{Y_{ii}}} \end{array}} \right)} {e^{2i{k_i}r}} + \\ \left( {\begin{array}{*{20}{c}} {{X_{12}}}\\ {{Y_{12}}} \end{array}} \right){e^{i({k_1} - {k_2})r}} + \left( {\begin{array}{*{20}{c}} {{X_{23}}}\\ {{Y_{23}}} \end{array}} \right){e^{i({k_2} - {k_3})r}} + \left( {\begin{array}{*{20}{c}} {{X_{31}}}\\ {{Y_{31}}} \end{array}} \right){e^{i({k_3} - {k_1})r}} + c.c, \end{array}$

(37)

the coefficients of (37) can be obtained by solving the set of linear equations about $e^0$, $e^{ik_{i}r}$, $e^{2ik_{i}r}$, $e^{2i(k_{i}-k_{j})r}$. We have

$\left(\begin{array}{ccc}X_{0}\\ Y_{0}\\ \end{array}\right)=\left(\begin{array}{ccc}C_{1}\\ C_{2}\\ \end{array}\right)(|W_{1}^2|+|W_{2}^2|+|W_{3}^2|),$

(38)

where

$C_{1}=-\frac{R(2v^{\ast}S_{1}^2+4u^{\ast}S_{1})(a_{12}^{\ast}+a_{22}^{\ast})}{a_{11}^{\ast}a_{22}^{\ast}-a_{12}^{\ast}a_{21}^{\ast}},\\ C_{2}=\frac{R(2v^{\ast}S_{1}^2+4u^{\ast}S_{1})(a_{11}^{\ast}+a_{21}^{\ast})}{a_{11}^{\ast}a_{22}^{\ast}-a_{12}^{\ast}a_{21}^{\ast}}.\\ X_{i}=S_{1}Y_{i}(i=1,2,3),$

(39)

$\left(\begin{array}{ccc}X_{ii}\\ Y_{ii}\\ \end{array}\right)=\left(\begin{array}{ccc}E_{1}\\ E_{2}\\ \end{array}\right)W_{i}^2,$

(40)

where

$E_{1}=-\frac{R(a_{12}^{\ast}+a_{22}^{\ast}-4d_{12}k_{T}^2-4d_{22}k_{T}^2)(v^{\ast}S_{1}^2+2u^{\ast}S_{1})}{(a_{11}^{\ast}-4d_{11}k_{T}^2)(a_{22}^{\ast}-4d_{22}k_{T}^2)-(a_{12}^{\ast}-4d_{12}k_{T}^2)(a_{21}^{\ast}-4d_{21}k_{T}^2)},\\ E_{2}=\frac{R(a_{11}^{\ast}+a_{21}^{\ast}-4d_{11}k_{T}^2-4d_{21}k_{T}^2)(v^{\ast}S_{1}^2+2u^{\ast}S_{1})}{(a_{11}^{\ast}-4d_{11}k_{T}^2)(a_{22}^{\ast}-4d_{22}k_{T}^2)-(a_{12}^{\ast}-4d_{12}k_{T}^2)(a_{21}^{\ast}-4d_{21}k_{T}^2)}.\\ \left(\begin{array}{ccc}X_{jk}\\ Y_{jk}\\ \end{array}\right)=\left(\begin{array}{ccc}F_{1}\\ F_{2}\\ \end{array}\right)W_{j}^2\bar{W}_{k}^2,$

(41)

where

$F_{1}=-\frac{R(a_{12}^{\ast}+a_{22}^{\ast}-3d_{12}k_{T}^2-3d_{22}k_{T}^2)(v^{\ast}S_{1}^2+2u^{\ast}S_{1})}{(a_{11}^{\ast}-3d_{11}k_{T}^2)(a_{22}^{\ast}-3d_{22}k_{T}^2)-(a_{12}^{\ast}-3d_{12}k_{T}^2)(a_{21}^{\ast}-3d_{21}k_{T}^2)},\\ F_{2}=\frac{R(a_{11}^{\ast}+a_{21}^{\ast}-3d_{11}k_{T}^2-3d_{21}k_{T}^2)(v^{\ast}S_{1}^2+2u^{\ast}S_{1})}{(a_{11}^{\ast}-3d_{11}k_{T}^2)(a_{22}^{\ast}-3d_{22}k_{T}^2)-(a_{12}^{\ast}-3d_{12}k_{T}^2)(a_{21}^{\ast}-3d_{21}k_{T}^2)}.\\ $

For the third order of $\varepsilon$, we have

$L_{T} \left(\begin{array}{ccc}u_{3}\\ v_{3}\\ \end{array}\right)=\frac{\partial}{\partial T_{2}}\left(\begin{array}{ccc}u_{1}\\ v_{1}\\ \end{array}\right)-b_{2}M\left(\begin{array}{ccc}u_{1}\\ v_{1}\\ \end{array}\right)-\left(\begin{array}{ccc}R(2v^{\ast}u_{1}u_{2}+2u^{\ast}u_{1}v_{2}+2u^{\ast}u_{2}v_{1}+u_{1}^2v_{1})\\ -R(2v^{\ast}u_{1}u_{2}+2u^{\ast}u_{1}v_{2}+2u^{\ast}u_{2}v_{1}+u_{1}^2v_{1})\\ \end{array}\right)\\ =\left(\begin{array}{ccc}F_{u}\\ F_{v}\\ \end{array}\right).$

(42)

According to the orthogonality condition

$(1,S_{2}) \left(\begin{array}{ccc}F_{u}^{i}\\ F_{v}^{i}\\ \end{array}\right)=0,$

(43)

we have

$\begin{array}{l} ({S_1} + {S_2})\frac{{\partial {W_1}}}{{\partial {T_2}}} = {b_2}({S_1}{b_{11}} + {S_1}{S_2}{b_{21}} + {b_{12}} + {S_2}{b_{22}}){W_1} + \\ 2R{S_1}(1 - {S_2})({v^*}{S_1} + 2{u^*})({{\bar Y}_2}{{\bar W}_3} + {{\bar Y}_3}{{\bar W}_2}) - ({G_1}|W_1^2| + {G_2}(|W_2^2| + |W_3^2|)){W_1}, \end{array}$

(44)

where

$G_{1}=R(S_{2}-1)(2v^{\ast}S_{1}(E_{1}+C_{1})+2u^{\ast}S_{1}(E_{2}+C_{2})+2u^{\ast}(E_{1}+C_{1})+3S_{1}^2),\\ G_{2}=R(S_{2}-1)(2v^{\ast}S_{1}(F_{1}+C_{1})+2u^{\ast}S_{1}(F_{2}+C_{2})+2u^{\ast}(F_{1}+C_{1})+6S_{1}^2).\\ $

We expand amplitude $A_{1}$ as the following form

$A_{1}=\varepsilon W_{1}+\varepsilon^2Y_{1}+\cdots.$

(45)

Substituting $A_{1}$ into (44), we can obtain the amplitude equation corresponding to $A_{1}$ as follows

$\tau_{0}\frac{\partial A_{1}}{\partial t}=\mu A_{1}+h\bar{A}_{2}\bar{A}_{3}-[g_{1}|A_{1}|^2+g_{2}(|A_{2}|^2+|A_{3}|^2)]A_{1},$

(46)

where

$\mu=\frac{b_{T}-b}{b_{T}}, \tau_{0}=\frac{(S_{1}+S_{2})}{b_{T}(S_{1}b_{11}+S_{1}S_{2}b_{21}+b_{12}+S_{2}b_{22})},\\ h=\frac{2RS_{1}(1-S_{2})(v^{\ast}S_{1}+2u^{\ast})}{b_{T}(S_{1}b_{11}+S_{1}S_{2}b_{21}+b_{12}+S_{2}b_{22})},\\ g_{1}=\frac{G_{1}}{b_{T}(S_{1}b_{11}+S_{1}S_{2}b_{21}+b_{12}+S_{2}b_{22})},g_{2}=\frac{G_{2}}{b_{T}(S_{1}b_{11}+S_{1}S_{2}b_{21}+b_{12}+S_{2}b_{22})}.$

Another two equations can be obtained through the transformation of the subscript of $A$, we have

$\left\{\begin{array}{lcr} \tau_{0}\frac{\partial A_{1}}{\partial t}=\mu A_{1}+h\bar{A}_{2}\bar{A}_{3}-[g_{1}|A_{1}|^2+g_{2}(|A_{2}|^2+|A_{3}|^2)]A_{1}\\ \tau_{0}\frac{\partial A_{2}}{\partial t}=\mu A_{2}+h\bar{A}_{3}\bar{A}_{1}-[g_{1}|A_{2}|^2+g_{2}(|A_{1}|^2+|A_{3}|^2)]A_{2}\\ \tau_{0}\frac{\partial A_{3}}{\partial t}=\mu A_{3}+h\bar{A}_{1}\bar{A}_{2}-[g_{1}|A_{3}|^2+g_{2}(|A_{1}|^2+|A_{2}|^2)]A_{3}.\\ \end{array}\right.$

(47)

---

Appendix B: Turing pattern analysis

(47) can be decomposed to mode $\rho_{i}=|A_{i}|$ and corresponding phase angle $\varphi_{i}$[15, 28]. Therefore, we will gain four differential equations of the real variables in the following form:

$\left\{\begin{array}{lcr} \tau_{0}\frac{\partial \varphi}{\partial t}=-\frac{h(\rho_{1}^2\rho_{2}^2+\rho_{1}^2\rho_{3}^2+\rho_{2}^2\rho_{3}^2)\sin \varphi}{\rho_{1}\rho_{2}\rho_{3}}\\ \tau_{0}\frac{\partial \rho_{1}}{\partial t}=\mu \rho_{1}+h\rho_{2}\rho_{3}\cos\varphi-g_{1}\rho_{1}^3-g_{2}(\rho_{2}^2+\rho_{3}^2)\rho_{1}\\ \tau_{0}\frac{\partial \rho_{2}}{\partial t}=\mu \rho_{2}+h\rho_{1}\rho_{3}\cos\varphi-g_{1}\rho_{2}^3-g_{2}(\rho_{1}^2+\rho_{3}^2)\rho_{2}\\ \tau_{0}\frac{\partial \rho_{3}}{\partial t}=\mu \rho_{3}+h\rho_{1}\rho_{2}\cos\varphi-g_{1}\rho_{3}^3-g_{2}(\rho_{1}^2+\rho_{2}^2)\rho_{3},\\ \end{array}\right.$

(48)

where $\varphi=\varphi_{1}+\varphi_{2}+\varphi_{3}$.

System (48) has four kinds of solutions.

(ⅰ) The stationary state $(0,0,0)$ is given by

$\rho_{1}=\rho_{2}=\rho_{3}=0.$

(49)

(ⅱ) Stripes (S) are given by

$\rho_{1}=\sqrt{\frac{\mu}{g_{1}}}, \rho_{2}=\rho_{3}=0.$

(50)

(ⅲ) Hexagons $(H_{0},H_{\pi})$ are given by

$\rho_{1}=\rho_{2}=\rho_{3}=\frac{|h|\pm\sqrt{h^2+4(g_{1}+2g_{2}\mu)}}{2(g_{1}+2g_{2})},$

(51)

and exist in the following condition

$\mu>\mu_{1}=\frac{-h^2}{4(g_{1}+2g_{2})},$

(52)

where $\varphi=0$ corresponds to $H_{0}$ Hexagon, $\varphi=\pi$ corresponds to $H_{\pi}$ Hexagon.

(ⅳ) The mixed states are given by

$\rho_{1}=\frac{|h|}{g_{2}-g_{1}},\rho_{2}=\rho_{3}=\sqrt{\frac{h-g_{1}\rho_{1}^2}{g_{1}+g_{2}}},$

(53)

with $g_{2}>g_{1}$.

In the following, we will give a discussion about the stability of the above four stationary solutions.

To stripes, we give a perturbation at stationary solution $(\rho_{0},0,0)$ for studying the stability of stationary solution (50), where $\rho_{0}=\sqrt{\frac{\mu}{g_{1}}}$. Let $\rho_{1}=\rho_{0}+\Delta\rho_{1}$, $\rho_{2}=\Delta\rho_{2}$, $\rho_{3}=\Delta\rho_{3}$, linearization of Eq.(48) can be written as

$\frac{\partial \rho}{\partial t}=L_{1}\cdot \rho,$

(54)

where

$L_{1}= \left(\begin{array}{ccc} \mu-3g_{1}\rho_{0}^2 & 0 & 0\\ 0 &\mu-g_{2}\rho_{0}^2 & |h|\rho_{0}\\ 0 & |h|\rho_{0} & \mu-g_{2}\rho_{0}^2\\ \end{array}\right),\rho=\left(\begin{array}{ccc} \Delta\rho_{1}\\ \Delta\rho_{2}\\ \Delta\rho_{3}\\ \end{array}\right).$

(55)

The characteristic equation of matrix $L_{1}$ can be obtained as

$\lambda^{3}+P_{1}\lambda^{2}+P_{2}\lambda+P_{3}=0,$

(56)

where

$P_{1}=(3g_{1}+2g_{2})\rho_{0}^2-3\mu$,

$P_{2}=(6g_{1}g_{2}+g_{2}^2)\rho_{0}^4-(4\mu g_{2}+h^2+6\mu g_{1})\rho_{0}^2+3\mu^2,$

$P_{3}=3g_{1}g_{2}^2\rho_{0}^6-(3g_{1}h^2+\mu g_{2}^2+6\mu g_{1}g_{2})\rho_{0}^4+2\mu^2g_{2}+3\mu^2g_{1}+\mu h^2\rho_{0}^2-\mu^3.$

The eigenvalues are

$\lambda_{1}=\mu-3g_{1}\rho_{0}^2, \lambda_{2}=\mu+h\rho_{0}-g_{2}\rho_{0}^2,\lambda_{3}=\mu-h\rho_{0}-g_{2}\rho_{0}^2.$

(57)

Substituting $\rho_{0}=\sqrt{\frac{\mu}{g_{1}}}$ into Eq.(57), we obtain

$\lambda_{1}=-2\mu, \lambda_{2}=\mu(1-\frac{g_{2}}{g_{1}})+|h|\sqrt{\frac{\mu}{g_{1}}}, \lambda_{3}=\mu(1-\frac{g_{2}}{g_{1}})-|h|\sqrt{\frac{\mu}{g_{1}}}.$

(58)

It is all known to us that Eq.(58) has stable solutions when the eigenvalues $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are all negative. since $\mu>0,\frac{g_{2}}{g_{1}}>1$, the three eigenvalues are negative if the following holds

$\mu>\mu_{3}=\frac{h^2g_{1}}{(g_{2}-g_{1})^2}.$

(59)

In the following, we discuss the stability of the hexagons, similar to the above process, we perturb Eq.(51) at the point $(\rho_{0},\rho_{0},\rho_{0})$ as follows

$\rho_{i}=\rho_{0}+\Delta\rho_{i},$

(60)

where $\rho_{0}=\frac{|h|\pm\sqrt{h^2+4(g_{1}+2g_{2})\mu}}{2(g_{1}+2g_{2})}$. Therefore, Eq.(48) can be linearized as

$\frac{\partial \rho}{\partial t}=L_{2}\cdot\rho,$

(61)

where

$L_{2}= \left(\begin{array}{ccc} \mu-(3g_{1}+2g_{2}) & |h|\rho_{0}-2g_{2}\rho_{0}^2 & |h|\rho_{0}-2g_{2}\rho_{0}^2\\ |h|\rho_{0}-2g_{2}\rho_{0}^2 &\mu-(3g_{1}+2g_{2}) & |h|\rho_{0}-2g_{2}\rho_{0}^2\\ |h|\rho_{0}-2g_{2}\rho_{0}^2 & |h|\rho_{0}-2g_{2}\rho_{0}^2 & \mu-(3g_{1}+2g_{2})\\ \end{array}\right),\rho=\left(\begin{array}{ccc} \Delta\rho_{1}\\ \Delta\rho_{2}\\ \Delta\rho_{3}\\ \end{array}\right).$

(62)

The characteristic equation of $L_{2}$ can be obtained as

$\lambda^{3}+Q_{1}\lambda^{2}+Q_{2}\lambda+Q_{3}=0,$

(63)

where

$Q_{1}=(9g_{1}+6g_{2})\rho_{0}^2-3\mu$,

$Q_{2}=(27g_{1}^2+36g_{1}g_{2})\rho_{0}^4+12|h|g_{2}\rho_{0}^3-(18\mu g_{1}+3h^2+12\mu g_{2})\rho_{0}^2+3\mu^2$,

$\begin{array}{l} {Q_3} = (54g_1^2{g_2} + 27g_1^3)\rho _0^6 + 36|h|{g_1}{g_2}\rho _0^5 + (6{h^2}{g_2} - 36\mu {g_1}{g_2} - 9{h^2}{g_1} - 27\mu g_1^2)\rho _0^4 - \\ (2|h{|^3} + 12\mu |h|{g_2})\rho _0^3 + (9{\mu ^2}{g_1} + 6{\mu ^2}{g_2} + 3\mu {h^2})\rho _0^2 - {\mu ^3}. \end{array}$

Solving the characteristic equation (63) and we have

$\lambda_{1}=\lambda_{2}=\mu-|h|\rho_{0}-3g_{1}\rho_{0}^2,\lambda_{3}=\mu-3g_{1}\rho_{0}^2-63g_{2}\rho_{0}^2+2|h|\rho_{0}.$

(64)

Substituting $\rho_{0}=\sqrt{\frac{\mu}{g_{1}}}$ into Eq.(64), we can obtain two cases of stability as follows.

For the stationary solution

$\rho^{-}=\frac{|h|-\sqrt{h^2+4(g_{1}+2g_{2}\mu)}}{2(g_{1}+2g_{2})},$

$\lambda_{1}$ and $\lambda_{2}$ are always positive, so the corresponding pattern is also always unstable.

For the stationary solution

$\rho^{+}=\frac{|h|+\sqrt{h^2+4(g_{1}+2g_{2}\mu)}}{2(g_{1}+2g_{2})}$

$\lambda_{i}(i=1,2,3)$ is negative when the parameter $\mu$ satisfies the following condition

$\mu<\mu_{4}=\frac{2g_{1}+g_{2}}{(g_{2}-g_{1})^2}h^2.$

(65)

Based on the above analysis, we can conclude

(Ⅰ) The stationary state (0, 0, 0) is stable for $\mu<\mu_{2}=0$ and unstable for $\mu>\mu_{2}$.

(Ⅱ) The stripe is stable when $\mu>\mu_{3}$.

(Ⅲ) The hexagon $H_{\pi}$ is stable only for $\mu<\mu_{4}$ and hexagon $H_{0}$ is unstable.

(Ⅳ) the mixed states can exist for $\mu>\mu_{3}$ and are always unstable.

---