General conditions for Turing and wave instabilities in reaction -diffusion systems

Necessary and sufficient conditions are provided for a diffusion-driven instability of a stable equilibrium of a reaction–diffusion system with n components and diagonal diffusion matrix. These can be either Turing or wave instabilities. Known necessary and sufficient conditions are reproduced for there to exist diffusion rates that cause a Turing bifurcation of a stable homogeneous state in the absence of diffusion. The method of proof here though, which is based on study of dispersion relations in the contrasting limits in which the wavenumber tends to zero and to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\infty $$\end{document}∞, gives a constructive method for choosing diffusion constants. The results are illustrated on a 3-component FitzHugh–Nagumo-like model proposed to study excitable wavetrains, and for two different coupled Brusselator systems with 4-components.


Introduction
The formation of patterns and structures in a chemically reacting medium due to differential rates of diffusion of an otherwise stable equilibrium was first proposed in the last published work of Turing (1952). It seems Turing was on the verge of The original online version of this article was revised as the article title was incorrect.
B Edgardo Villar-Sepúlveda edgardo.villar-sepulveda@bristol.ac.uk Alan R. Champneys a.r.champneys@bristol.ac.uk explaining how this instability underlies the formation of structures in plants such as the daisy; see Dawes (2016). The recent issue of Philosophical Transactions (Krause et al. 2021) highlights the state of the art in how the patterns that emerge from a Turing instability have been used to explain many phenomena in biology. There is also a parallel literature in the wider theory of pattern formation across many different length scales of physics from micromechanics to spatial ecology; see also the other recent volumes (Yochelis et al. 2021;Champneys 2021) and references therein. Nevertheless, although the conditions for having a Turing bifurcation have been generalized to systems with n components, sometimes it remains unclear whether a system has the opportunity to develop these spatial patterns. The most definitive results so far are due to Satnoianu et al. (2000). They used the definition of a so-called S-stable Jacobian matrix to determine necessary and sufficient conditions for a Turing bifurcation in a wide class of reaction-diffusion systems with an arbitrary number of components. Their approach, which relies on the Routh-Hurwitz criterion, does not provide precise details on the form of the dispersion relations upon adding diffusion, nor definitive conditions on which parameters might lead to a bifurcation.
We note also the work (Klika et al. 2012) that extended some of the results in Satnoianu et al. (2000) to systems that have components with zero diffusion coefficients, and Hoang and Hwang (2013) that shows that linear instability implies a nonlinear bifurcation under certain technical assumptions. Finally, we mention the work of Kuznetsov and Polezhaev (2020) who find new ways of obtaining Turing bifurcations by considering what happens to dispersion relations in the limit that some diffusion coefficients tend to zero. Their results, which are mostly restricted to 3-component systems, provide a key motivation to the method introduced in this paper.
In all of these works, though, it is implied that the critical eigenvalue at the Turing bifurcation is real. The case of a complex critical eigenvalue, leading to a so-called wave bifurcation (or finite wavenumber Hopf bifurcation) has received less attention. A number of authors, e.g. Korvasová et al. (2015) have shown that such an instability requires at least a three-species model. A notable work is that by Anma et al. (2012) that gives sufficient conditions for such a bifurcation for 3-component systems. See Sect. 3 below for a recapitulation of their result.
Given the recent focus on synthetic biology, several authors have sought to understand the design principles of interacting chemical components necessary to achieve Turing instability, see e.g. Diego et al. (2018), Hambric et al. (2021). Scholes et al. (2019) have performed a comprehensive search of the linear properties of all possible network topologies of two or three interacting chemical species to see how common and robust Turing instability is. They sample both the functional forms of interaction and the values of diffusion constants and other parameter values. Their conclusion is Turing patterns are common but not robust in some sense. Note that they do not consider conditions for bifurcations (e.g. upon increasing diffusion rates) nor the possibility of wave bifurcations.
A similar approach is taken by Haas and Goldstein (2021), except they use a sample of values of the Jacobian matrices of linearised kinetics in reaction-diffusion systems with n components up to n = 6. They argue that the minimum diffusion threshold required for instability decreases with n.
In the present paper, we develop a general approach by studying dispersion relations of n-dimensional reaction-diffusion systems linearised around a homogeneous steady state. We suppose that a researcher is studying a reaction-diffusion system and asks the question whether, for the given kinetic terms, there are values of the diffusion coefficients that lead to either a Turing or a wave bifurcation. In order to answer this question, we take various limits, which strictly-speaking are likely to go outside the bounds for which the model is strictly valid, for example that certain diffusion coefficients vanish or are sufficiently large. This leads to conditions on principal submatrices of the homogeneous system, for which we can make precise statements about the existence of either kind of bifurcation for a finite range of values of the diffusion rates.
In what follows, we consider a system of n reaction-diffusion equations in time, t, and space x ∈ R m , where each state variable u i (x, t) i = 1, . . . , n satisfies the system of partial differential equations (PDEs) (1) (a) (b) (c) Fig. 1 Sketches of dispersion relations, that is the locus of the real part of an eigenvalue λ(k), for different parameter values. Upon parameter variation we have: a absence of instability; b a bifurcation point giving onset to instability; c unstable patterns exist for wavenumbers k 1 < k < k 2 . For simplicity we depict only the dispersion curve of a mode for which max k (Re(λ i (k))) is largest Let P be a homogeneous equilibrium point of (1). That is, f(P) = 0. We will denote the Jacobian matrix of this system at P by This implies that the linear part of the full-system at P is given by J u f(P) + D ∇ 2 . Next, as we want to analyze Turing and wave instabilities, then we will be interested in the space of eigenfunctions of the operator ∇ 2 with negative eigenvalues, i.e., vector functions u that satisfy the following equation: where k = (k 1 , k 2 , . . . , k m ) . This implies that the Jacobian matrix of our system at P, restricted to this space, is given by J u f(P) − |k| 2 D.
In particular, considering |k| = k as a real number, a diffusion-driven instability arises if P is stable in the absence of diffusion, yet becomes unstable in the presence of a specific real wavenumber k and a diffusion matrix D. In particular, we will assume that all eigenvalues of J u f(P) have negative real part, then we shall seek a diffusion matrix D and wavenumber k ∈ R such that J u f(P) − k 2 D has eigenvalues in the righthalf of the complex plane. In particular, let λ(k) be an eigenvalue of J u f(P) − k 2 D for which Re(λ(k)) > 0 for some k > 0. If we allow the diffusion rates (entries of D) to act as continuous parameters, then we are interested in the transition depicted in Fig. 1. In panel (a), we have no diffusion-driven instability since the system is stable for every value of k. Panel (b) depicts a bifurcation point, which we define to occur at a critical value of a parameter (either diffusion coefficient or otherwise) at which there is a quadratic tangency of the dispersion curve Re(λ(k)) at Re(λ(k)) = 0. This could be a Turing bifurcation or a wave one depending on whether λ (k * ) is real or pure imaginary, respectively. Finally, panel (c) depicts that the instability patterns that exist correspond to perturbations to P that have wavenumbers k 1 < k < k 2 .
The rest of this paper is outlined as follows. In Sect. 2 we give some preliminary results that are used to prove the main theorems of the paper, ending with a simple proof of a simple part of the main result in Satnoianu et al. (2000), which will be used later. Section 3 considers conditions for wave instabilities, confirming that the minimal number of components for such instabilities to occur is n = 3, and gives sufficient conditions for this case. Section 4 then provides conditions for both Turing and wave instabilities to occur in general n-component models. Section 5 contains examples that illustrate the theory developed throughout the paper and Sect. 6 contains concluding remarks.

Preliminary results
In this section we introduce some general results that will form the core methodology of this work, see specifically Lemmas 5 and 6 below. We build up these results via some elementary lemmas that are largely already known from literature on matrix stability, see e.g. Cross (1978), Satnoianu et al. (2000), Gantmacher (1959) but will prove useful in establishing the setting and notation of our theory.
We begin with the following elementary result which will apply to characteristic polynomials of matrices.
Proof First, as λ 1 , λ 2 . . . , λ m are roots of a polynomial equation with real coefficients, then all complex numbers come in complex-conjugate pairs. With this, their product can be decomposed as where the first product on the right-hand side of (3) has all the complex non-real eigenvalues, whilst the second and third ones are the real positive and negative eigenvalues, respectively. The first two products are positive, so the sign of the whole expression is given by the sign of the third term, which depends on the parity of positive roots, as stated.
Next we shall introduce a notation for the principal submatrices and principal minors of a square matrix A. Note that the notation used in Cross (1978) uses the complement of the set of indices used here.
Definition 1 For each integer 1 ≤ k ≤ n, we define a word to be a k-tuple of nonrepeating integers (i 1 , i 2 , . . . , i k ) with 1 ≤ i j ≤ n for each j = 1, 2 . . . , k, and we shall call it an increasing word if i 1 < i 2 < · · · < i k . For each increasing word, we refer to the corresponding principal submatrix S i 1 i 2 ...i k of J u f(P) as the k × k matrix that takes into account only the rows and columns with indices i 1 , i 2 , . . . , i k . We also denote by M i 1 i 2 ...i k the principal minor of J u f(P) which is the determinant of S i 1 i 2 ...i k .
With this, note that Lemma 1 can be used directly to prove the following: Lemma 2 There exists an integer 1 ≤ k ≤ n and an increasing word ..i k has an odd number of positive real eigenvalues.
Proof Note that the condition (−1) k M i 1 i 2 ...i k < 0 is equivalent to saying that (−1) k times the k eigenvalues of S i 1 i 2 ...i k is negative. Hence, owing to Lemma 1, this is true if and only if there are an odd number of positive real eigenvalues of S i 1 i 2 ...i k .
Next, we state a result that will be implicitly assumed in what follows under the requirement that P is a stable equilibrium point in the absence of diffusion.
Lemma 3 Let P be an equilibrium point of (1) in the absence of diffusion. If P is stable, then for all k = 1, 2, . . . , n.
Proof Based on a standard result on the determinant of sums of matrices (Marcus 1990), we know that for any pair of complex-valued n × n matrices, the determinant of A + B can be written as where the inner sum is over all strictly increasing words α, β of length r ; A[α|β] is the r × r -square submatrix of A formed of the rows in α and columns in β; B(α|β) is the (n − r )-square submatrix of B formed of the rows complementary to α and columns complementary to β; and s(α) (respectively, s(β)) is the sum of the integers in α respectively, β. In particular, when r = 0 this summand is equal to det(B) and, when r = n, it is det(A). Next, by definition, the characteristic polynomial of J u f(P) is given by The eigenvalues of J u f(P) are given by setting this polynomial to zero, which can be multiplied by (−1) n to obtain n m=0 (−1) n−m λ m Now, by the Routh-Hurwitz criterion (Routh 1877; Hurwitz 1895), a necessary condition for all the roots of this equation to lie in the left-half plane is The following result provides a useful characterisation of the characteristic polynomial in the presence of diffusion.
b 0 = det(J u f(P)), and b n = (−1) n n j=1 D j . Proof Let us take A = J u f(P), and B = −μ D in Eq. (4). In particular, in the notation of that equation, note that if α = β, then a row or column of B(α|β) would be equal to zero since B is a diagonal matrix. Thus the only non-zero contributions to the sum are gotten when α = β. Thus, notice that for r = n − m, and β = (i 1 , i 2 , i 3 , . . . , i n−m ), we have that which is exactly what is required to conclude the proof.
The above lemma, enables us to establish the following result, which is a generalisation of Theorem 1 in Satnoianu et al. (2000). Indeed, we do not require the Jacobian matrix J u f(P) to be S-stable for this result to be true. Nevertheless, the following Theorem states something only about Turing bifurcations, not considering wave instabilities.
Theorem 1 Let P be a stable equilibrium point of (1) in the absence of diffusion. If (−1) k M i 1 i 2 ...i k ≥ 0 for all k = 1, 2, . . . , n − 1, and each increasing word (i 1 , i 2 , . . . , i k ), then P will not be able to undergo Turing bifurcation when we take into account the full system (1).
Proof Note that if P is a stable equilibrium point of (1) in the absence of diffusion, then we have (−1) n det(J u f(P)) > 0.
Then, observe that the characteristic polynomial of the Jacobian matrix of our system after the addition of diffusion is given by Since we want to determine conditions in which we can see a Turing instability pattern around P, we want to look for conditions in which there exists μ > 0 such that the characteristic equation has a root λ = 0. That is, we want to find μ > 0 such that: which is a polynomial equation in μ, whose coefficients are described in Lemma 4. Next, notice that if (−1) k M i 1 i 2 ...i k ≥ 0 for all 1 ≤ k ≤ n − 1 and every increasing word (i 1 , i 2 , . . . , i k ), then all the coefficients in the polynomial (6) have the same sign. Therefore, there cannot exist a positive real root of (6) for any choice of diffusion constants. Thus, there is no choice of diffusion matrix D such that P admits a Turing bifurcation.
The next two lemmas will be crucial to the method we use to find conditions under which either Turing or wave instability patterns can occur in (1).
Lemma 5 Let P be an equilibrium point of system (1) in the absence of diffusion. Assume that D = 0 for all ∈ {i 1 , i 2 , . . . , i k }, for 1 ≤ k ≤ n, while all the other diffusion rates are positive. Then, there exist k eigenvalues of J u f(P) − μ D that tend to the eigenvalues of S i 1 ,i 2 ,...,i k as μ → ∞, while the real part of the others tend to −∞ as μ → ∞.
Proof First, notice that if k = n, then all diffusion rates would be equal to 0 and the eigenvalues of S i 1 ,i 2 ,...,i k would be those of J u f(P) which are independent of μ, thus making the result trivially true in this case.
Next, if k < n, although we might not be able to find all the eigenvectors of J u f(P)−μ D, note that we can get all its eigenvalues by solving the following equation: where λ m is an eigenvalue of J u f(P) − μ D, and v m is a corresponding eigenvector, for each m = 1, 2, . . . , n. Let us denote the elements of J u f(P) and v m as follows: Therefore, the eigenvalue problem is equivalent to the following set of equations for each i = 1, 2, . . . , n. In particular, notice that for each Then, at least one of the following two possibilities occur: 1. There exists an eigenvector v p , with p ∈ {1, 2, . . . , n}, and q ∈ {i 1 , i 2 , . . . , i k } such that lim μ→∞ v p,q = 0. Moreover, by continuity, v p,q = 0 for μ sufficiently large, in which case (8) implies we can write Therefore, without loss of generality, by assuming v p = 1, then everything on the left-hand side of (10), except possibly μ, is finite. Therefore, if we take the limit μ → ∞, we must have that lim With this condition, for μ sufficiently large, we know that (8) can be simplified to which can be rewritten as This means that λ p tends to an eigenvalue of S i 1 ,i 2 ...,i k , as μ → ∞, andv p to an associated eigenvector of that matrix. Note that, this gives us information only about k eigenvalues, although due to the possibility of generalised eigenvectors, we cannot necessarily find k independent eigenvectors via this process. 2. There exists an eigenvector v p , with p ∈ {1, 2, . . . , n}, such that lim . . , i k } and μ sufficiently large, we obtain the set of equations This means that λ p is an eigenvalue of the submatrix of J u f(P) − μ D formed by the rows and columns that are complementary to i 1 , i 2 , . . . , i k . Therefore, because the diagonal elements of this matrix tend to −∞ as μ → ∞ and the other elements are fixed, then by the Gershgorin Circle Theorem, the eigenvalue λ p must tend to −∞ as μ → ∞. Again, the number of independent eigenvectors v p that fulfill the condition of this case depend on the number of generalized eigenvectors of the submatrix in question.
Finally, the first (respectively, second) possibility gives us information about only k (respectively, n − k) eigenvalues. Therefore, both possibilities must occur simultaneously to describe the complete set of n eigenvalues of J u f(P) − μ D.
Lemma 6 Let P be a stable equilibrium of system (1) in the absence of diffusion. If we assume that D l → ∞ for all ∈ {i 1 , i 2 , . . . , i k }, for 1 ≤ k ≤ n, while all the other diffusion rates are positive, but finite, then there exist k eigenvalues of J u f(P) − μ D having real parts that tend to −∞ for every μ > 0, while n − k eigenvalues of that matrix tend to the eigenvalues of the submatrix of J u f(P) formed out of rows and columns that are complementary to i 1 , i 2 , . . . , i k , as μ → 0 + . Proof Using the same notation as in the previous proof, we need to solve the equations (7) for each i = 1, 2, . . . , n. Again, in this case we have two possibilities (note that the proofs are analogous but have some key differences in the sets of indices): 1. There exists an eigenvector v p , with p ∈ {1, 2, . . . , n} and q / ∈ {i 1 , i 2 , . . . , i k } such that lim v p,q = 0 for a fixed 0 < μ 1 sufficiently small. In this case, (7) implies that we can write If we assume, without loss of generality, that v p = 1, then everything on the right-hand side of (12) is finite whilst everything on the left-hand side is finite, except possibly μD i . Now if we let D i → ∞, while μ, sufficiently small, remains fixed, then we must have that lim and 0 < μ 1 sufficiently small. With this condition, for 0 < μ 1, we know that, for m = p, (7) can be simplified to As before, this means that λ p tends to an eigenvalue of the submatrix of J u f(P) that takes into account the rows and columns that are complementary to i 1 , i 2 . . . , i k , for 0 < μ 1 sufficiently small. Note that this gives us information only about n − k eigenvalues. Again, we highlight that this does not mean that we will always be able to find n − k eigenvalues that fulfill the condition required for this case. That depends on the number of generalized eigenvectors of the submatrix of J u f(P) that takes into account only its rows and columns that are complementary to i 1 , i 2 , . . . , i k . 2. There exists an eigenvector v p , with p ∈ {1, 2, . . . , n} such that lim . . , i k }, and 0 < μ 1 sufficiently small. Therefore, the set of equations we need to solve for i ∈ {i 1 , i 2 , . . . , i k } is given by: This means that λ p tends to an eigenvalue of the submatrix of J u f(P)−μD formed out of rows and columns i 1 , i 2 , . . . , i k . Therefore, as the off-diagonal elements of that matrix are constant, while the diagonal entries tend to −∞, then by the Gerschgorin circle theorem, the real part of λ p tends to −∞ as D → ∞, for each ∈ {i 1 , i 2 , . . . , i k }, and 0 < μ 1. Furthermore, note that as some diffusion rates are infinitely large, then increasing the value of μ above 1 will not impede the use of the Gerschgorin circle theorem to say that, in this case, the real part of the eigenvalues will tend to −∞ for every μ > 0.
Finally, as the first (resp. second) possibility gives us information only about n − k (resp., k) eigenvalues, then they both must occur in order to describe the n eigenvalues of J u f(P) − μ D. Now we introduce a useful definition and an important general result that enables us to find specific ratios of diffusion constants that give rise to an instability.
Definition 2 For each ∈ Z + , let s be the sorting function where V ⊂ [1 . . . n] l is a subset of vectors with no repeated entries, and j σ (1) , . . . , j σ ( ) is an increasing word.

this result is invariant under multiplication of all the diffusion constants by a positive real number.
Proof First, let ε > 0 and, for each (5) will have coefficients of the form ε p î 1 ,î 2 ,...,î n− , where p î 1 ,î 2 , . . . ,î n− is the sum of the powers of ε determined by the diffusion rates Dˆi q , whereî q / ∈ {i 1 , i 2 , . . . , i } for each integer 1 ≤ q ≤ n − . Next, notice that for 1 ≤ m ≤ n, we have that Here, note that the lowest power of ε will be p ( j n−m+1 , j n−m+2 , . . . , j n ). That coefficient will be the sum of all the lowest exponents of ε. All the other terms are multiplied by the same number of diffusion rates, but have at least a difference in one index with the diffusion rates multiplying (−1) m M s( j 1 , j 2 ,..., j n−m ) . This means that the other terms will have a higher order in terms of ε. With this, taking ε > 0 sufficiently small, we will have that for every integer 1 ≤ k ≤ n.
Last but not least, note that if we pick any δ > 0 and perform the change which will not change the sign of these terms.

Wave instability for n ≤ 3
In this section, we shall try to understand the minimal ingredients for wave instabilities to occur. Essentially, this result was originally stated in Anma et al. (2012, Th. 1.1), but we noted that the theorem contains some slightly stricter conditions on the diffusion constants D i and also includes a claim (result (iv) in the theorem) which is false, in general. For example, in the notation of the present paper, the matrix can easily be shown to provide a counterexample when taking D 1 = 3 and D 2 , D 3 small. Finally, the general way we state our result is easily generalisable to arbitrary n ≥ 3, which forms the subject of the following Section.
First we state our general result.
Proposition 1 Let P be a stable equilibrium point of (1) in the absence of diffusion. Then the minimum number of components that the system needs to develop wave instabilities around P is 3. Furthermore, if n = 3 and there exist 1 ≤ i 1 < i 2 ≤ 3 such that M i 1 + M i 2 > 0 and M i 1 i 2 > 0, then we can choose diffusion rates D 1 , D 2 , D 3 ≥ 0 such that (1) develops wave instability patterns around P.
Before proceeding with a proof, we establish a few useful standard lemmas.
Lemma 8 Let n = 2. Suppose that P is a stable equilibrium point of (1) in the absence of diffusion. Then the system will not admit a wave instability pattern around P for any diffusion matrix D.
Proof Let us label the elements of J u f(P) as follows Given that all the eigenvalues of J u f(P) are in the left-half plane, then its trace must be negative. That is α + δ < 0. To have a wave instability around P, we need to find μ > 0 such that the eigenvalues of J u f(P) − μD are pure imaginary. We have that This means that the eigenvalues of J u f(P) − μ D are given by Finally, note that, in the case that λ 1,2 are complex conjugate, then (13) shows that (λ 1,2 ) < 0 for all μ > 0.
Next, let us consider the case n = 3, and use the notation for principal minors established in Definition 1. The following two lemmas establish a general condition for finding non-negative diffusion rates for a wave instability, and a sufficient condition under which it cannot occur.
Lemma 9 Let n = 3. Then the condition for the matrix J u f(P) − μD to have a pure imaginary pair of eigenvalues for some μ > 0 can be expressed as where Proof The characteristic polynomial of J u f(P) − μ D can be written as where Therefore, when we evaluate this characteristic polynomial at λ = ωi, we will have the following expression Setting this expression to zero we get two expressions for ω 2 that are well defined if and only if a 1 > 0, a 2 = 0 and a 1 a 2 − a 0 = 0. Upon substitution of the definitions of a 0 , a 1 , a 2 into this final equality, we arrive at the polynomial in μ given at the statement of the result.
and det(J u f(P)) ≤ 0, then the system will not be able to show a wave instability pattern around P.

Proof
The proof is clear because, under the hypotheses, the coefficients in Lemma 9 b 0 , b 1 , b 2 , b 3 will be non-negative, which means that (14) cannot have a positive real solution for μ.

Remark 2
The previous lemma says that if all 2 × 2 principal submatrices of the linearisation around a homogeneous steady state of a 3-component system have only eigenvalues with negative real part, then the system will not be able to undergo a wave instability. Proposition 1 states the converse, namely that if there is an unstable 2 × 2 principal submatrix having two eigenvalues with a positive real part, then there will exist diffusion rates for which there is a wave instability.
Proof Most of the statement is a trivial consequence of Lemma 3 for the case n = 3. The only thing that remains to be proved is that b 0 > 0. Upon writing the characteristic polynomial of J u f(P) as p(λ) = λ 3 + c 2 λ 2 + c 1 λ + c 0 . Then, from the Routh-Hurwitz criterion, P is stable if and only if c 0 , c 1 , c 2 > 0 and (c 2 c 1 − c 0 )/c 2 > 0. But this final condition can be rearranged to read Finally, we can prove Proposition 1.

Proof of Proposition 1
Lemma 8 shows that wave instability cannot occur for n = 2. Therefore, let n = 3. Lemma 10 shows that no wave instability can occur if J u f(P) has only 2 × 2 principal submatrices with a pair of negative real-part eigenvalues. The condition on principle submatrices stated in the Theorem is essentially to assume that there is one principal 2 × 2 submatrix with a pair of eigenvalues with a positive real part. Without loss of generality, let us assume that M 2 + M 3 > 0 and M 23 > 0. Then we can choose D 2 = D 3 = 0. With this, note that b 3 = 0 and coefficient a 1 defined in (15) is a decreasing function on μ, which is positive when μ = 0 and it is equal to zero when Therefore, using the notation of Lemma 9, if we define the function then we have that q(0) = b 0 which, by Lemma 11, is strictly positive, while is strictly negative. Therefore, thanks to the Intermediate Value Theorem, there must exist at least one solution to the equation q(μ) = 0 in the interval 0 < μ < μ * .

General conditions for Turing and wave bifurcations
So far, we have established sufficient conditions for the absence of a Turing instability in a system around P (Theorem 1) along with the minimum number of components, a sufficient condition for a wave instability to occur, and Proposition 1. In this section, we shall establish further general conditions for system (1) to be able to show a Turing or wave instability pattern around P for any number of components, and to show existence of Turing or wave bifurcations as parameters vary. First, we are going to state a general standard result that gives us some conditions in which we can have neither a Turing nor a wave instability around P.

Theorem 3 Let P be a stable equilibrium point of (1) in the absence of diffusion.
Suppose that all diffusion rates are equal, that is, there exists δ ≥ 0 such that D i = δ for every integer 1 ≤ i ≤ n. Then (1) can admit neither a Turing nor a wave instability around P.
Proof Under the hypothesis on equal diffusion constants, note that the eigenvalues of J u f(P) − μD satisfy the equation det(J u f(P) − (μδ + λ) I ) = 0, which implies that δμ + λ is an eigenvalue of J u f(P). Therefore, because P is stable, (μδ + λ) < 0 implies (λ) < 0 − μδ < 0. Hence the eigenvalues of J u f(P) − μ D cannot cross the imaginary axis for any μ > 0.
We are now ready to state our main theorem, which develops sufficient conditions for system (1) to admit a Turing or wave instability pattern around P; see also Fig. 2 for the main idea.
Theorem 4 Let P be a stable equilibrium point of (1) in the absence of diffusion. If there exists an integer 1 ≤ ≤ n − 1, and a word ( j 1 , j 2 , . . . , j n ) of non-repeated integers between 1 and n (not necessarily increasing), such that S s( j 1 , j 2 ,..., j ) has q > 0 eigenvalues with positive real part, then if M s( j 1 , j 2 ,..., j k ) = 0 for each integer l ≤ k ≤ n and (−1) k M s( j 1 , j 2 ,..., j k ) changes its sign p times as k increases from l to n, then:

If p is even (respectively, odd), then there exist a choice of non-negative diffusion
rates D i , i = 1, . . . n such that J u f(P) − μ D has an even (odd) number of zero eigenvalues, bounded above by p as μ varies from 0 to ∞; 3. Furthermore, if q > p, then J u f(P) − μ D has at least q − p 2 pairs of complex conjugate eigenvalues with zero real part as μ varies from 0 to ∞; 4. In particular, with those diffusion rates, the number of crossings to the axis (λ) = 0 will be at least q as μ varies from 0 to ∞.
Proof Let p > 0 and q > 0 be defined as in the statement of the theorem. Note that the definition of p and q is independent of the value of the diffusion constants. Therefore, we can start by setting D m = 0 for all m = j 1 , j 2 , . . . , j . Next, thanks to Lemma 7, we will be able to find values for the other n − diffusion constants such that the signs of b k are the same as (−1) k M s( j 1 , j 2 ,..., j k ) , for each + 1 ≤ k ≤ n. This implies that, at least, q eigenvalues will cross the real axis as μ → ∞. On the other hand, as the number of complex conjugate eigenvalues that cross this axis is even, then as the polynomial of μ in Eq. (5) changes sign p times, if we assume that p is even (resp. odd), then we will have an even (resp. odd) number of real zeros crossing the μ-axis, which implies that q is even (resp. odd).
Next, note that we will have only p changes of sign in the polynomial equation for μ, which, by Descartes' rule of signs, implies that we will have an even (resp. odd) number of values of μ > 0 such that J u f(P) − μ D has a zero eigenvalue if p is even (resp. odd), while (if q > p) at least q − p 2 values of μ > 0 such that J u f(P) − μ D has a pair of pure imaginary eigenvalues. Therefore, by a matter of continuity, we will be able to increase the zero-diffusion rates by sufficiently small amounts in order to keep the same number of positive eigenvalues and eigenvalues with a positive real part, even though the number of crossings may increase. This concludes the proof.

Remark 3
Theorem 4 effectively characterizes zeroes we can find in dispersion curves as μ = k 2 tends to ∞. Figure 2 shows one of the simplest cases that can occur in the dispersion relation of a system that has five components. In Fig. 2a, three diffusion rates are equal to zero while the others are positive. In that graph, we can recognize immediately that a submatrix S i 1 i 2 i 3 is unstable and has two eigenvalues with a positive real part. Furthermore, when setting D i 1 = D i 2 = D i 3 = 0, there is an even number of changes of sign in the polynomial equation in μ. Next, in Fig. 2b, all the diffusion rates are positive, so the eigenvalues that were converging to a fixed value as k → ∞ are now decaying. Finally, Fig. 2c shows what happens as we keep increasing the diffusion rates that were equal to zero up to a point in which one of the dispersion relation curves becomes tangent to the k-axis and we find a bifurcation point. It is worth highlighting that these are only depicting the real part of the eigenvalues and the figure could equally well apply to a higher-component system, where some of the dispersion curves represent complex conjugate eigenvalues. See also Fig. 3 below. A natural consequence of Theorem 4 is its use in showing instabilities that might not have been apparent if only looking for Turing instabilities. The following is a straightforward consequence of the theorem.

Corollary 1 Let P be a stable equilibrium point in the absence of diffusion. If P does not fulfill the conditions to show Turing instability patterns, but one of the principal submatrices S i 1 i 2 ,...i k has an even number of eigenvalues with a positive real part, then we can find non-negative diffusion rates such that P exhibits wave instability patterns.
Remark 4 To illustrate why the corollary is important, consider the following Jacobian matrix: The eigenvalues of this matrix are given by λ 1 = −31.8681, λ 2 = −2.13507, λ 3,4 = −0.0484039 ± 2.2506i. Thus, Theorem 1 implies the matrix J u f(P) will not be able to have a zero eigenvalue when adding the diffusion terms. Nevertheless, note that S 234 has the following eigenvalues:

Furthermore, this matrix satisfies
This implies, when choosing the diffusion rates D 2 = D 3 = D 4 = 0, there will exist a value of μ > 0 in which the Jacobian matrix J u f(P) − μ D goes unstable. Thus, J u f(P) presents a counterexample to the conjecture stated in Wang and Li (2001, p. 144) on necessary and sufficient conditions for stability.

Remark 5
Theorem 4 provides information about the kind of transitions we can find for the eigenvalues of J u f(P) − μ D from having a negative to positive real part, as μ ranges from 0 to ∞. If the transition occurs for a real (resp. complex) eigenvalue, we will be able to find a Turing (resp. wave) instability pattern. Nevertheless, Lemmas 5 and 6 tell us that, even though these transitions could involve real (resp. complex) eigenvalues, the values to which the eigenvalues of J u f(P) − μ D converge when setting some diffusion rates as zero could be complex (resp. real). Such a dispersion relation would give rise to both kinds of patterns for the same values of parameters and which has eigenvalues λ ± = −9.50278 ± 81.3221i and λ 3 = −0.29445. Moreover, S 23 has eigenvalues λ 1 = 1.94868 and λ 2 = 0.0513167. This implies that when choosing D 2 = D 3 = 0, two eigenvalues of J u f(P) − μ D will tend to λ 1 and λ 2 as μ → ∞. However, according to Proposition 1, the transition of the eigenvalues from having a negative to a positive real part occurs in a complex way. Figure 3 shows what happens to the two largest real-part eigenvalues of J u f(P) − k 2 D in this case. In particular, the crossing with the k axis occurs at k ≈ 18.3146, where the eigenvalues are λ ± = ±0.827364i, and λ 3 = −3373.55. Nevertheless, there is a value of k for which two real eigenvalues get separated from this single line. This implies that, when increasing the zero diffusion rates a bit, by continuity, we will have two kinds of diffusion driven instability patterns for different values of k.
Last but not least, until now we have only considered instabilities, corresponding to the existence of a zero real-part eigenvalue, and not bifurcation points at which a dispersion curve first crosses the real axis. Our final results give conditions under which bifurcations occur under the variation of parameters (including diffusion constants). We start with a lemma.

Lemma 12
Let P be a stable equilibrium point of (1) in the absence of diffusion. If for some diagonal matrix D with non-negative entries and μ > 0, J u f(P) − μD has at least one eigenvalue with a positive real part for some μ = μ * , then there exist a diffusion matrix D * with positive diffusion rates such that P goes through a Turing or wave bifurcation under a change of parameter that takes D to D * .
Proof Thanks to Lemma 6, we know that if all the diffusion rates tend to ∞, then the real parts of the eigenvalues of J u f(P) − μD tend to −∞ for all μ > 0. Therefore, as there exist a matrix D with non-negative entries and μ > 0 such that J u f(P) has an eigenvalue with a positive real part, then by continuity, if we increase all the diffusion rates sufficiently enough, we will be able to find a matrix D * such that J u f(P) − μD * has one eigenvalue that is tangent to the axis (λ) = 0 for some μ > 0, which will give rise to a bifurcation.
In particular, this lemma is helpful when proving one last Theorem: Theorem 5 Suppose that the function f in (1) depends on a parameter ε ∈ R. Also, suppose that P is a stable homogeneous steady-state for all ε sufficiently small. Assume that there exists only one increasing word (i 1 , i 2 , . . . , i k ) such that S i 1 i 2 ...i k has a single real eigenvalue (respectively, a pure imaginary pair of eigenvalues) that cross the imaginary axis as ε increases through zero. Then, for ε > 0 sufficiently small, there exists a positive diffusion matrix D * such that P goes through a Turing (respectively, wave) bifurcation under small variation as diffusion constants are generically varied from their values in D * .
Proof The proof is a straightforward consequence of the ideas we have developed up to now. First, let us choose D i 1 = D i 2 = · · · = D i k = 0, while leaving all the other diffusion rates positive. This implies, thanks to Lemma 2, that the eigenvalues of J u f(P) − μD will converge to numbers with a negative real part when ε < 0, while some of them will converge to numbers with a positive real part when ε > 0. Therefore, thanks to the previous lemma, we can increase the zero diffusion rates such that we find a bifurcation. The type of bifurcation is characterized by continuity. In fact, for a sufficiently small value of ε, one or two eigenvalues of J u f(P) − μ D will be converging to one or two eigenvalues with a positive real part, so increasing the diffusion rates just a little will not change the type of eigenvalue we are having as μ → ∞.

Remark 6
Note that, owing to linearity and the last part of Lemma 7, only the ratios between the diffusion rates matter. The wavenumber k that corresponds to an instability can be made arbitrarily small or large, by scaling all diffusion constants by the same factor δ > 0.
To close the section, we can state the following result, which is a straightforward consequence of Theorems 2 and 4.

Theorem 6 Let P be a stable homogeneous steady-state of (1). Then linear diffusiondriven instabilities can occur around P if and only if the system is stable about P in the absence of diffusion but there is a principal unstable subsystem. In particular,
• If J u f(P) is S-stable, then the system cannot show a Turing or wave instability around P; or • If a principal submatrix of J u f(P) is unstable, then there exist non-negative diffusion rates such that the full system shows a Turing or wave instability around P.
Furthermore, in the latter case, the diffusion rates can be chosen as stated in Theorem 4.

Applications
We consider several different examples to illustrate the theory, in each case a model for a biophysical process that has arisen in the literature. The final three examples each illustrate how Theorem 4 helps classifying the crossings (e: could you restate this?) of dispersion curves with the axis (λ) = 0 as being either real or complex. Then Lemmas 5 and 6 give an idea about the behaviour of the eigenvalues of J u f(P) − μ D as μ → ∞, when some diffusion rates are large or equal to zero. To see the importance of those results, see Remark 5. In particular, we highlight that there are two cases that generate wave instabilities automatically after choosing the diffusion rates in a sensible way. These are: First, one of the principal submatrices of J u f(P) has a complex-conjugate pair of eigenvalues with positive real part, which generates wave instabilities for large wavenumbers, or second, the number of eigenvalues with a positive real part grows faster than the number of changes of sign in the terms (−1) k M s( j 1 , j 2 ,..., j k ) stated in Theorem 4.

Lack of instability in a model for host-parasyte interaction
First though, we provide an example where our theory can show the non-existence of instability. We consider a realistic model for the spread of malaria, based on the classic Ross-MacDonald model (see Alonso et al. 2019, and references therein). This model treats the dynamics of the infection between humans and mosquitoes. In particular, let H (resp. I ) be the population density of healthy (resp. infected) people. Besides, let P be the density of population of infected mosquitoes. We set the following system to model the dispersion of malaria: Note that healthy mosquitoes do not appear in this model. Nevertheless, they are considered under the assumption that Q is the total mosquito population, so that Q − P is the population of healthy mosquitoes. To preserve the significance of this model, we explicitly state that all the variables and parameters are non-negative, b H > d H > 0 and P ≤ Q. This system has only two equilibria, 0 = (0, 0, 0) and P = (H * , I * , P * ), where , , paper, the authors found numerical evidence for both a Turing and wave bifurcations for particular values of the parameters. System (17) has, at most, three homogeneous steady-states. We shall focus on the one given by P : which exists provided L > 0. The Jacobian matrix of (17), evaluated at P, can be readily computed to be from which the principal minors can be easily computed Note that M 2 , M 3 < 0, M 23 > 0, whereas M 1 , M 12 M 13 and M 123 can take either sign, depending on the value of R 2 . From Lemma 3, the solution P is stable if the following conditions are satisfied which lead to a non-open set that can readily be written in terms of the problem parameters using (18). From here on, we will assume that conditions (19) hold so that P is a stable equilibrium in the absence of diffusion. Next, using the results proven in the previous sections, we can state the following.
Proposition 2 M 1 > 0 if and only if there exist non-negative diffusion rates D 1 , D 2 , D 3 such that (17) shows a Turing instability pattern around P.
Proof Note from (18) that M 13 = − w M 1 , so that M 1 > 0 implies M 13 < 0. Note also that the sequence (−1)M 1 , (−1) 2 M 13 , (−1) 3 M 132 changes its sign only once. Furthermore, as M 13 < 0, then we know that S 13 has only one positive eigenvalue. Thus, according to Theorem 4, we can find non-negative diffusion rates such that (17) admits a Turing instability pattern around P.
On the other hand, if M 1 ≤ 0, then M 12 ≥ 0 and M 13 ≥ 0. This implies that (−1) 3 M 123 , (−1) 2 M i j and (−1)M i are all strictly non-negative for each 1 < i < j < 3. Then, Theorem 1 implies that (17) will not be able to develop a Turing instability at P for any positive diffusion constants.
Proposition 3 If M 1 + M 2 > 0 and M 12 > 0, then we can choose non-negative diffusion rates D 1 , D 2 , D 3 such that (17) admits a wave instability around P.
Proof If M 1 + M 2 > 0 and M 12 > 0, then S 12 will have two eigenvalues with a positive real part. Furthermore, we will have no change of sign from −M 123 to M 12 . Then, according to Theorem 4, we will be able to find non-negative diffusion rates such that (17) admits a wave instability pattern around P.
For numerical illustration of these results, we shall consider the following parameter values extracted from Yochelis et al. (2008): which can easily be seen to lie inside the parameter region in which P is stable. Furthermore, at these parameter values we have M 1 = 0.151497 > 0, M 2 = −0.087 < 0, M 3 = −1 < 0, M 12 = 0.18682 > 0, M 13 = −0.151497 < 0 M 23 = 0.087 > 0 and M 123 = −0.0868198 < 0, from which we see that M 1 + M 2 > 0. Thus, the above two propositions show that there must exist a set of non-negative diffusion rates such that P undergoes a Turing bifurcation, and another set under which it undergoes a wave bifurcation. The eigenvalues of S 12 are given by 0.0322485 ± 0.431022i, which have positive real part, and give conditions for a wave bifurcation. Moreover, S 13 has eigenvalues given by −1 and 0.151497 which will give rise to a Turing instability. Consider first the wave bifurcation. Let us define a parameter δ > 0 and consider D 1 = D 2 = δ and D 3 = 0.1. The corresponding dispersion curves are plotted in Fig. 4a for δ ranging from 0 to 0.01 in uniform steps. Note that when δ = 0 one of the dispersion relation curves does indeed tend to a finite positive value (≈ 0.032249), as μ → ∞ as shown by Lemma 5. Therefore, for δ > 0 sufficiently small, there must exist an interval of values of k for which the dispersion curve is positive. As δ is increased through a critical value δ * ≈ 2.27312 × 10 −3 , there is a point of tangency between this dispersion curve and the zero axis, at a critical wavenumber k * ≈ 2.43287. To check that this is indeed a wave zero, we have computed the critical eigenvalues of J u f(P) − μD to be −1.55429 and ±0.355382i. Furthermore, Fig. 4b shows the results of a numerical computation in one spatial dimension, which was performed at a δ-value just beyond that at which the bifurcation occurs, which indeed confirms the onset of a spatio-temporal instability.
To find a Turing bifurcation, let us introduce diffusion rates D 2 = 5, while D 1 = D 3 = δ, where again, δ is a parameter. Figure 5a shows dispersion curves for this case as δ ranges from 0 to 0.5 in uniform steps. In addition, the curve for δ ≈ 0.28684 shows the bifurcation parameter value, for which the critical wavenumber is k ≈ 0.496617, and the eigenvalues of J u f(P) − μ D are 0 and −1.15507 ± 0.293739i. Thus, we have a Turing bifurcation for these parameter values.
This result is illustrated in Fig. 5 where we show both the details of the dispersion relation in this case, and also a computation in two spatial dimensions.  u,v,w, respectively, at the Turing instability pattern obtained after integrating the system in two spatial dimensions with δ = 0.03. Here we make a perturbation in all components equally corresponding to wave numbers (k x , k y ) = (3, 2) with an amplitude ε = 10 −4 (color figure online)

Symmetrically coupled Brusselator systems
In this section, let us consider a model of two symmetrically coupled Brusselator systems studied by Konishi and Hara (2018) on a finite domain given by Here, u(x, t), v(x, t), w(x, t), and z(x, t) are scalar fields, a, b, c, η > 0 are reaction parameters while D 1 , D 2 , D 3 , D 4 ≥ 0 are diffusion rates. In Konishi and Hara (2018) it is shown numerically that the symmetrically coupled system may have two distinct Turing bifurcations. We show here how this fits with our theory. The model has at most five homogeneous steady states. We are interested in the simplest one Next, note that Therefore, tr (S 13 ) > 0 if and only if b > c η + 1, and det (S 13 ) > 0 if and only if b > 2 c η + 1. Thus, if 2 c η + 1 < a 2 + 1, there will exist b > 0 such that 2 c η + 1 < b < a 2 + 1, making S 13 have two eigenvalues with a positive real part, while P remains stable in the absence of diffusion.
To show that these really are independent Turing bifurcations, we can find the critical eigenvector of the operators J u f(P) + D ∂ x x in each case; see Fig. 7.

Asymmetrically coupled Brusselator systems
As a final example, we shall show how a small change to the previous example, breaking the symmetry of the coupling, can lead to the existence of a wave bifurcation. In particular, let us consider model (21) with the addition of a term −2c v − b a to the right-hand side of the z-equation, so that it becomes: We assume, again, that a, b, c > 0 and D 1 , D 2 , D 3 , D 4 ≥ 0, but now η may take either sign. This model has at most nine homogeneous steady states but, again, we focus only on the simplest one (22), the Jacobian matrix about which is the same as (23) except a switch in the sign of the (4, 2) entry.
Note now that the matrix has the property that when tr (S 24 ) = 0, M 24 = c 2 > 0. Thus, this submatrix has a pair of complex conjugate eigenvalues that cross the imaginary axis when η = − c a 2 . Therefore, if the steady state is stable then, according to Theorem 5, we can find positive diffusion rates such that the modified system goes through a wave instability when η < − c a 2 .
we find the eigenvalues of J u f(P) to be approximately −5.42567 ± 1.1069i, and −0.808327±8.62189i, which means that P is a stable homogeneous steady state. Furthermore, we find that S 24 has eigenvalues that are approximately 0.546 ± 4.8i, which means that this matrix is unstable, while (−1) 3 M 124 ≈ 151.552 and (−1) 3 M 234 = 151.552 are both positive. This implies that (−1) 4 M 1234 > 0, (−1) 3 M 124 > 0, (−1) 3 M 234 > 0 and (−1) 2 M 24 > 0. Therefore, we have an unstable submatrix of the Jacobian matrix of the system and, upon setting D 2 = D 4 = 0, there is a change of sign in the polynomial (5) as μ varies from 0 to ∞. Thus, Theorem 4 shows the existence of wave instabilities for non-zero values of D 2 and D 4 in this case. Explicit parameter values for the wave bifurcation can be found if we choose (D 1 , D 2 , D 3 , D 4 ) = (10, δ, 10, δ). Figure 8a shows dispersion curves for different δ-values ranging from 0 to 0.025 in equal steps. The figure also shows that a wave bifurcation occurs at δ = 0.016939.
In particular, when δ = 0.016939, we find that the critical wavenumber is k = 3.97209, and the eigenvalues of the Jacobian matrix are, approximately, −169.084, −159.468, and ±4.80815i, confirming that the system indeed undergoes a wave bifurcation. Finally, Fig. 8b shows the result of numerical simulation for a δ-value within the wave unstable parameter region, which indeed indicates the presence of a spatiotemporal instability.

Discussion
This paper has established general conditions for the existence of both Turing and wave instabilities for reaction-diffusion systems with an arbitrary number of components, generalising and improving earlier partial results in the literature. Even though the fact that a Turing instability pattern is a special case of a wave, there were no previous studies mixing both of these phenomena in systems with n components and they need to be considered together in order to be able to find and characterize diffusion-driven instabilities you might find in this kind of system. Note that all our results were stated for reaction-diffusion systems with diagonal diffusion matrices. In principle, though, the results can readily be extended to the case that D is not diagonal, that is, where there are cross-diffusion terms.
One simple way to extend the theory to the case that the off-diagonal diffusion coefficients are small, is to first set them to zero. Then one can use the results developed in this paper to establish conditions for Turing and wave bifurcations. Next, one can use a homotopy parameter to increase the off-diagonal diffusion rates, and the theory will still hold provided no eigenvalues with positive real part cross the imaginary axis during the homotopy.
An alternative way to treat cases with cross-diffusion, assuming that D is positive semi-definite and has n distinct real eigenvectors, could be to undertake a linear change of co-ordinates of the dependent variables to diagonalize D. Further changes of parametrisation may be necessary to avoid diffusion constants appearing in the expressions for the local reaction terms f.
At a more general level, everything we developed throughout the paper was done assuming that we are working with a reaction-diffusion system in R m . However, all of the theory can in principle be generalized, up to a few technical assumptions, to the case of a system ∂ t u = f(u) + D L(u), where L is a linear self-adjoint operator in u with negative eigenvalues, and boundary conditions accepting homogeneous steady states. In that case, the instabilities around those states we have described would spatially be described by eigenfunctions of that operator, rather than simple sinusoidal functions. An obvious example would be to consider L(u) to be a constant matrix times the Laplacian, defined on a disk or sphere, see e.g. Madzvamuse et al. (2015). We shall leave the details of all such generalizations to future work.
As a final comment, we note that the existence of a Turing bifurcation does not necessarily imply the onset of a stable periodic pattern. It seems in many reaction-diffusion systems the bifurcation may be subcritical, which means the periodic patterned states are unstable. Instead, in such systems, particularly on long domains, stable spatially localised patterns are typically observed; see e.g. Al Saadi et al. (2021), Champneys et al. (2021) and references therein. The determination of whether a Turing (or wave) bifurcation is super-or sub-critical requires the computation of certain cubic coefficients that arise as part of the normal form. This can also help understanding the type of wave that persists after a wave bifurcation since, depending on boundary conditions, there may be two types of waves that arise when a homogeneous steady-state goes through that bifurcation, namely the traveling wave and the standing wave (see Takada et al. 2022;Cuiñas et al. 2008;Kaminaga et al. 2005, and references therein). Such nonlinear analysis is beyond the scope of this paper, which has focused exclusively on linear calculations, but it will form the subject of future work by the present authors.