Temporal oscillations in Becker-Doering equations with atomization

TL;DR

This paper demonstrates that a modified Becker-Doering coagulation-fragmentation system with atomization exhibits time-periodic solutions arising from Hopf bifurcation, modeling oscillatory behavior in physical chemistry gas evolution processes.

Contribution

It introduces a novel variant of the Becker-Doering equations incorporating atomization and proves the existence of oscillatory solutions via bifurcation analysis.

Findings

Time-periodic solutions are proven to exist.

Oscillations arise through Hopf bifurcation.

Model captures physical chemistry gas evolution oscillations.

Abstract

We prove that time-periodic solutions arise via Hopf bifurcation in a finite closed system of coagulation-fragmentation equations. The system we treat is a variant of the Becker-Doering equations, in which clusters grow or shrink by addition or deletion of monomers. To this is added a linear atomization reaction for clusters of maximum size. The structure of the system is motivated by models of gas evolution oscillators in physical chemistry, which exhibit temporal oscillations under certain input/output conditions.

Tables

Table 1. Table 1: Critical parameters for first Hopf bifurcation

$M$	$K$	${\kappa }_1=K\sqrt{M}$	$Im\lambda$
${10}^2$	0.39349	3.9349	0.021740
${10}^3$	0.075016	2.3722	3.6176e-4
${10}^4$	0.020376	2.0376	9.3596e-6
${10}^5$	0.0061392	1.9414	2.7777e-7
${10}^6$	0.0019118	1.9118	8.6091e-9

Equations

$$(j)+(k)\stackrel{a_{j,k}}{⟶}(j+k)\text{mbox}(aggregation),$$

$$(j)+(k)\stackrel{b_{j,k}}{⟵}(j+k)\text{mbox}(binarybreakup),$$

$$R_{j,k}=a_{j,k}{n}_j{n}_k-b_{j,k}{n}_{j+k}.$$

$${\partial }_t{n}_j=\frac{1}{2}\sum _{k=1}^{j-1}{R}_{j-k,k}-\sum _{k=1}^{\infty }{R}_{j,k},j=1,2,\dots$$

$$J_{\ell }=n_{\ell }{n}_1-n_{\ell +1},\text{mbox}for1\le \ell \le N,$$

$${\partial }_t{n}_{\ell }$$

$${\partial }_t{n}_M$$

$${\partial }_t{n}_1$$

$${\partial }_t\left(\sum _{\ell =1}^M\ell n_{\ell }\right)=0.$$

$$a_{i,j}=(i/j{)}^{\alpha }+(j/i{)}^{\alpha }$$

$${\partial }_t{n}_1$$

$${\partial }_t{n}_k$$

$${\partial }_t{N}_j=q_{j-1}{N}_{j-1}-(q_j+k_j)N_j,$$

$${\partial }_{t}x(t)=-ax(t-\tau )-bx(t),$$

$$\rho (\lambda ):=a{e}^{-\lambda \tau }+b+\lambda =0.$$

$$\rho \circ \gamma (s)=a{e}^{is\tau }+b-is,$$

$$\frac{d}{ds}\rho \circ \gamma (s)=i\tau (\rho \circ \gamma (s)-b+is)-i$$

$$M=25,\quad K=3,\quad n_{1}=4.2,\quad n_{\ell}=1+K=4\quad\mbox{ for $\ell\geq 2$.}$$

$${\partial }_t{n}_j=(j-1)n_1{n}_{j-1}-j{n}_1{n}_j-Kj{n}_j,\text{mbox}forj\ge 2.$$

$$n_j=\frac{e^{-\beta j}}{j},e^{-\beta }=\frac{n_1}{n_1+K}.$$

$${\partial }_t{n}_j=\sum _{k=1}^{j-1}(j-k)k{n}_{j-k}{n}_k-2\sum _{k=1}^{\infty }jk{n}_j{n}_k-Kj{n}_j,j\ge 2.$$

$${\bar{n}}_{\ell }=z(1-\alpha )+z^{\ell }\alpha ,\text{mbox}whereJ=(z^2-z)(1-\alpha ).$$

$$0=J-K{\bar{n}}_M=(z^2-z-Kz)(1-\alpha )-K{z}^M\alpha .$$

$$\alpha =\frac{z-1-K}{K{z}^N+z-1-K},$$

$${\bar{n}}_{\ell }=\frac{K{z}^M+z^{\ell }(z-1-K)}{K{z}^N+z-1-K},\ell =1,\dots ,M.$$

$${\bar{n}}_{\ell }=\frac{1+(M-\ell )K}{1+NK},\ell =1,\dots ,M.$$

$$m=\sum _{\ell =1}^M\ell n_{\ell }=\alpha {\mu }_M(z)+(1-\alpha )z{\mu }_M(1)$$

$${\mu }_M(z)=\sum _{\ell =1}^M\ell z^{\ell }=z\frac{d}{dz}\frac{1-z^{M+1}}{1-z}=\frac{z-z^{M+1}}{(1-z{)}^2}-\frac{M{z}^{M+1}}{1-z}.$$

$${\bar{n}}_{\ell }=A:=1+K,1\le \ell \le M.$$

$$n_k=A+v_k,1\le \ell \le M.$$

$$L_{\ell }=A{v}_1+A{v}_{\ell }-v_{\ell +1},1\le \ell \le N=M-1,$$

$${\partial }_t{v}_{\ell }$$

Full text

1-Department of Mathematics and Center for Nonlinear Analysis

Carnegie Mellon University, Pittsburgh, Pennsylvania, PA 12513, USA

email: rpego@cmu.edu

2- Institut für Angewandte Mathematik

Universität Bonn

Endenicher Allee 60

53115 Bonn, Germany

email: velazquez@iam.uni-bonn.de

*Keywords. *bubbling oscillator, shattering, bubbelator, time periodic solution

*Mathematics Subject Classification. * 34C23, 34A34, 82C99

1 Introduction

Coagulation-fragmentation equations are commonly used to model particle size distributions in a wide range of scientific and technological applications. These equations model binary reactions of clusters of size $j$ with clusters of size $k$ as indicated schematically by

[TABLE]

With rate coefficients $a_{j,k}$ for aggregation and $b_{j,k}$ for breakup, the net rate of this binary reaction is modeled by the law of mass action to be

[TABLE]

The coagulation-fragmentation equations accounting for the gain and loss rates for the number density $n_{j}(t)$ of groups of size $j$ then take the form

[TABLE]

To date, mathematical investigations of the dynamic behavior of solutions have largely focused on questions of convergence to equilibrium and the phenomenon of gelation, in which mass conservation fails (either in finite or infinite time) due to a flux to infinite size. We refer to classic work of Aizenman and Bak [1] who established an $H$-theorem for perhaps the simplest coagulation-fragmentation model with constant rate coefficients, and Ball, Carr and Penrose [2] for the first analysis of (infinite-time) gelation in the Becker-Döring equations. If fragmentation is weak, finite-time gelation can occur [12, 13, 19, 31] as it does for the case of pure coagulation about which there is now an extensive literature.

Regarding convergence to equilibrium, entropy methods have been effectively used to study general classes of coagulation-fragmentation equations that admit equilibria in detailed balance, meaning that $R_{j,k}=0$ for each individual reaction in the system, so the forward and backward reaction rates match. See work of Laurençot and Mischler [21] for the continuous-size case and Cañizo [7] for the discrete-size case. More recent studies of equilibration have examined rates of convergence and their relation to entropy-dissipation relations [18, 6, 27, 28, 5].

In the absence of detailed balance, however, one does not expect that an $H$-theorem always holds, and it is not clear whether the structure of coagulation-fragmentation reaction networks means that solutions necessarily always converge to some equilibrium. Sometimes it is indeed the case, as in cases when coagulation is weak [14] or for special systems that can be studied globally using transform methods, as in [9]. In [22], Laurençot and van Roessel analyzed a model with a critical balance of coagulation and fragmentation rates, and used transform methods to show that infinite-time gelation emerges through self-similar growth.

On the other hand, in studies of pure coagulation without fragmentation, the usual expectation of self-similar growth has sometimes been shown not to occur. For special rate kernels, solutions with fat tails are known to be capable of periodic and even chaotic behavior after rescaling [26]. Temporal oscillations can persist after rescaling without fat tails for Smoluchowski equations with diagonal rate kernel [20].

Our goal in the present work is to demonstrate that persistent oscillations in time are possible in a simple discrete-size coagulation-fragmentation model, by proving that Hopf bifurcations occur.

The particular system that we study is a modified system of Becker-Döring equations. (For a nice historical review of mathematical developments concerning the Becker-Döring equations, see [17].) As usual for Becker-Döring equations, we suppose that the coagulation of clusters of size $\ell$ with monomers proceeds at the rate $a_{\ell}n_{\ell}n_{1}$, and clusters of size $\ell+1$ lose monomers at the rate $b_{\ell+1}n_{\ell+1}$. We take these rates to apply only for a finite range of sizes $1\leq\ell\leq N$, however, and consider only the simplest case, always taking $a_{\ell}=b_{\ell+1}=1$. Thus the net flux of clusters from size $\ell$ to $\ell+1$ is $J_{\ell}=R_{1,\ell}$, as given by

[TABLE]

We suppose further that $M=N+1$ is the size of the largest clusters in the system, and these are also subject to a linear atomization reaction that converts an $M$-cluster into $M$ monomers and proceeds at rate $Kn_{M}$. Thus the governing equations take the following form:

[TABLE]

All solutions of the system (2)–(4) conserve mass, since

[TABLE]

A formal continuum analog of this system will be studied for illustrative purposes in Section 4.

2 Background and motivation

Model with nonlinear atomization. In the physical literature, recent work of Matveev et al. [25] and Brilliantov et al. [4] has identified a coagulation-fragmentation model with a different, nonlinear atomization mechanism that exhibits persistent temporal oscillations in numerical simulations. In this model, aggregation of clusters of size $i$ and $j$ proceeds at rate $a_{i,j}n_{i}n_{j}$ where

[TABLE]

and pairs of such clusters atomize upon collision into $i+j$ monomers at rate $\lambda a_{i,j}n_{i}n_{j}$. In total, the rate equations in [25] take the form

[TABLE]

This system has the feature that interactions between large clusters of similar size appear to be dominated by interactions between large clusters and small ones (for which either $i/j$ or $j/i$ is large). Oscillations are found for $\frac{1}{2}<\alpha\leq 1$ and small $\lambda>0$. Though the numerics is convincing, to our knowledge there is no proof yet that temporal oscillations persist in this system.

Bubbling oscillators. Our motivation for studying the system (2)–(4) comes from literature in physical chemistry concerning bubbling oscillators (often called ‘gas evolution oscillators’ in much of the literature). In these systems, dissolved gas (such as CO or CO2) is added slowly to a liquid solution, producing a super-saturated mixture. At some time, nucleation of gas bubbles occurs spontaneously and the bubbles grow rapidly and carry most of the dissolved gas out of the system. The first system of this kind was reported by J. S. Morgan in 1916, who found that a small concentration of formic acid mixed in sulphuric acid produced periodic bursts of carbon monoxide. Such systems were the subject of part of an extensive series of quantitative studies by R. M. Noyes and collaborators concerning chemical oscillators, including some of the original studies of chemical oscillators such as the BZ reaction and the Oregonator. Regarding gas evolution oscillators, we especially refer to [30, 32, 3]. The phenomenon of sudden outgassing of CO2 after slow buildup of supersaturation was responsible for the 1986 Lake Nyos disaster in Cameroon, which killed more than 1700 people.

In the work of Yuan, Ruoff and Noyes [32], this process was simulated numerically by grouping bubble sizes into a finite set corresponding to exponentially spaced radii $r_{j}$, and writing rate equations to model the number density $N_{j}$ of bubbles of size $r_{j}$. A key equation when $r_{j}$ greater than a critical value $r_{\rm eq}$ is

[TABLE]

where the coefficients $q_{j}$ are proportional to bubble growth rate and the $k_{j}$ are rate constants for escape. This resembles a linearized Becker-Döring equation or a discretized advection equation, and models the process of free bubble growth and escape. With $M=60$ size classes, numerical simulations in [32] exhibit temporal oscillations for a range of parameters designed to model experimental conditions.

Bar-Eli and Noyes [3] later devised a simplified, qualitative model for bubbling oscillators that involves a nonlinear differential-delay equation for the concentration of dissolved gas. When linearized about a constant steady-state, one obtains a constant-coefficient linear DDE of the form

[TABLE]

where the parameters $a$, $b$ and the delay time $\tau$ are positive constants. Whenever $a>b$, one finds there is an oscillatory transition from stability to instability as $\tau$ increases.

We sketch loosely how one can see this mathematically. (A detailed analysis of (8) can be found in work of Hadeler and Tomiuk [16].) Equation (8) has solution $e^{\lambda t}$ provided

[TABLE]

For $\tau=0$, naturally $\lambda<0$, and moreover $\lambda=0$ is never possible for any $\tau$. But for $a>b$ and $\tau$ sufficiently large, there are solutions with $\mathop{\rm Re}\nolimits\lambda>0$. To show this is so, one can consider the winding number around 0 of a curve $\rho\circ\gamma$, where $\gamma$ is a concatenation of a path $s\mapsto-is$ for $s\in[-R,R]$ and a path in the right half plane along the semicircle where $|\gamma|=R>a+b$. Along the semicircle, $\rho\circ\gamma$ can never cross the negative real axis $\mathbb{R}_{-}$. Along the imaginary axis, however,

[TABLE]

and this does cross $\mathbb{R}_{-}$ for $s$ between [math] and $2\pi/\tau$, if $\tau$ is large enough. Moreover, $\rho\circ\gamma(s)$ can only ever cross $\mathbb{R}_{-}$ going from the second quadrant to the third, since whenever $\rho\circ\gamma(s)<0$,

[TABLE]

and this has negative imaginary part. Consequently, the winding number of $\rho\circ\gamma$ around 0 is positive if $\tau$ is large enough, and this implies (9) has a root $\lambda$ inside $\gamma$.

Becker-Döring with linear atomization (our model). Now, the rough idea behind our model (2)–(4) is that the Becker-Döring equations involve a well-known advection mechanism that transports mass from small cluster sizes to large ones when the monomer concentration is supercritical. The atomization reaction added in (3) couples the advected wave back to the monomer concentration after a time delay that depends on the size of the system. Luckily enough, we find that for large $M$ there indeed is an oscillatory transition to instability as the parameter $K$ varies, in a certain parameter range where $K$ is small but $KM$ remains large. See Figure 1, where we plot the monomer concentration vs. time for a numerically computed solution of (2)–(4) with parameters and initial values given by

[TABLE]

Other models with linear atomization. Finally, we mention two other kinds of merging-splitting models involving a linear atomization reaction that have appeared in the literature. Alongside discussion of Niwa’s model [29] for animal group size, Ma et al. [24] described a “preferential attachment” model, which takes the form

[TABLE]

This model admits a simple logarithmic distribution in equilibrium, of the form

[TABLE]

(This is roughly similar to the distribution Niwa found to be a good description of empirical data on school size for pelagic fish.) A model of herd behavior by networks of colluding agents in financial markets introduced by Eguiluz & Zimmermann [11] takes the form

[TABLE]

D’Hulst and Rodgers [10] found a formula for equilibrium solutions of this model by use of generating functions. But as far as we are aware, no analysis of dynamics has been carried out for either of these models.

3 Equilibria, linearization, and main result

In this section, we find the general equilibrium solutions of the model (2)–(4), describe the special family of constant equilibria, and state our main rigorous result on the existence of Hopf bifurcations from this family, occurring at particular values of $K$, for large enough $M$.

3.1 General equilibria

We find the general equilibria as follows. In equilibrium, due to (2) the fluxes $J_{\ell}$ are all equal to the same value $J$ for $\ell=1,\ldots,N$, so the equilibrium number densities $\bar{n}_{\ell}$ satisfy the difference equation $\bar{n}_{\ell+1}=z\bar{n}_{\ell}-J$, when we require $\bar{n}_{1}=z$. For $z\neq 1$, the solution takes the form

[TABLE]

To obtain an equilibrium it remains to require that (3) hold, i.e.,

[TABLE]

Then it follows (recall $N=M-1$)

[TABLE]

and

[TABLE]

Note that (4) then holds also. For every $z>0$, $N>1$ and $K>0$, such an equilibrium exists and is positive. In case $z=1$, one finds directly that

[TABLE]

The total mass as a function of $z$ and $K$ is now

[TABLE]

where

[TABLE]

3.2 Linearization at constant equilibria

Particularly convenient for our analysis is the special family of equilibria that have constant densities, corresponding to $\alpha=0$. By (14) these take the form

[TABLE]

Corresponding to $K>0$ we require $A>1$. We will study the linearization of the system (2)–(4) about this equilibrium. We write:

[TABLE]

The linearized fluxes take the form

[TABLE]

and the linearized evolution equations are written as follows:

[TABLE]

Equivalently, after some computations, the system takes the more explicit form

[TABLE]

Equation (21) yields a combination of diffusion and transport. It is not able to yield oscillatory behavior by itself, but this will be generated through the ‘boundary conditions,’ or more precisely the equations with $\ell=M$ and $\ell=1.$

Looking for solutions of this system with the form

[TABLE]

leads to the eigenvalue problem (recalling $A=1+K$ and $M=N+1$)

[TABLE]

This system takes the form of an eigenvalue problem $BV=\lambda V$ for a vector $V=(V_{1},\ldots,V_{M})^{T}\in\mathbb{C}^{M}$, with $M\times M$ matrix $B$ having the structure

[TABLE]

Our goal is to understand the spectrum of $B$ in some detail, and eventually show that in a certain parameter range, some pair of complex eigenvalues of $B$ crosses the imaginary axis, and this forces the system (2)–(4) to undergo a Hopf bifurcation.

To begin, one can check that for any eigenvector $V=(V_{1},\ldots,V_{M})^{T}$ corresponding to a nonzero eigenvalue $\lambda\neq 0$, the mass conservation condition holds:

[TABLE]

This is due to the fact that $\bm{\ell}=(1,2,\ldots,M)$ is a left null vector of $B$. Thus $\lambda=0$ is an eigenvalue of $B$. If we represent the equations of equilibrium from (2)–(4) in vector form as ${\mathcal{F}}(\bar{n})=0$, then the matrix $B=\partial{\mathcal{F}}/\partial n$ evaluated at $\bar{n}$ with constant components $A$. Thus by differentiation it naturally follows from (15) that a right null vector $\bar{v}$ satisfying $B\bar{v}=0$ is given by

[TABLE]

In fact, we have the following.

Lemma 3.1**.**

For all $N>1$ and $K>0$, $\lambda=0$ is a simple eigenvalue of $B$.

Proof.

First, we show the null space of $B$ is one-dimensional. Whenever $BV=0$, the fluxes defined by

[TABLE]

must all take the same value due to (18), and for $V=\bar{v}$ this value is $K\bar{v}_{M}>0$ due to (19). If $BV=0$, then we can replace $V$ by a linear combination with $\bar{v}$ to make all fluxes $L_{\ell}=0$. But by (30) it follows $V_{2}=2AV_{1}$ and, by induction, $V_{\ell}=A_{\ell}V_{1}$ with $A_{\ell}>0$ for $\ell=2,\ldots,M$. Since $0=L_{N}=KV_{M}$, the only vector $V$ making all the fluxes vanish is $V=0$. It follows that $\bar{v}$ spans the null space of $B$.

Next, we claim there is no generalized eigenvector $V$ satisfying $BV=\bar{v}$. The reason is that, because $\bar{v}_{\ell}\geq 1>0$ for all $\ell$ and $\bm{\ell}$ is a left null vector of $B$, we would obtain a contradiction, via

[TABLE]

Thus the eigenvalue $\lambda=0$ has algebraic multiplicity one, so it is simple. ∎

3.3 Main results

If $\lambda$ is an eigenvalue of $B$, we say $\lambda$ is unstable if $\mathop{\rm Re}\nolimits\lambda>0$. We find that we can show the matrix $B=B(K,M)$ has unstable eigenvalues when $M$ is sufficiently large and $K$ is small but not too small, in a range proportional to $1/\sqrt{M}$. These eigenvalues are characterized as follows. It is convenient to state our results in terms of the parameter

[TABLE]

in place of $K=\kappa/\sqrt{M}$.

Theorem 3.2**.**

For each $k\in\mathbb{N}$ and $\beta_{0}\in(0,1)$, there exists $\beta_{k}>\beta_{0}$, and positive constants $M_{0,k}$, $\hat{C}_{k}$, such that for each $M>M_{0,k}$, the following hold:

If $\beta_{0}<\kappa<\beta_{k}$, then any unstable eigenvalue of $B$ is non-real and simple, and satisfies

[TABLE] 2. 2.

There are numbers $\kappa_{j}=\kappa_{j}(M)$ for $j=1,\ldots,k$, satisfying

[TABLE]

such that:

1. (a)

If $\kappa_{j}<\kappa<\kappa_{j+1}$ ($j=0,\ldots,k$), then $B$ has exactly $j$ complex-conjugate pairs $(\lambda,\bar{\lambda})$ of unstable eigenvalues. 2. (b)

There are analytic curves $\lambda_{j}:[\beta_{0},\beta_{k}]\to\mathbb{C}$, $j=1,\ldots,k$, such that $\lambda_{j}(\kappa)$ is an eigenvalue of $B$ that satisfies

[TABLE]

along with

[TABLE]

By the properties stated in part 2 of this theorem, the matrix $B$ has a unique pair of nonzero, purely imaginary eigenvalues $\pm\lambda_{j}(\kappa_{j})$ when $M$ is large and $\kappa=\kappa_{j}(M)$, and these cross transversely into the right half plane as $K$ increases.

The simple eigenvalue at zero, described in Lemma 3.1, is nominally an obstruction to applying the standard Hopf bifurcation theorem at this point. This eigenvalue is easily removed, however, by considering the dynamics of the nonlinear system (2)–(4) restricted to the invariant affine hyperplane determined by conservation of mass, i.e., the hyperplane where

[TABLE]

with $\bar{n}=(\bar{n}_{\ell})$ being the constant equilibrium state from (17). Within this hyperplane, the linearization of the system (2)–(4) is restricted to orthogonal complement of the left null vector $\bm{\ell}$ of $B$. In this subspace, the zero eigenvalue is removed, and the standard Hopf bifurcation theorem can be applied to yield the following result. (For a proof of the Hopf bifurcation theorem see [8, pp. 98–99]. For further discussion also see [15, pp. 150–152].)

Theorem 3.3**.**

Let $k\in\mathbb{N}$ and suppose $M>M_{0,k}$ as given by Theorem 3.2. Then for each $j=1,\ldots,k$, the system (2)–(4) admits a Hopf bifurcation as the bifurcation parameter $\kappa$ passes through $\kappa_{j}=\kappa_{j}(M)$. Thus a time-periodic solution exists for some value of $\kappa$ with $|\kappa-\kappa_{j}|$ small.

We have not managed to determine analytically whether the bifurcating solutions are stable (the supercritical case) or not. Many of our numerical computations, as in Fig. 1, are consistent with the presence of stable periodic solutions, however.

In Figure 2 we illustrate the location of the complex eigenvalues of $B$ computed numerically for the parameter values $M=100$ and $K=3$. The unstable eigenvalues shown correspond to values

[TABLE]

Besides the real eigenvalue $\lambda=1$ this matrix also has a large negative eigenvalue $\lambda\approx-410.94$. There are 49 complex-conjugate pairs of eigenvalues that lie close to an ellipse that we will describe formally in Remark 5.1 and Section 6. The presence of eigenvalues spaced closely along a smooth curve and a large isolated eigenvalue is reminiscent of the structure of the spectrum of differential-delay systems with large delay, see [23].

The real parts of eigenvectors for the first 3 complex eigenvalues closest to $\lambda=1$ are plotted in Figure 3. They appear to have a “smooth” structure except in a boundary layer near $\ell=M$.

In Table 1 we tabulate for various values of $M$ numerically computed critical values of $K$ that correspond to $\kappa_{1}$, the value at which the first pair of complex-conjugate eigenvalues crosses the imaginary axis. Eigenvalues were obtained by solving the equation (77) in Proposition 5.3 below using an iteration method. The first two rows were computed also by finding all eigenvalues of $B$ using the julia function eigen. The values of $\kappa_{1}$ in the third column can be compared to the value $\kappa_{\rm cr}\approx 1.89825$ described below in (92). This value is proved in Section 10 to be the limiting value of $\kappa_{1}$ as $M\to\infty$, see (152).

4 A formal continuum approximation

Here we describe and heuristically analyze a formal continuum approximation of the system (2)–(4) that we are studying. This will serve to preview how the proof of Theorem 3.2 will proceed, and also provides an approximate understanding of the origin of the oscillatory instability in the system and the smooth structure of the eigenvectors as shown in Figure 3.

We introduce a scaled atomization rate and scaled space and time variables via

[TABLE]

With these relations, we write a continuum approximation of (2)–(4) in terms of variables

[TABLE]

as follows: The evolution equations for $\partial_{t}n_{\ell}$, $\ell=2,\ldots,M$ in (2)–(3) are approximated by the PDE

[TABLE]

where the number flux relation $J_{\ell}=(n_{1}-1)n_{\ell}+n_{\ell}-n_{\ell+1}$ from (1) is approximated by

[TABLE]

Taking (2) to hold also for $\ell=M$ with $J_{M}=Kn_{M}$, the following boundary condition replaces (3):

[TABLE]

We replace the evolution equation for $\partial_{t}n_{1}$ by an equation equivalent to mass conservation, which says

[TABLE]

After integrating by parts, we require

[TABLE]

The system of equations (36)–(39) has the constant equilibrium $u(x)=A=1+K$ just like the discrete system. Using a superposed dot to denote differentiation with respect to a variational parameter leads to the linearized system

[TABLE]

Then looking for solutions of the form $e^{\Lambda\tau}v(x)$ leads to the eigenvalue problem

[TABLE]

The steps that we now take to analyze this continuum eigenvalue problem parallel the steps that we will take to analyze the discrete eigenvalue problem (24)–(26). First, we describe the solutions of (44) as linear combinations of solutions of the form $e^{-zx}$ where

[TABLE]

It is convenient to write a general solution of (44) in terms of the two solutions $z_{\pm}$ of (47) as

[TABLE]

Then a nonzero solution can satisfy the boundary conditions (45) and (46) if and only if

[TABLE]

[TABLE]

So far this is an exact treatment of the eigenvalue problem (44)–(46). But now we approximate, noting that the two solutions $z_{\pm}$ of (47) satisfy

[TABLE]

for small $\varepsilon$. For $\kappa>0$ and complex $z_{+}$ both of order $O(1)$, we neglect the exponentially small term $e^{z_{-}}=O(e^{-\kappa/\varepsilon})$ and keep only the leading order terms in the other entries of the matrix in (48), expressing $z_{-}$ in terms of $z_{+}$ using the relations in (50). Thus, writing $z=z_{+}$ we make the approximations

[TABLE]

Multiplying the second row of (48) by $z/\kappa$ and dividing the second column by $-\varepsilon\kappa$, we find that the equation determining eigenvalues approximately takes the form $Q(z;\kappa)=0$, where

[TABLE]

It is exactly this function $Q$ that we will find responsible for the appearance of unstable eigenvalues for the discrete problem (24)–(26) in the limit of large $M$. We analyze the complex roots of $Q$ in depth in Section 7. It turns out that a complex-conjugate pair of roots $z$ emerges into the right half plane $\mathop{\rm Re}\nolimits z>0$ as $\kappa$ increases past each value $\kappa=\kappa_{j}^{0}$ of an infinite sequence. The values $\kappa_{j}^{0}$ will be seen to be the limiting values of $\kappa_{j}(M)$ as appear in the statement of Theorem 3.2, in the limit $M\to\infty$.

In the remainder of this paper, we carry out the proof of Theorem 3.2 by performing an analogous analysis for the discrete eigenvalue problem, including rigorous estimates for all the error terms. For brevity’s sake, we forgo the formulation and rigorous demonstration of results analogous to Theorems 3.2 and 3.3 for the (parabolic) continuum model (36)–(39) described in this section. It should be evident from our analysis, though, that Hopf bifurcation occurs for this model in a similar way.

5 Reformulation of the eigenvalue equation

5.1 The difference equation

The eigenvalue equations (25) for $2\leq\ell\leq N$ comprise a family of second order difference equations. These difference equations have solutions of the form

[TABLE]

whenever

[TABLE]

which we can rewrite using $A=K+1$ as

[TABLE]

or as

[TABLE]

We take decreasing powers in (54) for reasons of scaling explained below.

We can then “connect” the values of $V_{1}$ and $V_{M}$ by means of a transition matrix depending on two constants (for each value of $\lambda$). More precisely, any solution of (25) takes the form

[TABLE]

whenever $\varphi_{1}$ and $\varphi_{2}$ are distinct roots of (57). Evidently the two roots are always related by $\varphi_{1}\varphi_{2}=1/A$, and for $\lambda=0$ the roots are $\varphi_{1}=1$ and $\varphi_{2}=1/A$.

The roots are distinct except when $\varphi_{1}=\pm A^{-1/2}$, which corresponds to

[TABLE]

For small $K>0$, we note that this becomes

[TABLE]

The roots $\varphi_{1}$, $\varphi_{2}$ are naturally functions of $\lambda$. However, it will be more convenient to recast the eigenvalue equations in terms of the variable $\varphi$ and regard $\lambda$ as a function of $\varphi$, given by the following equation equivalent to (56):

[TABLE]

Except when $\varphi=\pm A^{-1/2}$ (which will generate spurious roots below), corresponding to $M$ eigenvalues $\lambda$ there should exist $2M$ roots $\varphi$ of the relevant equations, which occur in pairs $\varphi$, $1/(A\varphi)$ that produce the same $\lambda$.

Remark 5.1*.*

We note that by (56), values of $\varphi$ on the unit circle, with $\varphi=e^{is}$ for $s$ real, produce values of $\lambda$ on an ellipse with

[TABLE]

This ellipse lies in the left half plane and passes through $\lambda=0$. In numerical computations such as those reported in Fig. 2, almost all the eigenvalues lie near this ellipse. By consequence we will expect to find most roots satisfying $|\varphi_{1}|\approx 1$ and $|\varphi_{2}|\approx 1/A<1$, with $|\varphi_{2}^{M}|$ extremely small. (This is the basic reason for the form we took in (54).) The possibility of transition to instability will depend upon the deviation of roots $\varphi_{1}$ from this ellipse in the vicinity $\varphi\approx 1$ where $\lambda\approx 0$.

5.2 Reduction to a $2\times 2$ determinant

We now use the expression (58) to write the “boundary conditions” for $V$, that correspond to the equations for $V_{M}$ and $V_{1}$ in (26) and (24) respectively. Using the fact that (56) holds for both $\varphi_{1}$ and $\varphi_{2}$, after some computation we find that these equations take the following form:

[TABLE]

Except in the degenerate cases when $\varphi_{1}=\varphi_{2}=\pm A^{-1/2}$ and (59) holds, the eigenvalue problem in (24)–(26) is therefore equivalent to the vanishing of a determinant:

[TABLE]

where $\varphi_{1}=\varphi$ and $\varphi_{2}=1/(A\varphi)$, and the functions $f$, $g$ are given by

[TABLE]

The function $\delta$ depends on $M$ and $K$, but this dependence will not be written explicitly for simplicity. We note the general root-exchange symmetry

[TABLE]

Because $\lambda=0$ is an eigenvalue we also know that $\delta$ has roots at $\varphi=1$ and $1/A$. Note that $\delta(\pm A^{-1/2})=0$ due to dependence of the columns, but these roots are spurious, unless double, as we now discuss.

The degenerate case. In the cases of (59) when the two roots of (57) coincide at $\varphi=\varphi_{1}=\pm A^{-1/2}=\varphi_{2}$, one checks that the difference equation (25) has the general solution

[TABLE]

by the expedient of replacing $c_{1}$, $c_{2}$ in (58) with

[TABLE]

and taking $\varphi_{1}\to\pm A^{-1/2}$. Doing the same with (63)–(64), we see that the eigenvalue condition (65) is replaced by the condition

[TABLE]

This is equivalent to the condition $\delta^{\prime}(\varphi)=0$ because one finds $\delta^{\prime}(\varphi)=-2\hat{\delta}(\varphi)$ at these points.

Remark 5.2*.*

In order to characterize Hopf bifurcation, we will use the fact that when $\varphi\neq\pm A^{-1/2}$, 1 or $A^{-1}$, $\varphi$ is a simple root of $\delta(\varphi)$ if and only if $\lambda=(A-\varphi^{-1})(\varphi-1)$ is a simple eigenvalue of $B$. See Lemma 5.4 and its proof in Section 10.

5.3 Sorting terms and removing singularities

For convenience in analysis, we sort the terms in (65) according to $M$th powers of $\varphi$ and $A$. Note that

[TABLE]

In order to remove singularities, we multiply (65) by $\varphi(\varphi-1)(1-A\varphi)$. Define

[TABLE]

where

[TABLE]

with the definitions

[TABLE]

By consequence we have the sorted representation

[TABLE]

where

[TABLE]

Observe that $F$ has a pole at $\varphi=0$ of order $M$, with $F(\varphi)\sim-A^{1-M}\varphi^{-M}$, because $R_{2}=H_{2}=-A$ at the origin. And for $|\varphi|\to\infty$ we find that

[TABLE]

Consequently $F$ must have exactly $2M+4$ zeros, counting multiplicities.

We may summarize the situation as follows.

Proposition 5.3**.**

A complex number $\lambda$ is an eigenvalue of $B$ if and only if (61) holds for some pair $\varphi$, $1/A\varphi$ satisfying

[TABLE]

except in the two cases $\lambda=-1-A\pm 2\sqrt{A}$. In these cases, $\lambda$ is an eigenvalue if and only if

[TABLE]

Of the $2M+4$ roots of $F$, four are spurious, counting $\varphi=\pm A^{-1/2}$, $\varphi=1$ and $\varphi=1/A$ once each, coming from the dependence of the columns in (65) and the factors used to remove singularities from $\delta$.

The polynomial $\varphi^{M}F(\varphi)$ of degree $2M+4$ is divisible by the factor

[TABLE]

and the remaining $2M$ roots of $\varphi^{M}F(\varphi)/S(\varphi)$ correspond in pairs $\varphi$, $1/(A\varphi)$ to the $M$ eigenvalues of $B$. The values $\varphi=1$ and $1/A$, are (at least) double roots of $F$ because they were already roots of $\delta$, and correspond to the simple eigenvalue $\lambda=0$. Concerning other roots of $F$, we have the following result whose proof we defer to Section 10.

Lemma 5.4**.**

Suppose $S(\varphi)\neq 0$ and $\lambda=(A-\varphi^{-1})(\varphi-1)$. Then $\varphi$ is a simple root of $F$ if and only if $\lambda$ is a simple eigenvalue of $B$.

6 Formal approximation

Before we begin a rigorous analysis of the zeros of $F(\varphi)$, we treat the problem approximately in the limit of large $M$ to gain insight. Numerical experimentation suggests that we can expect to find most solutions of (65) to satisfy $|\varphi_{1}|\approx 1$, and $|\varphi_{2}|\approx A^{-1}<1$, with $A^{-M}$ extremely small.

Thus we neglect the terms containing $A^{-M}$ in (76) and study the zeros of

[TABLE]

For any such zero, evidently

[TABLE]

unless both numerator and denominator vanish. The right-hand side is a ratio of polynomials of low degree, while for large $M$, the function $\varphi\mapsto\varphi^{M}$ expands a small region about any $M$th root of unity $e^{2\pi ik/M}$ to cover a large part of the complex plane. Roughly, then, we can expect (81) to have a solution near each $M$th root of unity. These should then provide eigenvalues spread out around the ellipse in (62).

We focus next on looking for imaginary roots $\varphi\approx 1$. We change variables from $\varphi$ to $z=M(\varphi-1)$, noting that

[TABLE]

With these relations we have

[TABLE]

and we find from (72)–(74) the exact expressions

[TABLE]

It turns out to be appropriate to require $K$ is small while $KM$ is large. Somewhat more precisely, we ask that

[TABLE]

Then we get the approximate relations

[TABLE]

By consequence, we find that

[TABLE]

If we suppose $K\sim\kappa/\sqrt{M}$ as $M\to\infty$, then

[TABLE]

where

[TABLE]

The complex roots of $Q$ provide an approximation for roots of $F_{0}(\varphi)$ when $M$ is large. These approximate eigenvalues $\lambda$ of (24)–(26) through (55), which may be written directly in terms of $z$ as

[TABLE]

Thus purely imaginary roots of $Q$ approximate eigenvalues $\lambda$ near the imaginary axis, and roots of $Q$ in the right half plane $\mathop{\rm Re}\nolimits z>0$ should approximate eigenvalues satisfying $\mathop{\rm Re}\nolimits\lambda>0$.

7 Analysis of roots of $Q$

In this section we establish basic properties of the roots $z$ of $Q(z;\kappa)$ as defined in (87). This will serve as the foundation to analyze the roots of $F_{0}$ and ultimately those of $F$, in subsequent sections.

Purely imaginary roots $z=it$ of $Q$ occur whenever

[TABLE]

Matching real parts demands that if $t\neq 0$ then $\sec t=1-\frac{t^{2}}{\kappa^{2}}<1$ hence $\sec t\leq-1$. Matching also the ratio of imaginary to real parts, one finds that (89) holds if and only if

[TABLE]

together with

[TABLE]

Each positive root of (90) provides a complex conjugate pair of imaginary roots $z=\pm it$ of $Q$. Let $t_{1}<t_{2}<\ldots$ denote the increasing sequence of all these positive roots of (90). The smallest occurs for $t=t_{1}\approx 4.4934095$ (less than $\frac{3}{2}\pi\approx 4.71238898$).  This corresponds to a critical value of $\kappa$ given by

[TABLE]

The roots $t_{k}$ approach $\frac{3}{2}\pi+2\pi k$ from below as $k\to\infty$. As $k$ increases, they correspond to larger values of $\kappa^{2}$, hence larger values of $K$ for a fixed $M$.

In the rest of this section, we shall prove that non-real roots of $Q$ are always simple, and purely imaginary roots must move into the right half plane $\mathop{\rm Re}\nolimits z>0$ as $\kappa$ increases, where they must remain in a bounded region. By this result and (88), when $\kappa>\kappa_{\rm cr}$ we can expect that for large enough $M$ with $K\sim\kappa/\sqrt{M}$, there will be some eigenvalue $\lambda$ of (24)–(26) in the right half plane, and when $0<\kappa<\kappa_{\rm cr}$ we can expect there will not.

Lemma 7.1**.**

For any $\kappa>0$, $Q$ has a double root $z=0$. All other complex roots are non-real and simple.

Proof.

Clearly $Q(0;\kappa)=0$, and for real $z\neq 0$ we have $Q>e^{z}-1-z>0$ by the convexity of $e^{z}$. In general we compute

[TABLE]

The root $z=0$ is double because $0=Q=\partial_{z}Q<\partial^{2}_{z}Q$ at 0. At a complex double root, on the other hand, necessarily $e^{z}=-\kappa^{2}/2$. This implies $z=r+i\pi k$ where $e^{r}=\kappa^{2}/2$ and $k$ is an odd integer. Then, however, it follows

[TABLE]

so $r=-\frac{1}{2}$ and we infer

[TABLE]

a contradiction. Hence the nonzero roots of $Q$ are all non-real and simple. ∎

For the next result, let $t_{0}=0$ and recall that $t_{1}<t_{2}<\ldots$ denotes the sequence of positive roots of (90).

Lemma 7.2**.**

The function $Q$ has exactly $k$ complex-conjugate pairs of roots $z$ in the right half plane $\mathop{\rm Re}\nolimits z>0$ if $\kappa^{2}=\sqrt{1+t^{2}}-1$ with $t\in(t_{k},t_{k+1}]$.

Proof.

First, we claim that the imaginary roots of $Q$ always cross into the right half plane $\mathop{\rm Re}\nolimits z>0$ as $\kappa$ increases. To see this, regard $w:=\kappa^{2}$ as a complex variable and note that $Q=0$ if and only if

[TABLE]

Because $\frac{d}{dz}(1+z)e^{-z}=-ze^{-z}$, we compute

[TABLE]

by using the identity $Q=0$ to eliminate $e^{-z}$. Multiplying by $\bar{z}|1+z|^{2}$, we find

[TABLE]

For $z=x+iy$ in the first quadrant, the imaginary part of this expression is negative, which implies

[TABLE]

Furthermore, provided $y^{2}>\kappa^{2}$ (which must be the case if $x=0$ by (91)), the real part of (95) is larger than $y^{4}-wy^{2}>0$, hence

[TABLE]

It follows from these computations that the roots $z=\pm it_{k}$ of $Q$ on the imaginary axis always pass into the right half plane as $\kappa$ increases, with derivative $dz/d\kappa$ remaining in the first quadrant. They can never escape to infinity, because any roots of $Q$ in the right half plane must lie in the bounded region where

[TABLE]

To finish the proof, we show that if $\kappa>0$ is small enough, then $Q$ has no roots with $\mathop{\rm Re}\nolimits z>0$. If $\kappa\in(0,1)$, any such root must satisfy

[TABLE]

and this implies $|z|<2\kappa$. Now it follows

[TABLE]

Therefore, for small enough $\kappa>0$, $Q$ does not vanish when $\mathop{\rm Re}\nolimits z>0$. ∎

Labeling the roots. Due to the results of the previous two lemmas, we may label all the non-real roots of $Q$ that cross the imaginary axis and lie in the upper half plane $\mathop{\rm Im}\nolimits z>0$ by analytic functions $z=z_{j}^{0}(\kappa)$, $j=1,2,\ldots$, defined for all $\kappa>0$ according to the property that

[TABLE]

Thus we can summarize as follows.

Lemma 7.3**.**

There are analytic curves $z_{j}^{0}\colon(0,\infty)\to\mathbb{C}$, $j=1,2,\ldots$, satisfying (98) and $\mathop{\rm Im}\nolimits z_{j}^{0}(\kappa)>0$ for all $\kappa>0$, such that when $\kappa\in(\kappa^{0}_{k},\kappa^{0}_{k+1}]$, the numbers $z_{1}^{0}(\kappa)\ldots,z_{k}^{0}(\kappa)$ comprise all the roots of $Q$ in the first quadrant. Moreover, $dz_{j}^{0}/d\kappa\neq 0$ for all $\kappa>0$, and

[TABLE]

Proof.

To show the curves $z_{j}^{0}$ are well defined and nondegenerate for all $\kappa>0$, we note that according to standard continuation theory for the ODE (94), a solution exists for real $w$ in a maximal interval $(w_{-},\infty)\subset(0,\infty)$ for which $dz/dw$ remains bounded. It is not possible that $w_{-}>0$, however, because the right-hand side of (94) cannot approach zero at the same time as (93) holds with $w\to w_{-}$, for the following reason: If (94) vanishes, then $0=2(1+z)+w+z^{2}$, hence $z=-1+i\tau$ with $\tau=\sqrt{1+w}>1$. But then (93) implies

[TABLE]

This implies $\tau^{2}=1-2e^{i\tau-1}$, so necessarily $\sin\tau=0$ but also $1<\tau^{2}<1+2/e$, and this is impossible. ∎

8 Analysis of roots of $F_{0}$

In this section we locate all the roots of the polynomial $F_{0}=\varphi^{M}P_{2}-P_{1}$ in (80) of degree $M+4$, to a rough approximation, provide bounds on roots that may correspond to unstable eigenvalues, and establish the convergence in (86) in a precise sense. Let $B(z,r)\subset\mathbb{C}$ denote the closed disk with center $z\in\mathbb{C}$ and radius $r>0$. We fix a constant $\gamma>2$. (Actually, $\gamma=3$ suffices.) Depending on some large $\beta>1$ (to be chosen in the proof of Theorem 3.2), we presume throughout that

[TABLE]

8.1 Rough locations of all roots

Locations of the $M+4$ roots of $F_{0}$ will be identified as follows. We recall that the four values $\varphi=1$, $A^{-1}$, $\pm A^{-1/2}$, which comprise the roots of the polynomial

[TABLE]

from (79), are already known to be roots of the function $F$ that $F_{0}$ approximates. Note that the three roots of $S(\varphi)$ with $\varphi\neq 1$ satisfy

[TABLE]

Proposition 8.1**.**

Fix $\gamma>2$. Then for any $\beta>1$ there exists $\alpha_{0}>0$ and $M_{0}>0$ such that whenever $M>M_{0}$ and (100) holds, the polynomial $F_{0}$ has exactly:

• (i)

one double root at $\varphi=1$.

• (ii)

one simple root in each of the following disks of radius $r_{K}=K/8$:

[TABLE]

• (iii)

one simple root in $B(-M,1)$.

• (iv)

$M-2$* roots in the punctured annulus*

[TABLE]

Proof.

Recall $F_{0}(\varphi)=\varphi^{M}P_{2}-P_{1}$. where we can write

[TABLE]

with $S(\varphi)$ as in (79), and

[TABLE]

Step 1. First we establish (i). Note that $F_{0}(1)=0$, since

[TABLE]

Furthermore, $F_{0}^{\prime}(1)=MK^{3}+P_{2}^{\prime}(1)-P_{1}^{\prime}(1)=0$ since

[TABLE]

Hence $\varphi=1$ is at least a double root. But one also checks

[TABLE]

(e.g., by computer algebra) so $F_{0}^{\prime\prime}(1)\neq 0$ when $KM\geq 1$.

Step 2. Next we claim that the only roots of $F_{0}$ in the disk $B(0,M^{-\gamma/M})$ are three as described in (ii). We can write

[TABLE]

where

[TABLE]

It suffices to show that for all $\varphi$ in $B(0,M^{-\gamma/M})$ outside the balls listed in (ii),

[TABLE]

for $M$ large enough. For then our claim follows from Rouché’s theorem, since each of the balls in (ii) contains one simple root of $S$.

Observe that $|P_{2}|\leq CM$ for $|\varphi|\leq 1$, therefore

[TABLE]

(Here and below $C$ denotes a generic constant which may depend on $\beta$ and $\gamma$ but is independent of $M$ and $K$, whose value may change from instance to instance.) To complete the proof of (107), we consider three sub-cases:

[TABLE]

In case (a), for each $\hat{\varphi}\in\{1,A^{-1},A^{-1/2}\}$ (i.e., for each positive root of $S$), we have $1-K<|\varphi-\hat{\varphi}|<2$, therefore

[TABLE]

Because $\gamma>2$ and $K^{2}M\geq\beta^{-2}$ it follows that for $|\varphi+A^{-1/2}|\geq r_{K}=K/8$, with $M$ large enough we have

[TABLE]

(We could replace $r_{K}$ by say $20/M$ here, but we have no need.)

In case (b), each positive root of $S$ satisfies $|1-\hat{\varphi}|\leq K$, hence

[TABLE]

Consequently

[TABLE]

and (for $K<\frac{1}{2}$)

[TABLE]

Therefore as in (110) we get

[TABLE]

for $M$ large enough depending on $\beta$ and $\gamma$.

In case (c), we have

[TABLE]

with min and max taken over positive roots of $S$. Therefore for $M$ large, when $|\varphi-\hat{\varphi}|\geq r_{K}=K/8$ (chosen to separate the roots) we find

[TABLE]

for some $c>0$ depending on $\beta$.

This finishes the proof of (107). The conclusion in (ii) now follows, and also the fact that $F_{0}$ has no other roots in $B(0,M^{-\gamma/M})$.

Step 3. Next we show that $F_{0}$ has no roots satisfying

[TABLE]

for large enough ${\alpha_{0}}$ depending on $\beta$, and deduce (iii) and (iv). The estimates in (113) imply

[TABLE]

Observe

[TABLE]

where

[TABLE]

(To get this last, expand $(A\varphi-1)^{2}=(K+A(\varphi-1))^{2}$ and cancel a term.)

We now show the ratios $S_{3}/S_{4}$ and $\varphi^{-M}P_{1}/S_{4}$ are uniformly small for $\varphi$ satisfying (113), by estimating six terms as follows:

(a) The first term of the ratio $S_{3}/S_{4}$ is bounded using (114) as follows:

[TABLE]

(b) To bound the next term in $S_{3}/S_{4}$, observe

[TABLE]

For $|\varphi|>M/2$, since $|M+\varphi|\geq 1$, for $M>4$ we have

[TABLE]

while for $|\varphi|\leq M/2$ we have $|M+\varphi|>M/2$ and infer from (114) that

[TABLE]

(c) The last term in the ratio $S_{3}/S_{4}$ satisfies the bound

[TABLE]

(d) The terms in $\varphi^{-M}P_{1}/S_{4}$ are estimated as follows. By (114),

[TABLE]

Further, $A\varphi^{2}-1=A(\varphi-A^{-1/2})(\varphi+A^{-1/2})$ and

[TABLE]

Therefore, since $|\varphi+A^{-1/2}|<2|\varphi|$ and recalling $|M+\varphi|^{-1}\leq 2|\varphi|/M$,

[TABLE]

Hence, since $KM+1<2KM$, the last term in $\varphi^{-M}P_{1}/S_{4}$ is bounded by

[TABLE]

(e) For the next term in $\varphi^{-M}P_{1}/S_{4}$, we have the bound

[TABLE]

(f) Lastly we have the bound

[TABLE]

Assembling the estimates in (a)-(f), we conclude that if $\alpha_{0}\geq\alpha_{0}(\beta)$ and $M\geq M_{0}(\beta)$, then

[TABLE]

for all $\varphi$ satisfying (113). Part (iii) now follows by Rouché’s theorem since $S_{4}$ has only one simple zero at $\varphi=-M$ inside $B(-M,1)$. Part (iv) follows since we have shown that $F_{0}$ has exactly 6 roots (counting multiplicity) in the complement of the punctured annulus $D_{a}$. ∎

We record here several estimates that follow from the proof above.

Corollary 8.2**.**

Under the conditions of Proposition 8.1, we have the following estimates, for some $c>0$ depending on $\beta$:

• (i)

$\displaystyle|F_{0}(\varphi)|\geq\frac{\beta^{-2}}{20}$*  if   $|\varphi+A^{-1/2}|=r_{K}$.*

• (ii)

$\displaystyle|F_{0}(\varphi)|\geq\frac{c}{M}$*    if   ${|\varphi-A^{-1}|=r_{K}}$ or   ${|\varphi-A^{-1/2}|=r_{K}}$.*

• (iii)

$\displaystyle|F_{0}(\varphi)|\geq\frac{c}{M}$*   if   ${|\varphi|=M^{-\gamma/M}}$.*

• (iv)

$\displaystyle|\varphi^{-M}F_{0}(\varphi)|\geq\frac{1}{2}M^{3}$*  if   ${|\varphi+M|=1}$.*

• (v)

$\displaystyle|F_{0}(\varphi)|\geq\frac{K\alpha_{0}^{2}}{2M^{2}}$*   if   $\displaystyle{|\varphi|^{-1}=1-\frac{\alpha_{0}}{M}}$.*

Proof.

Part (i) follows from (110) in case (a) of Step 2, because $|F_{0}(\varphi)|\geq\Delta_{0}$. Similarly, part (ii) follows from (112) in case (c) of Step 2, and part (iii) follows from all cases of Step 2. To infer part (iv), note that (127) of Step 3 implies that for $|\varphi+M|=1$ we have

[TABLE]

because $|A\varphi-1|\geq A|\varphi-1|-K\geq AM-K\geq M$. Part (v) follows similarly, since $|\varphi|\geq 1+\alpha_{0}/M$ and therefore $|S_{4}|\geq K(|\varphi|-1)^{2}\geq K\alpha_{0}^{2}/M^{2}$. ∎

8.2 Bounds for roots relevant to instability

Next we focus on roots of $F_{0}$ that may be related to eigenvalues $\lambda$ of the matrix $B$ having non-negative real part. It turns out these are roots $\varphi$ in the punctured annulus $D_{a}$ of Proposition 8.1 that are near $1$. Recall the relation (61) between eigenvalues of the matrix $B$ and roots $\varphi$ of $F$, namely

[TABLE]

Lemma 8.3**.**

Under the conditions of Proposition 8.1, if $M$ is large enough, then whenever (61) holds with $\varphi\in D_{a}$, then $\mathop{\rm Re}\nolimits\lambda\geq 0$ implies

[TABLE]

Proof.

By (61), $\lambda=A\varphi-A-1+\varphi^{-1}$, hence if $\mathop{\rm Re}\nolimits\varphi<0$ then $\mathop{\rm Re}\nolimits\lambda<-A-1$. Writing

[TABLE]

we then have $\mu\geq-1$ and

[TABLE]

For $\varphi\in D_{a}$ and $M$ large, we infer $|\varphi|^{-2}\leq M^{2\gamma/M}\leq 1+{4\gamma}M^{-1}\log M$, then

[TABLE]

Now because $|\varphi|^{-2}\leq(1+\mu)^{-2}$, we deduce from (130) that

[TABLE]

This entails $\mu\geq 0$, due to (131). Since $|\varphi|^{-1}\geq 1-\frac{\alpha_{0}}{M}$ implies $\mathop{\rm Re}\nolimits\varphi<1+\frac{2\alpha_{0}}{M}$, we have established the desired bounds on $\mathop{\rm Re}\nolimits\varphi$.

Now since $|\varphi|^{2}=(1+\mu)^{2}+\nu^{2}$ and $0\leq\mu<\frac{2\alpha_{0}}{M}$, we deduce from (130) that

[TABLE]

Since we may presume $\beta\leq\alpha_{0}$, therefore $|\nu|<2\alpha_{0}M^{-3/4}$ as claimed. ∎

Any roots of $F_{0}$ in the region where (129) holds actually satisfy a tighter bound, namely $|\varphi-1|=O(1/M)$, as we now show.

Proposition 8.4**.**

Under the conditions of Proposition 8.1, there exist positive constants $\alpha_{1}=\alpha_{1}(\beta)$ and $M_{1}=M_{1}(\beta)$ such that whenever $M>M_{1}$, any zeros $\varphi\in D_{a}$ of $F_{0}$ that satisfy the bounds in (129) must satisfy $|\varphi-1|\leq\frac{\alpha_{1}}{M}$. Moreover,

[TABLE]

for all $\varphi$ that satisfy

[TABLE]

Proof.

In the expression $F_{0}=\varphi^{M}P_{2}-P_{1}$ we seek to show that the first term dominates, provided (133) holds for some $\alpha_{1}$. Writing $\zeta=A(\varphi-1)$ for convenience, we have $A\varphi-1=K+\zeta$, so by (8.1),

[TABLE]

By (133) we have $\zeta/K=O(M^{-\frac{1}{4}})$ and $|\varphi|\leq 1+O(M^{-3/4})$, so

[TABLE]

Because

[TABLE]

for $\alpha_{1}>2\beta$ and large enough $M$ we infer that

[TABLE]

On the other hand, due to (122) we have

[TABLE]

therefore from (8.1) we obtain the upper bound

[TABLE]

if $\alpha_{1}>40\beta^{2}$, say, and $M$ is large enough. Since $|\varphi|\geq 1$ if (133) holds, the result follows. ∎

8.3 Convergence of $K^{-3}F_{0}(1+z/M)$

After the results of the previous subsection, to study unstable eigenvalues of $B$ we are motivated to make the change of variables

[TABLE]

as in Section 6. According to Proposition 8.4, for any zeros $\varphi\in D_{a}$ of $F_{0}$ that correspond to $\mathop{\rm Re}\nolimits\lambda\geq 0$, the quantity $z=M(\varphi-1)$ must satisfy

[TABLE]

As in (31), let us now define $\kappa=\kappa(K,M)=K\sqrt{M}$ and $\varepsilon=1/\sqrt{M}$. Then the formal approximations in Section 6 are rigorous, with errors that are uniform over the values of $(z,\kappa)\in\mathbb{C}\times\mathbb{C}$ such that

[TABLE]

where $\hat{\alpha}>\alpha_{1}$ is an arbitrary constant (to be chosen later), and

[TABLE]

for some small $\hat{\gamma}>0$. (We allow $\kappa$ to be complex with small argument here, to simplify derivative estimates later.) By consequence, the convergence in (86) holds, in the following sense.

Proposition 8.5**.**

Uniformly for $(z,\kappa)$ satisfying (137)–(138), with $M=1/\varepsilon^{2}$ and $K=\kappa{\varepsilon}$ we have that

[TABLE]

9 Analysis of roots of $F$

Recall from (76) we have

[TABLE]

where $R_{1}$, $R_{2}$ are low-degree polynomials that may be written in the form

[TABLE]

For large $M$, $A^{-M}$ is exponentially small, with the bound

[TABLE]

We now roughly characterize the location of the $2M+4$ roots of $F$.

Proposition 9.1**.**

Under the conditions of Proposition 8.1, there exists $M_{2}=M_{2}(\beta)$ such that whenever $M>M_{2}$, $F$ has (counting multiplicities):

• (i)

one double root at $\varphi=1$, and one double root at $\varphi=A^{-1}$.

• (ii)

$M-2$* roots in the punctured annulus $D_{a}$, and $M-2$ roots with $(A\varphi)^{-1}\in D_{a}$ which satisfy $|\varphi|<1-\frac{3}{4}K<M^{-\gamma/M}$.*

• (iii)

one simple real root in $B(-M,1)$, and one with $(A\varphi)^{-1}\in B(-M,1)$.

• (iv)

one simple real root at $\varphi=A^{-1/2}$ and one at $\varphi=-A^{-1/2}$.

Proof.

We note that due the root symmetry (66), the multiplicity of each root $\varphi$ of $F$ is the same as the multiplicity of $1/(A\varphi)$, unless $\varphi=\pm A^{-1/2}$. Also, all non-real roots of $F$ come in complex-conjugate pairs when $K$ is real.

For $|\varphi|\geq M^{-\gamma/M}$ we then have $|\varphi|^{-M}<M^{\gamma}$ and it follows $A^{-M}|F_{1}(\varphi)|$ is exponentially small.

Combining the lower bounds in parts (iii)–(v) of Corollary 8.2 with the count of roots of $F_{0}$ in parts (i), (iii) and (iv) of Proposition 8.1, we conclude from Rouché’s theorem that $F$ has a simple root inside the ball $B(-M,1)$, and $M$ roots inside the closed annulus $D_{a}\cup\{1\}$, the same as $F_{0}$.

By examining (140)–(141), we find $\varphi^{M}F_{1}=\varphi^{M}R_{1}+R_{2}$ has at least a double root at $\varphi=1$, due to the fact that the expression

[TABLE]

has a double root at $\varphi=1$. Then, because $A^{-M}F_{1}^{\prime\prime}(1)$ is exponentially small, it follows from (104) that $F^{\prime\prime}(1)\neq 0$. This proves (i).

Now (iii) follows and also (ii), due to the fact that for $\varphi\in D_{a}$ and $M$ large,

[TABLE]

To infer (iv) we can simply recall that we know $F(\pm A^{-1/2})=0$ due to the root symmetry relation (66). These roots must be simple, since we have accounted for all $2M+4$ roots of $F$. ∎

Next, we can characterize zeros of $F$ that may correspond to unstable eigenvalues of $B$ as follows.

Proposition 9.2**.**

Under the conditions of Propositions 9.1 and 8.4, there exists $M_{3}=M_{3}(\beta)$ such that whenever $M>M_{3}$ and $\lambda$ is an eigenvalue of $B$ with $\mathop{\rm Re}\nolimits\lambda\geq 0$, then $\lambda=(A-\varphi^{-1})(\varphi-1)$ for some root $\varphi$ of $F$ that satisfies

[TABLE]

Proof.

Under the correspondence between $\lambda$ and $\varphi$ in (61), the zeros of $F$ described in parts (iii) and (iv) of Proposition 9.1 correspond to negative real values of $\lambda$, and the roots in part (i) correspond to $\lambda=0$. So, given $M$ is large enough, for any nonzero eigenvalue $\lambda$ satisfying $\mathop{\rm Re}\nolimits\lambda\geq 0$, necessarily (61) holds for some $\varphi\in D_{a}$. This $\varphi$ must satisfy the bounds in (129), due to Lemma 8.3. For these values of $\varphi$, we have $|F_{1}(\varphi)|\leq CM^{1/4}$, so $A^{-M}|F_{1}|$ is exponentially small. Then we can conclude from Proposition 8.4 that

[TABLE]

for all $\varphi$ that satisfy (133). The conclusion follows. ∎

Further, the convergence in Proposition 8.5 holds with $F$ in place of $F_{0}$:

Proposition 9.3**.**

Let $\hat{\alpha}>\alpha_{1}$, and let $\hat{\gamma}>0$ be small. Uniformly for $(z,\kappa)$ satisfying (137)–(138), with $M=1/\varepsilon^{2}$ and $K=\kappa{\varepsilon}$ we have that

[TABLE]

Furthermore, for each pair of integers $j,k\geq 0$, the derivatives

[TABLE]

uniformly for all $z$ and $\kappa$ satisfying

[TABLE]

Proof.

For $|z|\leq\hat{\alpha}$ and $\varphi=1+\frac{z}{M}$, the factor $|\varphi|^{-M}$ is bounded by $e^{2\hat{\alpha}}$. Hence again $A^{-M}F_{1}$ is exponentially small, and the convergence of $Q^{\varepsilon}(z;\kappa)$ follows from Proposition 8.5.

The convergence of derivatives follows from the Cauchy integral formula representation for such derivatives, since $Q^{\varepsilon}(z;\kappa)$ is analytic for $z$ satisfying (137) and $\kappa$ satisfying (138). ∎

Curves of roots. Recall that the non-real roots $z$ of $Q=Q(z;\kappa)$ are simple and those that may satisfy $\mathop{\rm Re}\nolimits z\geq 0$ lay on the curves $z_{j}^{0}(\kappa)$ described by Lemma 7.3. Moreover, due to (98) and (99), only a finite number of these curves provide values that can satisfy (136), corresponding to values of $\varphi=1+\frac{z}{M}$ that satisfy (142). In particular, we note the following.

Corollary 9.4**.**

For $j\in\mathbb{N}$, if $\beta\geq\kappa^{0}_{j}$ and $\alpha_{1}$ is given by Proposition 8.4, then

[TABLE]

Proof.

Suppose $\alpha_{1}\leq|z_{j}^{0}(\kappa)|$ for some $\kappa\in[\kappa^{0}_{j},\beta]$. Recall $z=z_{j}^{0}(\kappa)$ satisfies $Q(z,\kappa)=0$, $\mathop{\rm Re}\nolimits z\geq 0$. Then for $M$ large enough, $\varphi=1+\frac{z}{M}$ satisfies (133), and

[TABLE]

due to (143). But this contradicts the convergence result in Proposition 9.3. ∎

Any finite number of the curves $z_{j}^{0}$ of simple zeros of $Q$ perturb to curves $z_{j}^{\varepsilon}$ of simple zeros of $Q^{\varepsilon}$ as a consequence of the implicit function theorem, as follows.

Proposition 9.5**.**

For $j\in\mathbb{N}$, suppose $\beta>\kappa^{0}_{j}$. Let $\alpha_{1}$ be given by Proposition 8.4, and suppose

[TABLE]

Then for sufficiently small $\varepsilon>0$, there is a curve $z_{j}^{\varepsilon}:[\beta^{-1},\beta]\to B(0,\hat{\alpha})$ that is real analytic, with the following properties:

• (i)

For each $\kappa\in[\beta^{-1},\beta]$, $z_{j}^{\varepsilon}(\kappa)$ is a simple root of $Q^{\varepsilon}(z,\kappa)$.

• (ii)

$z_{j}^{\varepsilon}(\kappa)\to z_{j}^{0}(\kappa)$* as $\varepsilon\to 0$, uniformly for $\kappa\in[\beta^{-1},\beta]$, together with any finite number of derivatives in $\kappa$.*

• (iii)

There exists $\zeta_{j}<\beta$ satisfying $\zeta_{j}^{\varepsilon}\to\kappa_{j}^{0}$ as $\varepsilon\to 0$, such that $\mathop{\rm Re}\nolimits z_{j}^{\varepsilon}(\kappa)\geq 0$ if and only if $\kappa\geq\zeta_{j}^{\varepsilon}$, and

[TABLE]

Proof.

The existence of the curve, its analyticity in $\kappa$, and properties (i), (ii) and (iii), follow from standard implicit function theorem arguments using the simplicity of the roots of $Q$, the convergence in Proposition 9.3, and Lemma 7.3. ∎

10 Analysis of eigenvalues of $B$

The $M$ eigenvalues $\lambda$ of $B$ are generated via the relation (61) by: one of the roots of $F$ at $\varphi=1$, the one near $-M$, and the $M-2$ roots in $D_{a}$. The roots $\pm A^{-1/2}$, one root at $1$, and one root at $A^{-1}$ are spurious, as discussed earlier. We have not characterized the multiplicity of all the eigenvalues or all the roots, but each eigenvalue must correspond to some root of $F$, and vice versa.

10.1 Curves of unstable eigenvalues

Recall that zeros $z$ of $Q^{\varepsilon}(z;\kappa)$ correspond to eigenvalues $\lambda$ of the matrix $B$ via the relation (88). We rescale this relation by defining

[TABLE]

Clearly $\Lambda(z;\kappa,\varepsilon)\to z$ as $\varepsilon\to 0$, together with derivatives, uniformly for $z$, $\kappa$ satisfying (146).

When $\varepsilon=0$, of course we have $\mathop{\rm Re}\nolimits\Lambda(z;\kappa,0)\geq 0$ if and only if $x=\mathop{\rm Re}\nolimits z\geq 0$, for any $\kappa>0$. By stardard implicit function theorem arguments, for small enough $\varepsilon>0$ there is a real analytic function $(y,\kappa)\mapsto\hat{x}(y,\kappa,\varepsilon)$ such that for $|z|\leq\hat{\alpha}$ and $\kappa\in[\beta^{-1},\beta]$,

[TABLE]

Let ${\mathcal{I}}_{\varepsilon}\subset B(0,\hat{\alpha})\times[\beta^{-1},\beta]$ denote the surface on which this holds, i.e., where $\mathop{\rm Re}\nolimits\Lambda=0$. When $\varepsilon=0$, the imaginary axis ${\mathcal{I}}_{0}$ meets each curve $z_{j}^{0}$ transversely due to the computation in (97). Therefore, for sufficiently small $\varepsilon>0$, the surface ${\mathcal{I}}_{\varepsilon}$ meets each curve $z_{j}^{\varepsilon}$ provided by Proposition 9.5 transversely. By consequence, each curve given by

[TABLE]

provides a curve of eigenvalues of $B$ that must cross the imaginary axis transversely as $\kappa$ increases, exactly once for $\kappa\in[\beta^{-1},\beta]$.

10.2 Proof of Theorem 3.2

Let $\beta_{0}\in(0,1)$ and $k\in\mathbb{N}$. Recalling that the curves $z_{j}^{0}(\kappa)$ and numbers $\kappa^{0}_{j}>1$ were defined in (98), we fix $\beta_{k}\in(\kappa^{0}_{k},\kappa^{0}_{k+1})$, and note

[TABLE]

Next, choose $\beta>\max(\beta_{k},\beta_{0}^{-1})$, let $\alpha_{1}=\alpha_{1}(\beta)$ be determined by Proposition 8.4, and choose $\hat{\alpha}>\alpha_{1}$ such that

[TABLE]

If $M$ is sufficiently large (i.e., $M>M_{0,k}$ for some $M_{0,k}$ depending on $\beta$) then analytic curves $z_{j}^{\varepsilon}(\kappa)$ are defined by Proposition 9.5 and $\lambda_{j}^{\varepsilon}(\kappa)$ by (150). Let

[TABLE]

Due to Propositions 9.2 and 9.5 and the discussion above, each curve $\lambda_{j}$ crosses the imaginary axis transversely at some point $\kappa_{j}=\kappa_{j}^{\varepsilon}\in[\zeta_{j}^{\varepsilon},\beta]$ that satisfies

[TABLE]

By consequence, for small enough $\varepsilon>0$ we have $\kappa_{j-1}^{\varepsilon}<\kappa_{j}^{\varepsilon}<\beta_{k}$ for $j=1,\ldots,k$, where we set $\kappa_{0}^{\varepsilon}=\beta_{0}$. Also we have the monotonicity relations in (33).

Since $|z_{j}^{\varepsilon}(\kappa)|\leq\hat{\alpha}$, the eigenvalues of $B$ given by $\lambda_{j}(\kappa)$, $j=1,\ldots,k$ satisfy the bound

[TABLE]

for $M$ large. Furthermore, due to Lemma 5.4 (proved below), every such eigenvalue $\lambda_{j}(\kappa)$ is a simple eigenvalue of $B$, since the roots $z=z_{j}^{\varepsilon}(\kappa)$ of $Q^{\varepsilon}(z;\kappa)$ are simple.

It remains to prove that for $\kappa\in[\beta_{0},\beta_{k}]$, if $\hat{\lambda}\neq 0$ is an eigenvalue of $B$ with $\mathop{\rm Re}\nolimits\hat{\lambda}\geq 0$, and $\mathop{\rm Im}\nolimits\lambda\geq 0$, then necessarily $\hat{\lambda}=\lambda_{j}(\kappa)$ for some $j\leq k$ with $\kappa\geq\kappa_{j}^{\varepsilon}$. According to Proposition 9.2, necessarily such an eigenvalue must satisfy

[TABLE]

where $Q^{\varepsilon}(\hat{z};\kappa)=0$, $\mathop{\rm Re}\nolimits\hat{z}\geq 0$ and $|\hat{z}|\leq\alpha_{1}$.

Now, for any $r>0$ sufficiently small, note that the balls $B(z_{j}^{0}(\kappa),r)$ do not overlap or contain [math] for any $\kappa$, and each must contain a simple root $z_{j}^{\varepsilon}(\kappa)$ of $Q^{\varepsilon}(z,\kappa)$. Fix some such $r>0$, and let $\Omega_{r}$ be the set of $(z,\kappa)$ such that

[TABLE]

and $\kappa\in[\beta_{0},\beta_{k}]$. Because $\beta_{k}<\kappa^{0}_{k+1}$, for sufficiently small $r>0$ we have

[TABLE]

From the convergence in Proposition 9.3 it follows

[TABLE]

if $\varepsilon>0$ is sufficiently small. Then it follows that $|\hat{z}-z_{j}^{0}|<r$ for some $j\leq k$, whence necessarily $\hat{z}=z_{j}^{\varepsilon}(\kappa)$. And $\kappa\geq\kappa^{\varepsilon}_{j}$ since $\mathop{\rm Re}\nolimits\hat{z}\geq 0$.

This completes the proof of Theorem 3.2.

10.3 Simplicity of eigenvalues

It remains to prove Lemma 5.4, which shows in particular that simple roots of $F$ provide simple eigenvalues of $B$.

Proof of Lemma 5.4.

First, we show that the kernel of $B-\lambda I$ is one-dimensional. Recall from Section 5 that whenever $(B-\lambda I)V=0$, then the components $V_{\ell}$ have the form (58) for some constants $c_{1}$, $c_{2}$. More generally, if $V=V(\varphi)$ has the form (58) with $\varphi_{1}=\varphi$, $\varphi_{2}=(A\varphi)^{-1}$, and if $\lambda(\varphi)=(A-\varphi^{-1})(\varphi-1)$, then equations (63)–(64) are equivalent to the equation

[TABLE]

where $e_{j}$ denotes the $j$th standard basis vector, and

[TABLE]

The value $\lambda$ is an eigenvalue if and only if ${\mathcal{D}}(\varphi)$ is singular. The matrix ${\mathcal{D}}(\varphi)$ does not vanish in this case, however, for the following reason. Since $S(\varphi)\neq 0$ and $A\varphi_{1}\varphi_{2}=1$, necessarily $\varphi_{1}$ and $\varphi_{2}$ are distinct and have the same sign. But the function $\varphi f(\varphi)=A\varphi^{M}+\varphi-1$ is convex and cannot have two distinct roots with the same sign. Hence it is not possible that $f(\varphi_{j})=0$ for both $j=1$ and $2$.

It follows that the kernel of $B-\lambda I$ is one dimensional, and the eigenspace is spanned by $V(\varphi)$, taking

[TABLE]

Next, we determine when $\lambda$ is simple, i.e., when it has algebraic multiplicity one. Since $(B-\lambda I)V=0$, this is the case if and only if the equation

[TABLE]

has no solution. Letting ′ denote differentiation with respect to $\varphi$, it follows by differentiating (154) (while keeping $c_{1}$, $c_{2}$ fixed), that

[TABLE]

Now, $\lambda^{\prime}=A-\varphi^{-2}\neq 0$ whenever $\varphi\neq\pm A^{-1/2}$, so it follows that a solution to (156) exists if and only if $\lambda^{\prime}U=V^{\prime}(\varphi)-\hat{U}$ where $\hat{U}$ is a solution to

[TABLE]

As in Section 5, necessarily $\hat{U}_{\ell}=\hat{c}_{1}\varphi_{1}^{M-\ell}+\hat{c}_{2}\varphi_{2}^{M-\ell}$ for some constants $\hat{c}_{1}$, $\hat{c}_{2}$ that satisfy

[TABLE]

Writing $f_{j}=f(\varphi_{j})$, $f_{j}^{\prime}=f^{\prime}(\varphi_{j})\varphi_{j}^{\prime}$ and similarly for $g_{j}$, $g_{j}^{\prime}$, the fact that ${\mathcal{D}}(\varphi)$ is singular means

[TABLE]

and a left null vector is given by $(g_{1},-f_{1})$ or $(g_{2},-f_{2})$ (since ${\mathcal{D}}(\varphi)\neq 0$). Supposing $f_{1}\neq 0$, applying the left null vector to (157) we find that a solution of (157) exists if and only if

[TABLE]

where we used (158) to replace $g_{1}f_{2}$ by $f_{1}g_{2}$. If $f_{2}\neq 0$ similarly the criterion is $0=f_{2}\delta^{\prime}(\varphi)$. Thus an eigenvalue $\lambda$ is simple if and only if $\delta^{\prime}(\varphi)\neq 0$, and this is equivalent to $F^{\prime}(\varphi)\neq 0$. ∎

Acknowledgements

The authors acknowledge support from the Hausdorff Center for Mathematics and the CRC 1060 on Mathematics of emergent effects, Universität Bonn. This material is based upon work supported by the National Science Foundation under grants DMS 1515400 and 1812609, partially supported by the Simons Foundation under grant 395796, and by the NSF Research Network Grant no. RNMS11-07444 (KI-Net).