DESTABILIZATION THRESHOLD CURVES FOR DIFFUSION SYSTEMS WITH EQUAL DIFFUSIVITY UNDER NON-DIAGONAL FLUX BOUNDARY CONDITIONS

. This article deals with destabilizations of Turing type for diﬀusive systems with equal diﬀusivity under non-diagonal ﬂux boundary conditions. Stability-instability threshold curves in the complex plane are described as the graph of a piecewise analytic function for simple m -dimensional domains ( m ≥ 1). Also analyzed are eﬀects caused by imposing homogeneous boundary conditions of Dirichlet or Neumann type on appropriate portions of the domain boundary.


1.
Introduction. Turing's diffusion-induced-instability [7,6] for reaction-diffusion systems requires the presence of a substantial difference in the diffusion rates between activator and inhibitor species. This "restriction" universally applies to reacting and diffusing systems in which no interactions are involved on the boundary of domain, i.e., under no flux boundary conditions. In many cell-biological systems, however, relevant reagents have nearly equal diffusivities, and hence, the traditional Turing-mechanism may not be applied to such systems in order to account for the emergence of spatial or temporal inhomogeneities out of uniform states.
Levine and Rappel may be the first to point out in [5] that a mechanism of diffusion induced destabilization is still operative in diffusive systems with equal diffusivity, provided that suitable interactions on the boundary (cell membrane) are taken into consideration. In [5] they proposed a hypothetical toy model of reactiondiffusion equations with boundary flux-interaction. Together with linear stability analysis, numerical simulations are performed on this model to display diffusioninduced-instabilities of steady and oscillatory types under the equal diffusivity of reagents.
This article tries to put the rather surprising phenomena found in [5] on a mathematical footing. 642 KUNIMOCHI SAKAMOTO 1.1. Problem. The problem to be investigated in this article is slightly different from, but related to, the model treated in [5].
Let Ω ⊂ R m (m ≥ 1) be a bounded domain with smooth boundary ∂Ω. We consider the following system of diffusion equations under the Robin type boundary conditions ∂ t u = D u in Ω, D∂ n u = Ju on ∂Ω, (1.1) for the vector of N -species u(t, x) = (u 1 (t, x), . . . , u N (t, x)), where • ∂ t = ∂/∂t is the partial derivative with respect to time t; • is the m-dimensional Laplace operator; • D = diag(d 1 , . . . , d N ) is the diagonal diffusion matrix with d j > 0; • n is the unit outward normal vector field on ∂Ω; • ∂ n = ∂/∂n is the derivative along n on the boundary; • D∂ n u stands for the vector of fluxes projected onto n at the boundary, and • J is a non-diagonal N × N real matrix called the mass transfer matrix. The system (1.1) is the linearization of the nonlinear problem in which f : R N → R N is a nonlinear mapping modeling the mass transfer mechanisms across the boundary. The initial value problem for the nonlinear system is known to be well-posed (see, [1,3]). If u * ∈ R N is a constant solution of the non-linear problem, i.e., f (u * ) = 0, then (1.1) is the linearization of the nonlinear problem around the constant solution, where J = ∂ u f (u * ) is the Jacobian matrix of f at u * . We exclusively deal with the linear problem with emphasis on identifying the thresholds between stability and instability of the constant solution. Weakly nonlinear analyses, such as normal form and bifurcation analyses, will be presented elsewhere.
In order to determine whether the trivial solution u ≡ 0 ∈ R N of (1.1) is stable or unstable, the relevant eigenvalue problem is The complex number λ ∈ C is called an eigenvalue when (1.2) has a nontrivial solution φ ≡ 0. The eigenvalue problem (1.2) is obtained from (1.1) by substituting the ansatz u(t, x) = e λt φ(x). Therefore, if the eigenvalues of (1.2) have negative real part then the trivial solution of (1.1) is asymptotically stable. If, on the other hand, (1.2) has an eigenvalue with Re λ > 0, then the system (1.1) is unstable. 4) If J is unstable (but not necessarily symmetric) with a positive eigenvalue and D is close to a scalar matrix (nearly equal diffusivities), then (1.1) is unstable. In 3) and 4) above, the condition that D must be close to a scalar matrix is in general indispensable. Counter examples for 3) and 4) without this condition are given in [2], in which these examples were interpreted as the manifestations of the mechanisms of Turing-type for 3) and of anti-Turing-type for 4).
When the diffusion matrix D is scalar, i.e., D is a positive constant multiple of the N × N identity matrix, the following theorem was established in [2]. [2]). There exists an open set U ⊂ C with the following properties.
(i) The set U has the reflection symmetry with respect to the real axis, and satisfies U ⊂ {z ∈ C | Re z > 0}.  (v) If α/d ∈ C for some eigenvalue α of J, then (1.2) has a purely imaginary eigenvalue λ.
The threshold set C plays the role of the imaginary axis which is the stability/instability threshold in the usual context of reaction-diffusion systems under no flux boundary conditions, while S (or, U) corresponds to the left (or, right) half plane. The main focus of this article is to give a concrete representation of C ⊂ C in terms of a piecewise smooth function, mapping from the imaginary axis to the real axis, when the domain Ω has a simple geometry.
2. Main results. As mentioned in Introduction, our interest is to give a concrete description of the critical set C for m-dimensional (m ≥ 1) domains Ω with special geometries. It turns out that C is either analytic or piecewise analytic for such domains.
2.1. When ∂Ω has two connected components. We first deal with the situation in which ∂Ω consists of two connected components.
where M is a smooth closed manifold with dim M ≥ 1. Then, the threshold set C in Theorem 1.1 is the graph of a piecewise analytic function g : R → R, namely, C = {z ∈ C | Re z = g(Im z)} (and U := {z ∈ C | Re z > g(Im z)}) .
Moreover, for the scalar diffusion matrix D = diag(d, . . . , d), an eigenvalue λ of (1.2) transversely crosses the imaginary axis as α/d transversely crosses the curve C, where α is an eigenvalue of J.
Let us give a more concrete form to the function g. We use two functions g 0 (s) and g 1 (s) to define g(s) := min[g 0 (s), g 1 (s)] for s ≥ 0. Then for s ≤ 0, g(s) is simply defined by reflection, i.e., g(s) = g(−s) for s ≤ 0. The functions g k (s) (k = 0, 1) for s ≥ 0 are implicitly defined by The right hand sides in (2.1-Im) -(2.2-Re) appear to be the functions of √ 2τ , but actually they are analytic functions of τ . The right hand sides of (2.1-Im) and (2.2-Im) are strictly increasing functions of τ , taking all values s ≥ 0. Solving these equations in τ as τ = τ 0 (s) and τ = τ 1 (s), and substituting these into (2.1-Re) and (2.2-Re), respectively, we arrive at the definitions of g 0 (s) and g 1 (s). It is elementary to show that g 0 (s) = 1 3 s 2 + O(s 4 ) and g 1 (s) = 1 + 1 5 By plotting the graphs of Re z = g k (Im z) for k = 0, 1 as in Figure 1, we observe that these two curves intersect many times (and actually, infinitely many times). This is the reason why the critical curve C in Theorem 2.1 must be piecewise analytic, but not analytic. From the discussion so far, an interesting corollary follows.

2.2.
When ∂Ω is connected. A typical example of Ω with ∂Ω being connected is the flat Möbius band For the scalar diffusion matrices D = diag (d, . . . , d), the stability properties of (1.1) are characterized as follows. with the scalar diffusion matrix D = diag(d, . . . , d) transversely cross the imaginary axis iR\{0} back and forth infinitely many times as d varies over all positive real numbers.

2.3.
Other Ω with connected ∂Ω. We denote by B 2 the two dimensional unit open disk with center at the origin.
Then, the threshold set C in Theorem 1.1 is the graph of a piecewise-analytic function g : R → R. The functionĝ is even and satisfieŝ Moreover, for the scalar diffusion matrix D = diag(d, . . . , d), an eigenvalue λ of (1.2) transversely crosses the imaginary axis as α/d transversely crosses the curve C, where α is an eigenvalue of J.
To define the functionĝ, we first define a family of functions {ĝ n } ∞ n=0 . The graph of the functionĝ 0 , {z ∈ C | Re z =ĝ 0 (Im z)}, is the image of the imaginary axis under the mapping while the graph of the functionĝ n (n ≥ 1) is the image of imaginary axis under the mapping We then defineĝ(s) := inf{ĝ n (s) | n = 0, 1, . . .}. Numerical evaluations suggest thatĝ(s) =ĝ 0 (s), in which caseĝ(s) is analytic, instead of being piecewise-analytic. To rigorously prove this would require delicate analyses involving modified Bessel functions of non-negative integer orders. The graph ofĝ is numerically depicted in  To draw Figures 2 and 3, approximations are made by replacing the upper limit ∞ of the summations by 25, and λ = i τ (τ ∈ R) is substituted with |τ | ≤ 50 for Figure 2 and |τ | ≤ 500 for Figure 3.
has a nontrivial solution φ ≡ 0, where α ranges over the eigenvalues of J. This means that one can enumerate the eigenvalues of (1.2) as the family of eigenvalues λ j (α) of (3.1), parametrized by the eigenvalues of J. Since J is a real matrix, its non real eigenvalues appear in complex conjugate pairs α, α. Therefore, if the triplet (λ, φ, α) is a nontrivial solution of (3.1), then so is its complex conjugate (λ, φ, α). This observation will imply that the critical set C, to be defined below, is symmetric with respect to the real axis. The scalar problem (3.1) is now reformulated in terms of the so called Dirichletto-Neumann map. Let us begin with recalling the definition and several properties of this map (see, [4]).
Consider Dirichlet boundary value problem for the C-valued function v(x): where p ∈ H 3/2 (∂Ω) is a given function (Dirichlet data) and λ ∈ C plays the role of a parameter. The Dirichlet-to-Neumann map is defined as follows.
The definition of T (λ) implies that (3.1) is equivalent to This allows us to approach the scalar problem (3.1) from another viewpoint. Namely, we now ask how the eigenvalues of T (λ) depend on λ ∈ C\σ( Dir ). Let us denote by E(λ) an eigenvalue of T (λ) for λ ∈ C\σ( Dir ). From the observation made above in the paragraph following (3.1), we find that E(λ) satisfies E(λ) = E(λ). It is expected that E(λ) depends analytically in λ. In fact, for the domain Ω under consideration in Theorems 2.1, 2.2 and 2.3 (as well as Theorems 4.1 and 4.2, in §4), we will show that the eigenvalues E(λ) are analytic in λ. This will be done in general context as follows. First of all, note that E(0) must be one of the Steklov eigenvalues {ν j } ∞,1 j=0 . We then analytically continue E(0) = ν j for all possible values λ ∈ C, and denote the continuation by E j (λ). It is shown in [2] that such analytic continuations are possible for all ν j that is a simple eigenvalue of T (0). In the situations of Theorems 2.1, 2.2 and 2.3, however, we can bypass these steps, since the Fourier expansions of H 2 (Ω)-functions allow us to directly obtain the analyticity of E(λ), even when some of the Steklov eigenvalues are not simple.
3.1. Proof of Theorem 2.1. We now compute the eigenvalues E 0 (λ) and E 1 (λ) for Ω = (−1, +1) ⊂ R. In this case, it is easy to show that the Dirichlet-to-Neumann map T (λ) is represented by the 2 × 2 matrix The eigenvalues E 0 (λ) and E 1 (λ) of the matrix are given by By expanding the numerator and denominator, one finds that these eigenvalues are analytic functions of λ. The domains of analyticity are, respectively, Moreover, the excluded sets together constitute σ( Dir ); Let us denote by C k the image of E k (λ) (k = 0, 1) as λ ranges over the imaginary axis {λ = iτ | τ ∈ R}. The parametric representations of these curves are respectively given by (2.1) and (2.2), or equivalently, Re z = g 0 (Im z) for C 0 and Re z = g 1 (Im z) for C 1 .
It is elementary to show that d dλ E(λ) λ=iτ = 0 for all τ ∈ R. To prove the transverse crossing of eigenvalues, we use the Cauchy-Riemann relation for k = 0, 1. Note that ∂E k (λ) ∂τ is the tangent vector along C k , while ∂ ∂ς means that λ transversely crosses the imaginary axis from left to right in the right angle π/2. Then the relation in (3.4) says that E k (λ) transversely crosses C k from left to right in the right angle π/2. By reversing the roles, this means that as E k (λ) transversely crosses the curve C k from left to right, λ transversely crosses the imaginary axis from left to right. Moreover, the left half plane minus σ( Dir ) is mapped by E k (·) into the domain left to the curve C k , while the right half plane is mapped into the domain right to the curve. To show that the curves C 0 and C 1 intersect infinitely many times, let τ > 0 be such that √ 2τ = n + 1 2 π for an even integer n > 0. Then, the corresponding values of s, g k (s) for (2.1) and (2.2) satisfy For n + 1 (odd integer) the relations are reversed as follows.
If τ is such that √ 2τ = nπ (n ∈ N), then s 1 ), i.e., C k intersect the diagonal line Re z = Im z infinitely many times. Therefore, the curves C k (k = 0, 1) interlace with each other around the diagonal Re z = Im z, and we conclude that they intersect infinitely many times. The proof of Theorem 2.1 is complete for Ω = (−1, +1) ⊂ R.
To prove the theorem for Ω = (−1, +1) × M , we only need to describe the eigenvalues of T (λ). Applying the variable separation, we find that the eigenvalues are given by E 0 (λ + κ j ) and j=0 is the set of the eigenvalues of M , the Laplacian on M . These eigenvalues satisfy Thanks to the arguments above, it is easy to see that the curves generated by E k (iτ + κ j ) for j ∈ N, k = 0, 1 are shifted to the right by the amount corresponding to κ j > 0. Therefore, the threshold curves are determined by E k (iτ + κ 0 ) = E k (iτ ) which are the same as the critical eigenvalue for the interval Ω = (−1, +1). This completes the proof of Theorem 2.1.

Proof of Theorem 2.2. On the flat Möbius band
our eigenvalue problem is λu = u xx + u θθ for (x, θ) ∈ (−1, +1) × (0, 2π), where u must satisfy u(x, 0) = u(−x, 2π). Therefore, in the variable separation u(x, θ) = A(x)B(θ), A is even, A(x) = A(−x), and B(θ) must be 2π-periodic in θ. Substituting the ansatz, and using the 2π-periodicity, we have a family of problems parametrized by non-negative integers k. The solutions are exhausted by the constant multiples of e √ λ+k 2 x + e − √ λ+k 2 x . Therefore, the eigenvalues of the Dirichletto-Neumann map in this case are given by This is nothing but E 0 (λ+k 2 ), where E 0 is the one that has appeared in §3.1. Therefore, the remaining part of the proof has been already supplied. This completes the proof of Theorem 2.2.
In terms of the modified Bessel function [8] of the first kind of order n, I n (z), the solution is expressed as and hence the corresponding eigenvalue E n (λ) of T (λ) is From the series expansion it follows that E n (λ) is analytic in λ: Let {j n,k | k ∈ N} be the set of the positive zeros of the Bessel function J n (z). Since, by definition, I n (z) = e −inπ/2 J n (iz), we find from (3.5) that E n (λ) is defined on C\{−(j n,k ) 2 | k ∈ N}. The excluded sets together constitute the Dirichleteigenvalues of the Laplacian on B 2 ; σ( Dir ) = ∞ n=0 {−(j n,k ) 2 | k ∈ N}. Clearly, the Steklov eigenvalues are given by {E n (0) = n} ∞ n=0 . For each non-negative integer n, let us define the curve C n := {E n (iτ ) | τ ∈ R}, the image of the imaginary axis under the analytic map E n ( · ). Evidently, C n intersects the real axis at n, one of the Steklov eigenvalues.

Fixed boundary conditions on part of the boundary.
It is technically difficult to study (1.2) for rectangular domains Ω. This is partly due to the fact that the boundary of such a domain is piecewise smooth but not smooth. For such domains, the Dirichlet-to-Neumann map tends to be singular. Moreover, the procedure of separating variables do not work for rectangular domains. However, by imposing fixed boundary conditions on appropriate portions of the boundary, and thus making smooth the free remaining boundary portions on which the Dirichlet-to-Neumann map is defined, our analysis easily extends to piecewise smooth domains. Let us, therefore, consider the following variant of (1.1).
the threshold curve is given by C = {z ∈ C | Re z = g(Im z)}, where g is the same function as appeared in §2.1.
(v) For the solid cylinder Ω = (−1, +1) × B 2 with the threshold curve is given by C = {z ∈ C | Re z = g 0 (Im z)}, where g 0 is the same function as appeared in §2.1.
We now apply the variable separation u(x, y) = A(x)B(y), and obtain the onedimensional problem for each non-negative integer k. The proof now follows from that of Theorem 2.1.
(ii) Arguing as in the proof for (i) above, we arrive at The proof is completed by following the same line of arguments in §3.
The proofs of (iii), (iv) and (v) are similar.