Pattern formation in reaction-diffusion systems in the presence of short-term memory

TL;DR

This paper investigates how short-term memory influences pattern formation in reaction-diffusion systems, extending traditional models beyond Markovian assumptions to better understand and potentially control spatio-temporal patterns.

Contribution

It introduces a non-Markovian reaction-diffusion model, derives a complex Ginzburg-Landau equation incorporating memory effects, and demonstrates how memory alters pattern properties and stability.

Findings

Memory significantly affects pattern properties.

Derived a memory-dependent complex Ginzburg-Landau equation.

Numerical simulations confirm the impact of memory on patterns.

Abstract

We study reaction-diffusion systems beyond the Markovian approximation to take into account the effect of memory on the formation of spatio-temporal patterns. Using a non-Markovian Brusselator model as a paradigmatic example, we show how to use reductive perturbation to investigate the formation and stability of patterns. Focusing in detail on the Hopf instability and short-term memory, we derive the corresponding complex Ginzburg-Landau equation that governs the amplitude of the critical mode and we establish the explicit dependence of its parameters on the memory properties. Numerical solution of this memory dependent complex Ginzburg-Landau equation as well as direct numerical simulation of the non-Markovian Brusselator model illustrate that memory changes the properties of the spatio-temporal patterns. Our results indicate that going beyond the Markovian approximation might be necessary to study the formation of spatio-temporal patterns even in systems with short-term memory. At the same time, our work opens up a new window into the control of these patterns using memory.

Equations

$$\frac{\partial P_{\xi }(t)}{\partial t}={\int }_0^{t}d{t}^{\prime }\sum _{\mu }[W_{\xi \mu }(t-t^{\prime })P_{\mu }(t^{\prime })-W_{\mu \xi }(t-t^{\prime })P_{\xi }(t^{\prime })],$$

$$W_{\mu \xi }(t^{\prime })=F_{\mu \xi }\delta (t^{\prime }),$$

$$\frac{\partial P_{\xi }(t)}{\partial t}=\sum _{\mu }[F_{\xi \mu }{P}_{\mu }(t)-F_{\mu \xi }{P}_{\xi }(t)].$$

$$\frac{\partial P(x,t)}{\partial t}=D\frac{{\partial }^{2}P(x,t)}{\partial x^2},$$

$$W_{\mu \xi }(t)=F_{\mu \xi }\mathcal{K}(t).$$

$$\frac{\partial P_{\xi }(t)}{\partial t}={\int }_0^{t}d{t}^{\prime }\mathcal{K}(t-t^{\prime })\sum _{\mu }[F_{\xi \mu }{P}_{\mu }(t^{\prime })-F_{\mu \xi }{P}_{\xi }(t^{\prime })].$$

$$\frac{\partial P(x,t)}{\partial t}={\int }_0^{t}d{t}^{\prime }\mathcal{K}(t-t^{\prime })D\frac{{\partial }^{2}P(x,t^{\prime })}{\partial x^2},$$

$$\left\{\begin{array}[]{cc}\frac{\partial U({\bf r},t)}{\partial t}=f_{1}(U,V;\mu)+\int_{0}^{t}dt^{\prime}{\cal K}_{U}(t-t^{\prime})D_{U}\nabla^{2}U({\bf r},t^{\prime})\\ \frac{\partial V({\bf r},t)}{\partial t}=f_{2}(U,V;\mu)+\int_{0}^{t}dt^{\prime}{\cal K}_{V}(t-t^{\prime})D_{V}\nabla^{2}V({\bf r},t^{\prime})\end{array}\right.,$$

$$\left\{\begin{array}[]{cc}\frac{\partial U({\bf r},t)}{\partial t}=&A-(B+1)U+U^{2}V+D_{U}\nabla^{2}U({\bf r},t)\\ \frac{\partial V({\bf r},t)}{\partial t}=&BU-U^{2}V+\int_{0}^{t}dt^{\prime}{\cal K}_{V}(t-t^{\prime})D_{V}\nabla^{2}V({\bf r},t^{\prime})\end{array}\right..$$

$${\int }_0^{t}dz{\mathcal{K}}_V(z)D_V{\nabla }^{2}V(\mathbf{r},t-z).$$

$$\left\{\begin{array}[]{cc}\frac{\partial U({\bf r},t)}{\partial t}=&A-(B+1)U+U^{2}V+D_{U}\nabla^{2}U({\bf r},t)\\ \frac{\partial V({\bf r},t)}{\partial t}=&BU-U^{2}V+D_{V}\nabla^{2}\big{(}\int_{0}^{t}{\cal K}_{V}(z)V({\bf r},t)dz\\ &-\int_{0}^{t}z{\cal K}_{V}(z)\frac{\partial V({\bf r},t)}{\partial t}dz\big{)}\end{array}\right..$$

$$\mathcal{K}(t^{\prime })=\gamma e^{-\gamma t^{\prime }},$$

$$\left\{\begin{array}[]{cc}\frac{\partial U({\bf r},t)}{\partial t}=&A-(B+1)U+U^{2}V+D_{U}\nabla^{2}U({\bf r},t)\\ \frac{\partial V({\bf r},t)}{\partial t}=&BU-U^{2}V+D^{\prime}_{V}(t)\nabla^{2}V({\bf r},t)\\ &-D^{\prime\prime}_{V}(t)\nabla^{2}\frac{\partial V({\bf r},t)}{\partial t}\end{array}\right..$$

$$\left\{\begin{array}[]{cc}\frac{\partial U}{\partial t}=&A-(B+1)U+U^{2}V+D_{U}\nabla^{2}U\\ \frac{\partial V}{\partial t}=&BU-U^{2}V+D^{\prime}_{V}(t)\nabla^{2}V-BD^{\prime\prime}_{V}(t)\nabla^{2}U\\ &+D^{\prime\prime}_{V}(t)\nabla^{2}(U^{2}V)-D^{\prime}_{V}(t)D^{\prime\prime}_{V}(t)\nabla^{2}(\nabla^{2}V)\end{array}\right..$$

$$\left\{\begin{array}[]{cc}\frac{\partial U}{\partial t}=&A-(B+1)U+U^{2}V+D_{U}\nabla^{2}U\\ \frac{\partial V}{\partial t}=&BU-U^{2}V+D_{V}\nabla^{2}V-BD_{V}\tau_{V}\nabla^{2}U\\ &+D_{V}\tau_{V}\nabla^{2}(U^{2}V-D_{V}\nabla^{2}V)\end{array}\right..$$

$$U(\mathbf{r},t)=U_0+u(\mathbf{r},t),V(\mathbf{r},t)=V_0+v(\mathbf{r},t).$$

$$\left({{\begin{array}[]{cc}u\\ v\\ \end{array}}}\right)=\left({{\begin{array}[]{cc}c_{1}\\ c_{2}\\ \end{array}}}\right)exp(\lambda t+i{\bf k}\cdot{\bf r}),$$

$$\left({{\begin{array}[]{cc}a_{11}&a_{12}\\ a_{21}&a_{22}\\ \end{array}}}\right)\left({{\begin{array}[]{cc}c_{1}\\ c_{2}\\ \end{array}}}\right)={\bf 0},$$

$$a_{11}$$

$$a_{12}$$

$$a_{21}$$

$$a_{22}$$

$${\lambda }^2+g(k,{\tau }_V)\lambda +h(k,{\tau }_V)=0,$$

$$g(k,\tau_{V})=1-B+A^{2}+\big{[}D_{U}+D_{V}+D_{V}\tau_{V}A^{2}\big{]}k^{2}+D_{V}^{2}\tau_{V}k^{4},$$

$$\begin{gathered} h(k,{\tau }_V)=A^{2}B(1+D_V{\tau }_V{k}^2)+ \\ (1-B+D_U{k}^2)(A^2+D_V{k}^2+D_V{\tau }_V{A}^2{k}^2+D_V^2{\tau }_V{k}^4). \end{gathered}$$

$$k^2=-\frac{D_U+D_V+D_V{\tau }_V{A}^2}{2D_V^2{\tau }_V}.$$

$$k_{cH}=0,B_{cH}=1+A^2,$$

$$\left({{\begin{array}[]{cc}\dot{u}\\ \dot{v}\\ \end{array}}}\right)={\bf L}\left({{\begin{array}[]{cc}u\\ v\\ \end{array}}}\right)+{\bf N}+{\bf D},$$

$${\bf L}=\left({{\begin{array}[]{cc}B-1&A^{2}\\ -B&-A^{2}\\ \end{array}}}\right)$$

$$D_1=$$

$$D_2=$$

$$\displaystyle 2A\nabla^{2}(uv)+\frac{B}{A}\nabla^{2}u^{2}+\nabla^{2}(u^{2}v)-D_{V}\nabla^{2}(\nabla^{2}v)\big{)}.$$

$$\mathbf{L}={\mathbf{L}}_0+\mu {\mathbf{L}}_1+{\mu }^2{\mathbf{L}}_2+...,$$

Full text

I Introduction

The term “non-Markovian process” covers all stochastic processes with the exception of the small minority that happens to have the Markov property van Kampen . Despite the fact that stochastic processes in nature are generally non-Markovian, they are usually treated in the context of a Markov approximation. In this approximation, it is assumed that the correlation time between the system and environment is infinitely short so that memory effects can be neglected. Actually, treating non-Markovian processes with a Markov approximation can be as poor as, for instance, approximating highly nonlinear dynamical systems by a harmonic oscillator van Kampen . Since realistic systems typically possess a finite correlation or scattering time, they are non-Markovian by nature and considering memory effects is often inevitable.

Specifically, non-Markovian processes appear in many different fields including quantum optics Breuer ; Gardiner ; Lambropoulos , solid state physics Lai , quantum chemistry Shao ; Pomyalov , quantum information processing Aharonov ; Maniscalco , and even in the description of biological systems Rebentrost . Strongly coupled systems are also non-Markovian where the collisions can not be considered as instantaneous Prigogine2 . Recently, attention has been paid to memory effects in non-Markovian systems in different areas Schumann ; Trimper ; Serva ; Fritzsche ; Harris ; Dalmaroni . For instance, taking into account long-range memory in extreme events by using fractional Lévy processes Watkins ; Moloney allows one to make a time-dependent hazard assessment of future events based on events observed in the past Schumann . The effect of memory-dependent transport on the survivability of a population is investigated in Harris . The extension of conventional diffusive transport, within the framework of evolution equations, is considered in Trimper to take into account memory. Memory effects in correlated anisotropic diffusion are studied in nanoporous crystalline solids in Fritzsche . In this paper we take into account memory in reaction-diffusion systems to investigate its effect on the formation of patterns out of equilibrium.

Reaction-diffusion systems are extensively used in the study of self-organized phenomena that occur in open systems out of equilibrium Nicolis ; Cross ; Desai ; Kerner ; Buceta ; Mimura ; Dividsen . They are useful in many fields such as biology Harrison ; Meinhardt ; Murray1 ; Murray2 ; Kariga ; Chaplain ; Sherratt1 ; Sherratt2 ; Gatenby ; Kondo ; Shin , chemistry Kapral ; Dividsen2 ; Beta ; Castets ; Ruan ; Lengyel ; Gafiychuk ; Horvath ; Horvath2 ; Porjai ; Showalter ; Malchow , medicine Cameron ; Fenton ; Echebarria ; Alonso , neuroscience Liang ; Lefevre ; Tuckwell ; Habib ; Wu , physics Astrov ; Arecchi ; Staliunas ; Nakoa and ecology Holmes ; Skellam ; Gowda ; Meron . Reaction-diffusion systems were first proposed by Alan Turing in the study of morphogenesis Turing . Actually, Turing noticed that adding a diffusion term to a reaction system can drive the system to instability and plays an important role in the formation of patterns out of equilibrium. The characteristic feature of most of the studied reaction-diffusion systems is that the diffusion is considered to be a normal one. From microscopic point of view, the normal diffusion is derived from the (Markovian) master equation Reichl ; Srokowski .

The master equation is one of the most important equations in statistical physics that governs the dynamics of stochastic processes. As a system of stochastic variables evolves in time, transitions occur between various states of the system. To formulate the master equation, it is assumed that the probability of each transition depends only on the preceding time step and not on any previous history, which is exactly the Markov approximation. Since “normal” transport theory employs the master equation as its point of departure, normal diffusion does not take into account memory effects. The significance of the diffusion term in the formation of patterns and the importance of understanding and controlling pattern formation in far from equilibrium systems — an area of significant importance Bertram ; Mikhailov ; Zhang2 ; Cui ; Hu ; Hata ; Li ; Zhang3 ; Hu — motivates us to study pattern formation in reaction-diffusion systems in the presence of memory. Using analytical and numerical analysis, we show that memory can change the properties of spatio-temporal patterns. This includes the stability of spiral waves, which occur in many chemical and biological systems such as the heart Bar ; Fenton2 ; FentonH ; Bernus ; Dividsen . This opens up a new mechanism to control pattern formation using memory instead of external perturbations such as applying periodic forcing Zhang2 ; Cui or time-delayed feedback Hu ; Li ; Zhang3 .

The paper is organized as follows. In section II we construct non-Markovian reaction-diffusion systems by starting from the generalized master equation to take into account memory effects. Linear stability analysis is used in section III to investigate the effect of memory on the instability in a non-Markovian Brusselator model as a paradigmatic example of non-Markovian reaction-diffusion systems. In section IV, reductive perturbation method is exploited to study the behavior of the system near a Hopf instability and investigate the effect of memory on spatio-temporal patterns. The results indicate that a complex Ginzburg-Landau equation that depends on the memory governs the amplitude equation of the critical mode for short-term memories. Section VI contains the numerical solution of the obtained complex Ginzburg-Landau equation as well as direct numerical simulations of the non-Markovian Brusselator model to illustrate the effect of memory on the formation of spatio-temporal patterns. Finally, a summary and discussion is presented in section VI.

II Non-Markovian reaction-diffusion systems

The normal diffusion equation, ${\partial n(x,t)}/{\partial t}=D{\nabla^{2}n(x,t)}$, is a macroscopic equation where $n(x,t)$ is the density of the substrate at point $x$ and time $t$. This equation can be obtained starting from the master equation Reichl . In this section we are going to establish a general equation governing non-Markovian reaction-diffusion systems. For this purpose, we first review how to introduce a non-Markovian diffusion equation by starting from the generalized master equation.

Various methods have been used to obtain the generalized master equation, however, the procedure of Zwanzig Zwanzig ; Zwanzig2 distinguishes itself by its elegance and economy of effort. His procedure was based on projection operator techniques Nakajima ; Zwanzig ; Zwanzig2 ; Prigogine , which are widely used in open systems and non-equilibrium statistical mechanics Haake ; Balescu ; Grabert ; Kubo . Zwanzig started from the Liouville-von Neumann equation for the density operator, and obtained the generalized master equation Zwanzig ; Zwanzig2 ; Kenkre ; Srokowski

[TABLE]

where $P_{\xi}$ are the diagonal elements of density matrix and denote the probability of finding the system in state $\xi$. $W_{\xi\mu}$ are the transition rates from state $\mu$ to $\xi$ that generally are time dependent. Integration over all previous times indicates the non-Markovian nature of the processes. Since this equation is hard to solve, physicists use Markov approximation to simplify this equation. By this approximation, they ignore the dependence of transition rates to previous times Kenkre and consider

[TABLE]

where $F_{\mu\xi}$ are the transition rates in Markovian case which are time independent. By neglecting the dependence on previous times they assume that there is no memory in the system. As a result, the equation (1) reduces to the well known master equation,

[TABLE]

Applying the nearest-neighbor approximation on the transition rates and taking the continuum limit Reichl ; Kenkre ; Srokowski , one can derive

[TABLE]

where $P(x,t)$ is the probability of finding a particle in $x$ at time $t$ and $D$ is the diffusion coefficient. Generalizing this equation to three dimensions and noting that the density of the substrate is proportional to this probability, one can obtain the normal diffusion equation.

Now, if we want to be more precise and consider realistic phenomena beyond the Markovian assumption, we should start from the generalized master equation (1). However, it is hard to work with this equation. A simplification may often be possible whereby the time dependence of $W_{\xi\mu}$ is assumed independent of the states $\xi$ and $\mu$ Kenkre :

[TABLE]

Using this assumption, (1) can be rewritten in a simpler form

[TABLE]

The nature of transport is thus decided by the memory ${\cal K}$ as well as by the rates $F_{\xi\mu}$. By applying the nearest-neighbor approximation on the transition rates $F_{\mu\xi}$ and taking the continuum limit, one obtains the non-Markovian Fokker-Planck equation Kenkre

[TABLE]

where ${\cal K}(t^{\prime})$ is the memory kernel. This equation expresses the fact that an interaction process is in general a non instantaneous event and thus nonlocal in time. The memory kernel describes the finiteness of the correlation or scattering time Manne ; Kenkre2 . Obviously, for the Markovian case, the memory kernel is ${\cal K}(t^{\prime})=\delta(t^{\prime})$ and (7) converts to (4). In the study of realistic systems one encounters different forms of the memory kernel. For instance, the excitation energy transfer between molecules indicates the presence of an exponential memory kernel Kenkre . However, the description of reactions in polymer systems needs more complicated kernels Cherayil ; Bandyopadhyay .

Using (7), we can write a general two component non-Markovian reaction-diffusion system as

[TABLE]

which is a non-Markovian activator-inhibitor system. Here, $f_{1}$ and $f_{2}$ are reaction terms. $U$ and $V$ are the densities of activator and inhibitor, respectively. $D_{U}$ and $D_{V}$ are their corresponding diffusion coefficients, and ${\cal K}_{U}$ and ${\cal K}_{V}$ are the memory kernels associated with them. $\mu$ is the bifurcation parameter and, as $\mu$ varies, the system might move from a steady state to an oscillating or a patterned state via a Hopf or Turing instability, respectivelyKuramoto .

III Linear stability analysis

In the previous section we constructed a general two component non-Markovian reaction-diffusion system. In this section we use linear stability analysis to see how the memory affects the critical behavior of the system. In order to study the problem analytically in more detail, we consider the Brusselator model with memory in the diffusion of the inhibitor. The Brusselator model is a two-component activator-inhibitor system and it is one of the paradigmatic models for nonlinear chemical systems Reichl . This system in the presence of memory, in the diffusion of the inhibitor, takes the form

[TABLE]

Here, $U$ and $V$ are chemical concentrations that can vary in space and time. $A$ and $B$ are constants and B is considered as the control parameter of the system to generate patterns. The steady state solution of the non-Markovian Brusselator model, (9), is $(U_{0},V_{0})=(A,B/A)$. As the system (9) is difficult to treat analytically in the presence of general memory and also many systems exhibit short-term memory, we focus on this latter case in order to present an analytical study of non-Markovian reaction-diffusion systems. To this end, we use a change of variable $z=t-t^{\prime}$ to rewrite the integral in (9) as

[TABLE]

In the presence of short-term memory only small $z$’s are important because the memory kernel, ${\cal K}_{V}(z)$ is almost zero for large $z$’s. So, we can expand $V({\bf r},t-z)$ up to the first order in $z$, $V({\bf r},t-z)=V({\bf r},t)-z\frac{\partial V({\bf r},t)}{\partial t}$, and rewrite (9) as

[TABLE]

To proceed further, we need to consider a specific form for the memory kernel. A widely used form for the memory kernel function is an exponential one Srokowski ; Kenkre

[TABLE]

where $\gamma$ measures the reciprocal of the characteristic scattering time or in other words the memory kernel decay time $(\tau=1/\gamma)$. Note that, the kernel becomes the delta function in the limit $\gamma\rightarrow\infty$ (the Markovian case). By considering this memory kernel, (10) takes the form

[TABLE]

Here, $D^{\prime}_{V}(t)=D_{V}(1-e^{-t/\tau_{V}})$, $D^{\prime\prime}_{V}(t)=D_{V}(\tau_{V}-(t+\tau_{V})e^{-t/\tau_{V}})$ and $\tau_{V}=1/\gamma_{V}$ is the characteristic memory kernel decay time for inhibitor. $D^{\prime}_{V}(t)$ and $D^{\prime\prime}_{V}(t)$ are monotonically increasing functions of time with $0<D^{\prime}_{V}<D_{V}$ and $0<D^{\prime\prime}_{V}<D_{V}\tau_{V}$. In the second relation in (12), $\frac{\partial V({\bf r},t)}{\partial t}$ has appeared on both sides. We replace the right one by the full expression for $\frac{\partial V({\bf r},t)}{\partial t}$. Since $0<D^{\prime\prime}_{V}<D_{V}\tau_{V}$ and $\tau_{V}$ is small, if $D_{V}$ is not so large, we can truncate the process to obtain

[TABLE]

We are interested in the long time behavior of the system to study the formation of patterns. On the other hand, in the case of short-term memory where $\tau_{V}$ is small, the exponential term in $D^{\prime}_{V}(t)$ and $D^{\prime\prime}_{V}(t)$ goes to zero rapidly. Therefore, equation (13) can be approximated by

[TABLE]

In order to perform a linear stability analysis, small perturbations, $u({\bf r},t)$ and $v({\bf r},t)$, about the steady state $(U_{0},V_{0})$ are considered

[TABLE]

Inserting the above in Eq. (14) and choosing the linear part, we obtain a linear equation governing the dynamics of perturbations. Using the normal mode ansatz

[TABLE]

in the linearized equations leads to

[TABLE]

where,

[TABLE]

$\lambda$ is the perturbation growth rate, ${\bf k}$ is the wave vector and $k=|\bf{k}|$. The characteristic equation in the presence of memory takes the form

[TABLE]

where

[TABLE]

and

[TABLE]

Using the characteristic equation, one can find the perturbation growth rate, $\lambda$, in terms of $k$. Wave vectors that result in a nonnegative real part of the growth rate $\lambda$ are critical. Note that this equation depends on the memory properties of the system and in the limit of no memory ($\tau_{V}\rightarrow 0$) we recover the characteristic equation for the Markovian Brusselator model Torabi .

The steady state is linearly stable if and only if both $g(k,\tau_{V})$ and $h(k,\tau_{V})$ are non-negative for all $k$. Clearly, this stability condition can be violated in either of the following two ways:

1. $h(k,\tau_{V})$ vanishes for some $k$, but $g(k,\tau_{V})$ and $h(k,\tau_{V})$ remain positive for all $k$. This condition together with ${\partial h(k,\tau_{V})}/{\partial k}|_{k=k_{cT}}=0$ determine the critical parameter, $B_{cT}$, and also $k_{cT}$ for a Turing bifurcation.

2. $g(k,\tau_{V})$ vanishes for some $k$, but $g(k,\tau_{V})$ and $h(k,\tau_{V})$ remain positive for all $k$. This condition together with ${\partial g(k,\tau_{V})}/{\partial k}|_{k=k_{cH}}=0$ leads to the critical values $k_{cH}$ and $B_{cH}$ for a Hopf bifurcation. The frequency of oscillation is $\omega=\pm\sqrt{h(B_{cH},k_{cH})}$.

Focusing on the latter, we find that ${\partial g(k,\tau_{V})}/{\partial k}=0$ leads to a solution of $k=0$ as well as

[TABLE]

Since the right hand of (19) is always negative, it is not an acceptable solution and we obtain that

[TABLE]

which means that the critical values for Hopf instability do not change when memory is present in the diffusion of inhibitor. On the other hand, since the characteristic function depends on the memory property of the system, the growth rates of the modes depend on memory as well. This affects the modulation of unstable modes near the Hopf bifurcation and, thus, influences the formation of spatio-temporal patterns. In fact, linear analysis only provides us with insight about the critical behavior of the system. However, in order to investigate the effect of memory on pattern formation, we need to go beyond the linear analysis. In the next section, we use a reductive perturbation method to investigate the effect of memory on the formation of spatio-temporal pattern near a Hopf instability.

IV Reductive perturbation method

Near a bifurcation the neighboring modes of the critical mode, which have decay times of the same order of magnitude as the critical mode, play an important role in the long-time behavior of the system that results in the generation of patterns. According to the slaving principle Kuramoto , near a bifurcation we are left with a couple of relevant dynamical variables whose time scales are distinguishably slower than the other dynamical variables so that the latter can be eliminated adiabatically using rescaled space-time coordinates. In this section, we investigate the non-Markovian Brusselator model analytically using the reductive perturbation method. The results indicate that memory changes the properties of the spatio-temporal patterns.

Starting from equation (14), considering small perturbations about the steady state and keeping nonlinear terms as well as the linear terms results in the equation governing the fluctuations as

[TABLE]

where

[TABLE]

is the Jacobian matrix, {\bf N}=\left({{\begin{array}[]{cc}\beta\\ -\beta\\ \end{array}}}\right) is the nonlinear part with $\beta=\frac{B}{A}u^{2}+2Auv+u^{2}v$ and {\bf D}=\left({{\begin{array}[]{cc}D_{1}\\ D_{2}\\ \end{array}}}\right) is the diffusion term with

[TABLE]

We introduce a new bifurcation parameter by considering $\mu=\frac{B-B_{cH}}{B_{cH}}$. By this definition the steady state is stable for $\mu<0$ and is unstable for $\mu>0$. Near the instability, ${\bf L}$, ${\bf N}$, ${\bf D}$, \left({{\begin{array}[]{cc}u\\ v\\ \end{array}}}\right) and the eigenvalues of the Jacobian matrix can be expanded in powers of $\mu$ as

[TABLE]

where $\lambda_{0}=i\omega_{0}=iA={\bf U_{L}}{\bf L_{0}}{\bf U_{R}}$ and $\lambda_{i}=\sigma_{i}+i\omega_{i}={\bf U_{L}}{\bf L_{1}}{\bf U_{R}}$. The vectors ${\bf U_{R}}$ and ${\bf U_{L}}$ are right and left eigenvectors of ${\bf L_{0}}$, respectively. ${\bf L_{0}}$ and ${\bf L_{1}}$ are given by

[TABLE]

and the right and left eigenvectors of ${\bf L_{0}}$ are

[TABLE]

We must get inside the neighborhood of the critical mode to consider just the important modes based on the slaving principle. We do this by rescaling the time and space variables as Kuramoto ; Torabi

[TABLE]

where $\varepsilon^{2}\chi\equiv\mu$ and $\chi=\rm sgn(\mu)$. Since we are interested in the supercritical Hopf bifuracation, $\chi$ is positive. Reductive perturbation method is a multi-scale method, so, we consider \left({{\begin{array}[]{cc}u\\ v\\ \end{array}}}\right) as a function of $t$, $T$, $R$. This means that we are dealing with the long-time, long-wavelength modes in their natural variables $T$ and $R$, and reserving $t$ for the overall periodic motion (limit cycle). Substitution of (IV) into (21) and equating coefficients of different powers of $\varepsilon$, yields to a set of equations in the form of

[TABLE]

where the first three $\bf B$’s are

[TABLE]

and ${\bf\hat{D}}$ is given by

[TABLE]

Note that, in the limit of $\tau_{V}\rightarrow 0$ this matrix reduces to \left({{\begin{array}[]{cc}D_{U}&0\\ 0&D_{V}\\ \end{array}}}\right), which is identical to ${\bf\hat{D}}$ in Markovian Brusselator model Torabi . The solution for $\nu=1$ is

[TABLE]

where $c.c.$ stands for the complex conjugate and $W(T,R)$ is a complex amplitude to be determined. Using this solution in the equation for $\nu=2$ leads to an expression for \left({{\begin{array}[]{cc}u_{2}\\ v_{2}\\ \end{array}}}\right). Putting these two solutions into the equation for $\nu=3$ together with the solvability condition for the set of equations (34) Kuramoto , results in the equation governing the amplitude $W$

[TABLE]

where $\lambda^{(1)}$, $d$ and $g$ are generally complex numbers and are given by

[TABLE]

[TABLE]

[TABLE]

Equation (35) is the well known complex Ginzburg-Landau equation, but note that the coefficient $d$ depends on memory kernel decay time. In fact, the presence of non-Markovian diffusion of inhibitor results in the change in real and imaginary part of the coefficient $d$. For $\tau_{V}\rightarrow 0$, we recover the complex Ginzburg-Landau equation in Markovian Brusselator model Kuramoto ; Torabi . With a redefinition as follows

[TABLE]

the complex Ginzburg-Landau equation above the criticality can be written in a more convenient form (after dropping the primes)

[TABLE]

where

[TABLE]

and

[TABLE]

Since the parameter $c_{1}$ depends on the memory property of the system (see (40)), one can change the stability of solutions of the complex-Ginzburg Landau equation by changing the memory property of the system and, thus, control the resulting spatio-temporal pattern. In fact, by looking at the phase diagram of the complex Ginzburg-Landau equation Aranson , one can see how the solutions of the complex Ginzburg-Landau and their stability change over the parameter space. By introducing memory or changing its properties (here by varying $\tau_{V}$), we can cross the boundaries of different regimes in any spatial spatial dimension, which enables us to change the selected spatio-temporal patterns. To explicitly illustrate how to achieve this and also to provide further support for our analytical results, we exploit numerical simulations in the next section using $\tau_{V}$ as our control parameter.

V Numerical simulations

In this section, we perform numerical simulations of the complex Ginzburg-Landau equation obtained analytically in section IV as well as direct numerical simulations of the non-Markovian Brusselator model (9), with the exponential memory kernel (11). For integrating the complex Ginzburg-Landau equation we use a pseudospectral method to perform numerical computations in Fourier space based on the method of exponential time differencing (ETD2) Cox . Small-amplitude random initial data about $W=0$ and periodic boundary conditions are used. For the direct numerical simulation of the non-Markovian Brusselator model, we also consider small amplitude random initial condition about the steady state and use periodic boundary condition. We chose to focus on two spatial dimensions, which displays a rich spectrum of dynamical behavior including spiral waves.

As a specific example, we focus on the stability of spiral waves. Fig. 1a shows rotating spiral waves that are generated in the absence of memory for the Markovian approximation of the Brusselator model, which is achieved by considering a delta function for the memory kernel in (9). For the selected parameters and no memory, equations (40)-(41) give rise to $(c_{1},c_{2})=(-0.33,0.96)$ for the coefficients of the corresponding complex Ginzburg-Landau equation (39). As expected based on our analytical calculations, for these parameters the complex Ginzburg-Landau equation gives basically identical results as can be seen in Fig. 1c. Note that the change in scales between Fig. 1a and Fig. 1c is a consequence of rescaling (38) in obtaining (39). Most importantly, $(c_{1},c_{2})=(-0.33,0.96)$ belong to the spiral bound states in parameter space Aranson , consistent with the observed behavior. On the other hand, if we introduce memory in the Brusselator model, direct numerical simulations of (9) indicate that the selected pattern changes and spiral waves break up (see Fig.1b). In fact, by introducing memory with $\tau_{V}=1/\gamma_{V}=1/3.5$, the coefficient $c_{1}$ of the complex Ginzburg-Landau equation changes and the new coefficients $(c_{1},c_{2})=(-0.97,0.96)$ belongs to the defect turbulence regime beyond the absolute instability line where spiral waves are unstable Aranson ; AransonA . This is directly confirmed by Fig. 1d providing further support for our analytic calculations above. Thus, $\tau_{V}$ can be used as a control parameter to select stable or turbulent spiral patterns.

Changing the memory property of the system ($\tau_{V}$) and consequently changing the coefficient $c_{1}$ in the complex Ginzburg-Landau equation, not only helps us to break up spiral waves but it also enables us to change the properties of spiral patterns as well. To this end, one can start from the monotonic range Aranson ; AransonK and cross the oscillatory range line Aranson ; AransonK by changing $c_{1}$ to end up in the bound state regime. Fig.2 illustrates the behavior of the complex Ginzburg-Landau equation, obtained in (39), for two sets of coefficients $(c_{1},c_{2})$ with different $c_{1}$. Fig. 2a shows $Re(W)$ for $(c_{1},c_{2})=(0.2,0.82)$ in the monotonic range where pairs of spirals form that drift and annihilate by interacting with each other. However, Fig. 2b corresponds to $(c_{1},c_{2})=(-0.33,0.82)$ in the spiral bound state regime of the oscillatory range. In order to illustrate the difference between the properties of the solutions, $|W|$ is shown in Fig. 2c and Fig. 2d. There are no shock lines in Fig. 2c for the monotonic range in contrast to Fig. 4d for the oscillatory range where cellular structures are formed.

The change in the selected spatio-temporal pattern, due to changing the memory property of the system, is not limited to spiral patterns. One can change the coefficient $c_{1}$ in complex Ginzburg-Landau equation to make a transition, for instance, from defect turbulence regime to phase turbulence regime as well. Also, by focusing on the solution of complex Ginzburg-Landau equation in 1D or 3D, we can study the effect of memory on plane waves or vortex filaments, respectively.

VI Summary and discussion

Realistic systems possess a finite correlation or scattering time and are generally non-Markovian. In spite of this fact, they are usually treated in the context of Markov approximation where memory is neglected. We studied the effect of memory on the instability and spatio-temporal pattern formation in reaction-diffusion systems and showed that considering memory near the instability might be inevitable. For this purpose, we constructed a general non-Markovian reaction-diffusion system by starting from the generalized master equation under the assumption (5), to take into account memory effect.

We considered the Brusselator model as a typical reaction-diffusion system in the presence of memory in the diffusion of the inhibitor. In the limit of short-term memory and for small diffusion coefficients of the inhibitor, the non-Markovian Brusselator model reduced to a set of equations that could be solved analytically. Linear stability analysis indicated that, in the presence of memory, the characteristic equation depends on the characteristic memory kernel decay time. Then, we used the reductive perturbation method near the Hopf instability to investigate the effect of memory on the formation of spatio-temporal patterns. We found that a memory dependent complex Ginzburg-Landau equation governs the amplitude of the critical mode. Since the derived complex-Ginzburg Landau equation was memory dependent, we can change the stability of its solutions by changing the memory property of the system and, thus, control the selected spatio-temporal pattern. These analytical findings were confirmed by numerical simulations of both the non-Markovian Brusselator model and the complex Ginzburg-Landau equation. The results indicate that going beyond the Markovian approximation might be necessary to study the formation of spatio-temporal patterns in many cases and also opened up a new window to the control of these patterns using memory.

We would like to point out that the analytical approach we established here is general. In particular, it is not limited to the non-Markovian Brusselator model. In fact, it can be used to study other reaction-diffusion systems beyond the Markovian approximation in any dimensions. While we showed that considering short-term memory in the diffusion of inhibitor can result in a significant impact on the selected spatio-temporal patterns, considering memory in the diffusion of the activator can lead to an even richer dynamical behavior. In the latter case, the right hand side of equation (19) changes and would not necessarily be negative. Therefore, for some parameters, the system can admit nonzero solutions for critical wave numbers ($k_{c}\neq 0$) where considering them in the system demands more investigations. In addition, one can not readily use reductive perturbation method to study the system in the presence of any nonzero modes. It is because reductive perturbation method is based on the coarse graining of the system to consider important long wavelength modes (with $k=0$ and its neighboring modes). Similarly, while we largely focused on the case of a Hopf bifurcation, the presented framework for studying reaction-diffusion systems in the limit of short-term memories directly gives the conditions for a Turing bifurcation (See equation (16) for the case of Brusselator model). This allows one to investigate the effect of short-term memories on Turing patterns using weakly nonlinear analysis Zhang ; Pena . The mathematical treatment of this case is challenging because of the complicated equations beyond the Markov approximation. Since the main goal of this paper is to provide a proof of concept of how memory can affect the formation of spatio-temporal patterns in analytically tractable way, the detailed investigation of memory in the diffusion of the activator as well the effect of memory on Turing patterns remains an interesting topic for the future.

JD was financially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). RT thanks the University of Calgary for their hospitality.