#  Instability of all regular stationary solutions
to reaction-diffusion-ODE systems

Szymon Cygan  ID Instytut Matematyczny, Uniwersytet Wrocławski, pl.
Grunwaldzki 2/4, 50-384 Wrocław, Poland
orcid.org/0000-0002-8601-829X<EMAIL_ADDRESS>http://www.math.uni.wroc.pl/ scygan ,  Anna Marciniak-Czochra  ID Institute of
Applied Mathematics, Interdisciplinary Center for Scientific Computing (IWR)
and BIOQUANT, University of Heidelberg, 69120 Heidelberg, Germany
orcid.org/0000-0002-5831-6505<EMAIL_ADDRESS>http://www.biostruct.uni-hd.de/ ,  Grzegorz Karch  ID Instytut Matematyczny,
Uniwersytet Wrocławski, pl. Grunwaldzki 2/4, 50-384 Wrocław, Poland
orcid.org/0000-0001-9390-5578<EMAIL_ADDRESS>http://www.math.uni.wroc.pl/ karch and  Kanako Suzuki ID Graduate School of
Science and Engineering, Ibaraki University, 2-1-1 Bunkyo, Mito 310-8512,
Japan
orcid.org/0000-0001-8018-9621<EMAIL_ADDRESS>
###### Abstract.

A general system of several ordinary differential equations coupled with a
reaction-diffusion equation in a bounded domain with zero-flux boundary
condition is studied in the context of pattern formation. These initial-
boundary value problems may have regular (i.e. sufficiently smooth) stationary
solutions. This class of close-to-equilibrium patterns includes stationary
solutions that emerge due to the Turing instability of a spatially constant
stationary solution. The main result of this work is instability of all
regular patterns. It suggests that stable stationary solutions arising in
models with non-diffusive components must be far-from-equilibrium exhibiting
singularities. Such discontinuous stationary solutions have been considered in
our parallel work [Stable discontinuous stationary solutions to reaction-
diffusion-ODE systems, preprint (2021)].

###### Key words and phrases:

Reaction-diffusion equations; stationary solutions, stability, close-to-
equilibrium patterns

###### 2010 Mathematics Subject Classification:

35K57; 35B35; 35B36; 92C15

The authors appreciated helpful comments on this work by Dr. Chris Kowall. S.
Cygan acknowledges a support by the Polish NCN grant 2016/23/B/ST1/00434. This
work was also supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 -
390900948 (the Heidelberg STRUCTURES Excellence Cluster and SFB1324 (B05), and
JSPS the Grant-in-Aid for Scientific Research (C) 18K03354 and 19K03557.

##  1\. Introduction

We study properties of solutions of a general system of $n$ ordinary
differential equations coupled with a scalar reaction-diffusion equation

(1.1) $\displaystyle\boldsymbol{u}_{t}$
$\displaystyle=\boldsymbol{f}(\boldsymbol{u},v),$ $\displaystyle
x\in\overline{\Omega},\quad t>0,$ $\displaystyle v_{t}$ $\displaystyle=\Delta
v+g(\boldsymbol{u},v),$ $\displaystyle x\in\Omega,\quad t>0,$

with an unknown vector field (denoted using a bold-face font) and a scalar
function:

(1.2)
$\displaystyle\boldsymbol{u}=\boldsymbol{u}(x,t)=\left(\begin{array}[]{c}u_{1}(x,t)\\\
\vdots\\\ u_{n}(x,t)\end{array}\right)\qquad\text{and}\qquad v=v(x,t).$

We consider arbitrary $C^{2}$-nonlinearities

(1.3)
$\displaystyle\boldsymbol{f}=\boldsymbol{f}(\boldsymbol{u},v)=\left(\begin{array}[]{c}f_{1}(\boldsymbol{u},v)\\\
\vdots\\\ f_{n}(\boldsymbol{u},v)\end{array}\right)\qquad\text{and}\qquad
g=g(\boldsymbol{u},v)$

and define system (1.1) in a bounded open domain $\Omega\subset\mathbb{R}^{N}$
for $N\geq 1$ with $C^{2}$-boundary $\partial\Omega$. Moreover, the reaction-
diffusion equation in (1.1) is supplemented by the homogeneous Neumann
boundary condition

(1.4) $\displaystyle\partial_{\nu}v=0\qquad\text{for}\quad
x\in\partial\Omega,\quad t>0,$

where $\partial_{\nu}=\nu\cdot\nabla$ with $\nu$ denoting the unit outer
normal vector to $\partial\Omega$. We also impose an initial condition

(1.5)
$\displaystyle\boldsymbol{u}(x,0)=\boldsymbol{u}_{0}(x)=\left(\begin{array}[]{c}u_{0,1}(x)\\\
\vdots\\\ u_{0,n}(x)\end{array}\right),\qquad v(x,0)=v_{0}(x).$

It is well-known that problem (1.1)-(LABEL:ini) has a unique local-in-time
solution corresponding to a continuous initial condition, see e.g. [20, 21].
Our goal is to study stability of stationary solutions
$(\boldsymbol{U},V)=\big{(}\boldsymbol{U}(x),V(x)\big{)}$ satisfying the
relation

(1.6) $\displaystyle\boldsymbol{f}\big{(}\boldsymbol{U}(x),V(x)\big{)}=0,\quad
x\in\Omega$

and the boundary value problem

(1.7) $\displaystyle\Delta V+g(\boldsymbol{U},V)=0,$ $\displaystyle
x\in\Omega,$ $\displaystyle\partial_{\nu}V=0,$ $\displaystyle
x\in\partial\Omega.$

Problem (1.6)-(LABEL:eq:ImFun2) may have different types of solutions. Here,
we generalize results from works [16, 15, 21, 9, 22, 8, 19, 20, 6, 10, 24, 36,
37, 39] where the following two types of stationary solutions were considered.

  1. (1)

Stationary solutions of the first type are called regular. In this case, the
function $\boldsymbol{U}(x)$ is obtained from $V(x)$ by solving equation (1.6)
in such a way that $\boldsymbol{U}(x)=\boldsymbol{k}\big{(}V(x)\big{)}$ for
some $C^{2}$-function $\boldsymbol{k}=(k_{1},...,k_{n})$.

  2. (2)

The second class of stationary solutions is called jump-discontinuous, because
$\boldsymbol{U}(x)$ is obtained from $V(x)$ in a discontinuous way by choosing
different branches of solutions to equation (1.6).

In this work, we deal with the first class of solutions and we show that all
regular stationary solutions are unstable – see Theorems 2.7 and 2.8 below for
precise statements of our results. In our second work [3], we study
discontinuous stationary solutions to problem (1.1)-(LABEL:ini) and we find
sufficient conditions for their existence and stability. Other types of
stationary solutions, different from those mentioned above, may also exist and
we comment them in Remark 2.10 below.

It is well known that a single autonomous reaction diffusion equation
$v_{t}=\Delta v+h(v)$ in a bounded convex domain and with the Neumann boundary
condition does not support interesting patterns and only constant solutions
can be stable, see Casten-Holland [1] and Matano [23]. Similarly, it was shown
by Kishimoto and Weinberger [12] that there is no stable non-constant
equilibrium solutions of homogeneous cooperative-diffusion systems in convex
domains with no-flux boundary conditions. Result from [12] were recently
extended by Wang [38] on reaction-diffusion-ODE cooperative systems.

In this paper, we contribute to the theory by showing that all close-to-
equilibrium (regular) patterns of a general reaction-diffusion-ODE problem
(1.1)-(LABEL:ini) are unstable, regardless of the particular structure
assumption on nonlinearities. In particular, it implies that such models
cannot exhibit stable Turing patterns and the only possible stable stationary
solutions of such models have to be somewhat singular or discontinuous, except
a certain marginal class of degenerate models discussed below in Remark 2.10.
In this way, we generalize a series of previous works where either an
existence or stability of either regular or discontinuous stationary solutions
have been analysed in the case of specific reaction-diffusion-ODEs from
mathematical biology [20, 21, 13, 16, 9]. Finally, we refer the reader to our
work [3] for a discussion of discontinuous stationary solutions for some
particular models considered e.g. in papers [36, 8, 9, 13, 18, 15, 16, 20, 21,
22, 37]. Related results on stationary solutions of reaction-diffusion-ODE
systems can be also found in [2, 25, 31, 39, 10, 24, 27].

In Section 2, we present and discuss main results of this work. In Proposition
2.5, we construct a non-constant family of regular stationary solutions to
problem (1.1)-(LABEL:ini) using the bifurcation theory. The proof of that
proposition is postponed to Section 3. Instability results for regular
stationary solutions are stated in Theorem 2.7 and Theorem 2.8 whose proofs
are contained in Section 5. In Section 4, we characterize spectra of
linearised operators corresponding to problem (1.1)-(LABEL:ini).

### Notation.

By the bold font, e.g. $\boldsymbol{A},\boldsymbol{u},$ we denote either
matrices or vector valued functions in order to distinguish them from scalar
quantities. For $\boldsymbol{f}$ and $g$ defined in (1.3) we set

(1.8) $\displaystyle\boldsymbol{f}_{\boldsymbol{u}}$
$\displaystyle=\begin{pmatrix}\frac{\partial f_{1}}{\partial
u_{1}}&\cdots&\frac{\partial f_{1}}{\partial u_{n}}\\\ \vdots&&\vdots\\\
\frac{\partial f_{n}}{\partial u_{1}}&\cdots&\frac{\partial f_{n}}{\partial
u_{n}}\\\ \end{pmatrix},$ $\displaystyle\boldsymbol{f}_{v}$
$\displaystyle=\begin{pmatrix}\frac{\partial f_{1}}{\partial v}\\\ \vdots\\\
\frac{\partial f_{n}}{\partial v}\\\ \end{pmatrix},\;g_{\boldsymbol{u}}$
$\displaystyle=\begin{pmatrix}\frac{\partial g}{\partial
u_{1}}&\cdots&\frac{\partial g}{\partial u_{n}}\\\ \end{pmatrix},$
$\displaystyle\;g_{v}$ $\displaystyle=\frac{\partial g}{\partial v}.$

$Y^{n}={Y\times\cdots\times Y}$ ($n$-times) is the product of a given space
$Y$ with the norm denoted by $\|\boldsymbol{y}\|_{Y}$ (instead of
$\|\boldsymbol{y}\|_{Y^{n}}$). The symbol $\sigma(L)$ means the spectrum of a
linear operator $\big{(}L,D(L)\big{)}$ and $s(L)=\sup\\{{\rm
Re\,}\lambda:\lambda\in\sigma(L)\\}$ is the spectral bound of an operator
$\big{(}L,D(L)\big{)}$. The usual Sobolev space $W^{1,2}(\Omega)$ has the
scalar product

(1.9) $\displaystyle\langle u,v\rangle_{1,2}=\int_{\Omega}\nabla u\cdot\nabla
v\,dx+\int_{\Omega}uv\,\text{d}x.$

We use the symbol $\Delta_{\nu}$ for the Laplace operator in a bounded domain
$\Omega$ with the Neumann boundary condition. It is defined in the usual way
via a bilinear form on $W^{1,2}(\Omega)$, and $-\Delta_{\nu}$ has the
eigenvalues $\mu_{k}$ satisfying
$0=\mu_{0}<\mu_{1}\leqslant\cdots\leqslant\mu_{n}\to\infty$. Constants in the
estimates below are denoted by the same letter $C$, even if they vary from
line to line. Sometimes we shall emphasize dependence of such constants on
parameters used in our calculations.

##  2\. Main results

The goal of this work is to construct certain non-constant stationary
solutions of problem (1.1)-(LABEL:ini) and to study their stability. Thus, we
deal with a solution $(\boldsymbol{U},V)=\big{(}\boldsymbol{U}(x),V(x)\big{)}$
to the boundary value problem

(2.1) $\displaystyle\boldsymbol{f}(\boldsymbol{U},V)$ $\displaystyle=0,$
$\displaystyle x\in\overline{\Omega},$
$\displaystyle\Delta_{\nu}V+g(\boldsymbol{U},V)$ $\displaystyle=0,$
$\displaystyle x\in\Omega,$

with arbitrary $C^{2}$-functions $\boldsymbol{f}$ and $g$ of the form (1.3)
and in a bounded domain $\Omega\subset\mathbb{R}^{N}$ with $C^{2}$-boundary.

######  Definition 2.1.

A pair $(\boldsymbol{U},V)=\big{(}\boldsymbol{U}(x),V(x)\big{)}$ is a weak
solution of problem (LABEL:seq1) if

  * •

$\boldsymbol{U}$ is measurable,

  * •

$V\in{W^{1,2}(\Omega)}$,

  * •

$g(\boldsymbol{U},V)\in\big{(}{W^{1,2}(\Omega)}\big{)}^{\star}$ $($the dual of
the space ${W^{1,2}(\Omega)})$,

  * •

the equation $\boldsymbol{f}\big{(}\boldsymbol{U}(x),V(x)\big{)}=0$ is
satisfied for almost all $x\in\Omega$,

  * •

the equality

(2.2) $\displaystyle-\int_{\Omega}\nabla
V(x)\cdot\nabla\varphi(x)\,\text{d}x+\int_{\Omega}g\big{(}\boldsymbol{U}(x),V(x)\big{)}\varphi(x)\,\text{d}x=0$

holds true for all test functions $\varphi\in{W^{1,2}(\Omega)}$.

######  Definition 2.2.

A weak solution in the sense of Definition 2.1 to problem (LABEL:seq1) is
called regular if $\boldsymbol{U},V\in{L^{\infty}(\Omega)}$ and, moreover,
there exists a $C^{2}$-function $\boldsymbol{k}:\mathbb{R}\to\mathbb{R}^{n}$
such that $\boldsymbol{U}(x)=\boldsymbol{k}(V(x))$ for all $x\in\Omega$.

######  Remark 2.3.

Notice that every regular solution of problem (LABEL:seq1) satisfies the
equation

(2.3)
$\displaystyle\boldsymbol{f}\big{(}\boldsymbol{U}(x),V(x)\big{)}=\boldsymbol{f}\big{(}\boldsymbol{k}(V(x)),V(x)\big{)}=0\quad\text{for
all}\quad x\in\Omega,$

where $V=V(x)$ is a solution of the elliptic Neumann problem

(2.4) $\displaystyle\Delta_{\nu}V+h(V)=0\qquad\text{for}\quad x\in\Omega$

with $h(V)=g\big{(}\boldsymbol{k}(V),V\big{)}$. Since we require
$V\in{L^{\infty}(\Omega)}$ and since $\partial\Omega$ is $C^{2}$, by a
standard elliptic regularity, we have $V\in W^{2,p}(\Omega)$ for each
$p\in(1,\,\infty)$ and consequently $\boldsymbol{U}\in W^{2,p}(\Omega)^{n}$
for every $p\in(1,\infty)$. In particular, every regular stationary solution
satisfies $(\boldsymbol{U},V)\in C(\overline{\Omega})^{n+1}$.

######  Remark 2.4.

Assume that $\big{(}\boldsymbol{U}(x),V(x)\big{)}$ is a non-constant regular
solution to problem (2.1). Then, there exists $x_{0}\in\Omega$ such that a
vector
$(\overline{\boldsymbol{U}},\overline{V})=\big{(}\boldsymbol{U}(x_{0}),V(x_{0})\big{)}$
is a constant solution to this problem. To prove this fact, it suffices to
integrate the second equation in (2.1) over $\Omega$ and to use Neumann
boundary condition to obtain
$\int_{\Omega}g\big{(}\boldsymbol{U}(x),V(x)\big{)}dx=0$. Since
$\boldsymbol{U}$ and $V$ are continuous there exists $x_{0}$ such that
$g\big{(}\boldsymbol{U}(x_{0}),V(x_{0})\big{)}=0$ and by the first equation in
(2.1) we obtain
$\boldsymbol{f}\big{(}\boldsymbol{U}(x_{0}),V(x_{0})\big{)}=0$.

Solutions to the general Neumann elliptic problem (2.4) has been constructed
in several works. For example, classical bifurcation methods has been used in
the works [33, 32, 15, 16] and the phase portrait method in [20]. We refer the
reader also to the review by Ni [26] and to the references therein. Here, we
recall one possible construction of solutions to problem (2.4) using a
variational approach to bifurcation methods.

######  Proposition 2.5.

Let $N\leqslant 6$. Consider a constant stationary solution of problem
(LABEL:seq1), namely, the constant vector
$\big{(}\overline{\boldsymbol{U}},\overline{V}\big{)}\in\mathbb{R}^{n+1}$ such
that

(2.5)
$\displaystyle\boldsymbol{f}(\overline{\boldsymbol{U}},\overline{V})=0\quad\text{and}\quad
g(\overline{\boldsymbol{U}},\overline{V})=0.$

Define the following matrices

(2.6)
$\displaystyle\boldsymbol{A}_{0}=D_{\boldsymbol{u}}\boldsymbol{f}(\overline{\boldsymbol{U}},\overline{V}),\quad\boldsymbol{B}_{0}=D_{v}\boldsymbol{f}(\overline{\boldsymbol{U}},\overline{V}),\quad\boldsymbol{C}_{0}=D_{\boldsymbol{u}}g(\overline{\boldsymbol{U}},\overline{V}),\quad
d_{0}=D_{v}g(\overline{\boldsymbol{U}},\overline{V})$

and assume that

(2.7) $\displaystyle\det\boldsymbol{A}_{0}\neq
0\quad\text{and}\quad\frac{1}{\det\boldsymbol{A}_{0}}\det\begin{pmatrix}\boldsymbol{A}_{0}&\boldsymbol{B}_{0}\\\
\boldsymbol{C}_{0}&d_{0}\end{pmatrix}=\mu_{k}>0,$

where $\mu_{k}$ is one of the eigenvalues of $-\Delta_{\nu}$. Then, there
exists a sequence of real numbers $d_{\ell}\to 1$ such that the following
“perturbed” problem

(2.8) $\displaystyle\boldsymbol{f}(\boldsymbol{U},V)=0,$ $\displaystyle
x\in\overline{\Omega},$ $\displaystyle
d_{\ell}\Delta_{\nu}V+(1-d_{\ell})(V-\overline{V})+g(\boldsymbol{U},V)=0,$
$\displaystyle x\in\Omega$

has a non-constant regular solution.

######  Remark 2.6.

Let us illustrate Proposition 2.5 in the case of the linear problem

(2.9) $\displaystyle\boldsymbol{A}_{0}\boldsymbol{U}+\boldsymbol{B}_{0}V$
$\displaystyle=0,$ $\displaystyle x\in\overline{\Omega},$
$\displaystyle\Delta_{\nu}V+\boldsymbol{C}_{0}\boldsymbol{U}+d_{0}V$
$\displaystyle=0,$ $\displaystyle x\in\Omega,$

with arbitrary constant coefficient matrices (not necessary as those in (2.6))

(2.10) $\displaystyle\boldsymbol{A}_{0}$
$\displaystyle=\begin{pmatrix}a_{11}&\dots&a_{1n}\\\ \vdots&\ddots&\vdots\\\
a_{n1}&\dots&a_{nn}\end{pmatrix},$ $\displaystyle\boldsymbol{B}_{0}$
$\displaystyle=\begin{pmatrix}b_{1}\\\ \vdots\\\ b_{n}\\\ \end{pmatrix},$
$\displaystyle\boldsymbol{C}_{0}$
$\displaystyle=\begin{pmatrix}c_{1}&\dots&c_{n}\end{pmatrix},$ $\displaystyle
d_{0}$ $\displaystyle=d.$

Under the assumption

(2.11) $\displaystyle\det\boldsymbol{A}_{0}\neq
0\quad\text{and}\quad-\boldsymbol{C}_{0}\boldsymbol{A}^{-1}_{0}\boldsymbol{B}_{0}+d_{0}=\mu_{k}>0$

problem (2.9) has non-constant regular stationary solutions of the form

(2.12) $\displaystyle\begin{pmatrix}\boldsymbol{U}\\\
V\end{pmatrix}=\begin{pmatrix}-\boldsymbol{A}_{0}^{-1}\boldsymbol{B}_{0}\Phi_{k}\\\
\Phi_{k}\end{pmatrix},$

where $\Phi_{k}$ is an eigenfunction of $-\Delta_{\nu}$ corresponding to the
eigenvalue $\mu_{k}$. By relation (3.17) below, assumptions (2.11) take the
form

(2.13) $\displaystyle\det\boldsymbol{A}_{0}\neq
0\quad\text{and}\quad\frac{1}{\det\boldsymbol{A}_{0}}\det\begin{pmatrix}\boldsymbol{A}_{0}&\boldsymbol{B}_{0}\\\
\boldsymbol{C}_{0}&d_{0}\end{pmatrix}=\mu_{k}>0.$

Proposition 2.5 shows how to extend this construction to nonlinear systems by
using a bifurcation argument.

In order to formulate our instability results, we associate with a matrix

(2.14) $\displaystyle\boldsymbol{A}$
$\displaystyle=\boldsymbol{A}(x)=\begin{pmatrix}a_{11}(x)&\dots&a_{1n}(x)\\\
\vdots&\ddots&\vdots\\\ a_{n1}(x)&\dots&a_{nn}(x)\end{pmatrix},$

where

(2.15) $\displaystyle a_{ij}\,\in
L^{\infty}(\Omega),\quad\text{for}\;i,j\in\\{1,\cdots,n\\},$

the corresponding multiplication operator on ${L(\Omega)}^{n}$, with
$p\in[1,\,\infty)$, given by

(2.16)
$\displaystyle\boldsymbol{\varphi}\in{L(\Omega)}^{n}\mapsto\boldsymbol{A}\boldsymbol{\varphi}\in{L(\Omega)}^{n}.$

We denote by $\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}$ the spectrum of this
mapping and by $s\big{(}\boldsymbol{A}(\cdot)\big{)}$ the spectral bound.
Below in Lemma 4.4, we show that in the case of $a_{i,j}\in
C(\overline{\Omega})$ we have

(2.17) $\displaystyle s\big{(}\boldsymbol{A}(\cdot)\big{)}=\sup\Big{\\{}{\rm
Re\,}\lambda(x):\,\text{where }\lambda(x)\text{ is an eigenvalue of
}\boldsymbol{A}(x)\text{ for some }x\in\overline{\Omega}\Big{\\}}$

We are in a position to state main results of this work. We prove in Theorems
2.7 and 2.8 below that all (i.e. not only those from Proposition 2.5) regular
stationary solutions to problem (1.1)-(LABEL:ini) are unstable except one
degenerate case (see assumption (2.20)) which we discuss below in Remark 2.10.

######  Theorem 2.7.

Let $N\geqslant 1$ and $\Omega\subset\mathbb{R}^{N}$ be bounded and open with
a smooth boundary. Let $(\boldsymbol{U},V)$ be a weak regular stationary
solution to problem (1.1)-(LABEL:ini). If

(2.18) $\displaystyle
s\big{(}f_{\boldsymbol{u}}(\boldsymbol{U}(\cdot),\,V(\cdot))\big{)}>0$

then $(\boldsymbol{U},V)$ is unstable in $C(\overline{\Omega})^{n+1}$.

######  Theorem 2.8.

Let $N\geqslant 1$ and $\Omega\subset\mathbb{R}^{N}$ be bounded, open, convex
and with a smooth boundary. Let
$(\boldsymbol{U},V)=\big{(}\boldsymbol{U}(x),V(x)\big{)}$ be a non-constant
regular stationary solution of problem (1.1)-(LABEL:ini) such that

(2.19) $\displaystyle
s\big{(}f_{\boldsymbol{u}}(\boldsymbol{U}(\cdot),\,V(\cdot))\big{)}\leqslant
0$

and

(2.20) $\displaystyle\det
f_{\boldsymbol{u}}\big{(}\boldsymbol{U}(x),\,V(x)\big{)}\neq 0\quad\text{for
all}\quad x\in\overline{\Omega}.$

Then $(\boldsymbol{U},V)$ is unstable in $C(\overline{\Omega})^{n+1}$.

######  Remark 2.9.

Let us illustrate both instability Theorems 2.7 and 2.8 by applying them to
the constant coefficients linear problem

(2.21) $\displaystyle\boldsymbol{u}_{t}$
$\displaystyle=\boldsymbol{A}_{0}\boldsymbol{u}+\boldsymbol{B}_{0}v,$
$\displaystyle x\in\overline{\Omega},\quad t>0,$ $\displaystyle v_{t}$
$\displaystyle=\Delta_{\nu}v+\boldsymbol{C}_{0}\boldsymbol{u}+d_{0}v,$
$\displaystyle x\in\Omega,\quad t>0,$

which under condition (2.13) has a non-constant stationary solution (2.12).
This solution is unstable. Indeed, the second expression in (2.13) is the
quotient of $n+1$ eigenvalues of the matrix
$\bigl{(}\begin{smallmatrix}\boldsymbol{A}_{0}&\boldsymbol{B}_{0}\\\
\boldsymbol{C}_{0}&d_{0}\end{smallmatrix}\bigr{)}$ and of $n$ eigenvalues of
the matrix $\boldsymbol{A}_{0}$. Obviously, all these eigenvalues are either
real or pairwise conjugate. Thus, the second inequality in (2.13) implies that
there exists at least one real eigenvalue $\lambda>0$. If $\lambda$ is an
eigenvalue of matrix $\boldsymbol{A}_{0}$ then solution (2.12) is unstable by
Theorem 2.7. Otherwise, $\lambda$ has to be an eigenvalue of matrix
$\bigl{(}\begin{smallmatrix}\boldsymbol{A}_{0}&\boldsymbol{B}_{0}\\\
\boldsymbol{C}_{0}&d_{0}\end{smallmatrix}\bigr{)}$ and solution (2.12) is
unstable by Theorem 2.8, provided $\Omega$ is convex. In fact, this second
instability result can be also obtained directly from the fact that zero is an
unstable solution of system (2.9) without the diffusion and the convexity of
the domain is not needed.

######  Remark 2.10.

Condition (2.20) has been imposed because some particular versions of problem
(1.1)-(LABEL:ini) may have smooth stationary solutions which are either stable
or unstable if $\det
f_{\boldsymbol{u}}\big{(}\boldsymbol{U}(x),\,V(x)\big{)}=0$ for some
$x\in\overline{\Omega}$. As the simplest example, for $\boldsymbol{f}\equiv
0$, $g\equiv 0$, this problem has a solution of the form
$(\boldsymbol{U},V)=(\boldsymbol{U},C)$ where
$\boldsymbol{U}=\boldsymbol{U}(x)$ is an arbitrary smooth vector field and
$C\in\mathbb{R}$ is a constant. Obviously this is a stable solution, but not
regular in the sense of Definition 2.2. As another example, we consider the
problem

(2.22) $\displaystyle u_{t}=0,$ $\displaystyle x\in\overline{\Omega},\quad
t>0,$ $\displaystyle v_{t}=\Delta_{\nu}v+(\mu_{1}+1)u-v,$ $\displaystyle
x\in\Omega,\quad t>0,$

which has the regular stationary solution $u=v=\Phi_{1}$ (recall that
$\mu_{1}>0$ is the eigenvalue of $-\Delta_{\nu}$ with the eigenfunction
$\Phi_{1}$) which is stable but not not asymptotically stable. On the other
hand, the following modified problem

$\displaystyle u_{t}=0,$ $\displaystyle x\in\overline{\Omega},\quad t>0,$
$\displaystyle v_{t}=\Delta_{\nu}v+(\mu_{1}-1)u+v,$ $\displaystyle
x\in\Omega,\quad t>0,$

has the same regular stationary solution $u=v=\Phi_{1}$ which is now unstable.

######  Remark 2.11.

For one ODE coupled with one PDE, assumption (2.18) reduces to the following
inequality
$f_{\boldsymbol{u}}\big{(}\boldsymbol{U}(x_{0}),\,V(x_{0})\big{)}>0$, for some
$x_{0}\in\overline{\Omega}$, which is called as an autocatalysis condition in
the work [21]. In this case, an instability of all regular stationary
solutions has been shown in [21, Thm. 2.1].

##  3\. Existence of regular stationary solutions

First, we prove a general result on the existence of non-constant solutions to
Neumann boundary value problem (2.4). This is rather standard result obtained
for example in [33, Sec. 3] or [32, Sec. 2] in the case of simple eigenvalues.
We recall one possible construction of such solutions using variational
methods in the following lemma which is a generalization of [30, Theorem
11.32].

######  Lemma 3.1.

Let $N\leq 6$. Assume that $r\in C_{b}^{2}(\mathbb{R})$ satisfies
$r(0)=r^{\prime}(0)=0$. There exists a sequence of numbers $d_{\ell}\to 1$ and
a sequence of non-constant functions $v_{\ell}\in W^{1,2}(\Omega)$ such that
$\|v_{\ell}\|_{W^{1,2}}\to 0$ and which satisfy the boundary value problem

(3.1) $\displaystyle
d_{\ell}\Delta_{\nu}v_{\ell}+(\mu_{k}+1-d_{\ell})v_{\ell}+r(v_{\ell})$
$\displaystyle=0\qquad\text{for}\quad x\in\Omega.$

###### Proof.

We prove this lemma by using the Rabinowitz Bifurcation Theorem [29, Thm. 0.2]
for equations with a variational formulation. It is assumed in that approach
that

  * •

$E$ is a real Hilbert space,

  * •

$I\in C^{2}(E,\mathbb{R})$ with $I^{\prime}(u)=Lu+H(u)$,

  * •

$L$ is linear and $H(u)=o(\|u\|)$ at $u=0$,

  * •

$\mu$ is an isolated eigenvalue of $L$ of a finite multiplicity.

Under these assumptions, by [29, Thm. 0.2], the couple
$(\mu,0)\in\mathbb{R}\times E$ is a bifurcation point of the equation

(3.2) $\displaystyle G(\lambda,v)\equiv Lv+H(v)-\lambda v=0$

which means that each neighborhood of $(\mu,0)$ contains a solution
$(\lambda,v)$ with $\|v\|\neq 0$ of equation (3.2). In our case, we use the
the usual Sobolev space $E=W^{1,2}(\Omega)$ with the equivalent scalar product

(3.3) $\displaystyle\langle u,v\rangle_{W^{1,2}(\Omega)}=\int_{\Omega}\nabla
u\cdot\nabla v\,\text{d}x+\int_{\Omega}uv\,\text{d}x$

as well as the functional

(3.4) $\displaystyle
I(v)=\frac{\mu_{k}+1}{2}\int_{\Omega}v^{2}\,\text{d}x+\int_{\Omega}R(v)\,\text{d}x$

with $R(v)=\int_{0}^{v}r(s)\,ds$. It is an elementary calculation to show that

  * •

$I\in C\big{(}W^{1,2}(\Omega),\mathbb{R}\big{)}$,

  * •

it is differentiable in the Fréchet sense and for each $v\in W^{1,2}(\Omega)$

(3.5) $\displaystyle
DI(v)\varphi=(\mu_{k}+1)\int_{\Omega}v\varphi\,\text{d}x+\int_{\Omega}r(v)\varphi\,\text{d}x,$

  * •

$DI\in
C\Big{(}W^{1,2}(\Omega),\;$Lin$\big{(}W^{1,2}(\Omega),\mathbb{R}\big{)}\Big{)}$.

The second Fréchet derivative at the point $v\in W^{1,2}(\Omega)$ is
represented by the bilinear form

(3.6)
$\displaystyle\big{\langle}D^{2}I(v)\varphi,\psi\big{\rangle}=(\mu_{k}+1)\int_{\Omega}\varphi\psi\,\text{d}x+\int_{\Omega}r^{\prime}(v)\varphi\psi\,\text{d}x.$

Let us show that $D^{2}I\in C\Big{(}W^{1,2}(\Omega),$
Lin$\big{(}W^{1,2}(\Omega),$ Lin$(W^{1,2}(\Omega),\mathbb{R})\big{)}\Big{)}$.
For $v_{n}\to v$ in $W^{1,2}(\Omega)$ and $\psi,\varphi\in W^{1,2}(\Omega)$ we
estimate

(3.7)
$\displaystyle\left|\big{\langle}\big{(}D^{2}I(v_{n})-D^{2}I(v)\big{)}\varphi,\psi\big{\rangle}\right|$
$\displaystyle\leqslant\int_{\Omega}|r^{\prime}(v_{n})-r^{\prime}(v)||\varphi||\psi|\,\text{d}x$
$\displaystyle\leqslant\|r^{\prime\prime}\|_{\infty}\int_{\Omega}|v_{n}-v||\varphi||\psi|\,\text{d}x$
$\displaystyle\leqslant\|r^{\prime\prime}\|_{\infty}\|v_{n}-v\|_{3}\|\varphi\|_{3}\|\psi||_{3}$
$\displaystyle\leqslant\|r^{\prime\prime}\|_{\infty}\|v_{n}-v\|_{W^{1,2}}\|\varphi\|_{W^{1,2}}||\psi||_{W^{1,2}}.$

The last inequality follows from the Sobolev embeddings with the assumption
$N\leqslant 6$.

In particular, for each test function $\varphi\in W^{1,2}(\Omega)$, we obtain

(3.8) $\displaystyle
I^{\prime}(v)(\varphi)=(\mu_{k}+1)\int_{\Omega}v\varphi\,\text{d}x+\int_{\Omega}r(v)\varphi\,\text{d}x\equiv
Lv(\varphi)+H(v)(\varphi),$

where, by the assumption on $r=r(v)$, we obtain immediately that
$H(v)=o\big{(}\|v\|_{W^{1,2}}\big{)}$ as $\|v\|_{W^{1,2}}\to 0$.

Notice now that $\lambda=1$ is an isolated eigenvalue of finite multiplicity
of the operator $L$. Indeed, this claim is equivalent to the equality

(3.9) $\displaystyle Lv(\varphi)=\langle
v,\varphi\rangle_{W^{1,2}(\Omega)}\qquad\text{for all}\quad\varphi\in
W^{1,2}(\Omega),$

that is, to the equation

(3.10) $\displaystyle(\mu_{k}+1)\int_{\Omega}v\varphi\,dx=\int_{\Omega}\nabla
v\cdot\nabla\varphi\,\text{d}x+\int_{\Omega}v\varphi\,\text{d}x\qquad\text{for
all}\quad\varphi\in W^{1,2}(\Omega),$

which obviously reduces to the eigenvalue problem for $\Delta_{\nu}$. Now, we
apply the fact that $\mu_{k}$ is an isolated eigenvalue of finite
multiplicity.

Thus, by the Rabinowitz Theorem [29, Thm. 0.2], the couple $(1,0)$ is a
bifurcation point of equation (3.2) which means that there exists a sequence
of numbers $d_{\ell}\to 1$ and nonzero $\\{v_{\ell}\\}\subset W^{1,2}(\Omega)$
such that $\|v_{\ell}\|_{W^{1,2}}\to 0$, satisfying the following equation

(3.11) $\displaystyle Lv_{\ell}(\varphi)+H(v_{\ell})(\varphi)-d_{\ell}\langle
v,\varphi\rangle_{1,2}=0\qquad\text{for all}\quad\varphi\in W^{1,2}(\Omega)$

which, in our setting, is equivalent to the equation satisfied by the weak
solutions $v_{\ell}\in W^{1,2}(\Omega)$ to problem (3.1):

(3.12) $\displaystyle-d_{\ell}\int_{\Omega}\nabla
v_{\ell}\cdot\nabla\varphi\,\text{d}x+(\mu_{k}+1-d_{\ell})\int_{\Omega}v_{\ell}\varphi\,\text{d}x+\int_{\Omega}r(v_{\ell})\varphi\,\text{d}x=0$

for all $\varphi\in W^{1,2}(\Omega)$. ∎

###### Proof of Proposition 2.5..

We apply Lemma 3.1 to construct non-constant solutions to the reaction-
diffusion-ODE system. In the following, we denote by
$B_{\varepsilon}(\overline{V})$ an open ball of radius $\varepsilon>0$ centred
at $\overline{V}$. First let us prove that, in Lemma 3.1, only a finite number
of $v_{\ell}$ can be constant. Indeed, if there exists a subsequence of
constants $\\{v_{\ell_{n}}\\}$ satisfying equation (3.1) such that
$v_{\ell_{n}}\to 0$ then passing to the limit we obtain that
$r^{\prime}(0)=-\mu_{k}$ which is a contradiction.

Next, since $\det
D_{\boldsymbol{u}}\boldsymbol{f}(\overline{\boldsymbol{U}},\overline{V})\neq
0$, by Implicit Function Theorem, there exists $\varepsilon>0$ and a function
$\boldsymbol{k}\in C^{2}\big{(}B_{\varepsilon}(\overline{V})\big{)}$ such that
$\boldsymbol{k}(\overline{V})=\overline{\boldsymbol{U}}$ and
$\boldsymbol{f}(\boldsymbol{k}(v),v)=0$ for all
$v\in\big{(}B_{\varepsilon}(\overline{V})\big{)}$. First, we show that the
function $h(v)\equiv g\big{(}\boldsymbol{k}(v),v\big{)}$, defined for all
$v\in\big{(}B_{\varepsilon}(\overline{V})\big{)}$, satisfies the relations

(3.13) $\displaystyle h(\overline{V})=0\qquad\text{and}\qquad
h^{\prime}(\overline{V})=\gamma\mu_{k}>0.$

Indeed, the first one results from the equations
$h(\overline{V})=g(\boldsymbol{k}(\overline{V}),\overline{V})=g(\overline{\boldsymbol{U}},\overline{V})=0$.
For the second one, differentiating the function
$h(v)=g\big{(}\boldsymbol{k}(v),v\big{)}$ yields

(3.14) $\displaystyle
h^{\prime}(v)=D_{\boldsymbol{u}}g(\boldsymbol{k}(v),v)\boldsymbol{k}^{\prime}(v)+D_{v}g(\boldsymbol{k}(v),v).$

On the other hand, we differentiate both sides of the equation
$\boldsymbol{f}\big{(}\boldsymbol{k}(v),v\big{)}=0$ with respect to $v$ to
obtain the equality
$D_{\boldsymbol{u}}\boldsymbol{f}\big{(}\boldsymbol{k}(v),v\big{)}\boldsymbol{k}^{\prime}(v)+D_{v}\boldsymbol{f}\big{(}\boldsymbol{k}(v),v\big{)}=0$,
or, equivalently,

(3.15)
$\displaystyle\boldsymbol{k}^{\prime}(v)=-D^{-1}_{\boldsymbol{u}}\boldsymbol{f}(\boldsymbol{k}(v),v)D_{v}\boldsymbol{f}(\boldsymbol{k}(v),v).$

Finally, substituting equation (3.15) into equation (3.14), choosing
$v=\overline{V}$, and applying the notation from the statement of Theorem 2.5
we obtain the equations

(3.16) $\displaystyle h^{\prime}(\overline{V})$
$\displaystyle=-D_{\boldsymbol{u}}g(k(\overline{V}),\overline{V})D^{-1}_{\boldsymbol{u}}\boldsymbol{f}(\boldsymbol{k}(\overline{V}),\overline{V})D_{v}\boldsymbol{f}(\boldsymbol{k}(\overline{V}),\overline{V})+D_{v}g(k(\overline{V}),\overline{V})$
$\displaystyle=-\boldsymbol{C}_{0}\boldsymbol{A}^{-1}_{0}\boldsymbol{B}_{0}+d_{0}.$

Notice that the right-hand side satisfies the following relation

(3.17) $\displaystyle
d_{0}-\boldsymbol{C}_{0}\boldsymbol{A}^{-1}_{0}\boldsymbol{B}_{0}$
$\displaystyle=\det\begin{pmatrix}\boldsymbol{I}&\boldsymbol{0}\\\
\boldsymbol{0}&d_{0}-\boldsymbol{C}_{0}\boldsymbol{A}_{0}^{-1}\boldsymbol{B}_{0}\end{pmatrix}$
$\displaystyle=\frac{1}{\det\boldsymbol{A}_{0}}\det\begin{pmatrix}\boldsymbol{A}_{0}&\boldsymbol{0}\\\
\boldsymbol{C}_{0}&d_{0}-\boldsymbol{C}_{0}\boldsymbol{A}_{0}^{-1}\boldsymbol{B}_{0}\end{pmatrix}$
$\displaystyle=\frac{1}{\det\boldsymbol{A}_{0}}\det\left(\begin{pmatrix}\boldsymbol{A}_{0}&\boldsymbol{B}_{0}\\\
\boldsymbol{C}_{0}&d_{0}\end{pmatrix}\begin{pmatrix}\boldsymbol{I}&-\boldsymbol{A}_{0}^{-1}\boldsymbol{B}_{0}\\\
\textbf{0}&1\end{pmatrix}\right)$
$\displaystyle=\frac{1}{\det\boldsymbol{A}_{0}}\det\begin{pmatrix}\boldsymbol{A}_{0}&\boldsymbol{B}_{0}\\\
\boldsymbol{C}_{0}&d_{0}\end{pmatrix}.$

Hence, $h^{\prime}(\overline{V})=\mu_{k}>0$.

Next, we denote by $\bar{h}$ an arbitrary extension of the function $h$ to the
whole line $\mathbb{R}$ which satisfies

(3.18) $\displaystyle\bar{h}\in
C_{b}^{2}(\mathbb{R})\quad\text{and}\quad\bar{h}(v)=h(v)\quad\text{for
all}\quad v\in B_{\varepsilon}(\overline{V}).$

By Lemma 3.1, we obtain a sequence $d_{\ell}\to 1$ such that the Neumann
boundary value problem

(3.19) $\displaystyle
d_{\ell}\Delta_{\nu}V_{\ell}+(1-d_{\ell})(V_{\ell}-\overline{V})+\bar{h}(V_{\ell})$
$\displaystyle=0\quad\text{for}\quad x\in\Omega$

has a non-constant solution $V_{\ell}\in W^{1,2}(\Omega)$. Indeed, it suffices
to look for these solutions in the form $V_{\ell}=\overline{V}+z_{\ell}$,
where $z_{\ell}$ solves the following problem

(3.20) $\displaystyle
d_{\ell}\Delta_{\nu}z_{\ell}+\big{(}\bar{h}^{\prime}(\overline{V})+1-d_{\ell}\big{)}z_{\ell}+r(z_{\ell})$
$\displaystyle=0\quad\text{in}\quad\Omega,$

with $\bar{h}^{\prime}(\overline{V})=\mu_{k}$ and
$r(z)=\bar{h}(\overline{V}+z)-\bar{h}^{\prime}(\overline{V})z$ satisfying
$r\in C_{b}^{2}(\mathbb{R})$, $r(0)=0$, and $r^{\prime}(0)=0$.

Next, by Lemma 3.1, the solutions $z_{\ell}$ to problem (3.20) satisfy
$\|z_{\ell}\|_{W^{1,2}(\Omega)}\to 0$, hence, by the standard elliptic theory,
we also have $\|z_{\ell}\|_{W^{2,2}(\Omega)}\to 0$. By the bootstrap arguments
utilizing the elliptic $L^{p}$-estimates and the Sobolev embedding theorem, we
conclude that $\|z_{\ell}\|_{W^{2,q}(\Omega)}\to 0$ for $q>N/2$ and hence
$\|z_{\ell}\|_{L^{\infty}(\Omega)}\to 0$. In particular, if
$\|z_{\ell}\|_{\infty}\leq\varepsilon$, by (3.18), we obtain

(3.21)
$\displaystyle\bar{h}(V_{\ell})=\bar{h}(\overline{V}+z_{\ell})=h(\overline{V}+z_{\ell})=h(V_{\ell}).$

Thus, $V_{\ell}=\overline{V}+z_{\ell}$ is a nontrivial solution of problem
(3.19) with $\bar{h}$ replaced by $h$.

Finally, we define $U_{\ell}=\boldsymbol{k}(V_{\ell})$ in order to obtain a
nontrivial regular solution of problem (2.8). ∎

##  4\. Linear equation and spectral analysis

###  4.1. Linearised operator

We analyse the stability of stationary solution to problem (1.1)-(LABEL:ini)
via the usual linearisation procedure. In order to apply that approach in this
section, we discuss stability properties of the zero solution to the following
linear reaction diffusion-ODE system

(4.1) $\displaystyle\boldsymbol{\varphi}_{t}$
$\displaystyle=\boldsymbol{A}\boldsymbol{\varphi}+\boldsymbol{B}\psi,$
$\displaystyle x\in\overline{\Omega},\quad t>0,$ $\displaystyle\psi_{t}$
$\displaystyle=\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+d\psi,$
$\displaystyle x\in\Omega,\quad t>0,$

with matrices

(4.2) $\displaystyle\boldsymbol{A}$
$\displaystyle=\boldsymbol{A}(x)=\begin{pmatrix}a_{11}(x)&\dots&a_{1n}(x)\\\
\vdots&\ddots&\vdots\\\ a_{n1}(x)&\dots&a_{nn}(x)\end{pmatrix},$
$\displaystyle\boldsymbol{B}$
$\displaystyle=\boldsymbol{B}(x)=\begin{pmatrix}b_{1}(x)\\\ \vdots\\\
b_{n}(x)\\\ \end{pmatrix},$ $\displaystyle\boldsymbol{C}$
$\displaystyle=\boldsymbol{C}(x)=\begin{pmatrix}c_{1}(x)&\dots&c_{n}(x)\end{pmatrix},$
$\displaystyle d$ $\displaystyle=d(x),$

where

(4.3) $\displaystyle a_{ij},\;b_{i},\;c_{i},\;d\in
L^{\infty}(\Omega),\quad\text{for}\;i,j\in\\{1,\cdots,n\\}.$

In the following, we study properties of the operator
$\big{(}\mathcal{L}_{p},D(\mathcal{L}_{p})\big{)}$ defined by the formula

(4.4) $\displaystyle\mathcal{L}_{p}\begin{pmatrix}\boldsymbol{\varphi}\\\
\psi\end{pmatrix}\equiv\begin{pmatrix}\boldsymbol{A}\boldsymbol{\varphi}+\boldsymbol{B}\psi\\\
\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+d\psi\end{pmatrix}\quad\text{with}\quad
D(\mathcal{L}_{p})=L^{p}(\Omega)^{n}\times W^{2,p}_{\nu}(\Omega)$

for each $p\in(1,\,\infty)$, where

(4.5) $\displaystyle W^{2,p}_{\nu}(\Omega)=\\{v\in
W^{2,p}(\Omega):\,\partial_{\nu}v=0\text{ on }\partial\Omega\\}$

is the Sobolev space $W^{2,p}(\Omega)$ supplemented with the Neumann boundary
condition (see e.g. [40, Ch. 2, sec 2.4]).

######  Proposition 4.1.

For each $p\in(1,\infty)$ the operator
$\big{(}\mathcal{L}_{p},D(\mathcal{L}_{p})\big{)}$ generates an analytic
semigroup of linear operators on $L^{p}(\Omega)^{n+1}$.

###### Proof.

This fact is well-known in the case of the Laplace operator with the Neumann
boundary condition, see e.g. [5, Ch. II, Sec. 5]. Here, we consider bounded
perturbations of this operator. ∎

In order to study a nonlinear stability, we need a counterpart of the
semigroup from Proposition 4.1 acting on a space of bounded continuous
functions. First, we recall that the operator
$\big{(}\Delta_{\nu},D(\Delta_{\nu})\big{)}$ with the domain

(4.6) $\displaystyle D(\Delta_{\nu})=\left\\{u\in\bigcap_{p\geqslant
1}W^{2,p}_{loc}(\Omega):\quad u,\,\Delta u\in
C(\overline{\Omega}),\quad\partial_{\nu}u=0\text{ on }\partial\Omega\right\\}$

is a sectorial operator on the space $C(\overline{\Omega})$, see [17,
Corollary 3.1.24].

######  Proposition 4.2.

If $a_{ij},\;b_{i},\;c_{i},\;d\in C(\overline{\Omega})$ for
$i,j\in\\{1,\cdots,n\\}$ then the operator given by the formula

(4.7) $\displaystyle\mathcal{L}_{\infty}\begin{pmatrix}\boldsymbol{\varphi}\\\
\psi\end{pmatrix}\equiv\begin{pmatrix}\boldsymbol{A}\boldsymbol{\varphi}+\boldsymbol{B}\psi\\\
\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+d\psi\end{pmatrix}\quad\text{with}\quad
D(\mathcal{L}_{\infty})=C(\overline{\Omega})^{n}\times D(\Delta_{\nu})$

generates an analytic semigroup of linear operators on
$C(\overline{\Omega})^{n+1}$.

One proves this proposition in the same way as Proposition 4.1.

###  4.2. Spectra of multiplication operators

As usual, we denote by $\sigma(L)$ the spectrum of a linear operator
$\big{(}L,D(L)\big{)}$ and $s(L)=\sup\\{{\rm
Re\,}\lambda:\lambda\in\sigma(L)\\}$ its spectral bound. In particular, for a
square matrix $\boldsymbol{M}$ with constant elements

(4.8) $\displaystyle s(\boldsymbol{M})=\max\\{{\rm
Re\,}\lambda:\lambda\in\mathbb{C}\text{ is an eigenvalue of
}\boldsymbol{M}\\}.$

Now, we discuss spectra of matrix multiplication operators induced by matrices
with $x$-dependent and bounded coefficients. The following results can be
found either in the work [4] or in the PhD dissertation [14] and we recall
them together with the proofs in slightly less general versions.

First we show that the spectrum of the multiplication operator consist only of
approximative eigenvalues.

######  Lemma 4.3.

For $\boldsymbol{M}(x)=\big{(}m_{i,j}(x)\big{)}_{i,j=1}^{n}$ with $m_{i,j}\in
L^{\infty}(\Omega)$, we define the multiplication operator on
$L^{p}(\Omega)^{n}$ with $p\in[1,\infty)$ by the formula
$\big{(}\boldsymbol{M}(\cdot)\boldsymbol{\xi}\big{)}(x)=\boldsymbol{M}(x)\boldsymbol{\xi}(x)$.
Then its spectrum has the form

(4.9)
$\displaystyle\sigma\big{(}\boldsymbol{M}(\cdot)\big{)}=\left\\{\lambda\in\mathbb{C}:\begin{array}[]{l}\forall{\varepsilon>0}\quad\exists{\Omega_{\varepsilon}\subset\Omega}\
open\ set\quad\exists{\boldsymbol{\xi}_{\varepsilon}\in\mathbb{R}^{n}}\\\
such\ that\
\|\boldsymbol{M}(\cdot)\boldsymbol{\xi}_{\varepsilon}-\lambda\boldsymbol{\xi}_{\varepsilon}\|_{L^{\infty}(\Omega_{\varepsilon})}\leqslant\varepsilon|\boldsymbol{\xi}_{\varepsilon}|\end{array}\right\\}.$

###### Proof.

Denote by $\boldsymbol{M}_{ess}$ the set on the right-hand side of equality
(4.9) and let $\lambda\in\boldsymbol{M}_{ess}$. Suppose that
$\lambda\notin\sigma(\boldsymbol{M}(\cdot))$. Hence, there a exists a constant
$C>0$ such that

(4.10) $\displaystyle\|\boldsymbol{\xi}\|_{p}\leqslant
C\|(\boldsymbol{M}(\cdot)-\lambda\boldsymbol{I})\boldsymbol{\xi}\|_{L(\Omega)}\quad\text{for
all}\quad\boldsymbol{\xi}\in L^{p}(\Omega)^{n}.$

Now by the definition of $\boldsymbol{M}_{ess}$, for each $\varepsilon>0$
there exist open set $\Omega_{\varepsilon}\subset\Omega$ and
$\boldsymbol{\xi}_{\varepsilon}$ such that
$\|\boldsymbol{M}(\cdot)\boldsymbol{\xi}_{\varepsilon}-\lambda\boldsymbol{\xi}_{\varepsilon}\|_{L^{\infty}(\Omega_{\varepsilon})}\leqslant\varepsilon|\boldsymbol{\xi}_{\varepsilon}|$.
Let
$\boldsymbol{\xi}=\frac{1}{|\Omega_{\varepsilon}|}\mathds{1}_{\Omega_{\varepsilon}}\boldsymbol{\xi}_{\varepsilon}$
in inequality (4.10). Then,

(4.11) $\displaystyle\|\boldsymbol{\xi}\|_{L^{p}(\Omega)}$
$\displaystyle\leqslant
C\|\big{(}\boldsymbol{M}(\cdot)-\lambda\boldsymbol{I}\big{)}\boldsymbol{\xi}\|_{L(\Omega)}$
$\displaystyle\leqslant
C\|\big{(}\boldsymbol{M}(\cdot)-\lambda\boldsymbol{I}\big{)}\|_{L^{\infty}(\Omega_{\varepsilon})}\|\boldsymbol{\xi}\|_{L^{p}(\Omega_{\varepsilon})}\leqslant
C\varepsilon\|\boldsymbol{\xi}\|_{L(\Omega)}.$

For sufficiently small $\varepsilon>0$, we have
$\|\boldsymbol{\xi}\|_{L(\Omega)}=0$ which is contradiction.

Conversely, let $\lambda\notin\boldsymbol{M}_{ess}$. By the definition of
$\boldsymbol{M}_{ess}$, there exists $\varepsilon>0$ such that, for all
$\Omega_{\varepsilon}\subset\Omega$ and for all vectors
$\boldsymbol{\xi}_{\varepsilon}\in\mathbb{R}^{n}$ with
$|\boldsymbol{\xi}_{\varepsilon}|=1$ we have
$\|\boldsymbol{M}(\cdot)\boldsymbol{\xi}_{\varepsilon}-\lambda\boldsymbol{\xi}_{\varepsilon}\|_{L^{\infty}(\Omega_{\varepsilon})}\geqslant\varepsilon$.
It implies that the inequality
$|\boldsymbol{M}(x)\boldsymbol{\xi}_{\varepsilon}-\lambda\boldsymbol{\xi}_{\varepsilon}|\geqslant\varepsilon|\boldsymbol{\xi}_{\varepsilon}|$
holds true for almost all $x\in\Omega$. Thus,

(4.12)
$\displaystyle\|\big{(}\boldsymbol{M}(\cdot)-\lambda\boldsymbol{I}\big{)}\boldsymbol{\xi}\|_{L(\Omega)}=\left(\int_{\Omega}\big{|}\big{(}\boldsymbol{M}(x)-\lambda\boldsymbol{I}\big{)}\boldsymbol{\xi}(x)\big{|}^{p}\,\text{d}x\right)^{\frac{1}{p}}\geqslant\varepsilon\|\boldsymbol{\xi}\|_{L(\Omega)}$

and the operator $\boldsymbol{M}(\cdot)-\lambda\boldsymbol{I}$ is invertible,
hence $\lambda\notin\sigma\big{(}\boldsymbol{M}(\cdot)\big{)}$. ∎

The following lemma explains a relation between
$\sigma\big{(}\boldsymbol{M}(\cdot)\big{)}$ and the eigenvalues of the
matrices $\boldsymbol{M}(x)$ for each $x\in\Omega$.

######  Lemma 4.4.

Assume that $\boldsymbol{M}(x)=\big{(}m_{i,j}(x)\big{)}_{i,j=1}^{n}$ with
$m_{i,j}\in L^{\infty}(\Omega)$. Assume that there exists a closed set
$\Omega^{\prime}\subset\Omega$ such that $|\Omega^{\prime}|=0$ and
$m_{i,j}\big{|}_{\Omega\setminus\Omega^{\prime}}$ is continuous. Then

(4.13)
$\displaystyle\sigma\big{(}\boldsymbol{M}(\cdot)\big{)}=\overline{\bigcup_{x\in\Omega\setminus\Omega^{\prime}}\sigma\big{(}\boldsymbol{M}(x)\big{)}}.$

###### Proof.

Let
$\lambda\in\overline{\bigcup_{x\in\Omega\setminus\Omega^{\prime}}\sigma\big{(}\boldsymbol{M}(x)\big{)}}$.
Hence, there exists a sequence
$\\{x_{n}\\}_{n=1}^{\infty}\subset\Omega\setminus\Omega^{\prime}$ such that
the matrices $\boldsymbol{M}(x_{n})$ have eigenvalues satisfying
$\lambda_{n}\to\lambda$. Denote by $\boldsymbol{\xi}_{\lambda_{n}}$ the
corresponding eigenvector, i.e.
$(\boldsymbol{M}(x)-\lambda_{n})\boldsymbol{\xi}_{\lambda_{n}}=0$. By the
assumption, matrix elements $m_{i,j}(x)$ are continuous functions on
$\Omega\setminus\Omega^{\prime}$. Since eigenvalues of the matrix
$\boldsymbol{M}(x)$ are also continues functions of $x$ (see Remark 4.5,
below), for each $\varepsilon>0$ there exists a ball $B_{\delta}(x_{n})$ such
that for all $y\in B_{\delta}(x_{n})$

(4.14)
$\displaystyle|(\boldsymbol{M}(y)-\lambda_{n})\boldsymbol{\xi}_{\lambda_{n}}|\leqslant\varepsilon|\boldsymbol{\xi}_{\lambda_{n}}|.$

We take $n$ sufficiently large and apply the triangle inequality to obtain

(4.15)
$\displaystyle|(\boldsymbol{M}(y)-\lambda)\boldsymbol{\xi}_{\lambda_{n}}|\leqslant|(\boldsymbol{M}(y)-\lambda_{n})\boldsymbol{\xi}_{\lambda_{n}}|+|(\lambda_{n}-\lambda)\boldsymbol{\xi}_{\lambda_{n}}|\leqslant
2\varepsilon|\boldsymbol{\xi}_{\lambda_{n}}|.$

Thus, by Lemma 4.3, we have $\lambda\in\sigma(\boldsymbol{M}(\cdot))$.

Conversely, assume that
$\lambda\notin\overline{\bigcup_{x\in\Omega\setminus\Omega^{\prime}}\sigma\big{(}\boldsymbol{M}(x)\big{)}}$.
Thus, for each $x\in\Omega\setminus\Omega^{\prime}$, the matrix
$\big{(}\boldsymbol{M}(x)-\lambda\boldsymbol{I}\big{)}$ is invertible and
there exists $\varepsilon_{x}>0$ such that for all vectors
$\boldsymbol{\xi}_{\varepsilon}\in\mathbb{R}^{n}$ we have

(4.16)
$\displaystyle\big{|}(\boldsymbol{M}(x)-\lambda\boldsymbol{I})\boldsymbol{\xi}_{\varepsilon}\big{|}\geqslant\varepsilon_{x}|\boldsymbol{\xi}_{\varepsilon}|.$

We show that there exists $\varepsilon>0$ such that
$\varepsilon_{x}\geqslant\varepsilon>0$ uniformly for all
$x\in\Omega\setminus\Omega^{\prime}$. Indeed, suppose a contrario that there
exists a sequence $x_{n}\in\Omega\setminus\Omega^{\prime}$ satisfying
$\varepsilon_{x_{n}}\to 0$. We take the sequence of vectors
$\\{\boldsymbol{\xi}_{\varepsilon_{x_{n}}}\\}_{n=1}^{\infty}\subset\mathbb{R}^{n}$,
with $|\boldsymbol{\xi}_{\varepsilon_{x_{n}}}|=1$ satisfying

(4.17)
$\displaystyle\big{|}(\boldsymbol{M}(x_{n})-\lambda\boldsymbol{I})\boldsymbol{\xi}_{\varepsilon_{x_{n}}}\big{|}=\varepsilon_{x_{n}}\to
0.$

Letting $n\to\infty$ we have that $\lambda$ is arbitrary close to an
eigenvalue of matrix $\boldsymbol{M}(x_{n})$, hence
$\lambda\in\overline{\bigcup_{x\in\Omega\setminus\Omega^{\prime}}\sigma\big{(}\boldsymbol{M}(x)\big{)}}$.
Therefore, there exists $\varepsilon>0$ such that for all
$\boldsymbol{\xi}_{\varepsilon}\in\mathbb{R}^{n}$ we have

(4.18)
$\displaystyle\big{|}(\boldsymbol{M}(x)-\lambda\boldsymbol{I})\boldsymbol{\xi}_{\varepsilon}\big{|}\geqslant\varepsilon|\boldsymbol{\xi}_{\varepsilon}|,\quad\text{for
all}\quad x\in\Omega\setminus\Omega^{\prime}.$

Lemma 4.3 implies again that $\lambda\notin\sigma(\boldsymbol{M}(\cdot))$. ∎

######  Remark 4.5.

Notice that eigenvalues of a matrix depend continuously on the matrix
elements, see eg. [7]. Such a continuous dependence on a matrix elements may
be false for eigenvectors, see eg. [11, Example 5.3 on p. 111].

###  4.3. Spectrum of the reaction-diffusion-ODE operator

######  Theorem 4.6.

Let $\big{(}\mathcal{L}_{p},D(\mathcal{L}_{p})\big{)}$, with $p\in(1,\infty)$,
be the operator given by formula (4.4) with matrices (4.2)-(4.3). Denote by
$\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}$ be the spectrum of multiplication
operator $\boldsymbol{A}(\cdot)$ on $L^{p}(\Omega)^{n}$. Then, the spectrum of
$\big{(}\mathcal{L}_{p},D(\mathcal{L}_{p})\big{)}$ can be characterized as

(4.19)
$\displaystyle\sigma(\mathcal{L}_{p})=\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}\cup\\{\text{eigenvalues
of }\;\mathcal{L}_{p}\\}.$

In particular, the spectrum $\sigma(\mathcal{L}_{p})$ is independent of $p$.

###### Proof.

We improve the reasoning form [21, Thm 4.5], where one ODE coupled with one
PDE were considered. First, we prove that
$\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}\subset\sigma(\mathcal{L}_{p})$ by
showing that for each $\lambda\in\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}$
the operator

(4.20)
$\displaystyle\mathcal{L}_{p}-\lambda\boldsymbol{I}:L^{p}(\Omega)^{n}\times{W^{2,p}_{\nu}(\Omega)}\rightarrow
L^{p}(\Omega)^{n}\times L^{p}(\Omega)$

defined by formula

(4.21)
$\displaystyle(\mathcal{L}_{p}-\lambda\boldsymbol{I})\begin{pmatrix}\boldsymbol{\varphi}\\\
\psi\end{pmatrix}=\begin{pmatrix}(\boldsymbol{A}-\lambda)\boldsymbol{\varphi}+\boldsymbol{B}\psi\\\
\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+(d-\lambda)\psi\end{pmatrix}$

cannot have a bounded inverse. Indeed, suppose a contrario that
$(\mathcal{L}_{p}-\lambda\boldsymbol{I})^{-1}$ exists and is a bounded
operator from $(L^{p}(\Omega))^{n+1}$ into $D(\mathcal{L}_{p})$ equipped with
the graph norm. The graph norm of $\mathcal{L}_{p}$ is equivalent to the norm
on the product space $(L^{p}(\Omega))^{n}\times({W^{2,p}_{\nu}(\Omega)})$ by
[40, Section 2.4]. Thus, there exists a number $K>0$ such that

(4.22)
$\displaystyle\|\boldsymbol{\varphi}\|_{L(\Omega)}+\|\psi\|_{W^{2,p}(\Omega)}\leqslant
K\big{(}$
$\displaystyle\|(\boldsymbol{A}-\lambda\boldsymbol{I})\boldsymbol{\varphi}+\boldsymbol{B}\psi\|_{L(\Omega)}$
$\displaystyle+$
$\displaystyle\|\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+(d-\lambda)\psi\|_{L(\Omega)}\big{)}$

for all $(\boldsymbol{\varphi},\psi)\in
L^{p}(\Omega)^{n}\times{W^{2,p}_{\nu}(\Omega)}$. Let $\varepsilon>0$ be
arbitrary (to be chosen later on). By Lemma 4.3 there exists
$\boldsymbol{\xi}\in\mathbb{R}^{n}$ with $|\boldsymbol{\xi}|=1$ and a ball
$B_{\varepsilon}$ such that

(4.23)
$\displaystyle\|(\boldsymbol{A}(\cdot)-\lambda\boldsymbol{I})\boldsymbol{\xi}\|_{L^{\infty}(B_{\varepsilon})}\leqslant\varepsilon.$

Let us show that, for arbitrary $\psi\in C_{c}^{\infty}(\Omega)$ with
$\text{supp}\;\psi\subset B_{\varepsilon}$, we may find
$\boldsymbol{\varphi}\in L^{p}(\Omega)^{n}$ with the following properties

  * •

$\boldsymbol{\varphi}=\rho\boldsymbol{\xi}$ for some $\rho\in{L(\Omega)}$ with
$\text{supp}\,\rho\subset B_{\varepsilon}$,

  * •

$\zeta=\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+(d-\lambda)\psi$
satisfies
$\|\zeta\|_{L(\Omega)}\leqslant\varepsilon\|\boldsymbol{\varphi}\|_{L(\Omega)}$.

Indeed, we cut the function $\boldsymbol{C}$ at the level $\varepsilon$ in the
following way

(4.24)
$\displaystyle\boldsymbol{C}^{\varepsilon}=\boldsymbol{C}^{\varepsilon}(x)\equiv\begin{cases}\boldsymbol{C}(x)&\text{if}\quad\big{|}\boldsymbol{C}(x)\boldsymbol{\xi}\big{|}>\varepsilon,\\\
(\varepsilon,\cdots,\varepsilon)&\text{if}\quad\big{|}\boldsymbol{C}(x)\boldsymbol{\xi}\big{|}\leqslant\varepsilon.\end{cases}$

Thus we obtain

(4.25)
$\displaystyle\zeta=\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+(d-\lambda)\psi=\Delta_{\nu}\psi+\boldsymbol{C}^{\varepsilon}\boldsymbol{\varphi}+(d-\lambda)\psi+(\boldsymbol{C}-\boldsymbol{C}^{\varepsilon})\boldsymbol{\varphi}$

for

(4.26)
$\displaystyle\boldsymbol{\varphi}=\frac{-(\Delta_{\nu}\psi+(d-\lambda)\psi)}{\boldsymbol{C}^{\varepsilon}\boldsymbol{\xi}}\boldsymbol{\xi}=\rho\boldsymbol{\xi}\in
L^{p}(\Omega)^{n}\quad\text{and}\quad\zeta=(\boldsymbol{C}-\boldsymbol{C}^{\varepsilon})\boldsymbol{\varphi}\in
L^{p}(\Omega)$

with
$\|(\boldsymbol{C}-\boldsymbol{C}^{\varepsilon})\boldsymbol{\xi}\|_{L^{\infty}(\Omega)}\leqslant
n\varepsilon$.

Now, noting that supp $\boldsymbol{\varphi}\subset B_{\varepsilon}$ and
$\boldsymbol{\varphi}(x)=\rho(x)\boldsymbol{\xi}$, we obtain the inequality

(4.27)
$\displaystyle\|(\boldsymbol{A}-\lambda\boldsymbol{I})\boldsymbol{\varphi}\|_{L(\Omega)}\leqslant\|(\boldsymbol{A}-\lambda\boldsymbol{I})\boldsymbol{\xi}\|_{L^{\infty}(B_{\varepsilon})}\|\rho\|_{L(\Omega)}\leqslant\varepsilon\|\boldsymbol{\varphi}\|_{{L(\Omega)}}.$

Thus, substituting the functions $\boldsymbol{\varphi},\psi,\boldsymbol{\xi}$
into inequality (4.22), we obtain the estimate

(4.28) $\displaystyle\|\boldsymbol{\varphi}\|_{{L(\Omega)}}$
$\displaystyle+\|\psi\|_{W^{2,p}(\Omega)}$ $\displaystyle\leqslant
K\big{(}\|(\boldsymbol{A}-\lambda\boldsymbol{I})\boldsymbol{\varphi}\|_{{L(\Omega)}}+\|\boldsymbol{B}\psi\|_{{L(\Omega)}}+\|\zeta\|_{L(\Omega)}\big{)}$
$\displaystyle\leqslant
K\big{(}(n+1)\varepsilon\|\boldsymbol{\varphi}\|_{{L(\Omega)}}+\|\boldsymbol{B}\|_{L^{\infty}(\Omega)}\|\psi\|_{L(\Omega)}\big{)}.$

Hence, choosing $\varepsilon>0$ small enough we compensate the term
$K\varepsilon(n+1)\|\boldsymbol{\varphi}\|_{L(\Omega)}$ by its counterpart on
the left-hand side of inequality (4.28) and obtain estimates

(4.29)
$\displaystyle\|\psi\|_{W^{2,p}(\Omega)}\leqslant\big{(}1-K(n+1)\varepsilon\big{)}\|\boldsymbol{\varphi}\|_{L(\Omega)}+\|\psi\|_{W^{2,p}_{\nu}(\Omega)}\leqslant
K\|\boldsymbol{B}\|_{L^{\infty}(\Omega)}\|\psi\|_{L(\Omega)}$

which obviously can not hold for all $\psi\in C_{c}^{\infty}(\Omega)$ such
that supp $\psi\subset B_{\varepsilon}$. This completes the proof that
$\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}\subset\sigma(\mathcal{L}_{p})$.
Moreover, by Lemma 4.4, the spectrum
$\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}$ is independent of $p$.

Next, we prove that the remainder of the spectrum
$\big{(}\mathcal{L}_{p},D(\mathcal{L}_{p})\big{)}$ consists of a discrete set
of eigenvalues by analysing the corresponding resolvent equations

(4.30)
$\displaystyle(\boldsymbol{A}-\lambda\boldsymbol{I})\boldsymbol{\varphi}+\boldsymbol{B}\psi$
$\displaystyle=\boldsymbol{F},\quad x\in\overline{\Omega},$
$\displaystyle\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+(d-\lambda)\psi$
$\displaystyle=G,\quad x\in\Omega$

with arbitrary $(\boldsymbol{F},G)\in L^{p}(\Omega)^{n}\times L^{p}(\Omega)$.
Since, $\lambda\notin\sigma(\boldsymbol{A}(\cdot))$ we can solve the first
equation with respect to $\boldsymbol{\varphi}$ and substitute into the second
equation to obtain the boundary value problem

(4.31) $\displaystyle\Delta_{\nu}\psi+q(\lambda)\psi=r(\lambda)$

where

(4.32) $\displaystyle q(\lambda)$
$\displaystyle=q(x,\lambda)=-\boldsymbol{C}(x)\big{(}\boldsymbol{A}(x)-\lambda\boldsymbol{I}\big{)}^{-1}\boldsymbol{B}(x)+d(x)-\lambda,$
$\displaystyle r(\lambda)$
$\displaystyle=r(x,\lambda)=-\boldsymbol{C}(x)\big{(}\boldsymbol{A}(x)-\lambda\boldsymbol{I}\big{)}^{-1}\boldsymbol{F}(x)+G(x).$

By the Analytic Fredholm Theorem for Banach spaces from [35], used in the same
way as in [21, proof of Thm 4.5] we conclude that the boundary value problem
(4.31) has no unique solution only for a countable set of
$\lambda\in\mathbb{C}$. Eigenfunctions corresponding to those eigenvalues
satisfy elliptic equation (4.31), hence, by standard elliptic regularity, they
belong to ${W^{2,p}_{\nu}(\Omega)}$ for each $p\in(1,\infty)$. This means that
the set of eigenvalues of $\mathcal{L}_{p}$ is independent of $p$. ∎

######  Remark 4.7.

By Theorem 4.6, the spectrum of the operator
$\big{(}\mathcal{L}_{\infty},D(\mathcal{L}_{\infty})\big{)}$ given by formula
(4.7) with the matrix coefficients

(4.33) $\displaystyle a_{ij},\;b_{i},\;c_{i},\;d\in
C(\overline{\Omega}),\quad\text{for}\;i,j\in\\{1,\cdots,n\\}$

satisfies

(4.34)
$\displaystyle\sigma\left(\mathcal{L}_{p}\right)=\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}\cup\\{\text{eigenvalues
of }\;\mathcal{L}_{p}\\}\subset\sigma(\mathcal{L}_{\infty}).$

Indeed, by the elliptic regularity, eigenfunctions of $-\Delta_{\nu}$ belong
to ${W^{2,p}_{\nu}(\Omega)}$ for each $p\in(1,\infty)$, thus they are
continuous up to the boundary. Using analogous estimates as those in the proof
of Theorem 4.6 with $L^{p}$-norm replaced by the $L^{\infty}$-norm on
$C(\overline{\Omega})$, we obtain that
$\sigma\big{(}\boldsymbol{A}(\cdot)\big{)}\subset\sigma(\mathcal{L}_{\infty}).$

##  5\. Instability of regular steady states

We are in a position to prove the main results of this work.

###### Proof of Theorem 2.7.

Let $(\boldsymbol{U},V)$ be a regular stationary solution (either constant or
non-constant). We consider the operator
$\big{(}\mathcal{L}_{p},D(\mathcal{L}_{p})\big{)}$ given by formula (4.4) and
the matrices

(5.1) $\displaystyle\boldsymbol{A}(x)$
$\displaystyle=D_{\boldsymbol{u}}\boldsymbol{f}({\boldsymbol{U}(x)},{V(x)}),$
$\displaystyle\boldsymbol{B}(x)$
$\displaystyle=D_{v}\boldsymbol{f}({\boldsymbol{U}(x)},{V(x)}),$
$\displaystyle\boldsymbol{C}(x)$
$\displaystyle=D_{\boldsymbol{u}}g({\boldsymbol{U}(x)},{V(x)}),$
$\displaystyle d(x)$ $\displaystyle=D_{v}g({\boldsymbol{U}(x)},{V(x)}).$

Under the assumption
$s\big{(}D_{\boldsymbol{u}}\boldsymbol{f}({\boldsymbol{U}(x)},{V(x)})\big{)}>0$
we obtain that $s(\mathcal{L}_{p})>0$ by Theorem 4.6. Since
$(\boldsymbol{U},V)$ is a regular solution, the matrix coefficients belong to
the space $C(\overline{\Omega})$ (see Remark 2.3). Consequently, Remark 4.7
provides $s(\mathcal{L}_{\infty})>0$ which implies that $(\boldsymbol{U},V)$
is linearly unstable in $C(\overline{\Omega})$. The nonlinear instability in
the Lyapunov sense of this solution is an immediate consequence of Theorem A.2
from Appendix. ∎

###### Proof of Theorem 2.8.

Now, we consider non-constant regular solution
$\big{(}\boldsymbol{U}(x),V(x)\big{)}$ such that
$s\big{(}f_{\boldsymbol{u}}(\boldsymbol{U}(\cdot),\,V(\cdot))\big{)}\leqslant
0$ and
$\det\boldsymbol{f}_{\boldsymbol{u}}\big{(}\boldsymbol{U}(x),\,V(x)\big{)}\neq
0$ for all $x\in\overline{\Omega}$. Moreover, we assume that $\Omega$ is
convex. Here, we are inspired by methods from the papers [1, 23] where
instability of non-constant solutions to one reaction-diffusion equation was
shown. We consider the operator
$\big{(}\mathcal{L}_{2},D(\mathcal{L}_{2})\big{)}$ given by formula (4.4) with
the matrices (5.1) and we study the eigenvalue problem

(5.2) $\displaystyle\boldsymbol{A}\boldsymbol{\varphi}+\boldsymbol{B}\psi$
$\displaystyle=\lambda\boldsymbol{\varphi},$
$\displaystyle\Delta_{\nu}\psi+\boldsymbol{C}\boldsymbol{\varphi}+d\psi$
$\displaystyle=\lambda\psi.$

It follows from assumptions (2.19) and (2.20) combined with Lemma 4.4 that
$\det\big{(}\boldsymbol{A}(x)-\lambda\boldsymbol{I}\big{)}\neq 0$ for all
$\lambda\geqslant 0$ and all $x\in\overline{\Omega}$, hence, the matrix
$(\boldsymbol{A}(x)-\lambda\boldsymbol{I})$ is invertible. Thus, we substitute
$\boldsymbol{\varphi}$ computed from the first equation in system (5.2) into
the second equation which takes the following form

(5.3)
$\displaystyle\Delta_{\nu}\psi+\Big{(}-\boldsymbol{C}(\boldsymbol{A}-\lambda\boldsymbol{I})^{-1}\boldsymbol{B}+d\Big{)}\psi=\lambda\psi.$

For all $\lambda\geqslant 0$ we define the family of the operators

(5.4) $\displaystyle G(\lambda)$
$\displaystyle:\psi\mapsto\Delta_{\nu}\psi+\Big{(}-\boldsymbol{C}(\boldsymbol{A}-\lambda\boldsymbol{I})^{-1}\boldsymbol{B}+d\Big{)}\psi$

with the domain $D\big{(}G(\lambda)\big{)}={W^{2,2}_{\nu}(\Omega)}$. We study
their eigenvalues, namely, the numbers $\eta\in\mathbb{C}$ for which the
equation

(5.5) $\displaystyle G(\lambda)\psi=\eta\psi$

has a nonzero solution and we are going to show that $\eta=\lambda$ for some
$\lambda>0$. This will give a positive eigenvalue of system (5.2) because,
under our assumption, system (5.2) is equivalent to the equation
$G(\lambda)\psi=\lambda\psi$.

First we show that

(5.6) $\displaystyle\sup_{\begin{subarray}{c}x\in\overline{\Omega}\\\
\lambda\geqslant
0\end{subarray}}\left|-\boldsymbol{C}(x)(\boldsymbol{A}(x)-\lambda\boldsymbol{I})^{-1}\boldsymbol{B}(x)+d(x)\right|<\infty.$

Let $\lambda\geqslant 0$. By assumption (2.20), the continuity of the
determinant, and regularity of solutions, there exists
$\overline{\varepsilon}>0$ such that every
$\rho\in\sigma\big{(}\boldsymbol{A}(\cdot)-\lambda\boldsymbol{I}\big{)}$
satisfies $|\rho|\geqslant\overline{\varepsilon}$. Consequently,
$|\det(\boldsymbol{A}(x)-\lambda\boldsymbol{I})|\geqslant\left(\max\\{\overline{\varepsilon},\lambda\\}\right)^{n}$
for every $x\in\overline{\Omega}$ and hence

(5.7) $\displaystyle\|(\boldsymbol{A}(x)-\lambda\boldsymbol{I})^{-1}\|$
$\displaystyle=\frac{1}{|\det(\boldsymbol{A}(x)-\lambda\boldsymbol{I})|}\left\|{\text{adj}(\boldsymbol{A}(x)-\lambda\boldsymbol{I})}\right\|$
$\displaystyle\leqslant\frac{{\|\text{adj}(\boldsymbol{A}(\cdot)-\lambda\boldsymbol{I})\|_{\infty}}}{\left(\max\\{\overline{\varepsilon},\lambda\\}\right)^{n}}\leqslant\frac{C(\lambda^{n}+1)}{\left(\max\\{\overline{\varepsilon},\lambda\\}\right)^{n}}\leqslant
C.$

Now, denote by $\eta_{0}(\lambda)$ the largest eigenvalue of problem
$G(\lambda)\psi=\eta_{0}(\lambda)\psi.$ The operator $\big{(}G,D(G)\big{)}$ is
self-adjoint, hence, the eigenvalue $\eta_{0}(\lambda)$ is given by Rayleigh
quotient

(5.8) $\displaystyle\eta_{0}(\lambda)$
$\displaystyle=\sup_{\begin{subarray}{c}\psi\in{W^{1,2}(\Omega)}\\\
\|\psi\|_{L^{2}(\Omega)}=1\end{subarray}}\big{\langle}G(\lambda)\psi,\psi\big{\rangle}$
$\displaystyle=\sup_{\begin{subarray}{c}\psi\in{W^{1,2}(\Omega)}\\\
\|\psi\|_{L^{2}(\Omega)}=1\end{subarray}}\left(-\int_{\Omega}|\nabla\psi|^{2}\,\text{d}x+\int_{\Omega}\Big{(}-\boldsymbol{C}(\boldsymbol{A}-\lambda\boldsymbol{I})^{-1}\boldsymbol{B}+d\Big{)}\psi^{2}\,\text{d}x\right).$

By inequality (5.6), we immediately obtain that $\sup_{\lambda\geqslant
0}\eta_{0}(\lambda)<\infty$. In the next step of this proof, we show that
$\eta_{0}(0)>0$.

In the one dimensional case $\Omega=(a,b)\subset\mathbb{R}$, differentiating
equations (LABEL:seq1) we obtain

(5.9)
$\displaystyle\boldsymbol{A}\boldsymbol{U}_{x}+\boldsymbol{B}V_{x}=0\quad\text{and}\quad(V_{x})_{xx}+\boldsymbol{C}\boldsymbol{U}_{x}+dV_{x}=0$

which yields

(5.10)
$\displaystyle(V_{x})_{xx}+(-\boldsymbol{C}\boldsymbol{A}^{-1}\boldsymbol{B}+d)V_{x}=0.$

Choosing $\psi=V_{x}/\|V_{x}\|_{L^{2}(\Omega)}$ in the functional in the right
hand side of equation (5.8), integrating by parts, and using relations (5.10)
we obtain the inequality

(5.11)
$\displaystyle\eta_{0}(0)\geqslant\frac{1}{\|V_{x}\|_{L^{2}(\Omega)}}\int_{\Omega}-V_{xx}V_{xx}+(-\boldsymbol{C}\boldsymbol{A}^{-1}\boldsymbol{B}+d)V_{x}^{2}\,\text{d}x=0.$

The hypothesis $\eta_{0}(0)=0$ implies that $V_{x}$ is a maximizer of the
functional in the right hand side of equation (5.8) with $\lambda=0$, hence,
$V_{x}$ satisfies equation (5.10) with Neumann boundary conditions, namely,
$V_{xx}(x)=0$ for $x\in\\{a,b\\}=\partial\Omega$. This ordinary differential
equation for $V_{x}$ supplemented with the boundary conditions

(5.12) $\displaystyle V_{x}(x)=0\quad\text{and}\quad
V_{xx}(x)=0\quad\text{for}\quad x\in\\{a,b\\}=\partial\Omega$

has the unique solution $V_{x}\equiv 0$. Hence, $V$ (and so $U$) is a constant
function which is a contradiction. Thus, $\eta_{0}(0)>0$.

The proof in higher dimensions is analogous. Differentiating equations
(LABEL:seq1) with respect to $x_{j}$ for $j\in\\{1,...,n\\}$, we obtain

(5.13)
$\displaystyle\boldsymbol{A}\boldsymbol{U}_{x_{j}}+\boldsymbol{B}V_{x_{j}}=0\quad\text{and}\quad\Delta
V_{x_{j}}+\boldsymbol{C}\boldsymbol{U}_{x_{j}}+dV_{x_{j}}=0.$

Thus, by direct calculations, analogous to those in (5.11) and by the
divergence theorem, we obtain

(5.14) $\displaystyle\sum_{j=1}^{N}\langle
G(0)V_{x_{j}},V_{x_{j}}\rangle=-\frac{1}{2}\int_{\partial\Omega}\partial_{\nu}|\nabla
V|^{2}\,d\sigma.$

Now, we use the use the inequality $\partial_{\nu}|\nabla V|^{2}\leqslant 0$
on $\partial\Omega$ which is valid for every $V\in C^{2}(\overline{\Omega})$
and for every convex $\Omega\subset\mathbb{R}^{n}$. Detailed calculations
leading to this inequality can be found either in [1, p. 269] or in [23, Lemma
5.3]. Thus, at least one term on the left hand side of equality (5.14), is
nonnegative, which implies that $\eta_{0}(0)\geqslant 0$ in the variational
formula (5.8). If $\eta_{0}(0)=0$, by a variational argument analogously as in
the one dimensional case, equations (5.13) have to be satisfied by the
function $V_{x_{j_{0}}}$ together with the Neumann boundary conditions for
$V_{x_{j_{0}}}$. Applying the Hopf maximum principle (see e.g. [28, Theorems 5
and 7 in Sec. 3]), we obtain immediately that $V_{x_{j_{0}}}\equiv 0$ which is
a contradiction because we deal with the non-constant solution $V$. Hence, we
have got $\eta_{0}(0)>0$.

Therefore, since $\eta_{0}(\lambda)$ is a continuous and bounded function on
$[0,\infty)$ satisfying $\eta_{0}(0)>0$, there exists $\overline{\lambda}>0$
such that $\eta_{0}(\overline{\lambda})=\overline{\lambda}$. The number
$\overline{\lambda}$ is a positive eigenvalue of system (5.2) and hence
$s(\mathcal{L}_{2})>0$. Since $(\boldsymbol{U},V)$ is a regular solution, we
conclude by Remark 4.7 that the number $\overline{\lambda}$ is also an
eigenvalue of $\mathcal{L}_{\infty}$. Since $s(\mathcal{L}_{\infty})>0$, the
nonlinear instability follows immediately from Theorem A.2 in Appendix. ∎

##  Appendix A Linearisation principle

In order to show that a linear instability of regular stationary solutions of
problem (1.1)-(LABEL:ini) implies their nonlinear instability in the Lyapunov
sense, we recall the reasoning from [21, Thm 2.1]. That approach concerns
general evolution equation

(A.1) $w_{t}=\mathcal{L}w+\mathcal{N}(w),\qquad w(0)=w_{0},$

where $\mathcal{L}$ is the generator of a $C_{0}$-semigroup of linear
operators $\\{e^{t\mathcal{L}}\\}_{t\geq 0}$ on a Banach space $X$ and
$\mathcal{N}$ is a nonlinear operator such that
${\mathcal{N}(w)}=o\left(\|{w}\|\right)$ as $w\to 0$. Here, we recall one of
the possible instability results, where the existence of a positive part of
the spectrum of the linear operator $\mathcal{L}$ is sufficient to show that
the zero solution of equation (A.1) is unstable.

######  Theorem A.1 ([34, Thm 1]).

Consider an abstract problem (A.1), where

  1. (1)

the linear operator $\mathcal{L}$ generates a strongly continuous semigroup of
linear operators on a Banach space $X$,

  2. (2)

the intersection of the spectrum of $\mathcal{L}$ with the right half-plane
$\\{\lambda\in\mathbb{C}:\;{\rm Re\,}\lambda>0\\}$ is nonempty.

  3. (3)

$\mathcal{N}:X\to X$ is continuous and there exist constants $\rho>0$,
$\eta>0$, and $C>0$ such that $\|\mathcal{N}(w)\|_{X}\leq C\|w\|_{X}^{1+\eta}$
for all $\|w\|_{X}<\rho$.

Then the zero solution of this equation is (nonlinearly) unstable.

######  Theorem A.2.

Let $(\boldsymbol{U},V)$ be a regular stationary solution of problem
(1.1)-(LABEL:ini). Assume that the corresponding linearisation operator
$\mathcal{L}_{\infty}$, defined in Proposition 4.2 and with matrices (5.1),
satisfies $s(\mathcal{L}_{\infty})>0$. Then $(\boldsymbol{U},V)$ is
nonlinearly unstable in the space $C(\overline{\Omega})^{n+1}$.

###### Proof.

Applying the usual linearisation procedure we obtain that the perturbation of
a stationary solution $(\boldsymbol{U},V)$ satisfies the semi-linear problem

(A.2) $\displaystyle\frac{\partial}{\partial
t}\begin{pmatrix}\boldsymbol{\varphi}\\\
\psi\end{pmatrix}=\mathcal{L}_{\infty}\begin{pmatrix}\boldsymbol{\varphi}\\\
\psi\end{pmatrix}+\begin{pmatrix}\boldsymbol{R}_{1}(\boldsymbol{\varphi},\psi)\\\
R_{2}(\boldsymbol{\varphi},\psi)\end{pmatrix}$

with the operator $\mathcal{L}_{\infty}$ given by formula (4.4) with matrices
(5.1) and where the remainders $\boldsymbol{R}_{1}$ and $R_{2}$ are obtained
from nonlinearities $\boldsymbol{f},\,g$ via the Taylor expansion. We apply
Theorem A.1 with the Banach space $X=C(\overline{\Omega})^{n+1}$. The
remainders $\boldsymbol{R}_{1},\,R_{2}$ in system (A.2) satisfy the condition
(3) of Theorem A.1 with $\eta=1$. ∎

## References

  * [1] R. G. Casten and C. J. Holland, Instability results for reaction diffusion equations with Neumann boundary conditions, J. Differential Equations, 27 (1978), pp. 266–273.
  * [2] L. H. Chuan, T. Tsujikawa, and A. Yagi, Asymptotic behavior of solutions for forest kinematic model, Funkcial. Ekvac., 49 (2006), pp. 427–449.
  * [3] S. Cygan, A. Marciniak-Czochra, G. Karch, and K. Suzuki, Stable discontinuous stationary solutions to reaction-diffusion-ODE systems, preprint in arXiv, (2021).
  * [4] K.-J. Engel, Operator matrices and systems of evolution equations, Surikaisekikenkyusho Kokyuroku, (1996), pp. 61–80.  Nonlinear evolution equations and their applications (Japanese) (Kyoto, 1995).
  * [5] K.-J. Engel and R. Nagel, One-parameter semigroups for linear evolution equations, vol. 194 of Graduate Texts in Mathematics, Springer-Verlag, New York, 2000.  With contributions by S. Brendle, M. Campiti, T. Hahn, G. Metafune, G. Nickel, D. Pallara, C. Perazzoli, A. Rhandi, S. Romanelli and R. Schnaubelt.
  * [6] Y. Golovaty, A. Marciniak-Czochra, and M. Ptashnyk, Stability of nonconstant stationary solutions in a reaction-diffusion equation coupled to the system of ordinary differential equations, Commun. Pure Appl. Anal., 11 (2012), pp. 229–241.
  * [7] G. Harris and C. Martin, The roots of a polynomial vary continuously as a function of the coefficients, Proc. Amer. Math. Soc., 100 (1987), pp. 390–392.
  * [8] S. Härting and A. Marciniak-Czochra, Spike patterns in a reaction-diffusion ODE model with Turing instability, Math. Methods Appl. Sci., 37 (2014), pp. 1377–1391.
  * [9] S. Härting, A. Marciniak-Czochra, and I. Takagi, Stable patterns with jump discontinuity in systems with Turing instability and hysteresis, Discrete Contin. Dyn. Syst., 37 (2017), pp. 757–800.
  * [10] Y. Hosono and M. Mimura, Singular perturbations for pairs of two-point boundary value problems of Neumann type, in Mathematical analysis on structures in nonlinear phenomena (Tokyo, 1978), vol. 2 of Lecture Notes Numer. Appl. Anal., Kinokuniya Book Store, Tokyo, 1980, pp. 79–138.
  * [11] T. Kato, Perturbation theory for linear operators, Classics in Mathematics, Springer-Verlag, Berlin, 1995.  Reprint of the 1980 edition.
  * [12] K. Kishimoto and H. F. Weinberger, The spatial homogeneity of stable equilibria of some reaction-diffusion systems on convex domains, J. Differential Equations, 58 (1985), pp. 15–21.
  * [13] A. Köthe, A. Marciniak-Czochra, and I. Takagi, Hysteresis-driven pattern formation in reaction-diffusion-ODE systems, Discrete Contin. Dyn. Syst., 40 (2020), pp. 3595–3627.
  * [14] C. Kowall, Uniform Shadow Limit Reduction for Reaction-Diffusion-ODE Systems, PhD thesis, University of Heidelberg, 01 2021.
  * [15] Y. Li, A. Marciniak-Czochra, I. Takagi, and B. Wu, Bifurcation analysis of a diffusion-ODE model with Turing instability and hysteresis, Hiroshima Math. J., 47 (2017), pp. 217–247.
  * [16] , Steady states of FitzHugh-Nagumo system with non-diffusive activator and diffusive inhibitor, Tohoku Math. J. (2), 71 (2019), pp. 243–279.
  * [17] A. Lunardi, Analytic semigroups and optimal regularity in parabolic problems, Modern Birkhäuser Classics, Birkhäuser/Springer Basel AG, Basel, 1995.  [2013 reprint of the 1995 original] [MR1329547].
  * [18] A. Marciniak-Czochra, Receptor-based models with hysteresis for pattern formation in hydra, Math. Biosci., 199 (2006), pp. 97–119.
  * [19] , Reaction-diffusion models of pattern formation in developmental biology, in Mathematics and life sciences, vol. 1 of De Gruyter Ser. Math. Life Sci., De Gruyter, Berlin, 2013, pp. 191–212.
  * [20] A. Marciniak-Czochra, G. Karch, and K. Suzuki, Unstable patterns in reaction-diffusion model of early carcinogenesis, J. Math. Pures Appl. (9), 99 (2013), pp. 509–543.
  * [21] , Instability of Turing patterns in reaction-diffusion-ODE systems, J. Math. Biol., 74 (2017), pp. 583–618.
  * [22] A. Marciniak-Czochra, M. Nakayama, and I. Takagi, Pattern formation in a diffusion-ODE model with hysteresis, Differential Integral Equations, 28 (2015), pp. 655–694.
  * [23] H. Matano, Asymptotic behavior and stability of solutions of semilinear diffusion equations, Publ. Res. Inst. Math. Sci., 15 (1979), pp. 401–454.
  * [24] M. Mimura, M. Tabata, and Y. Hosono, Multiple solutions of two-point boundary value problems of Neumann type with a small parameter, SIAM J. Math. Anal., 11 (1980), pp. 613–631.
  * [25] G. Mulone and V. A. Solonnikov, Linearization principle for a system of equations of mixed type, Nonlinear Anal., 71 (2009), pp. 1019–1031.
  * [26] W.-M. Ni, Qualitative properties of solutions to elliptic problems, in Stationary partial differential equations. Vol. I, Handb. Differ. Equ., North-Holland, Amsterdam, 2004, pp. 157–233.
  * [27] B. Perthame and J. Skrzeczkowski, Fast reaction limit with nonmonotone reaction function, arXiv:2008.11086, (2020).
  * [28] M. H. Protter and H. F. Weinberger, Maximum principles in differential equations, Springer-Verlag, New York, 1984.  Corrected reprint of the 1967 original.
  * [29] P. H. Rabinowitz, A bifurcation theorem for potential operators, J. Functional Analysis, 25 (1977), pp. 412–424.
  * [30] , Minimax methods in critical point theory with applications to differential equations, vol. 65 of CBMS Regional Conference Series in Mathematics, Published for the Conference Board of the Mathematical Sciences, Washington, DC; by the American Mathematical Society, Providence, RI, 1986.
  * [31] K. Sakamoto, Construction and stability analysis of transition layer solutions in reaction-diffusion systems, Tohoku Math. J. (2), 42 (1990), pp. 17–44.
  * [32] R. Schaaf, Global behaviour of solution branches for some Neumann problems depending on one or several parameters, J. Reine Angew. Math., 346 (1984), pp. 1–31.
  * [33] , Stationary solutions of chemotaxis systems, Trans. Amer. Math. Soc., 292 (1985), pp. 531–556.
  * [34] J. Shatah and W. Strauss, Spectral condition for instability, in Nonlinear PDE’s, dynamics and continuum physics (South Hadley, MA, 1998), vol. 255 of Contemp. Math., Amer. Math. Soc., Providence, RI, 2000, pp. 189–198.
  * [35] S. Steinberg, Meromorphic families of compact operators, Arch. Rational Mech. Anal., 31 (1968/69), pp. 372–379.
  * [36] I. Takagi and C. Zhang, Existence and stability of patterns in a reaction-diffusion-ode system with hysteresis in non-uniform media, Discrete Contin. Dyn. Syst., 41 (2021), pp. 3109–3140.
  * [37] , Pattern formation in a reaction-diffusion-ODE model with hysteresis in spatially heterogeneous environments, J. Differential Equations, 280 (2021), pp. 928–966.
  * [38] J. Wang, The stability of equilibria for a reaction-diffusion-ODE system on convex domains, Appl. Math. Lett., 93 (2019), pp. 147–152.
  * [39] H. F. Weinberger, A simple system with a continuum of stable inhomogeneous steady states, in Nonlinear partial differential equations in applied science (Tokyo, 1982), vol. 81 of North-Holland Math. Stud., North-Holland, Amsterdam, 1983, pp. 345–359.
  * [40] A. Yagi, Abstract parabolic evolution equations and their applications, Springer Monographs in Mathematics, Springer-Verlag, Berlin, 2010\.