Pattern analysis in a benthic bacteria-nutrient system

We study pattern formation in a reaction-diffusion system for a benthic bacteria-nutrient model in a marine sediment, which originally contains some spatially varying coefficients and with these shows some layering of patterns. Using the Landau reduction for the system with homogeneous coefficients, we locally analyze Turing bifurcations in 1D and 2D. Moreover, we use the software pde2path to compute the corresponding branches globally and find a number of snaking branches of patterns over patterns. This shows that spatially varying patterns are not necessarily due to spatially varying coefficients. We find a type of hexagon patches on a homogeneous background with no prior mention in the literature. We show the first numerically calculated solution-branch, which connects two different types of hexagons in parameter space. We call states on this branch rectangles. We check numerically, if the stability changes for hexagons and stripes, which are continued homogeneously into the third dimension. We find that stripes and one type of hexagons have the same stable range over bounded 2D and 3D domains, while the other type of hexagons becomes unstable earlier. Here we find a snaking branch of solutions, which are spatial connects between hexagonal prisms and a genuine 3D pattern (balls).


Introduction
A reaction-diffusion system for a benthic bacteria-nutrient model in a marine sediment is set up and analyzed in [2]. In dimensionless form this system is described by Here u = u(t,x,ỹ) denotes the population density of bacteria and v = v(t,x,ỹ) the concentration of a nutrient, wherex andỹ are the horizontal and vertical spatial coordinates in the sediment, respectively. t is the time. δ u and δ v are diffusion constants of bacteria and nutrient, respectively. γ describes the strength of active stimulation in a Holling II nonlinearity, k the associated half saturation constant, m the mortality of bacteria, and ε the rate of bacteria inflow, while v 0 is the nutrient concentration in the sea water above the sediment, i.e., atỹ = 0. The functionσ describes the bioirrigation.
Rescaling (1) with x =x √ δu and y =ỹ √ δu , yields where δ = δv δu = 50 and σ(y) :=σ(ỹ). In order to simplify the notation we set w = (u, v), D = 1 0 0 δ , and The function f is called the reaction term or kinetic of the full system (2). With this simplifications we can write (2) as To understand bioirrigation, we consider the system for a homogeneous σ. The unique solution of (5) is given by which converges to v 0 for t → ∞. Thus the bioirrigation σ determines the speed at which the nutrient concentration in the sediment adapts to v 0 . Here and in [2] the parameter set k = 1, v 0 = 4.125, ε = 0.005, m = 0.3175, δ u = 2 · 10 −5 , δ v = 10 −3 (6) is used and spot patterns are classified in cold and hot, which means that they have a minimum and maximum in the center of every spot, respectively. Because of the predator prey structure of (1) we have a hot-spot pattern for u, when we have a cold-spot pattern for v and vice versa. Thus we always present only the pattern of u, and when we are saying that a solution of (1) is hot or cold, this means that this is the case for u.
Quasi-stationary layering of patterns involving stripes and hot spots are observed for µ = 10, γ = 0.25, and α = 0.1 (see Fig.11 of [2]). By changing the diffusion constant δ x to 7 · 10 −6 and using µ = 3, γ = 0.25, and σ = 0.2, a quasi-stationary solution is found, which involves three different types of patterns. These are homogeneous, cold-spot, and stripe patterns (see Fig.12 of [2]). Here we describe the bioirrigation by σ(y) = 0.128 1 + e 0.011(y−480) (7) and obtain a quasi-stationary solution involving all five pattern types shown in Fig.10 of [2] (see Fig.1). Using the Landau analysis described below, we predict a transition of stability between cold spots and stripes and between stripes and hot spots in the ranges y ∈ [265, 296] and y ∈ [380, 418], respectively. One can see in Fig.1 that this is a sensible prediction. The goal of this paper is to continue the investigations of [2] by understanding the bifurcation scenarios of (1) for homogeneous bioirrigations and to find and investigate more stationary patterns via bifurcation analysis. Mostly we see γ as a given and fixed parameter, while we use σ as a bifurcation parameter, i.e., we examine how solutions and their types of stability change and what kind of new solutions bifurcate by varying σ.
It is shown in [2] that the problem of determining the homogeneous solutions can be reduced to computing the zeros of a cubic polynomial. The coefficients are given with some typos. The first γ of a 1 must be σ. The zeros of this polynomial are not derived analytically in [2] because of their complex terms. In Section 2 we show the reduction to the cubic polynomial again and sort the coefficients with respect to σ and γ, because we investigate (1) in the following for different values of σ and γ, while we use the set (6) for all other parameters of (1). Furthermore, we apply analytical formulas for the three zeros of the cubic polynomial to determine the regions of a bounded domain of the σ-γ-plain, where the homogeneous solutions are real and to see for : Density plots of u at times t = 1200 and 6200, which we found by using time-integration methods for a small random perturbation of (u, v) ≡ (1, 1). Here we use (2) for γ = 0.25 and δ = 50. The bioirrigation is described by (7).
which combinations of σ and γ exist only one or rather three homogeneous states. It follows an introduction of Turing-instabilities and bifurcations. By using a fine discretization of the σ-γ-domain, we determine Turing unstable ranges. With this informations we are able to predict the bifurcation locations and wavenumbers of periodic 1D and 2D solutions.
In Section 3 we recall how to reduce a general two species reaction-diffusion system to the Landau amplitude equation system on a hexagonal lattice. The Landau reduction can be used to approximate solutions and develop bifurcation diagrams, but this only works close to the bifurcation from the homogeneous solution. Hence, we use the continuation and bifurcation software pde2path [23] to get a global picture.
In Section 4 we study non trivial patterns over 1D domains and show global bifurcation diagrams for γ = 0.3, which we developed by using pde2path. We find that some stripe solutions bifurcate from the homogeneous branch and terminate in a different bifurcation on the homogeneous branch by holding their wavelength. Moreover, we find stripe branches of the same wavelength which are not connected, but change their wavelength and connect to other stripe branches. Using the Landau reduction, we choose a γ-value where stripes bifurcate subcritically, to develop bifurcation diagrams with snaking branches of localized stripes on homogeneous backgrounds.
In Section 5 we investigate 2D patterns and develop some global bifurcation diagrams by using pde2path. Furthermore, we show solutions and solution branches of localized hexagonal spots with a planar interface to striped backgrounds and patches of localized hexagons on homogeneous backgrounds. The patches which we find first have a flat overlying front such that the snake does not show a strong snaking behavior. Hence, we use the Landau reduction to find γ-values where the front is steeper and thus the snaking behavior is more pronounced. We end in Section 6 with a discussion and an outlook.
From this we see that a layering of patterns is not necessarily an effect of space dependent parameters but can occur for homogeneous bioirrigations. A bifurcation analysis for depth dependent functions will be presented elsewhere.

Homogeneous solutions and Turing instabilities
Homogeneous steady states are solutions of f (u, v) = 0 which are space and time independent. For a nonzero rate of bacteria inflow ε and nonzero product of the bioirrigation σ and the nutrient concentration v 0 in the sea water it holds for homogeneous solutions of (2). Adding the first to the second equation of (3), yields the linear relationship From (8), (9), and the fact that the population density of the bacteria u and the concentration of the nutrient v cannot be smaller than zero follows Substituting (9) into f (u, v) = 0, reduces the problem to The zeros of (10) are given by The polynomial (10) can have one or three real zeros for fixed σ and γ. The regions of the σ-γ-plain, where a homogeneous solution is real, can be calculated analytically. Determining analytically, where the population density of the bacteria and the nutrient concentration of such a solution is additionally positive, seems not so trivial. In Fig.2(a) we show the regions on a bounded domain of the σ-γ-plain, where homogeneous solutions are real. They are also positive for the regions shown in Fig.2(b). Comparing (a) and (b), we see that the positivity is not fulfilled for all homogeneous solutions (u, v) of (10). Here (u 2 , u 3 ) and (v 1 , v 2 ) are negative in the upper horizontal gray band and the lower left gray region of Fig.2a), respectively.  We verified that the conditions hold for all discretization points used to generate Fig.2, if all three homogeneous states are real. Thus (u 1 , v 1 ) is the state with the highest population density of bacteria and lowest concentration of nutrient of these three homogeneous states, while this is the opposite case for (u 3 , v 3 ). The equilibrium (u 2 , v 2 ) lies in the middle.
The Jacobian of f is given by The linearization of (4) in a homogeneous state w * is given by the linear operator L(∆) = J f + D∆. It holds L(∆)e i(x,y)·k =L(|k|)e i(x,y)·k withL(|k|) = J f − D |k| 2 and k ∈ R 2 .
ii) w * is linear unstable in the full system (4), i.e., it exists a k such that It can be shown (see [20]) that a homogeneous state of (4) is • We call a solution a Turing endpoint and Turing bifurcation point, if it lies at a transition from Turing-unstable to unstable and from Turing-unstable to stable, respectively. Turing bifurcation points occur, if detL(|k|) vanishes. We call the corresponding σ and |k| critical bioirrigation σ c and critical wavenumber k c , respectively. In Fig.3 we illustrate the stabilities of the three different homogeneous states in the σ-γdomain of Fig.2. First of all we can see in Fig.3(c) that (u 3 , v 3 ) is always stable, if it is real in our chosen region. By using the strength of bioirrigation σ as bifurcation parameter and the activity stimulation γ as a fixed parameter, we can classify the bifurcation scenarios into five different types. We do this by partitioning the γ-interval into the following five sections  The least interesting interval is I 5 . Here the state (u 1 , v 1 ) is always real, positive, and stable, while the other two are not real or not positive. For all other intervals we have Turing-unstable ranges with Turing bifurcation points at the right boundaries. The different behaviors at the left boundaries and what happens beyond is as follows.
, stops to be real, and a fold occurs. The branch continues on (u 2 , v 2 ) and is unstable. I 2 : (u 1 , v 1 ) becomes unstable (Turing endpoint) before following the behavior of I 1 . For I 3 and I 4 there is always only one real homogeneous state, which is also positive. By increasing σ from 0 to 0.25, the stability of this state changes for the different intervals as follows: Thus for I 4 we have one Turing unstable range, which ends up with Turing bifurcation points, while the Turing unstable range is non-contiguous for I 3 . Here we have two Turing bifurcation points and two Turing end points.
Turing patterns branch from Turing bifurcation points as discussed below. It is unclear for us which kind of solutions bifurcate from Turing endpoints. Here we have Re[µ ± (0)] = 0, while Im[µ ± (0)] = 0 such that we have Hopf bifurcations at the Turing endpoints. Currently pde2path [23] does not handle Hopf bifurcations and thus we postpone this analysis.
Example bifurcation diagrams for I 1 and I 3 can be seen in Fig.3 (a) and (d), respectively.

Landau reduction
First we recall how the Landau reduction works in general for systems of the form (4).
Let w * = (u * , v * ) be a homogeneous solution of f and i, j two non-negative integers. We . We use an analog notation for g and h.
Taylor-expanding f We stop at order three, because we also use the Landau reduction only up to third order. Substituting the expansion above into (2), the system becomes where L(∆) = J f (w * )+ ∆ 0 0 d∆ . B and C are symmetric bilinear and trilinear forms, respectively. For p, q, r ∈ R 2 they have the form: The eigenvalue µ − is negative in Turing-unstable ranges, while the eigenvalue µ + has positive parts. Thus we always consider µ + in the following and write µ for simplicity. The spot patterns shown in Fig.1 and [2] have a hexagonal structure, because every spot has six direct neighbors. Thus we start our Landau reduction on hexagonal lattices with the following ansatz to reduce (16) to a system of space-independent amplitudes. Here Φ = φ(k) is the eigenvector ofL(k), e j = e i(x,y)·k j , A j = A j (t) ∈ C for j = 1, 2, 3, Figure 5: Sketch of the vectors k1, k2, k3.
The length (wavenumber k) of the wave vectors k 1 , k 2 , k 3 is chosen such that there is a bifurcation parameter value σ k for which the curve of eigenvalues has a zero in k. It holds To eliminate quadratic terms from the residual Res(w) = −∂ t w + L(∆)w + B(w, w) + C(w, w, w) which do not correspond to modes e 1 , e 2 , or e 3 , we extend the ansatz (17) to Notice that φ ii and φ ij are independent of i and j. Using this notation, one can see that φ ii and φ ij correspond to A 2 i and A i A j , respectively. Substituting (18) into (16) and sorting with respect to e n m , yields To remove terms of order e i from the residual, we extend the ansatz (18) tow = w + 3 i=1 φ 3i e i . Substitutingw into (16) and sorting with respect to e 1 , e 2 , e 3 , yields The summands R 1 , R 2 , R 3 represent all higher order terms, e.g., A 1 |A 2 | 4 is a term of R 1 . By the Fredholm alternative there exists a solution for (19), iff the r.h.s. is an element of ker(L(k) H ) ⊥ . Let Φ * be the adjoint eigenvector ofL(k) to the eigenvalue µ(k) evaluated in σ k , i.e.,L(k) H Φ * = µ(k)Φ * , and let Φ * be normalized such that Φ, Φ * = 1. Multiplying (19) with Φ * and setting R 1 = R 2 = R 3 = 0, yields One can see in [22] that evaluating them all in σ can give better approximations. However, in the following we use the Landau reduction to predict the existence and bifurcation directions of stationary states which branch from homogeneous solutions. To approximate solutions and follow their branches we use pde2path. This method works also and yields the same system (20) and Landau coefficients, if we use phase shifted space coordinates in the ansatz, i.e., using (x + ψ, y + ψ) with ψ ∈ (0, 2π) instead of (x, y).

1D Patterns
First we consider solutions over one dimensional domains. Here the modes e 2 and e 3 do not exist. Hence, the system (20) reduces to Stationary amplitudes solve Clearly A = 0 solves (22). Setting this solution into the ansatz (18), yields the homogeneous solution. More interesting are the second type of solutions which we obtain from (22). They fulfill and generate periodic solutions by substituting (A 1 , A 2 , A 3 ) = (A, 0, 0) into (18). This type of solutions exists also in 2D and we call these solutions stripes because of their 2D-density plot. When we use numerical methods to determine the stripes in the following, we use Neumann boundary conditions. Stripes which fulfill Neumann boundary conditions over a domain (−lπ/k, lπ/k) with l ∈ N correspond to the amplitudes If the bacteria density u has its maximum (minimum) in x = 0 for such a stripe solution, we call it hot (cold) stripes. Notice that we do not have automatically hot and cold stripes for S + and S − , respectively. All other amplitudes which fulfill (23) generate space shifts of the hot resp. cold stripes. Let σ k be a bioirrigation value and k ∈ R + a wavenumber such that the curve of eigenvalues µ for σ k has a single zero in k. Let σ s and σ l be two bioirrigation values, which are infinitesimally smaller and larger than σ k , respectively. It holds that the curve of eigenvalues µ is positive in k for σ s or σ l , while it is negative for the other. The stripes can only exist, where (23) is fulfilled such that the algebraic sign of c 3 evaluated in σ k tells us in which direction the stripes bifurcate.

Changing wavelength
We start a global analysis of 1D patterns in (2) for γ = 0.3. From our analysis above we already know that we still have one homogeneous solution for all σ ∈ (0, 0.25) with an unstable range bounded by Turing bifurcation points. Furthermore, we know that there are two Turing endpoints in this unstable range. We use our data set of Fig.3 to find out that the right and left Turing bifurcation points are given by σ r c ≈ 0.11 and σ l c ≈ 0.025 with corresponding critical wave numbers k r c ≈ 0.187 and k l c ≈ 0.067, respectively. The curve of eigenvalues µ ± (see (12)) has a zero in k, when detL(k) = 0. This is the case for Clearly in the Turing unstable range it holds k ± ∈ R. We proved that this is also the case between the Turing endpoints. c3 and c f = c 2 2 /(4(c3 + 2c4) 2 ) evaluated in σc as function of γ are shown in (b) and (c) and used to predict the strength of the subcriticality of stripes and hexagons, respectively. Hexagon patterns will be introduced in Section 5.
The Landau formalism above predicts that periodic solutions of the type bifurcate from (u * , v * ) at σ r c with k = k r c , if we consider the problem over the 1D domain Ω r = (−4π/k r c , 4π/k r c ). We are able to prove analytically that there are bioirrigations σ 35 , σ 45 , σ 3 with σ r c > σ 35 > σ 45 > σ 3 , where branches R35, R45, and R3 of periodic solutions of the type (24) bifurcate with k = 3.5k r c /4, 4.5k r c /4, and 3k r c /4, respectively. Let s ∈ {σ r c , σ 35 , σ 45 , σ 3 } and κ be the corresponding wavenumber. In Fig.6(a) we see that a bioirrigation value s 2 = s in the unstable range exists for which the eigenvalue curve µ ± has a real zero in κ such that a stripe solution with wavenumber κ also branches in s 2 . Following the branches with pde2path [23] which bifurcate in s, we see that they terminate in s 2 (see Fig.7).
One may conjecture that stripe branches which correspond to the same wavenumber are connected, but this is not always the case. We also computed the branch L4 which bifurcates in σ l c with the critical wavenumber k l c over the domain Ω l = (−4π/k l c , 4π/k l c ). In Fig.6(a) we see that there is σ = σ l c , where stripe solutions of the wavelength k l c bifurcate. We call the corresponding solution branch L4 * . In Fig.8 we see that L4 and L4 * are not connected, but they connect to bifurcations on L8 and L16, which are stripe solutions of 4 and 16 periods, respectively. One might guess that this depends on the Turing endpoints. However, we also compute the branches for γ = 0.4 (not shown). Here we have no Turing endpoints, but the same effects.

Localized patterns and snaking
The Landau coefficient c 3 evaluated in σ c are positive for γ ∈ (0, 0.209) (see Fig.6(b)) and thus stripes bifurcate subcritically for these active stimulations. We calculate the branch s of hot-stripes (u s , v s ) by using pde2path for γ ≈ 0.004 over the domain Ω = (−24π/k c , 24π/k c ) with k c = 0.212. s bifurcates subcritically from the homogeneous solution (u * , v * ) = (u 1 , v 1 ) at σ c ≈ 0.196, as predicted by the Turing and Landau analysis above. A fold occurs at σ ≈ 0.1985 and the stripes become stable such that there is a bistable range between (u * , v * ) and (u s , v s ).
There are 10 bifurcation points on s on the way from its bifurcation to the fold. The first and 10th, second and 9th, third and 8th, 4th and 7th, 5th and 6th are connected pairwise by branches of stationary fronts between stripes and homogeneous states, which we call l1, l2, l3, l4, and l5, respectively (see Fig.9, 10, and 11).  Solutions of l1 are oscillations of the type with monotonously spatially growing amplitudes A and B on Ω. For simplicity we still describe the behavior of u in the following. Clearly the amplitude A s for u s is constant. By following l1 from the first bifurcation point on s to the first fold, which lies near the 171st solution, u transforms to a front between u * and u s . After this l1 "snakes" back and forth by changing its stability. Along the snake the front moves from the right to the left, and the front position shifts by π/(2k c ) between two successive folds. Beyond the last fold A starts to grow at the left boundary and l1 returns to s.
For every solution of l1 also exists a reflection. We call the corresponding branch l1'. Clearly the illustration of l1' is congruent with l1 and both branches together generate a loop from one bifurcation point to the other and back again.
Splitting the domain Ω into Ω l = (−24π/k c , 0) and Ω r = (0, 24π/k c ), we can describe solutions on l2 as two heteroclinics on Ω l and Ω r , where the one on Ω r is a reflection of the one on Ω l . Here the snake is shorter with fewer wiggles, because the lengths of the domains Ω l and Ω r are shorter than the length of Ω such that the fronts move faster to the boundaries of Ω l and Ω r . We can split the domain Ω into 3, 4, and 5 parts and use heteroclinics like above to describe solutions of l3, l4, and l5, respectively. Here the branches do not show any snaking behavior, because the partitions of Ω are too small. The same branches for connections between cold stripes and the homogeneous solutions can be found on the cold-stripe branch. Additional localized stripes of the form exist for periodic boundary condition. The illustrations of the corresponding branches d1, d2, d3, d4, and d5 (not shown) are not congruent with l1, l2, l3, l4, and l5, respectively. The folds of l1 and d1,..., l5 and d5 are connected pairwise by branches, which are called rungs. Using Dirichlet or Neumann boundary conditions, we can only find solutions of the form (26) or (25), respectively. Currently periodic boundary conditions are not implemented in pde2path but will be added. For more details on localized patterns and snaking over bounded dimensional domains see [4,11,17,16]. Seminal results for localized patterns over unbounded domains can be found in [6,7,3]. For a detailed analysis by using the Ginzburg-Landau formalism and beyond all order asymptotics see [10,13].

2D Patterns
The stripe patterns also exist over two dimensional domains. Here we present additionally some genuine 2D patterns, which can be analyzed via the Landau system (20). Setting A 1 = A 2 = A 3 =: A, the system (20) reduces to Stationary amplitudes fulfill Two types of non homogeneous solutions can be found from this equation. These are hexagon and triangle patterns with A ∈ R and A ∈ C \ R, respectively. The hexagon amplitudes are given by The triangle amplitudes are given by a 1 + a 2 i solving (20) one can find mixed mode solutions. We classify them into bean and rectangle patterns, which fulfill |A| > |B| and |A| < |B|, respectively. Fig.12 shows some bifurcation diagrams over a small 2D domain. In Fig.12(a) we see that when the hot stripes lose their stability, a branch of hot beans bifurcates subcritically, which terminates in a bifurcation on the hot hexagon branch in the point where these become stable. An analog scenario occurs in the bistable range between cold hexagons and cold stripes (see Fig.12(b)). Over this small domain the hot-bean branch has 4 bifurcation points (not indicated). By increasing the horizontal direction by a factor of 4, the number of bifurcation points on the hotbean branch also increases by a factor of 4 (see Fig.13) with bifurcations of stationary connections between stripes and hot hexagons. Here everything is analog to the 1D case described above.
The bifurcation points are connected in the same way. From the first one bifurcates a snaking branch of solutions, which we see again as front between hot stripes and hot hexagons. From the second bifurcates a branch of two interfaced heteroclinics, and this holds analogously for the other bifurcating solutions of the hot-bean branch. In contrast to the 1D case the snake does not snake around a vertical line but in a slanted manner. This is a finite size effect; cf. [5,4,16,12]. Such localized patterns can also be found on cold bean branches (see Fig.14(a)). It can be seen in (27) that one hexagon type (hot or cold) bifurcates sub-and the other supercritically such that there is a bistable range between the homogeneous and hexagon solution.
We use a triangular domain to calculate branches of localized hexagons on homogeneous backgrounds for γ = 0.25. The solutions which bifurcate from the first and second bifurcation points of the cold hexagon branch are a single localized patch and multi localized patches of hexagons, respectively. Here only the single patch exhibits a snaking behavior with only one wiggle. By increasing the domain size, the number of wiggles should increase for both.
Under the assumption that the Landau coefficients c 2 , c 3 , c 4 change much slower than c 1 when varying σ away from σ c , we can see in (27) that the subcriticality increases, if c f := c 2 2 /(4(c 3 + 2c 4 ) 2 ) increases. c f evaluated in σ c as function of γ is shown in Fig.6(c). For γ ≈ 0.08 it holds c 3 + 2c 4 = 0 such that we assume that we can increase the strength of subcriticality and with that the front's steepness and the width of snaking ranges by choosing a γ-value closer to 0.08. A similar prediction of the strength of subcriticality and width of the snaking branch for bean branches on a Landau level can be found in [22]. In Fig.15(a) we see that the branch of single patches already shows a snaking behavior with more than one wiggle for γ = 0.12 over a domain which is smaller than the domain which we used for γ = 0.25. The front of the 574th solution becomes so steep that it looks like one single spot. From this point bifurcates a branch of a single hexagon patch which is rotated by π/6. Between the 574th and 900th solutions the branches in Fig.15(a) and (b) seem to be congruent. Behind the 900th solution the boundary affects the rotated patches and the branch moves to a branch of stretched hexagons. formalism and beyond all order asymptotics. There are not so many investigations of localized states over 2D domains. It is already pointed out in [21] on a Ginzburg-Landau level that all the localized 2D structures mentioned above should exist in reaction-diffusion systems. Some localized patches of hexagons and stationary fronts between hexagons and stripes are found for the quadratic-cubic Swift-Hohenberg equation in [15] via time integrations. Branches of localized single hexagon patches on a homogeneous background for this equation are first studied numerically and analytically in [19]. In the present paper we show the first numerical calculation of a branch of multi hexagon patches. Branches of localized hexagons on striped backgrounds are first investigated in [22]. We restrict ourselves to planar interfaces between hexagons and stripes. Localized hexagon patches on striped backgrounds are shown in the outlook of [22]. Here the hexagons are stretched and the authors were not able to find such patches for regular hexagons and neither the corresponding vertical fronts, and it is still an open problem whether patches of regular hexagons on striped backgrounds exist.
We also developed numerically some 1D global bifurcation diagrams. We found two stripe branches of the same wavelength, which are not connected, but connect to branches of a different wavelength (see Fig.8), and branches of the same wavelength, which are connected (see Fig.7). Understanding these behaviors is an interesting study for a further work.
From these studies we learn that a layering of two patterns is not necessarily an effect of space dependent bifurcation parameters but can occur for homogeneous bioirrigations. Thus we postpone a bifurcation analysis of patterns for depth dependent bioirrigations. Furthermore, we postpone the following issues. For γ = 0.25 exist rectangle patterns, when the stripes are stable and the hexagons are unstable. Their branch connects the branches of cold and hot hexagons. We computed branches over bounded domains of the form (−2lπ/k c , 2lπ/kc) × (−2π/( √ 3k c ), 2π/( √ 3k c )) and found stable rectangles, which bifurcate from the hot and cold hexagons and become unstable in the middle. Increasing l, increases the unstable part. Thus we suspect that the rectangle branch is always unstable over unbounded domains. However, over bounded domains of the form like above there is a bistable range between rectangles and stripes. When the rectangles lose their stability, a snaking branch (with stable parts) bifurcates for which the solutions look like fronts between rectangles and stripes, but we can not find where the branch terminates.
We did not find bistable ranges of hot and cold hexagons, which exist for instance for the Rayleigh-Bénard convection (see Fig.10 of [14]). In this bistable range a branch of triangle patterns exists which connects both hexagon branches. This is analogously to bean branches and we suspect that branches of fronts between cold and hot hexagons bifurcate from this triangle branch. In this paper we did not study any triangle patterns, because they fulfill periodic but not Neumann boundary conditions. The current version of pde2path does not handle periodic boundary conditions, but an add-on is already under construction.
It can already be seen by considering the Jacobians of (20) and (21) that the stripes which correspond to S + bifurcate stable in 1D and unstable in 2D for negative c 3 and positive c 2 and c 4 . We expect that some of the stable solutions presented here, trivially extended in z direction, are not stable over 3D domains. Since the sediment is 3D, an investigation of the bacteria-nutrient model (1) over 3D domains is appropriate. A 3D version of pde2path is under construction, and preliminary numerical studies of Turing bifurcations in 3D show that the many more different patterns make an analysis much more difficult, as expected. See also [18,1] for experimental and numerical investigations of 3D patterns. Classifications of 3D Turing patterns by symmetries and Landau reductions can be found in [8,9].