Spatiotemporal patterns in the Turing-Takens-Bogdanov scenario

(1)

Spatiotemporal patterns in the Turing-Takens- Bogdanov scenario

Pablo Moreno Spiegelberg

Master’s Thesis

Master’s degree in Complex systems at the

UNIVERSITAT DE LES ILLES BALEARS

Academic year 2018-2019

28/10/2019

Master’s Thesis Supervisor: Damià Agustí Gomila Villalonga

(2)

Abstract

Some spatial dinamical systems exhibit, for close values of the parameter, diffusion drive instability (Turing bifurcation) and a Homoclinic bifurcation of the homogeneous solution. However, the interaction between these bifurcations has not been studied in detail in the literature. In this thesis we explore the interaction between a Turing and a Homoclinic bifurcation in a Reaction-Diffusion system. For this purpose we incorporate a diffusion term to the normal form for the Cusp Takens-Bogdanov codimension-3 point, in such a way that a Turing instability might occur. We analyse the spatio-temporal bifurcations and their interactions. These bifurcation curves converge in a new high codimension point, that we call Turin-Takens-Bogdanov point. The system shows a wide variety of stable solutions such as steady patterns, homogeneous oscilatory states , different more complex spatio-temporal periodic solution, pseudo-periodic states and turbulent regimes.

(3)

Acknowledgements

I would like to thank my advisor Dr. Dami`a Gomila for his guidance, advice and patience. This work would not have been possible without him. I also would like to thank Dr. Manuel Matias and Dr. Andreu Arinyo for their help, collaboration and valuable comments during the development of this thesis. I would also like to acknowledge to all the IFISC members their teachings during this last year. And last but not least, I am grateful to my family, friends and classmates for all the love and support they have given me.

(4)

1 Introduction

In 1952 Allan Turing first noticed an anti-intuitive mathematical behaviour of reaction-diffusion models [1]. In some cases, adding diffusion to a dynamical system with a stable homogeneous solution, instead of acting as a stabilizing term, destabilizes the homogeneous state generating a inhomogeneous solution.

This process is usually know as Turing instability. The conditions for this phe- nomenon to occur were difficult to replicate in a laboratory, and it was not until 1990 that the inhomogeneous structures where found in a chemical system [2].

Turing-like models are used nowadays as a theoretical framework to understand pattern formation in chemistry system, morphogenesis [3], ecological systems [4], and in the development of technological devices [5].

In some cases a system displaying a Turing instability shows also oscillatory instabilities (Hopf) of the homogeneous solution. The Turing-Hopf codimension- 2 point leads to complex spatiotemporal patterns such as bistability regions of steady spatial patterns and homogeneous oscillations, or hybrid non-homogeneous oscillatory solutions. General analytic approaches to Turing-Hopf scenario have been studied in [6][7][8] and the predicted spatiotemporal patterns have been found experimentally [9].

On the other hand the hand, oscillatory temporal systems can also undergo a Homoclinic bifurcation as a result of a more general scenario. The Homoclinic bifurcation, also known as saddle loop, is a global bifurcation in which a limit cycle is destroyed when it collides with a saddle fixed point. This process has been observed in systems such as the forced Van der Pol oscillator [10] and plays a main role in the transition to chaos in the Lorenz system [11]. Homoclinic bifurcation is also used in neural models to describe type I excitability [12].

Both, Turing and Homoclinic bifurcations are present in systems such as spatial predator-pray models [13] or models for Posidonia Oceanica meadows [14]. ThePosidonia Oceanica is a marine clonal plant endemic of the Mediter- ranean Sea. It forms extensive meadows that define the dominant ecosystem on the Balearic coasts. In [14] the meadows dynamics are modeled using a system of partial diferential equations. This proposed model presents a Hopf and a Ho- moclinic bifurcation, responsible for the creation and destruction of oscillatory solution in the meadow; and a Turing instability, responsible for the creation of spatial patterns, (see Fig. 1). The role of the interaction between the Turing and the Homoclinic bifurcation in the creation of more complex spatiotemporal structures has not been studied in detail in the literature.

Motivated by these systems we will study the interaction between a Tur- ing instability and a Homoclinic bifurcation of the homogeneous solution in a Reaction-Diffusion system. We will start from the Cusp Takens-Bogdanov codimension 3 point normal form given in [15]. Adding an appropriate diffusion term to the system we will reach a general system where both, a Homoclinic bifurcation and a Turing instability, coexist. For proper parameters, the starting point of both bifurcation curves collide in a high codimension point, the Turing-Takens-Bogdanov point, where a SN, Turing, Hopf and Homoclinic bifurcation curves converge. In this thesis we will explore the dynamics of the system around this point.

(6)

Figure 1: Phase diagram of the equation proposed in [14], using the adimensional sensitivity of the plant in the y-axis and the mortality in the x-axis. The figure shows the bifurcation curves of the system, including Hopf, Homoclinic and Turing bifurcations.

2 Temporal system

To study the interaction of Turing instability with an homoclinic global bifurcation we will use the system of ODEs given in [15] for the focus case to describe the local dynamics of the different spieces in the Reaction-Diffusion model:

˙ u=v

˙

v=−u³+µ2u+µ1+v(ν+bu−u²) (1) where u(t) and v(t) are two real time dependent variables;µ₁, µ₂, b andν are real parameters; and the dot is the total derivative with respect to time.

This system of ODEs is the normal form of the globally stable codimension 3 point Cusp Takens-Bogdanov. It gives therefore a globally stable scenario with a very rich parameter space, including global bifurcation such as the Homoclinic bifurcation, (see Fig.2).

For simplicity, as the system has 4 different parameters we will work fixing b and ν and working onµ1, µ2planes.

2.1 Bifurcation analysis

The complete bifurcation diagram of the system given by Eq. (1) is really com- plicate, (see Fig. 2). A complete rigorous analysis of it can be consulted on [15].

In this chapter I will give a characterization of the bifurcations present in the parameter region of study and that will interact with the spatial instabilities.

2.1.1 Fixed points and Cusp Bifurcation

The creation and destruction of fixed points in the proposed system occur at two Saddle Node bifurcations which collide in a Cusp codimension 2 point. The fixed points of the system can be obtained setting to zero the Eq. (1) getting:

(7)

Figure 2: Complete bifurcation diagram of Eq. (1) [15].

v^∗= 0 (2)

−u^∗3+µ₂u^∗+µ₁= 0 (3) Using the discriminant of Eq. (3) the two SN branches and the cusp codimension 2 point where they meet can been obtained:

SN :µ1=±2 rµ³₂

27 (4)

Cusp:µ1=µ2= 0 (5)

This bifurcation divide the parameter space in two different regions:

In the region with|µ1|<2 qµ³₂

27 Eq. (3) has 3 different roots and, therefore, the system has 3 different fixed points. The fixed points with higher and smaller value ofu^∗ are stable while the other one is a saddle fixed point.

On the other hand, for |µ1| < 2 qµ³₂

27 or µ₂ < 0, the system has a single stable fixed point.

The linear stability of the fixed points is given by the eigenvalues of the Jacobian matrix. The Jacobian matrix, applying the condition (2) is given by:

J|v=0=

0 1

−3u^∗2+µ₂ ν+bu^∗−u^∗2

(6) whose eigenvalues are given by:

λ_1,2=τ±√

τ²−4∆

2 (7)

(8)

Figure 3: This figure shows the Hopfs bifurcations for on the{µ1, µ2}parameter spaces for b= 2.05 and ν =−1 on the left side; andb= 0.5 and ν = 2 on the right side.

whereτ and ∆ are respectively the trace and the determinant of the Jaco- bian:

τ=ν+bu^∗−u^∗2 (8)

∆ = 3u^∗2−µ2 (9)

The values ofu^∗ in these equations must fulfil the condition of fixed points given by (3). It is also important to note that, on the SN bifurcation, ∆ = 0 for the involved fixed point, and therefore one of its eigenvalues is zero.

2.1.2 Hopf Bifurcation

A Hopf bifurcation occurs when a fixed point change its stability while a pair of complex conjugate eigenvalues cross the imaginary axis. At this moment a limit cycle is created or destroyed.

These conditions are reached when, for a given fixed point, τ changes its sign while ∆>0. Asτ will be zero on the Hopf bifurcation we can obtain from Eq. (8) an expression for the value ofu^∗ at the bifurcation:

uH± =b±√ b²+ 4ν

2 b²>−4ν (10)

As the fixed point should also fulfil Eq. (3) the following expression for the Hopf bifurcation is deduced:

µ1H=u³_H−uHµ2H; µ2H<3u²_H (11) where the condition forµ_2H arise imposing ∆>0 on Eq. (9). As it can be seen on Fig. 3, the two Hopf bifurcations will emerge from the same SN branch whenν <0 (Fig. 3 left) and from different branches whenν >0 (Fig. 3 right)

In the region of interest in this thesis, b = 2.4,ν =−1, µ₁ ∈ (−1,1) and µ2∈(−2,2), the present Hopf bifurcation (from now onHopf1) is supercritical, which means that a stable limit cycle is created when the fixed involved point loses stability.

(9)

Figure 4: This figure shows the closer to Cusp Homoclinic bifurcation forb= 2.4 andν =−1.

2.1.3 Takens-Bogdanov codimension-2 point

The Hopf bifurcation line ends in one of its edges on the Saddle Node bifurcation in a Takens-Bogdanov codimension 2 point, also call Double zero bifurcation.

At this point the two eigenvalues become zero.

The location of the TB point on the{µ1, µ2}parameter space can be derived from previous equations. Using the equality ∆ = 0 and the equations (9) and (10) we obtainµ2T B= 3u²_H. Replacing it on Eq. (11) we obtainµ1T B=−2u³_H. In the neighbourhood of the TB codimension 2 point not only the Hopf and the SN lines emerge, but also a Homoclinic bifurcation where the limit cycle created at the Hopf is destroyed [16]. The complete unfolding of this point, for the parameter region of study, is shown at Fig. 4 [16].

2.1.4 Homoclinic bifurcation

A Homoclinic bifurcation is a global bifurcation in which a limit cycle touches a saddle fixed point creating a homoclinic orbit with infinite period. After the bifurcation the limit cycle no longer exist [17].

In our system there is not an easy way to obtain an expression for the Homo- clinic bifurcation. Nevertheless the curve can be obtain by means of numerical simulations. The obtained curve, shown in Fig. 4, starts at one side from the TB point and ends colliding with the SN in a Saddle-node Separatrix-loop. From this point on, the Saddle Node acts as a SNIC bifurcation [15].

3 Spatial System

By adding a spatial dependence to Eq. (1) we attempt to reach a system of equations where the homogeneus solutions (fixed points of Eq. (1)) lose the stability via a Turing instability. This effect can be obtained by adding a diffusion matrix to Eq. (1). In general this matrix is not diagonal in the basis {u,v}.

We choose to use a diagonal matrix because most part of models in chemistry, biology and physics have a diagonal diffusion.

In principle a diagonal matrix does not guaranty the observation of a Turing instability but, fortunately, as we will see in Section 3.2.1, just by doing linear

(10)

transformation of variables we can obtain a basis where a diagonal diffusion allows a Turing instability.

3.1 Diffusion driven instability.

Two components Reaction-Diffusion systems are given by a PDE of the form:

∂u(x, t)

∂t =f_µ(u) +Db· ∇²u (12)

where uis a 2-vector field which gives the concentration of the 2 different reagents in space and time, the local reaction between reagents is given by f_µ :R²⇒R²,∇²is the Laplace operator andDis a diagonal diffusion operator of the form:

Db =

D₁ 0 0 D₂

(13) where D₁ and D₂ are the diffusion coefficients of the reagents. These coefficients can be transformed to dimensionless units by rescaling the space as x⁰=√

D₁x. Removing the tildes we obtain:

Db = 1 0

0 r

(14) wherer= ^D_D²

1 is a positive variable which gives the ratio between the diffusion coefficients of the species.

Assumingf_µ(u^∗) = 0 then u^∗ is an stationary uniform solution for a given set of parameters µ. Taking u(x, t) = up(x, t) +u^∗ where up(x, t) is a small perturbation around the uniform solutionu^∗, the evolution of the perturbation takes the form:

δup(x, t)

δt =

Jµ(u^∗) +D∇²

up+O(u²_p) =Lup+O(u²_p) (15) whereJµ is the Jacobian matrix of (12).

The operatorLcan be diagonalized by blocks on the basisn e^iqx

0

, 0

e^iqx o

whereq∈_2πn

L , n∈Z when working on a system of size L with periodic boundary conditions. On this basis the linear evolution operator can be written as L=diag({Lq}) with Lq =

Jµ^∗+

−q² 0 0 −rq²

. Doing thisup(x, t) can be written on the linear approximation:

u_p(x, t) =X

q

hC_q1ϕ_q1e^σ¹^(q)t+C_q2ϕ_q2e^σ²^(q)ti

e^iqx (16)

Whereϕq1,2 are the eigenvectors ofLq with associated eigenvaluesσ1,2(q) andCq1,2 are given by the initial conditions of the fields.

IfRe[σi(q)]<0∀i, q any perturbation will decrease in time and the homogeneous solution will be linearly stable. On the other hand, if for some q and i Re[σi(q)]>0 the homogeneous solution will be unstable. Studying the changes in sign ofRe[σ_i(q)] as function of the parameters we can find bifurcations that involve periodic pattern formation.

(11)

Figure 5: Growth rate for a Turing instability.Re[σi(q)] is shown as function of q aroundqc before, on and after the instability. [18]

In this thesis we are interested in a particular case of instability, know as Turing instability or Type-I instability [18]. In this instabilityRe[σi(q)] become positive at a critical wave number|qc|>0 while increasing a control parameter p, having Re[σi(q)] a local maximum at qc when p = pc. For pslightly less than pc Re[σi(q)] < 0 at the neighbour of qc, for p= pc Re[σi(qc)] = 0 and forpslightly larger thanpc theRe[σi(q)]>0 for qin the neighbour ofqc, (see Fig. 5).

In systems with stabilizing higher order terms the unstable modes saturate and a spatial periodic steady state is created.

In a two components Reaction-Diffusion system, described by an equation of the form (12) and with diagonal positive diffusion matrix, a Turing instability might happen only if the signs of the Jacobian matrix have the form [6]:

sign(J) =

+ +

− −

or

+ −

(17)

3.2 Requirements for a Turing Instability. Change of vari- ables.

The Jacobian of the temporal system proposed during the first part of this thesis does not fulfil the condition set in (17). Anyway, by a change of variables it is possible to transform the ODE in such a way that its Jacobian has the configuration given in (17). For simplicity we will consider only the change of thev variables given by:

v=w+u (18)

This change of variables will not affect any of the equations introduced during the the analysis of the temporal system with the exception of the Jacobian matrix (6). Implementing the change of variables on Eq. (1) and adding a diffusion of the form Eq. (14) the following set of equations is obtained:

∂u

∂t =w+u+∇²u

∂w

∂t =−u³+µ2u+µ1+ (w+u)(ν−1 +bu−u²) +r∇²w (19)

(12)

With homogeneous solutions given by Eq. (3) andw^∗=−u^∗. The Jacobian matrix of the homogeneous solutions is given by:

J=

1 1

−∆ +τ−1 τ−1

(20) Where τ and ∆ are still given by Eq. (8) and (9). According to (17) the system might present a Turing instability whenτ−1<0.

3.3 Linear stability. Turing instability and revised Hopf bifurcation.

In this section we study the bifurcations of the homogeneous solutions of the spatial system (given by (19)). We study the grown rateσ(q) of perturbations with wave number q. This growth rate can be obtained as the eigenvalues of the matrix L_q:

L_q=

1−q² 1

−∆ +τ−1 τ−1−rq²

(21) which are given by:

σ1,2(q) =

τ_q±q

τ_q²−4∆_q

2 (22)

whereτ_q and ∆_q are respectively the trace and determinant ofN_q:

τq=τ−(r+ 1)q² (23)

∆q = ∆−(r+τ−1)q²+rq⁴ (24) There are two kinds of bifurcations that can change the stability under perturbations with wavenumber q.

On one hand we have the case where both eigenvaluesσ1,2(q) are real and one of them changes its sign. This bifurcation is a prolongation of the Turing instability and it happens when ∆q changes its sign (∆q= 0).

On the other hand we have the case where the two eigenvalues σ_1,2(q) are complex conjugate and they change the sign of their real parts. This bifurcation is a prolongation of the Hopf bifurcation studied on the temporal system and it happens whenτ_q changes its sign (τ_q = 0) while ∆_q >0.

Turing instability

The Turing instability takes place when, for aqm6= 0 minimum of ∆q, ∆q_m= 0.

Taking the derivative of Eq. (24) and setting to zero we obtain the position of qm:

qm=±

rr+τ−1

2r (25)

Which exist only whenr >1−τ. It can be seen in Fig. 6 that the homoge- noeus solution will only present a Turing instability if ∆>0. Replacing on the

(13)

Figure 6: Determinant of the Jacobian matrix ∆q as function of q. When r < 1−τ the function has a single minimum at qm = 0. In the other hand, whenr >1−τ, ∆q has two minimums atqm=±q

r+τ−1

2r and a maximum at q= 0

equation of ∆q and setting to zero we obtain the critical value of r at which the Turing instability takes place:

rc= 2∆ + 1−τ+ 2p

∆²+ ∆(1−τ) (26)

which always fulfil the conditionrc >1−τ. As shown in Fig. 7 whenr > rc

there is a region around qm where the growth rates have different signs. For r < rc on the other hand, growth rates aroundqm have all the same sign.

As the system is globally stable, the modes destabilised by the Turing instability are saturated by higher order terms, giving stable (at least in the mode direction) steady periodic states. Depending of the nonlinearities Sub and Sup- critical Turing scenarios can be observed.

As it can be seen in Fig. 8 left, in the Supercritical Turing scenario a stable steady periodic solution is created atr_cand the transition from the homogeneous solution to the pattern state is continous and smooth.

On the other hand, at the Subcritical Turing scenario (Fig. 8, right), an unstable state created atrcfolds into the stable state for ar < rc. This scenario ,unlike the supercritical, presents hysteresis and discontinuous jumps of states.

In the studied system both scenarios are present. In particular between the Saddle-Node bifurcations, near the cusp point, the Turing becomes subcritical.

This interaction between the Turing instability and the Cusp bifurcation has already been studied in optical system and is known as ”nascent bistability”

[19]. This change of criticality makes the study of the interaction of the Turing with the rest of bifurcations more difficult as we will see in more detail below.

Revised Hopf bifurcation

The function τq given by Eq. (23) is a parabola with its maximum atq = 0, with value τq(q= 0) =τ. This means that, when τ = 0 while ∆>0 the two eigenvalues for theq= 0 mode will change their stability as complex conjugates.

(14)

Figure 7: Growth rate (left) and ∆_q as function of q before, on and after the Turing instability. The curves correspond to the homogeneous solution at µ₁ = 0,µ₂ =−7 and ν =−4 forr > r_c (r = 77.33, dashed curve), r=r_c (r= 37.33, solid curve) andr < r_c (r= 17.33, dotted curve).

Figure 8: Super-critical (left) and sub-critical (right) scenarios for the Turing instability. The solid (dashed) curves represent the stable (unstable) steady solutions.

(15)

Figure 9: Growth rate (left) andτq as function of q before on and after the Hopf bifurcation. The curves correspond to the homogeneous solution for µ2 = 0, b= 2.4,ν=−1,r= 4.2 andµ1= 1.123 (caseτ >0, dashed curve),µ1= 0.423 (case τ = 0, continuous curve) andµ1 = 0.223 (case τ <0, dotted curve). ∆q

is positive for all q in all the cases.

This is the equivalent to the Hopf bifurcation studied in section 2.1.2, and it has the same expression.

Above the Hopf bifurcation, whenτ >0 both eigenvalues of the modes with

|q|<q _τ

r+1 have positive real part, if ∆_q >0 in this region. This implies that there is a a group of unstable modes aroundq= 0 (Fig. 9).

The linear temporal frequencies (natural frequency) given byω =Im[σ(q)]

of the unstable modes depend on q. A typical frequency profile can be seen in Fig. 10.

3.4 Other spatial bifurcations and high codimension points

The interaction between the previous bifurcations generates a big number of codimension 2 scenarios. Some of them can be crucial to study the different scenarios and solutions present in the system. In this thesis we will study three of them.

3.4.1 Quadruple-zero codimension 2 point.

This codimension 2 point is the place in the parameter space where the qc of the Turing instability tends towards zero and, therefore, the Turing instability collide with the Saddle Node bifurcation. This point is especially relevant in the study of the steady solutions of the system.

To obtain steady solutions of the system is useful to look at the spatial dynamics of the system. The stationary solutions are the solutions of the four order dynamic system in x obtained when impossing ^∂u_∂t =^∂w_∂t = 0 on Eq. (19).

At the quadruple zero (QZ) the fixed point of the spatial system associated with the steady Homogeneous solution will have a high degenerate Jacobian with 4 zero eigenvalues [20] [21].

The quadruple zero point is located where ∆ = 0 at the same time that r =rc. Using both conditions we reach to the relation r = 1−τ. Using the expression for τ given by (8) the value of u at which the quadruple-zero takes

(16)

Figure 10: Typical natural frequency profile above the Hopf bifurcation. The figure shows the natural frecuency as function of q of the homogeneous solution forb= 2.4,ν=−1,r= 4.2,µ₁= 1.123 and µ₂= 0. This profile correspond to the dashed growth rate shown in Fig. 9.

place (uQZ) is obtained. Replacing it on the the expression of ∆, Eq. (9), and the condition of the fixed points, Eq. (3), the location of the position of the quadruple zero point in the{µ1, µ₂}parameter plane is obtained:

µ_1QZ=−2u³_QZ µ_2QZ= 3u²_QZ

u_QZ= b±p

b²+ 4(ν+r−1)

2 (27)

The two solutions are the starting and ending points of the Turing instability curve in the parameter plane {µ₁, µ₂} (Fig. 11).

(17)

Figure 11: Parameter space for b= 2.4,ν =−1 andr= 10. The figure shows the two QZ points (purple and pink diamonds) where the Turing instability curve (black dash-doted curve) starts and ends.The two Hopf bifurcation lines (red and green lines) with their respective TB points (red and green dots); and the SN bifurcation curves (blue curves) are also shown.

Looking at these expressions some conclusions about the shape of the Turing instability can be inferred. Whenν+r−1>0 each of the two QZ points will be in a different SN branch of the Cusp (Fig. 11). This mean that the Turing instability will trace a loop around the Cusp and the Turing instability will affect the two ∆>0 fixed points (Fig. 14a). Whenν+r−1 = 0 one of the QZ points will coincide with the Cusp in a codimension 3 point destroying the loop. After this, for ^−b₄² < ν+r−1<0 both QZ points are in the same SN branch. For ν+r−1<^−b₄² there is not a Turing instability for any homogeneous solution.

3.4.2 Generalised Takens-Bogdanov.

The double-zero is a codimension 2 point at which the linearised equations around a homogeneous solution of the system has two zero eigenvalues for some finite wavenumber (Fig. 12). According to Eq. (22) this occurs whenτq = 0 and

∆q = 0 for the same q=qT B. Applying this conditions to the equations ofτq

and ∆q, Eq. (23) and (24), the expressions for qT B and the condition for the bifurcation are:

q²_{T B}=1−r±p

(r−1)²+ 4∆

2 τ= (r+ 1)q²_{T B} (28)

The obtained relations are really complex and obtaining an analytic equation for the bifurcation curve is difficult. Despite this some conclusion can be obtained from (28). Asq_{T B} must be real and r is defined as positive, the parameter space can be divide in two different regions where the TB bifurcations curves have substantial differences.

• r ≥ 1. In this region the minus solution of qT B is meaningless, giving not a real value for any ∆. On the other hand, the plus solution ofqT B

(18)

is real only when ∆>0 and, therefore, only exist for the upper branch homogeneous solutions.

The generalized TB curve start and ends at the homogeneous Takens- Bogdanov codimension 2 points.

• 0< r <1. Both expressions forqT B, plus and minus, can give solutions with physical meaning in this region.

For the plus case we obtain two different curves, one with ∆>0, which affects the stable homogeneous solutions, and another one with−^(r−1)₂ ² <

∆<0, which affect the saddle homogeneous solution. Both curves start and end at two codimension 2 points with ∆ = 0 andτ= 1−r².

For the minus case theqT B is real when−^(r−1)₂ ² <∆<0, therefore this bifurcation branch affects only to the saddle homogeneous solution. This solution starts and ends at the TB codimension 2 points.

Figure 12: Typical growth rate at the generalised Takens-Bogdanov bifurcation.

The figure shows the growth rate of the homogeneous solution as function of q forb= 4.54,ν =−1,r= 10,µ1= 1 andµ2= 0.

3.4.3 Turing-Hopf codimension 2 point

The unstable modes generated from a Turing instability, and the unstable modes from the Hopf bifurcation interact generating spatio-temporal periodic solutions.

Close to the codimension 2 point where both bifurcation curves (Turing and Hopf) cross the behaviour of the system is very rich. The Turing-Hopf (T-H) scenario has been studied in the literature when both, the Turing and the Hopf instabilities are supercritical [7][13][18]. At this codimension 2 point there are two complex conjugate eigenvalues with null real part for q=0 and a single eigenvalue is zero for a finite wavelengthq_c (Fig. 13).

(19)

(a) (b)

Figure 14: Parameter space for b = 2.4, ν = −1 andr = 10. The figure (a) shows a zoom in of Fig. 11 around the Cusp point, where can be appreciated the Saddle Node (blue), Hopf (red) and Homoclinic (purple) bifurcation curves.

Also the Turing instability curve (black dotted-dashed curve) is shown, which has a loop around the Cusp point and cut with the Hopf bifurcation in the Turing-Hopf codimension 2 point (black dot). Figure (b) presents the vicinity of the TH codimension 2 point, with the Hopf, Turing, T1 (purple dashed) and T2 (blue dashed) bifurcation lines and the different regions in which the parameter space is divided.

Figure 13: Typical growth rate at the Turing-Hopf bifurcation. The figure shows the growth rate of the homogeneous solution as function of q forb= 2.4, ν =−1,r= 10,µ1= 0.778 and µ2=−1.161.

The normal form of the T-H codimension 2 points, for supercritical instabilities, shows two extra bifurcation curves (T₁andT₂) arising from the T-H point dividing the neighbourhood of the T-H point in 6 different regions, as shown in Fig. 14b.

• Region I. In this region the homogeneos solution is stable for any q.

(20)

• Region II. Crossing the Turing instability to region II the homogeneous solution became unstable for a mode with q = q_c and of stable spatial periodic solution is created (Fig. 15).

• Region III. Crossing the Hopf bifurcation the homogeneous solution loses the stability for q= 0 modes and a time-periodic homogeneous solution is created. However, this solution is unstable against perturbation with a finite wavenumber and a steady spatial-periodic solution is the final stable solution (Fig. 16).

• Region IV. Crossing the T2 bifurcation to region IV the time-periodic solution become stable in any direction and a hybrid solution, inhomogeneous and oscillatory, is generated from the oscillatory solution. This hybrid solution is unstable and can only be observed as a transient state that will decay to spatial patterns or oscillatory solution. In this region the system presents bi-stability of spatial patterns and oscillatory homogeneous solution (Fig. 17).

• Region V. At T1 the hybrid state is annihilated and the pattern solution loses its stability under homogeneous perturbations. In region V the oscillatory homogeneous solution is the only stable solution of the system (Fig. 18).

• Region VI. Crossing the Turing instability the steady pattern solution is annihilated and the steady homogeneous solution recover the stability against perturbation with wavenumberqc. Despite this it is still unstable under homogeneous perturbations and the oscillatory solution is the only stable solution (Fig. 19).

• Finally crossing the Hopf bifurcation and returning to region I the oscillatory solution is destroyed and the homogeneous steady solution recover completely its stability.

The scenario in which the positions of T₁ and T₂ are interchanged is also contemplated in the literature. In this scenario, in the region between T1 and T2the hybrid solution is stable while both, the oscillatory homogeneous and the spatial pattern solutions, are unstable. However, this scenario has not be found in our system. To clarify the possible existence of this scenario and get a deeper knowledge of the T-H point in this system it would be necessary to perform a multi-scale analysis of the system as in [8].

(21)

(a) (b)

Figure 15: Non-homogeneous solutions in region II (Fig. 14b). In this region there is only a spatial pattern solution, because the homogeneous steady solution is unstable only around qc. A space-time plot of the evolution of the u and w field are shown in (a) and (b) for this solution.

(22)

(a) (b)

(c) (d)

Figure 16: Non-steady-homogeneous solution and pattern solution in region III (Fig. 14b). In this region the homogeneous oscillatory solution is observed only as a transient, as it is unstable. The asymptotic stable state is a steady periodic pattern. A space-time plot of the evolution of the u and w field are shown at (a) and (b) for the oscillatory solution; and at (c) and (d) for the pattern solution.

(23)

(a) (b)

(c) (d)

(e) (f)

Figure 17: Non-steady-homogeneous solutions in region IV (Fig.14b). In this region there are two stable solutions, the spatial pattern solution and the homogeneous oscillatory solution. There is also an unstable hybrid solution. This solution can only be reach as a transient state. A space-time plot of the evolution of the u and w field for this patterns solution are shown in (a) and (b); for the hybrid solution in (c) and (d); and for the oscillatory solution in (e) and (f).

(24)

(a) (b)

(c) (d)

Figure 18: Non-steady-homogeneous solutions in region V (Fig.14b). In this region there is one stable solution only, the homogeneous oscillatory state. There is also an unstable steady spatial pattern solution, this solution can only be reached as a transient state. A space-time plot of the evolution of the u and w field are shown at (a) and (b) oscillatory solution; and at (c) and (d) for the for the pattern solution.

(25)

(a) (b)

Figure 19: Oscillatory-homogeneous solutions at region VI (Fig.14b). At this region there is only a homogeneous oscillatory, because the homogeneous steady solution is unstable only around q = 0. The u and w field for this oscillatory solution are shown in (a) and (b).

3.4.4 Turing-Takens-Bogdanov codimension 3 point

By doingr= 1 the QZ and the TB points collide in the Turing-Takens-Bogdanov (TTB) codimension 3 point. On this point the Turing, Hopf, Homoclinic, Saddle node and Generalised Takens-Bogdanov bifurcation curves start. At this point the linear analysis of the associated homogeneous solution have a quadruple zero in the four dimensional spatial space and a double zero in the temporal system.

Applying the conditionr= 1 to (27) the following expression for the TTB points is reached:

µ1T T B=−2u³_H µ2T T B = 3u²_H r= 1 (29)

whereuH is given by (10).

All the interaction dynamics between the Turing and Homoclinic bifurcation are expected to develop around this high codimension point. In this thesis we will only be able to do a superficial study of the full scenario.

3.5 Turing-Hopf scenario close the Homoclinic bifurcation

In this section we will present an exploration of the parameter space between the SN branches for the parameters b=2.4,ν =−1 and r= 4.2 (Fig. 20). We will focus the study on the attractors with a basin that includes part of the neighbourhood of the upper branch homogeneous solution emerging from the µ₁<0 SN branch.

In this region the Turing bifurcation becomes sub-critical [19]. Crossing the Turing instability the system jumps to a large amplitude stationary pattern state (Fig. 21). This state is very robust and once formed it survives much before the Turing instability. This scenario does not allow observing hybrid states or interactions with oscillatory modes in the vicinity of the Turing-Hopf codimension 2 point because, after the Turing instability, all initial conditions

(26)

Figure 20: Parameter space forb= 2.4,ν =−1 andr= 4.2 around the Homo- clinic bifurcation. In the figure are shown the SN, Cusp, Hopf, Homoclinic and Turing bifurcation curves as well as the points where computational simulation have been performed.

around the unstable homogeneous solution tend to the stationary pattern solution. Anyway, some interaction between the unstable pattern solution created on the Turing and the oscillatory modes near the T-H is expected, even if it has not been observed on this thesis because such interaction leads probably to unstable states.

Crossing the Hopf bifurcation stable homogeneous oscillatory solutions are created. Separating from the Hopf bifurcation the oscillatory solution start to develop spatial patterns with period 2 oscillations (Fig. 22). Further from the Hopf bifurcation there are windows where the solution loses the temporal periodicity (Fig. 20). It has not been studied if this state shows a pseudo- periodicity or a long timescale periodicity (Fig. 23). The transition from time periodic solutions to nonperiodic ones needs further investigation for a complete understanding of the dynamics. The non-periodic solutions are followed by a turbulent regime, where the solutions do not have any spatial or temporal periodicity, (Fig. 24). When the system is close to the Turing instability, this turbulent solution is only reached as a transient state and, eventually, it falls to the stationary spatial solution, which is the only attractor close to the upper branch homogeneous solution. Further from the Turing instability, this turbulent solutions survives for long integration times (at least 3·10⁴time units) even once the Homoclinic bifurcation is crossed.

(27)

Figure 21: Fully developed spatial periodic solution reached slightly above the Turing instability. In the figure can be seen the pattern in the spatial space (up) and its Fourier modes (down) for the u (left) and w (right) fields. The blue, red and green lines represent the homogeneous solutions present in the system for this parameters.

(a) (b)

Figure 22: Oscillatory pattern solution. Notice that, unlike the hybrid solution near the Turing-Hopf codimension 2 point, this oscillatory solution has period 2. The figure shows the spatio-temporal evolution of the u field (a) and the w field (b).

(28)

(a) (b)

Figure 23: Non periodic pattern solution. The figure shows the spatio-temporal evolution of the u field (a) and the w field (b).

(a) (b)

Figure 24: Turbulent solution. The figure shows the spatio-temporal evolution of the u field (a) and the w field (b).

4 Conclusions

In this thesis we have performed a bifurcation analysis of the temporal system given by (1). In this manner we have found the principal bifurcations and codimension 2 points. These results are consistent with [15].

By adding diffusion to the temporal equations, we have created a system with a Turing instability. The Turing instability curve has been analytically characterised in relation to the position of the QZ codimension 2 points and numerically evaluated using a Newton-Raphson method. The obtained Turing instability shows a change from super-critical to sub-critical when approaching the Cusp bifurcation. This behaviour supports the results from [19].

The interaction between the Turing and the Hopf bifurcation has been stud-

(29)

ied around the T-H codimension 2 points, obtaining results consistent with the behaviour described in the literature ([7][13][18]) in the unstable hybrid state scenario for the case withb= 2.4,ν =−1 andr= 10.

We have also found two unexplored high codimension bifurcations in the spatial system: the generalised Takens-Bogdanov bifurcation and the Turing- Takens-Bogdanov codimension 3 point. In the former, the Jacobian matrix presents a double zero for a non-homogeneous mode. This bifurcation is hy- pothesized to be the starting point of a non-homogeneous Homoclinic bifurcation which destroys the oscillations generated in the Hopf bifurcation. In the latter (the TTB point), the Turing, Hopf, Homoclinic, Saddle node and Gen- eralised Takens-Bogdanov bifurcation curves have been found to collide. The different interaction scenarios between Turing, Hopf and Homoclinic bifurcation are expected to develop around this point.

Finally, in this thesis we have studied the parameter space with numerical simulations. As a result we have found spatio-temporally periodic solutions and turbulent regimes.

Due to the complexity of the studied system, there are many open questions for future research.

First of all, the system has been built imposing a basis for the diagonal diffusion, given by the condition (18). The implications of a specific choice of basis should be studied to reach a universal description of the Turing-Homoclinic interaction.

The change from super-critical to sub-critical in the Turing bifurcation has not been deeply studied on this thesis. To understand the exact reason for this change and its relation with the proximity to the Cusp it would be necessary to do a multiscale analysis at this bifurcation. On the other hand, we have not found the Turing-Hopf scenario with stable hybrid states documented in [18], neither have we studied the T-H neighbourhood in the case a sub-critical Turing bifurcation. These issues can be addressed by performing a two parameter multiscale analysis of the system around the T-H point.

To conclude, in further research it would be interesting to make a deeper study of the generalized Takens-Bogdanov bifurcation and around the Turing- Takens-Bogdanov point, as well as a more detailed analysis of the transition to a turbulent regime presented in section 3.5.

(30)

Appendix: Computational methods

Pseudospectral integration method

The partial differential equation which describe the system, Eq. (19), has been integrated with periodic boundary conditions using a pseudo-spectral integration method adapted from [22]. In this section the integration method will be described.

Any reaction diffusion equation can be written as:

∂_tu_i=L_iu_i+N_i(u(x, t)) (30) Where L are the diagonal linear terms and N(u) the non-linear plus the non-diagonal linear terms. This equation can be writing in Fourier space as:

∂tu˜i(q, t) =−αi(q)˜ui(q, t) + Φi(˜u(q, t)) (31) Where ˜ui(q, t) =F[ui(x, t)],−αi(q) is the Fourier transform of the operator Liand Φi(q,t) is the Fourier transform ofNi(u(x, t)).

The pseudospectral method consist of two consecutive steps after which the field u(x, t0+ 2δt) is calculated with errorO(δt³). The steps are given by the expressions:

˜

u_i(q, t+δt) =e^−αⁱ^(q)δtu(q, t) +˜ 1−e^−αⁱ^(q)δt

α_i(q) Φ_i(q, t) +O(δt²) (32)

˜

ui(q, t+ 2δt) =e^−2αⁱ^(q)δtu(q, t) +˜ 1−e^−2αⁱ^(q)δt

αi(q) Φi(q, t+δt) +O(δt³) (33) At each step the nonlinear term is compute as Φ_i(q, t) =F[N(F⁻¹[˜u_i(q⁰, t)])].

The Fourier transformations are done using the Fast Fourier Transform from In- tel math kernel libraries.

In the system under study the methods functions and parameters are the following: u1 = u, u2 = w, α1 = −1 +q², α2 =−ν+ 1 +rq², Φ1 =w and Φ2=−u³+ (µ2+ν−1)u+µ1+ (w+u)(bu−u²).

Newton Raphson

To obtain a numerically the Turing instability curve in Fig. 14 we have use a multivariable Newton-Raphson generalization. In this section I will explain the method used and its implementation in computing the Turing instability curve.

The using method attempts to find a root of a given function f(µ) with µ= (µ1, µ2), creating a sequence ofµ(n) starting from an initial guess µ(0).

The method is expected to converge to a given valueµ^∗ that fulfilsf(µ^∗) = 0.

Each iteration is obtain using the following equations.

µi(n) =µi(n−1)−γ f(µ(n−1))

∂µ_if(µ(n−1)) (34)

Whereγ∈(0,1) is a parameter of the method which control the size of the steps, and it is used to optimize the convergence of the system.

(31)

To find the Turing instability we look for the roots off(µ) =r_c(u^∗(µ),µ)−r where r_c is given by Eq. (26). The derivative of f(µ) used in the method is given by:

∂µ_if = 2∂µ_i∆−∂µ_iτ+(2∆ + 1−τ)∂µ_i∆−∆∂µ_iτ

p∆²+ ∆(1−τ) (35)

where:

∂_µ₁∆ = 6u^∗∂_µ₁u^∗ ∂_µ₁τ= (b−2u^∗)∂_µ₁u^∗ ∂_µ₁u^∗= ∆⁻¹ (36)

∂_µ₂∆ = 6u^∗∂_µ₂u^∗−1 ∂_µ₂τ= (b−2u^∗)∂_µ₂u^∗ ∂_µ₂u^∗= u^∗

∆ (37)

which are obtained taking the partial derivative of Eq. (3), (8) and (9).

(32)

References

[1] A. M. Turing, “The chemical basis of morphogenesis,”Philosophical Trans- actions of the Royal Society of London. Series B, Biological Sciences, vol. 237, no. 641, pp. 37–72, 1952.

[2] V. Castets, E. Dulos, J. Boissonade, and P. De Kepper, “Experimental evidence of a sustained standing turing-type nonequilibrium chemical pattern,”Phys. Rev. Lett., vol. 64, pp. 2953–2956, Jun 1990.

[3] A. Nakamasu, G. Takahashi, A. Kanbe, and S. Kondo, “Interactions between zebrafish pigment cells responsible for the generation of turing patterns,”Proceedings of the National Academy of Sciences, vol. 106, no. 21, pp. 8429–8434, 2009.

[4] M. Rietkerk and J. van de Koppel, “Regular pattern formation in real ecosystems,” Trends in Ecology and Evolution, vol. 23, no. 3, pp. 169 – 175, 2008.

[5] Z. Tan, S. Chen, X. Peng, L. Zhang, and C. Gao, “Polyamide membranes with nanoscale turing structures for water purification,”Science, vol. 360, no. 6388, pp. 518–521, 2018.

[6] D. Walgraef, Spatio-Temporal Pattern Formation With Examples from Physics, Chemistry, and Materials Science. Springer Science, 1997.

[7] A. Rovinsky and M. Menzinger, “Interaction of turing and hopf bifurcations in chemical systems,”Phys. Rev. A, vol. 46, pp. 6315–6322, Nov 1992.

[8] Y. Song, T. Zhang, and Y. Peng, “Turing–hopf bifurcation in the reaction–diffusion equations and its applications,”Communications in Nonlin- ear Science and Numerical Simulation, vol. 33, pp. 229 – 258, 2016.

[9] D. G. M´ıguez, S. Alonso, A. P. Mu˜nuzuri, and F. Sagu´es, “Experimen- tal evidence of localized oscillations in the photosensitive chlorine dioxide- iodine-malonic acid reaction,” Phys. Rev. Lett., vol. 97, p. 178301, Oct 2006.

[10] J. Guckenheimer, K. Hoffman, and W. Weckesser, “The forced van der pol equation i: The slow flow and its bifurcations,”SIAM Journal on Applied Dynamical Systems, vol. 2, no. 1, pp. 1–35, 2003.

[11] C. Sparrow, The Lorenz equations : bifurcations, chaos, and strange attractors. Springer-Verlag New York, 1982.

[12] E. Izhikevich, “Neural excitability, spiking and bursting,” International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, vol. 10, pp. 1171–1266, 06 2000.

[13] M. Baurmann, T. Gross, and U. Feudel, “Instabilities in spatially extended predator–prey systems: Spatio-temporal patterns in the neighborhood of turing–hopf bifurcations,” Journal of Theoretical Biology, vol. 245, no. 2, pp. 220 – 229, 2007.

(33)

[14] D. Ruiz-Reyn´es,Dinamics of Posidonia Oceanica Meadows. Unpublished doctoral thesis, Universitat de les Illes Balears and IFSC, 2019.

[15] F. Dumortier, R. Roussarie, J. Sotomayo, and H. ˙Zo ladek, Bifurca- tions of Planar Vector Fields Nilpotent Singularities and Abelian Integrals.

Springer-Verlag, 1991.

[16] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Sys- tems, and Bifurcations of Vector Fields. Springer New York, 2002.

[17] S. H. Strogatz, Nonlinear Dynamics and Chaos. Perseus Books, 1994.

[18] M. Cross and H. Greenside, Pattern formation and dynamics in nonequilibrium systems. Cambridge University Press, 2009.

[19] M. Tlidi, M. Georgiou, and P. Mandel, “Transverse patterns in nascent optical bistability,”Phys. Rev. A, vol. 48, pp. 4605–4609, Dec 1993.

[20] P. Parra-Rivas, D. Gomila, L. Gelens, and E. Knobloch, “Bifurcation struc- ture of localized states in the lugiato-lefever equation with anomalous dis- persion,”Phys. Rev. E, vol. 97, p. 042204, Apr 2018.

[21] P. Coullet, “Localized Patterns and Fronts in Nonequilibrium Systems,”

International Journal of Bifurcation and Chaos, vol. 12, pp. 2445–2457, Nov. 2002.

[22] R. Montagne, E. Hern´andez-Garc´ıa, A. Amengual, and M. San Miguel,

“Wound-up phase turbulence in the complex ginzburg-landau equation,”

Phys. Rev. E, vol. 56, pp. 151–167, Jul 1997.

Spatiotemporal patterns in the Turing-Takens-Bogdanov scenario