Applications of a class of Herglotz operator pencils
Andrea K. Barreiro Jared C. Bronski Zoi Rapti








Applications of a class of  Herglotz operator pencils.
    Andrea K. Barreiro [Department of Mathematics, Southern
    Methodist University] 

 Jared C. Bronski[Department of Mathematics, University of Illinois] 
 Zoi Rapti[Department of Mathematics, University of Illinois]  
    December 30, 2023
=============================================================================================================================================================================================================================




We identify a  class of operator pencils, arising in a number of
 applications, which have only real
eigenvalues. In the one-dimensional case we prove a novel version of
the Sturm oscillation theorem: if the dependence on the eigenvalue
parameter is of degree k then the real axis can be partitioned into
a union of k disjoint intervals, each of which enjoys a Sturm
oscillation theorem: on each interval there is an increasing sequence
of eigenvalues that are indexed by the number of roots of the
associated eigenfunction.  One consequence of this is that it
guarantees that 
the spectrum of these
operator pencils has finite accumulation points, implying that the
operators 
do not have compact resolvents.  
As an application we apply this theory to an epidemic model and several species dispersal models arising in biology.




§ INTRODUCTION



There are a number of contexts in which  eigenvalue pencils 
arise in applied mathematics <cit.>. 
When studying the stability of standing waves and other coherent structures 
it is frequently the case that one must consider a non-self-adjoint
eigenvalue pencil. In the context of integrable systems the scattering
problem that diagonalizes an integrable  nonlinear flow frequently
takes the form of a eigenvalue pencil; take, for instance the
Sine-Gordon equation in the laboratory frame <cit.>  or the derivative
nonlinear Schrödinger equation <cit.>. Additionally we
will present a number of other problems arising in the analysis of
biological problems where rational operator pencils arise.     
The difficulty with these problems, from the point of view of 
both analysis and numerical experimentation, is that one 
has no apriori information as to where the 
spectrum of the problem might lie. The purpose of this note is to
point out a class of pencils, arising in a number of applications,
which only admit real eigenvalues, namely those that depend on the
spectral parameter in a Herglotz way. We show that a
number of the properties of self-adjoint eigenvalue problems carry
over to Herglotz operator pencils: the eigenvalues are real and
semi-simple, and for second order rational Herglotz pencils we have 
an interesting generalization of the Sturm oscillation theorem. We begin 
by reminding the reader of the definition of a Herglotz function
<cit.>. These functions are also sometimes referred to as
Nevanlinna or Pick functions, but in the spectral theory literature
they are most often referred to as Herglotz functions, so that is the
usage that we will adopt here. 


A meromorphic function f(λ) is Herglotz if Im(λ)>0 (resp. Im(λ)<0) implies that Im(f(λ))>0 (resp. Im(f(λ))<0))
     
The following is standard:

If f(λ) is Herglotz then 


  * The zeroes and poles of f(λ) lie on the real axis.

  * If λ is real and not a pole then
f'(λ)>0. 

  * The zeroes and poles of f(λ) alternate on the real axis. 



For the purposes of this paper the canonical example of a Herglotz function 
is the rational function 

    f(λ) = ∑_i=1^N A_i/α_i - λ + B + C λ

with A_i,B,C real and positive and α_i real, as all of the
problems we consider lead to rational operator pencils.   More
generally one has a representation formula for any Herglotz function 

    f(λ) = B + C λ + ∫_ℝ( 1/x-λ -
    λ/1+λ^2) dμ(x)
 
for some measure dμ(x) but we are unaware of any applications that require
this general representation. 



§ MOTIVATING EXAMPLES


In this section we present some examples intended to motivate  the later
discussion.  As a first example we consider a reaction-diffusion equation of the form 

    u_t = u_xx + f(u,v) 
        v_t = g(u,v)

together with appropriate boundary conditions. Assuming that we can
find a steady state solution 

    u(x,t) = u̅(x) 
       v(x,t) = v̅(x)

we consider the linearization of the evolution around the steady state
(u̅, v̅). This leads to the following equation for the
tangent flow 

    p_t = p_xx + ∂ f/∂ u(u̅,v̅) p +∂ f/∂ v(u̅,v̅) q  
        q_t = ∂ g/∂ u(u̅,v̅) p +∂ g/∂ v(u̅,v̅) q

or, equivalently the linearized stability problem

    λ p = p_xx + ∂ f/∂ u(u̅,v̅) p
      +∂ f/∂ v(u̅,v̅) q  
      
       λ q = ∂ g/∂ u(u̅,v̅) p +∂ g/∂ v(u̅,v̅) q

This is an ordinary, though non-selfadjoint eigenvalue problem, but
after one eliminates the function q it becomes a rational eigenvalue
pencil. 

    λ p = p_xx + ∂ f/∂ u(u̅,v̅) p
    +∂ f/∂ v(u̅,v̅)  ∂
      g/∂ u(u̅,v̅) 1/λ - ∂ g/∂ v(u̅,v̅)


The main observation is that under the sign condition f_v g_u(u̅,v̅)≥ 0 the right hand side of the operator pencil is (for x
fixed), a Herglotz function which implies some very nice behavior in
the eigenvalue problem, including (as we shall see) reality of the spectrum and
semi-simplicity of the eigenvalues. This generalizes the case when
f_v=g_u, for which case the flow is a gradient and the second variation is
obviously self-adjoint.  

One could object that the original system (<ref>,
<ref>) is, in fact, self-adjoint under the modified inner
product ‖ (p,q)‖_2^2 = ∫ v^2 + f_w/g_vw^2 when
the sign condition f_w g_v>0 is met. In some sense this must happen,
since it is an old theorem of Drazin and  Haynsworth <cit.> that any matrix
having all real eigenvalues and n linearly  independent eigenvectors
is self-adjoint under some appropriate inner product.  However, in
general it is unclear how to check whether there exists an
inner product rendering a given operator self-adjoint. Checking that a
given function is Herglotz, on the other hand, is relatively
straightforward, and amounts (for a rational pencil) to checking a
finite
number of reality/sign conditions.  

One biological model that falls into this class is the following system modeling
the spread of a human disease through infective propagules <cit.> 

    u_t    =    d u_xx - a_11 u+ a_12 v = d u_xx + f(u, v) 
    
    v_t    =   g̃(u)- a_22 v = g(u, v).

Here, u denotes the spatial density of the infectious propagules, v denotes the 
density of the human infective class; the parameters d, a_ij >0; the function
g̃: ℝ_+ →ℝ_+ and 
g̃'(u)>0 for u>0.
Then, the corresponding eigenvalue problem reads 

    d p_xx -a_11 p = λ p - a_12g̃' 1/λ + a_22 p.

Since a_12g̃'>0 the condition for the reality of the spectrum follows. 

Equations of similar form also arise in combustion, typically solid
combustion. In this situation heat is free to diffuse but the solid
reactant does not diffuse, at least on combustion time-scales. For
work on this problem we refer the reader to the article of Ghazaryan,
Latushkin and Schechter <cit.>, but we do note that in
such combustion problems the sign condition is typically not
  met: in the combustion context one expects f_w g_u<0.

Similarly  we can consider a system with three reacting species, one of
which is allowed to diffuse 

    u_t = u_xx + f(u,v,w) 
           v_t = g(u,v,w) 
           w_t = h(u,v,w).

Linearizing about a steady state solution yields the following eigenvalue
problem

    λ p = p_xx + f_u p +f_v q + f_w s 
          λ q = g_u p + g_v q + g_w s
          λ s = h_u p + h_v q + h_w s






Solving the last two for q, s in terms of  p results in

    q = g_u(λ - h_w) + h_u g_w /(λ - g_v)(λ - h_w) - g_w h_v p,  
    s = (λ - g_v) h_u + h_v g_u /(λ - g_v)(λ - h_w) - g_w h_v p

and substituting into the first equation yields

    p_xx + f_u p = ( λ - 
    f_v g_u(λ - h_w)+ h_u g_w/(λ - g_v)(λ - h_w) - g_w h_v - 
    f_w h_u(λ - g_v)+ h_v g_u/(λ - g_v)(λ - h_w)-g_w h_v) p 
    ≡ H(λ) p 
           p_xx + f_u p = ( λ - α(x) λ + β(x)
       /λ^2 + γ(x) λ + δ(x))

The function λ - α(x) λ + β(x)
   /λ^2 + γ(x) λ + δ is (for fixed x) a
   Herglotz function of λ if and only if the following
conditions are met 


  * The roots of λ^2 + γ(x) λ + δ(x) are
  real.

  * The residues of the function α(x) λ + β(x)
   /λ^2 + γ(x) λ + δ at the roots are
   positive.

 The above two conditions translate into the following sign conditions
 on the coefficients.


  * βγ - αδ >0

  * β^2 - α (βγ - αδ)<0

Equivalently, reality of  the spectrum follows if 

    | [ f_v f_w; g_v g_w ]| | [ g_u g_v; h_u h_v ]| + | [ f_v f_w; h_v h_w ]| | [ g_u g_w; h_u h_w ]| >0   
    ( h_u | [ g_u g_w; h_u h_w ]| + g_u | [ g_u g_v; h_u h_v ]| ) (  f_v | [ f_v f_w; g_v g_w ]| + f_w | [ f_v f_w; h_v h_w ]|)>0
 





A system in the form  of (<ref>) is considered in <cit.> 
as a model of the dynamics of two plants v, w and one herbivore u

    u_t    =    d u_xx +(b_1 s_1 v + 
    b_2 s_2 w (1 - s_2 w/4 n (1+h_1 s_1 v+ h_2 s_2 w)) - m ) 
    u/1+h_1 s_1 v+ h_2 s_2 w
       =    d u_xx + f(u, v, w) 
    
    v_t    =     r_1 v ( 1 -(v+a_12 w) ) - 
    s_1 u v/1+h_1 s_1 v+ h_2 s_2 w  = g(u, v, w), 
    
    w_t    =     r_2 w ( 1 -(w+a_21 v) ) - 
    s_2 u w/1+h_1 s_1 v+ h_2 s_2 w(1 - s_2 w/4 n (1+h_1 s_1 v+ h_2 s_2 w))  = h(u, v, w),

where d, a_ij >0 and all other parameters are positive. 



To show the Herglotz property for H(λ), we will first use
well-known properties of  Herglotz functions, namely, 


  * the sum of two Herglotz functions is Herglotz,


  * if m(λ) is a Herglotz function, then so is  - 1/m(λ). 

We consider f_v ≠ 0 and f_w ≠ 0, since otherwise, the problem is 
straightforward.
Assume g_u ≠ 0. Then,

    H_1(λ)    =    - 
    f_v g_u(λ - h_w) + h_u g_w /(λ - g_v)(λ - h_w) - g_w h_v
       =    - f_v g_u λ - h_w + h_u g_w/g_u/(λ - g_v)(λ - h_w) - g_w h_v
       =   -f_v g_u /λ - g_v - h_u g_w/g_u + 
    ( h_u g_w/g_u -h_w ) (h_u g_w/g_u +g_v ) 
    + g_v h_w - g_w h_v/λ - h_w + 
    h_u g_w/g_u.

H_1(λ) has the form 

    H_1(λ) = -A/λ +B + C/λ+D ,

so  a sufficient condition for it to be Herglotz is that
A>0 and C<0 or equivalently

    f_v g_u >0    ( h_u g_w/g_u -h_w ) 
    (h_u g_w/g_u +g_v ) + g_v h_w - g_w h_v<0.
 
More compactly: 

    f_v g_u >0     - 1/g_u^2| [ g_u g_w; h_u h_w ]| 
    (h_u g_w+g_u g_v ) + | [ g_v g_w; h_v h_w ]| <0.
 
Similarly, if h_u  ≠ 0, then 

    H_2(λ)    =    - f_w 
    h_u(λ - g_v) + h_v g_u /(λ - g_v)(λ - h_w) - g_w h_v
       =    - f_w h_u λ - g_v + h_v g_u/h_u/(λ - g_v)(λ - h_w) - g_w h_v
       =   -f_w h_u /λ - h_w - h_v g_u/h_u + 
    ( g_u h_v/h_u -g_v ) (g_u h_v/h_u +h_w ) 
    + g_v h_w - g_w h_v/λ - g_v + 
    h_v g_u/h_u.

A sufficient condition for it to be Herglotz is that

    f_w h_u >0    ( g_u h_v/h_u -g_v ) (g_u h_v/h_u +h_w )   
    + g_v h_w - g_w h_v<0.
 
or equivalently: 

    f_w h_u >0    1/h_u^2| [ g_u g_v; h_u h_v ]| 
    (g_u h_v+h_u h_w ) + | [ g_v g_w; h_v h_w ]| <0.
 

Now, if g_u =0 = h_u the problem is straightforward again. Hence, we assume
not both are zero, say g_u=0 and h_u ≠ 0. Then, we define 

    H_3(λ)    =    - 
    f_v g_w h_u + f_w h_u (λ- g_v)/(λ - g_v)(λ - h_w) - g_w h_v
       =    - f_w h_u λ - g_v + f_v g_w/f_w/(λ - g_v)(λ - h_w) - g_w h_v
       =   -f_w h_u /λ - h_w - f_v g_w/f_w + 
    ( f_v g_w/f_w -g_v ) (f_v g_w/f_w + h_w )   
    + g_v h_w - g_w h_v .

A sufficient condition for it to be Herglotz is that

    f_w h_u >0    ( f_v g_w/f_w -g_v ) (f_v g_w/f_w + h_w )   
    + g_v h_w - g_w h_v<0.
 
or equivalently: 

    f_w h_u >0    1/f_w^2| [ f_v f_w; g_v g_w ]| 
    (g_w f_v+h_w f_w ) + | [ g_v g_w; h_v h_w ]| <0.
 
Finally, if g_u ≠ 0 but h_u=0, we define

    H_4(λ)    =    - 
    f_w g_u h_v + f_u g_u (λ- h_w)/(λ - g_v)(λ - h_w) - g_w h_v
       =    - f_u g_u λ - h_w + h_v f_w/f_u/(λ - g_v)(λ - h_w) - g_w h_v
       =   -f_u g_u /λ - g_v - h_v f_w/f_u + 
    ( h_v f_w/f_u -h_w ) (h_v f_w/f_u + g_v )   
    + g_v h_w - g_w h_v .

A sufficient condition for it to be Herglotz is that

    f_u g_u >0    ( h_v f_w/f_u -h_w ) (h_v f_w/f_u + g_v )   
    + g_v h_w - g_w h_v <0.
 
or equivalently: 

    f_u g_u >0      -1/f_u^2| [ f_v f_v; h_u h_v ]| 
    (h_v f_w+g_v f_u ) + | [ g_v g_w; h_v h_w ]| <0.
 
Added 2/2 Monotonicity
Monotonicity argument shows a x + b/c x^2 + d x + e is
Herglotz if ac>0 and c(b^2 c - a b d + a^2 f)>0. 
Plugging in to
above gives the conditions

    f_w h_u - f_v g_u    > 0 
    
    -4 (f_v f_w g_u g_v^2 h_u - f_v^2 g_u g_v g_w h_u - 2 f_w^2 g_v^2 h_u^2 + 
       3 f_v f_w g_v g_w h_u^2 - f_v^2 g_w^2 h_u^2 - f_v f_w g_u^2 g_v h_v + 
       f_v^2 g_u^2 g_w h_v + 3 f_w^2 g_u g_v h_u h_v - 4 f_v f_w g_u g_w h_u h_v + 
       f_w^2 g_w h_u^2 h_v - f_w^2 g_u^2 h_v^2 + f_v^2 g_u g_w h_u h_w - 
       2 f_w^2 g_v h_u^2 h_w + f_v f_w g_w h_u^2 h_w + f_v f_w g_u^2 h_v h_w + 
       f_w^2 g_u h_u h_v h_w - f_v f_w g_u h_u h_w^2)    > 0



  
In <cit.> three models of morphogen concentrations were analyzed. 
Nonlinear reaction-diffusion equations  were used to model
morphogen  gradients that arise in the patterning of Drosophila wings.
In a series of works, the authors showed that the linear
stability analysis of steady state solutions yields a nonlinear eigenvalue
problem. They also demonstrated that the spectrum is real and stable.
In this section, we will offer an alternative proof  of these facts using
our main result.

Following the notation in <cit.>, the dimensionless normalized 
concentration of the diffusing morphogen species is denoted with a(x, t) and 
that of the morphogens that are bound to receptors is denoted by b(x, t).
Then, the system reads

    ∂ a/∂ t = ∂^2 a/∂ x^2 -h_0 a(1-b) +f_0 b,
       x ∈ (0, 1),   t>0,   
          ∂ a(0,t)/∂ t = ν_0 -h_o a(0,t) (1-b(0,t)) + f_0 b(0,t),   a(1,t)=0 
          ∂ b/∂ t = h_0 a(1-b)  - (f_0 +g_0) b,
       x ∈ [0, 1],   t>0.

Here, the constants h_0, f_0, g_0, ν_0 are positive and correspond to 
normalized rates.




















If  one  linearizes about the time-independent steady-state solution 
( a̅(x),  b̅(x) ) considering solutions of  the form 

    a(x,t) = a̅(x) + e^λ tâ(x),   
    b(x,t) = b̅(x) +e^λ tb̂(x,t),






























and after some standard algebra to eliminate the algebraic equation for the bound 
morphogen b, the following nonlinear eigenvalue problem is obtained

    â_xx - ( λ + q(x, λ) ) â =0,     
    q(x, λ) = h_0 (f_0 +g_0)/f_0 + g_0 + h_0 a̅λ + g_0/λ + f_0 + g_0 + h_0 a̅,
 
together with appropriate boundary conditions.




One easily notices that q(x, λ) is a Herglotz function, since 

    f_0+g_0+h_0 a̅ -g_0 = f_0 + h_0 a̅ >0.
 Therefore,

    Q(x, λ) = λ + q(x, λ)
 is a Herglotz function as well.
Using our main result it then follows  that
the spectrum is real. 

As these examples illustrate, a wealth of biological models have similar structure.
Their common theme is that some degrees of freedom permit diffusion or mobility, while others do not.


In many compartmental systems arising in biological applications, 
some of the components are mobile (such as herbivores and pollinators) while others 
are not (such as plants) <cit.>. 
In epidemic compartmental models, such distinctions between mobile and 
motionless compartments also arise <cit.> due 
to the fact that infected animals are mainly responsible for spreading the disease, 
or when infective propagules disperse randomly while the infected  population is 
characterized by small  mobility.  Another setting where such models arise is in pattern formation, where diffusing
morphogens interact with others that are bound to cell receptors <cit.>.

The final example is more theoretical, and is motivated by work of Kapitula and Promislow <cit.> on
the so-called Krein matrix method. Consider a self-adjoint eigenvalue
problem that can be written in the following block-partitioned form

    ( H-λ) v = ([ A-λ   B; B^t C-λ ])
                                                                   v = 0,

with A =  A^t and C =  C^t. 
The main idea is that, if we algebraically eliminate certain degrees
of freedom, we are naturally led to a Herglotz pencil. To see this
note that, if λ is in the resolvent set of C then we
have the following identity.

    U H U^t    =    
    ([            I - B( C-λ)^-1;            0            I ])
    ([   A   B; B^t   C ])
      ([              I              0; -( C-λ)^-1 B^t              I ])
       =
     ([ A- λ - B^t ( C-λ)^-1 B                      0;                      0                 ( C-λ) ])

Therefore for λ in the resolvent set of C the operator
U is bounded and invertible and the
eigenvalues of the full operator are eigenvalues of the Herglotz
operator pencil A - λ -  B^t ( C-λ)^-1 B. If λ is in the spectrum of C (the complement
of the resolvent set) one must do
additional work to determine whether or not λ  is an eigenvalue of the full
problem – see the work of Kapitula and Promislow for details. This
example shows that Schur reduction of a self-adjoint eigenvalue
problem leads naturally to a Herglotz operator pencil, and further motivates
the idea that Herglotz pencils should behave
like self-adjoint eigenvalue problems. 



§ MAIN RESULTS



Our basic observation  is that the
standard elementary argument that the spectrum of a self-adjoint
operator is real carries over in a natural way to operator pencils
that depend on the eigenvalue in a Herglotz way.  The main claim of this paper is 

Suppose that we have an eigenvalue pencil of the form 

    Hψ = ( A_0 λ - ∑_k=1A_i/λ - α_i) ψ

where H is a self-adjoint operator on some domain ( H), and A_i are
self-adjoint positive definite operators on  ( H).  Then the eigenvalue pencil has only real
eigenvalues. Furthermore the eigenvalues are semi-simple. 


Suppose that λ is an eigenvalue with eigenfunction ψ. 
Taking the inner product of both sides with ψ we have that 

    ⟨ψ,  Hψ⟩ = ⟨ψ,  A_0
    ψ⟩λ - ∑_k=1⟨ψ, A_k ψ⟩/λ - α_i

From the self-adjointness of H we have that the left-hand side is 
purely real; the righthand side is a Herglotz function. 
Thus, if ψ(x) is an eigenfunction the λ is a root of some
Hergoltz function g̃_ψ(λ) and therefore must be real. 

To see the semi-simplicity of the eigenvalues note that the
condition for the existence of a Jordan chain for an operator pencil
is the existence of a solution to 

    M(λ_0) ψ_0  = 0   
        M(λ_0) ψ_1 +   d M/dλ(λ_0) ψ_0 = 0

with ψ_0 non-zero. Taking the inner product of the second
equation  with ψ_0 gives 

    ⟨ψ_0, d M/dλ(λ_0) ψ_0 ⟩ =0,

but the Herglotz nature of the operator implies that the operator
d M/dλ(λ_0) is positive definite. 



Next we specialize to the rational Sturm-Liouville problem 

    -(D(x) p_x)_x = λ W_0(x) p - ∑_i=1^N W_i(x)/λ - α_i
    - V(x) p = g(x, λ) p;    p_x(0) = 0 = p_x(L)


As motivation for the next result we note that in the special case in
which V(x),W_i(x) and D(x) are all constant we have that the
solution is given by 

    p(x) = cos√(%s/%s)g(λ)Dx
 
with the eigenvalue condition g(λ)=k^2 π^2 D. Since
g(λ) is Herglotz it is monotone and takes every real value for
λ∈ (α_i-1,α_i). Thus we have that for each such
interval (α_i-1,α_i) there is a sequence of
eigenfunctions p_k(x) which are indexed by k, the number of roots
in (α_i-1,α_i).   


Equation (<ref>) will typically not, of course, be
solvable in closed form but in the WKB approximation, which is valid
when g(x, λ) is large (which occurs for λ=α_I^-)
we have the approximate eigenfunction
 
    p_k(x,λ) ≈cos( ∫_0^x
      √(%s/%s)g(x, λ)D(x))/√(D(x)) (g(x, λ))^1/4

where the WKB quantization condition is that 

    ∫_0^L
    √(g(x, λ)/D(x))dx = k π.

 Since the
quantity  ∫_0^L
√(g(x, λ)/D(x))dx is a monotone
function of λ whenever g(x, λ) is positive it follows
that there is an increasing sequence of (approximate) eigenvalues
asymptotic to α_i^- in each interval (α_i-1,α_i).   

These heuristics motivate the following version of the
classical Sturm theorem.

  
Consider the rational Sturm Liouville pencil subject to Neumann boundary conditions

    -(D(x) p_x)_x = λ W_0(x) p - ∑_i=1^N W_i(x)/λ - α_i
    - V(x) p = g(x, λ) p,     
    p_x(a) = 0 = p_x(b)

with D(x),V(x),W_i(x) continuous on [a,b] and D(x),W_i(x) > 0 on
[a,b]. The quantities {α_i }_i=1^N will be
referred to as the singular values. For ease of notation we will let α_0=-∞ and
α_N+1=+∞, and will let I_j,   j ∈ 0… N denote the interval I_j =
(α_j,α_j+1). From the previous result the eigenvalues
are real and semi-simple. Additionally the following oscillation results hold.


  * Sturm Theorem  For each integer k≥ 0 there exists one
  eigenvalue λ_k^j in each interval I_j= (α_j,α_j+1)   j∈ 0 … N
  such that the corresponding eigenfunction  p_k^j(x) has exactly k zeroes
  in the open interval (a,b). 

  * By definition we have that λ_k^j >
  λ_k^j' if j>j'. Further we have that λ_k^j  >
  λ_k'^j if k>k'. 

  * Each α_i except α_0 is an accumulation point of
  eigenvalues, so α_i,   i ∈ 1 … N are in the
  essential spectrum.  



In other words we have a Sturm oscillation type theorem that holds in
every interval (α_i,α_i+1): the eigenvalues in each such
interval are increasing and  indexed by the number of roots of the
corresponding eigenfunction in (a,b).


We first note that the
fact that the pencil has finite accumulation points of eigenvalues
means that the original operator is neither compact nor does it have compact
resolvent. This illustrates that the problem in which some
species do not diffuse is a singular perturbation of the problem in which
those same species have a small diffusion coefficient, since in that case it
is straightforward to see that the operator has compact resolvent. 

We also note that the more general operator 

    -d/dx( D(x) u_x)  = λ u - ∑_i=1^NW_i(x)/λ - α_i(x)
has only real eigenvalues by the argument given previously but has essential spectrum for λ∈∪_i=1^N (α_i(x)). This type of operator has been
analyzed by Adamjan, Langer and Langer
<cit.>  for N=1, including construction
of the spectral measure corresponding to the essential spectrum. 
 
 

The proof hinges on the fact that the (eigenvalue dependent) potential
g(x, λ) has the property that 


  * g(x,λ) → -∞ as λ→α_i^+

  * g(x,λ) → +∞ as λ→α_i+1^-.

  * ∂ g/∂λ > 0,

which are essentially the conditions needed to apply the classical
Sturm oscillation theorem for a one dimensional second order
operator. For the convenience of the reader a proof is given in Appendix 
<ref>. 



 
To illustrate this theorem we consider the rational Sturm-Liouville
pencil

    -u_xx + sin(x) u = λ u - .2 + cos^2(x)/λ-2 u,
         u(0)=0=u(π)
 
which is equivalent to the system 

    -u_xx + sin(x) u + v = λ u;      u(0)=0=u(π) 
       α u + (.2 + cos^2(x)) u = λ v



There is one singular value, α_1=2, and two open intervals on
the real line for which H(λ) is defined,
(-∞,2) and (2,+∞). We solved the system
(<ref>) numerically as described in Appendix <ref>.

 Using the discretization described with n_x=100 (see Appendix <ref> for details), we will
 have 99 eigenvalues in the interval I_0 = (-∞,2), denoted
 λ^0_k, and 99 eigenvalues in the interval I_1 = (2,∞),
 denoted λ^1_k. The first set (λ^0_k) will have an
 accumulation point at α_1=2, the second set (λ^1_k) at
 α_2= ∞; this is shown in
 Fig. <ref>. The corresponding eigenfunctions we
 will index in the same way: i.e. as u^j_k for the kth eigenvalue
 in interval I_j. We compare these eigenvalues with the WKB
 prediction 1/π∫_0^π√(λ - sin(x) -.2
     + cos^2(x)/λ-2) dx = k, which we also evaluated
   numerically. 

The numerical values of
eigenvalues for k∈ 0 … 4 are:



λ Interval     k=0     k=1     k=2     k=3     k=4 
(-∞,2) (Numerical)     1.22     1.75     1.92     1.95     1.97 
(-∞,2) (WKB)     1.029     1.773     1.916     1.956     1.973 
(2,∞)  (Numerical)     2.59     4.88     9.65     16.53     25.41
(2,∞)  (WKB)     2.68     4.88     9.73     16.69     25.67









The corresponding eigenfunctions,


u^0_k and u^1_k, are shown in Fig.  <ref>. The left column shows the primary (u^0_k) and auxiliary eigenfunctions (v^0_k) for the interval I_0, while the right column shows u^1_k and v^1_k. Note that u^j_k has exactly k zeros, as expected, and that the nodal sets of u^j_k and v^j_k coincide for each j,k. The functions u^0_k and v^0_k are of opposite sign, since λ^0_k < 2, while u^1_k and v^1_k are same-signed.
























We next use
Eqn. (<ref>), (<ref>) to make the following prediction for the asymptotics of the
eigenvalues accumulating near λ=2:

    λ_k^0 ≈ 2- (1/π∫_0^π√(.2 +
          cos^2(x)))^2/k^2+ O(k^-3) ≈ 2-.649/k^2,

while the other branch of eigenvalues has the asymptotics

    λ_k^1 ≈ k^2 + O(1)

The plots in Figure <ref> depicts a graph of (λ_k^1-2)k^2
vs k for the eigenvalues in the region  (-∞,2) together with a plot
of  λ^2_k k^-2 vs. k for the eigenvalues in the region
(2, ∞). These are expected to asymptote to .649 and 1
respectively. This is repeated for three different values of Δ
x, the discretization size, in order to assess the role of
discretization error. In each case we see that the eigenvalue curve
hugs the asymptote for a while before falling away due to
discretization error, and that as the size of the discretization is
decreased the curve hugs the presumed asymptotic value for longer, as
one would expect. 








§ APPLICATION: SPATIAL MODELS IN EPIDEMIOLOGY


We now consider an example with detailed numerics. 





	Variable     Unit  

	S susceptible host density     foxes/km^2 

E exposed host density     foxes/km^2 

I infected hosts density     foxes/km^2 

t time     year 

Parameter     Value  

	a : average birth rate               1 per year 
 
b : average intrinsic death rate     0.5 per year 

1/α(x) : average duration of clinical disease     ≈ 5 days 

1/σ : average incubation period     28 days 

β(x) : disease transmission coefficient     ≈ 80 km^2 per year 
 
K : carrying capacity                      0.25 to 4 foxes per km^2 

D : diffusion coefficient     50  to 330 km^2 per year

	





In <cit.> the following one-dimensional model for the spread of rabies among foxes was considered.

    ∂ E/∂ t = β(x) I S - σ E - ( b + (a-b)N/K) E   
       ∂ I/∂ t=∂/∂ x( D(x) ∂ I/∂ x) + σ E - α(x) I -(b + (a-b)N/K)I   
       ∂ S/∂ t = (a-b) (1-N/K) S-β(x) I S   
       ∂ I/∂η|_∂Ω=0,

where S, E and I are the population densities of the susceptible, exposed and infectious foxes,
respectively and N=S+E+I is the total fox population. 
When susceptible hosts contact infectious ones, they become exposed at rate β(x).
Exposed hosts remain in that class for an average period of σ^-1 before they
transition into the infectious class. Hosts remain infectious for an average period 
of α(x)^-1 before they die. 
While infectious, hosts experience random wandering; this movement is modeled as diffusion with coefficient D(x).  
All host classes experience average birth rate a
and death rate b. They also compete for resources and so it is assumed that the
environmental carrying capacity is K.       
Typical parameter values are given in the table above. 



The authors <cit.> assumed that α,β, D are positive continuous functions 
depending on the spatial variable x.  In this case there is the disease-free
steady state (E=0, I=0, S=K). Linearizing about the steady state gives the 
eigenvalue problem 

    -(σ+a) E + β(x) K I = λ E   
       ∂/∂ x( D(x) ∂ I/∂ x) + σ E - (α(x)+a)I = λ I    
       ∂ I/∂ x  |_∂Ω = 0.

Algebraically eliminating E gives the eigenvalue pencil 

    ∂/∂ x(D(x) ∂ I/∂ x) -(α(x)+a) I = ( λ - σβ(x) K /λ + σ + a ) I = g(λ,x) I     ∂ I/∂ x |_∂Ω = 0.




Our first observation is that the function g(λ, x) = λ - σβ(x) K/λ + σ + a is a Herglotz function for all x as 
long as σ,β(x), K, a>0, and thus the eigenvalues of this problem are all real.

This model was later considered by Wang and Zhao<cit.>, who used the 
next-generation operator approach <cit.> to give a recipe for 
computing the reproduction number, R_0, and understanding its
dependence on parameters in the spatially heterogeneous problem. 
The basic reproductive number R_0 is an important epidemiological threshold,
since it determines whether the disease can invade the population
(with R_0>1 indicating that the disease-free steady state is unstable)
or not (with R_0<1 indicating that the disease-free steady state is locally asymptotically stable).  



Defining λ̃=-λ, we see that the eigenvalue problem directly above can be written as:

    -d/dx( D(x) d I/dx)  +(α(x)+a) I + σβ(x) K /λ̃ - (σ + a)     =   λ̃I ;   d I/dx|_∂Ω = 0

 to which we can immediately apply Theorem <ref> to see that this has a smallest real eigenvalue which is simple and which will govern stability of Eqn. (<ref>) (this is Lemma 4.1 from <cit.>). 
 
Wang and Zhao analyze R_0 through a related eigenvalue problem,



    -d/dx( D(x) dϕ/dx) + (α+a)ϕ = μσ K β(x)/σ+aϕ;

they show that the reproduction number (for Eq. (<ref>))
is the inverse of the principal eigenvalue of
Eqn. (<ref>): i.e. R_0 = 1/μ_1.
However, this requires us to solve a second eigenvalue problem, which is not obviously identical to the first.  

We next show that there is a relationship between λ and
R_0, which eliminates our need to solve the second eigenvalue
problem. The basic idea is to note that Eq (<ref>) with
λ̃=0 is the same as Eq (<ref>) with
μ=1. It follows from the Sturm oscillation theorem that if μ=1
is the k^th eigenvalue of (<ref>)  then
λ=0 is the  k^th eigenvalue of (<ref>), because the (identical) eigenfunctions must have the same number of zeros. Moreover, μ=1
lies between the k^th and (k+1)^st eigenvalues of
(<ref>) if and only if λ̃=0 lies between
the k^th and (k+1)^st eigenvalues of
(<ref>). In conclusion, solutions to (<ref>) are stable if the
lowest eigenvalue of (<ref>)  satisfies μ_1>1. 



We next reproduce some of Wang and Zhao's numerical results. One such finding concerns the impact of spatial heterogeneity; 
in particular that μ_1 ≤μ̅_1, where μ_1 is the principal eigenvalue of their associated eigenvalue problem, and μ̅_1 is the principal eigenvalue where β and α are replaced by their spatial averages (this is Lemma 4.4 in <cit.>). Since R_0 = 1/μ_1, this means that R_0 is minimized for the spatially homogeneous problem. 
We used the constant parameters: a = 0.0027, α=0.2, σ=0.0357, K=0.98. Then β(x) = 0.2192(1+c_1 cos(π x)).
Fig. <ref>A shows that R_0 is minimized when c_1=0. Similarly, restoring β = 0.2192 and varying α(x) = 0.2(1+c_2cos(π x)), R_0 is minimized for the spatially homogenous case.


We next look at the effect of varying D_0, the baseline diffusion coefficient, when spatial heterogeneity in β and/or α are present. 

To get a fair comparison for the role of heterogeneity, we used β(x) = 0.2192 (1+ c_1 cos(π x)) and α(x) = 2/π(1+ c_2 (sin(π x)-2/π)), so that [α] is preserved as c_2 changes; the variable diffusion is: D(x) = D_0(1 + cos(π x)).  Remarkably, D_0 has no effect on R_0 when both β and α are constant.

We can also vary the modulation magnitude of the diffusion, i.e. allow D(x) = 0.0685(1 + c_3 cos(π x)), where c_3 is varied. Wang found that again R_0 is constant when β and α are spatially homogeneous. We confirm this in Fig. <ref>C, again including heterogeneity in either β, or α, both, or neither.  

Wang and Zhao next considered the problem of designing a vaccine strategy, based on the hypothesis that the impact of the vaccine could be modeled by changes to disease transmission β(x).  They suppose that

    β(x) = 6x(1-x)

and consider how R_0 would change if we modified 

    β→β1/1+v_0(x)

where

    v_0 (x)  = {c_0/L    0 ≤ a_0 ≤ x < a_0+L ≤ 1
    
    						0     otherwise.

The authors interpret v_0/1+v_0 = 1-1/1+v_0 as the vaccine efficacy.  The spatial function can be motivated by the idea of distributing a total quantity of c_0 over an area of size L, beginning at a_0. 



To understand the impact of a vaccine strategy, we solved Eqn. (<ref>) by sweeping over the parameters c_0, a_0 and L. We surveyed 0 ≤ c_0, L ≤ 1; by symmetry of β(x) it suffices to consider 0 ≤ a_0 ≤ 0.5, and we limited our consideration to a_0 + L ≤ 1.  
Other parameters were as used previously: a = 0.0027, b=a/2, σ=0.0357, K=1.5, α(x) = α_0=0.2, D(x) = D_0 = 0.1371. 

First, we can ask how much vaccine is required to prevent an epidemic; i.e. how large does c_0 need to be for λ_0 < 0.  Not surprisingly, λ_0 decreases with c_0; the minimum value of c_0 for which stability was observed was 0.44.  
However, stability also depends strongly on where the vaccine is used, as well as its spatial distribution: in Fig. <ref>
we show λ_0 as a function of a_0 and L, for selected values of c_0: for c_0 ≥ 0.44 (three rightmost panels), the boundary of the stable region is shown as bold. Note that even for the largest vaccine quantity tested (c_0=0.95), a suboptimal distribution strategy can easily fail to contain the infection (i.e. λ_0 > 0).

Finally, we can examine the optimal vaccine strategy — the configuration that minimizes λ_0, given a total quantity c_0. 
For the values of c_0 shown in Fig. <ref>, we found the (a_0,L) values at which λ_0 was minimized, and show the corresponding vaccine strategy and its effect on the disease transmission rate (i.e. the modified β(x)) in Fig. <ref>. In all cases, the optimal strategy is symmetric; the vaccine should always be distributed on the area where β is largest. 





§ ACKNOWLEDGMENTS

The authors are grateful for support from the University of  Illinois Research 
Board through grant RB 17060 to Z.R., from the National Science
Foundation through DMS-1615418 to J.C.B.,  





§ PROOF OF THE STURM OSCILLATION THEOREM FOR RATIONAL HERGLOTZ
  PENCILS



We give a short proof of the Sturm property for the rational Herglotz
pencil problem 

    -∂_x(D(x)p_x) = ( -V(x) + ∑_i=1^N β_i(x)/λ - α_i
    + λ) p = g(x,λ) p      p_x(a) = 0 = p_x(b)
 
We take the convention that α_0=-∞ and
α_N+1=+∞. Then in each interval
(α_j-1,α_j),   j∈{1… N} there exists a strictly increasing sequence of (real)
eigenvalues {λ^j_k}_k=0^∞∈
(α_j-1,α_j) whose only accumulation point is α_j
such that the eigenfunction ϕ_k^j(x) has exactly k (simple)
roots in the interval (a,b). The exact boundary conditions used are 
not particularly important – the same basic arguments work with
Dirichlet, Robin or (with some additional modifications) periodic
boundary conditions – but we work with Neumann boundary conditions
since these are most appropriate for the models considered here. 

The easiest way to prove the Sturm theorem is via the Prüfer angle, a change
of variable for the second order linear differential equation defined
by 

    p(x,λ) = R(x,λ) cos(Θ(x,λ))
    
        p_x(x,λ) = R(x,λ) sin(Θ(x,λ))
    .
   
It is well-known (at least to those that well-know it) that if 
p(x,λ) satisfies the second order equation 

    -∂_x (D(x) p_x) = g(x,λ) p
 then 
the
Prüfer angle Θ(x,λ) satisfies the equation 

    Θ_x = - sin^2 Θ - g(x,λ) cos^2 Θ -D_x
      R/DsinΘcosΘ

while the radial function R(x,λ) satisfies the equation 

    R_x= (1-g/D) R sin(Θ) cos(Θ)-D_x/D R sin^2
    Θ.

We will mostly work with the angle Θ(x,λ): the radial
function R(x,λ) is important only insofar as equation (<ref>)
shows that R(x,λ) cannot vanish unless it
is identically zero, implying that the Prüfer coordinate change is
well-defined for all x. 

The basic idea of the proof is that the winding of the Prüfer angle
is a good measure of the winding of the eigenfunctions, and that this
can be understood via equations (<ref>, <ref>). 
Since R(x,λ) cannot vanish, the roots of p correspond to places where
Θ(x,λ)=k π while critical points of p correspond to places where 
Θ(x,λ)=(k+1/2)π. We will look in particular at the
solution to (<ref>) satisfying the initial condition
  Θ(a,λ)=0 (and, for concreteness sake, R(a,λ)=1), so 
that the corresponding p satisfies p(a,λ)=1,
p_x(a,λ)=0.
In this setting eigenvalues will correspond to values of λ=λ^*
such that Θ(b,λ^*)=k π. 

We note several facts. First note that, for fixed values of x the
function Θ(x,λ) is a decreasing function of λ. 
This follows from differentiating equation <ref> with
respect to λ, using the fact that d g/dλ>0 and
applying the Gronwall inequality. Next we note that, for λ
sufficiently close to α_j-1^+, Θ(x,λ) must
satisfy 0≤Θ(x,λ) < π/2, with
Θ(x,λ)=0 if and only if x=0. To see this note that we
can choose λ close enough to α_j-1 to guarantee that
g(x,λ)<0. In this case we have that Θ_x|_Θ=0>0
and Θ_x|_Θ=π/2<0. This guarantees that the
interval Θ∈ (0,π/2) is invariant in positive
x. Finally note that, if Θ = -(k+1/2)π then
Θ_x=-1, so the curve Θ(x,λ) always crosses the
line -(k+1/2)π transversely and in the same direction. 
   

The proof now follows from considering a homotopy in λ: we
begin with the curve Θ(x,λ) with λ sufficiently
close to α_j-1 so that Θ(x,λ) lies in the
interval (0,π/2). Thus no value of λ in this range
can be an eigenvalue since Θ(b,λ) is strictly between 0
and π/2, so we are below the spectrum.  

Next we consider increasing λ, which will decrease
Θ(b,λ). The first eigenvalue occurs when
Θ(b,λ)=0.
Call this value λ_0. While we do not know much about the shape
of the curve Θ(x,λ_0) we do know that it cannot cross the
line Θ=-π/2 since at any point where
Θ(x,λ)=-(k+1/2)π we have that Θ_x=-1. 
Hence as λ is increased the first time that the curve
Θ(x,λ) can cross the line Θ=-(k+1/2)π must be
at the endpoint x=b. 

This pattern continues as λ increases:
Θ(b,λ)
decreases and alternately crosses the lines Θ=-kπ and
Θ=-(k+1/2)π: the former correspond to eigenvalues and at
the latter p(x,λ) picks up a root. Since R(x,λ) is
non-vanishing and Θ and Θ_x cannot vanish simultaneously
it follows that the k^th eigenfunction has exactly k roots in
(a,b). It follows from standard estimate that away from poles of
g(x,λ) we have analytic dependence of Θ(b,λ) on
λ and so solutions to Θ(b,λ)=-(k+1/2)π 
must be isolated. It is also clear that by choosing λ
sufficiently close to α_j we can guarantee that Θ_x is
uniformly large and negative, and thus that α_j is an
accumulation point of eigenvalues. 



§ NUMERICAL METHODS



We use standard methods to solve the eigenvalue problem:

    H u + q(x) u + ∑_i=1^Nβ_i (x) u/λ - α_i    =    λ u;   0 ≤ x ≤ L

where H is a Sturm-Liouville operator, such as the second derivative (i.e. Hu = -u_xx).
We assume that β_i(x) > 0 on the domain [0, L] and that boundary conditions are of the typical kind:

    b^0_0 u(0) + b^1_0 ∂ u/∂ x(0)     =     0 
    
    b^0_L u(L) + b^1_L ∂ u/∂ x(L)     =     0

where b^0_0 ≠ 0 and/or b^1_0 ≠ 0, and  b^0_L ≠ 0 and/or b^1_L ≠ 0.
Introducing auxiliary functions v^(i) = β_i(x) u/λ-α_i, we find 
that Eqn. (<ref>) is equivalent to

    H u + q(x) u + ∑_i=1^N v^(i)    =    λ u;
    β_i (x) u + α_i v^(i)    =    λ v^(i),   i=1,...,N

We used standard second-order discretizations to evaluate H u at n_x-1 interior points of [0,L]; e.g. centered differencing for Hu = -u_xx. For the spatially variable diffusion in <ref>, we used the stencil

    (D(x) u_x)_x     ≈   1/(Δ x)^2 (D_i-1/2 u_i-1 - (D_i-1/2+D_i+1/2) u_i + D_i+1/2 u_i+1)

where x_i = i Δ x, u_i = u(x_i), D_i-1/2 = D(x_i-1/2), etc.; i.e. the diffusion function is evaluated at half-intervals.


Finally, the boundary conditions are incorporated with the second-order stencils 

    ∂ u/∂ x (0) ≈-3 u_0 + 4 u_1 - u_2/2 Δ x;     ∂ u/∂ x (L) ≈3 u_n_x - 4 u_n_x-1 + u_n_x-2/2 Δ x;

This results in a linear system of size (N+1) (n_x-1), representing the unknown functions u and v_i at each of n_x-1 interior points; we numerically approximate its eigenvalues using standard methods (we used the  function from Matlab R2017b).

Here, we list parameters used for specific cases in the paper.


  * In Example <ref>, we used L=π and n_x = 100; i.e. Δ x = π/n_x and the number of interior points is n_x - 1=99. For Figure <ref>, we repeated the computations with n_x = 200 and n_x=400.

  * For all computations presented in <ref>, L=1 and n_x = 200.
 
 








10

Adamjan.Langer.Langer.2001
V. Adamjan, H. Langer, and M. Langer.
A spectral theory for a λ-rational Sturm-Liouville
  problem.
J. Differential Equations, 171(2):315–345, 2001.

capasso1997
V. Capasso and R. E. Wilson.
Analysis of a reaction-diffusion system modeling man-environment-man
  epidemics.
SIAM Journal on Applied Mathematics, 57(2):327 – 346, 1997.

Haynsworth1962
M. P. Drazin and E. V. Haynsworth.
Criteria for the reality of matrix eigenvalues.
Math. Z., 78:449–452, 1962.

feng2013
Z. Feng, W. Huang, and D. L. DeAngelis.
Spatially heterogeneous invasion of toxic plant mediated by
  herbivory.
Mathematical Biosciences and Engineering, 10:1519 – 1538,
  2013.

Ghazaryan2010a
A. Ghazaryan, Y. Latushkin, and S. Schecter.
Stability of travelling waves for a class of reaction-diffusion
  systems that arise in chemical reaction models.
SIAM Journal on Mathematical Analysis, 42:2434–2472, 2010.

Ghazaryan2010b
Y. Ghazaryan, A. andLatushkin, S. Schecter, and A. J. de Souza.
Stability of gasless combustion fronts in one-dimensional solids.
Archive for Rational Mechanics and Analysis, 198(3):981–1030,
  Dec 2010.

kapitula2012
T. Kapitula and K. Promislow.
Stability indices for constrained self-adjoint operators.
Proceedings of the American Mathematical Society, 140:865–880,
  2012.

KaupNewell
D. J. Kaup and A. C. Newell.
An exact solution for a derivative nonlinear schrödinger equation.
Journal of Mathematical Physics, 19(4):798–801, 1978.

kollar2014
R. Kollar and P. D. Miller.
Graphical krein signature theory and evans-krein functions.
SIAM Review, 56:73–123, 2014.

lewis1994
M.A. Lewis.
Spatial coupling of plant and herbivore dynamics: the contribution of
  herbivore dispersal to transient and persistent waves of damage.
Theoretical Population Biology, 45:277–312, 1994.

lou2004
Y. Lou, Q. Nie, and F. Y. M. Wan.
Nonlinear eigenvalue problems in the stability analysis of morphogen
  gradients.
Studies in Applied Mathematics, 113:183–215, 2004.

morris1997
W. F. Morris and G. Dwyer.
Population consequences of constitutive and inducible plant
  resistance: herbivore spatial spread.
The American Naturalist, 149:1071–1090, 1997.

murray1986
J. D. Murray, E. A. Stanley, and D. L. Brown.
On the spatial spread of rabies among foxes.
Proceedings of the Royal Society of London B, 229:111–159,
  1986.

sanchez2011
F. Sánchez-Garduno and V. F. Brena-Medina.
Seaching for spatial patterns in a pollinator-plant-herbivore
  mathematical model.
Bulletin of Mathematical Biology, 73:1118–1153, 2011.

Simon
B. Simon.
Orthogonal polynomials on the unit circle.
American Mathematical Soc., 2005.

Takhtajan
L. A. Takhtajan and L. D. Faddeev.
Essentially nonlinear one-dimensional model of classical field
  theory.
Theoretical and Mathematical Physics, 21(2):1046–1057, Nov
  1974.

vanden2002
P. van den Driessche and J. Watmough.
Reproduction numbers and sub-threshold endemic equilibria for
  compartmental model of disease transmission.
Mathematical Biosciences, 180:29–48, 2002.

wang2012
W. Wang and X.-Q. Zhao.
Basic reproduction numbers for reaction-diffusion epidemic models.
SIAM Journal of Applied Dynamical Systems, 11(4):1652–1973,
  2012.



