A Connection Between Orthogonal Polynomials and Shear Instabilities in the Quasi-geostrophic Shallow Water Equations

TL;DR

This paper links orthogonal polynomial roots to shear instability in quasi-geostrophic shallow water equations, proving unique unstable modes for each wave number and estimating their growth rates.

Contribution

It establishes a novel connection between orthogonal polynomials and fluid instability analysis, providing rigorous proofs of unstable modes in QG equations.

Findings

Existence of a unique unstable mode for each wave number 0<k<1

Mathematically rigorous estimates of growth rates

Connection between polynomial roots and shear instability

Abstract

In this paper we demonstrate a connection between the roots of a certain sequence of orthogonal polynomials on the real line and the linear instability of a $x$-directionally homogeneous background velocity profile $u^b(x,y) = \cos(y)$ in the quasi-geostrophic shallow water (QG) equation in a domain with periodic boundaries in the $y$-direction. Using the relationship we establish, we then prove that there exists a unique unstable mode for each horizontal wave number $0<k<1$ and provide mathematically rigorous estimates of the associated growth rate.

Equations

$$f^{\prime \prime }(y)-\left(k^2+\frac{u^b(y{)}^{\prime \prime }-c/\text{Bu}}{u^b(y)-c}\right)f(y)=0.$$

$$q_t+{\psi }_x{q}_y-{\psi }_y{q}_x=0$$

$$\left(\Delta -\frac{1}{\text{Bu}}\right){\psi }_t^p+({\psi }^b{)}^{\prime \prime \prime }{\psi }_x^p-({\psi }^b{)}^{\prime }\Delta {\psi }_x^p=0.$$

$$L_c[f]:=f^{\prime \prime }(y)-\left(\frac{u^b(y{)}^{\prime \prime }-c/\text{Bu}}{u^b(y)-c}\right)f(y).$$

$$\mathcal{D}(L_c):=\left\{f\in \mathcal{H}:\text{f′′(y) exists and }{f}^{\prime \prime }(y)-\left(\frac{u^b(y{)}^{\prime \prime }-c/\text{Bu}}{u^b(y)-c}\right)f(y)\in \mathcal{H}\right\}.$$

$${\sigma }_d(L)$$

$${\sigma }_c(L)$$

$${\sigma }_s(L)$$

$$z_j=\frac{1}{2}{\left(\frac{(j^2+k^2-1)((j+1{)}^2+k^2-1)}{(j^2+k^2+1/\text{Bu})((j+1{)}^2+k^2+1/\text{Bu})}\right)}^{1/2},$$

$$x{p}_n(x)=a_n{p}_{n+1}(x)+b_n{p}_n(x)+a_{n-1}{p}_{n-1}(x).$$

$$u^b(y_0{)}^{\prime \prime }-u^b(y_0)/\text{Bu}=0.$$

$$c_i{\int }_0^{2\pi }\frac{u^b(y{)}^{\prime \prime }-u^b(y)/\text{Bu}}{(u^b(y)-c_r{)}^2+c_i^2}|f(y){|}^{2}dy=0.$$

$$((u^b-c{)}^2{g}^{\prime }(y){)}^{\prime }-\left(k^2(u^b-c{)}^2-\frac{c(u^b-c)}{\text{Bu}}\right)g(y)=0.$$

$$-{\int }_0^{2\pi }(u^b-c{)}^2(|g^{\prime }(y){|}^2+k^2|g(y){|}^2)dy+\frac{c}{\text{Bu}}{\int }_0^{2\pi }(u^b-c)|g(y){|}^{2}dy=0.$$

$${\int }_0^{2\pi }{u}^b(y)Q(y)dy=c_r{\int }_0^{2\pi }Q(y)dy+\frac{1}{2\text{Bu}}{\int }_0^{2\pi }{u}^b|g(y){|}^{2}dy$$

$${\int }_0^{2\pi }{u}^b(y{)}^{2}Q(y)dy=(c_r^2+c_i^2){\int }_0^{2\pi }Q(y)dy+\frac{1}{\text{Bu}}{\int }_0^{2\pi }{u}^b(y{)}^2|g(y){|}^{2}dy$$

$$(r^2-|c{|}^2){\int }_0^{2\pi }Q(y)dy\ge \frac{1}{\text{Bu}}{\int }_0^{2\pi }{u}^b(y{)}^2|g(y){|}^{2}dy$$

$$\left[\begin{array}[]{ccccc}b_{0}&a_{0}&0&0&\dots\\ a_{0}&b_{1}&a_{1}&0&\dots\\ 0&a_{1}&b_{2}&a_{2}&\dots\\ 0&0&a_{2}&b_{2}&\dots\\ \vdots&\vdots&\vdots&\vdots&\ddots\\ \end{array}\right]$$

$$x{p}_n(x)=a_n{p}_{n+1}(x)+b_n{p}_n(x)+a_{n-1}{p}_n(x),n\ge 1.$$

$${\mathcal{E}}_A({\lambda }_i):=\text{span}\{[p_0({\lambda }_j),p_1({\lambda }_j),\dots ,p_{n-1}({\lambda }_j){]}^T\}.$$

$${\sigma }_d(A)=\left\{\lambda :\sum _{n=1}^{\infty }|p_n(\lambda ){|}^2<\infty \right\}.$$

$$\int |x{|}^{n}d\mu (x)<\infty$$

$$\int p_j(x)p_k(x)d\mu (x)=0\text{for all j,k with j=k.}$$

$$m(z):={\int }_{\sigma (A)}\frac{1}{x-z}d\mu (x),z\notin \sigma (A).$$

$$m(z)=-\sum _j\frac{m_j}{z^{j+1}},\text{where}{m}_j=\int x^{j}d\mu (x).$$

$$m(z)=\frac{-1}{z-b_1+a_1^2\left[\frac{-1}{z-b_2+a_2^2[\dots ]}\right]}.$$

$$⟨p(x),q(x){⟩}_{\mu }:={\int }_{\mathbb{R}}p(x)q(x)d\mu (x).$$

$$r_{(n+1),1}<r_{n,1}<r_{(n+1),2}<r_{n,2}<\cdots <r_{n,n}<r_{(n+1),n}.$$

$$\sigma (A{)}_{\ell }\cap [r_{n,j},r_{n,j+1}]\ne 0,$$

$$[r_{n,j},r_{n,j+1}]\cap \bigcup _{n\ge 1}\{\lambda :p_n(\lambda )=0\}.$$

$$\underset{n\to \infty }{\lim }{a}_n=r,\underset{n\to \infty }{\lim }{b}_n=x_0.$$

Full text

1 Introduction

In this paper, we obtain for each wave number $k>0$ rigorous bounds for the eigenvalues of the (modified) Rayleigh stability equation

[TABLE]

with periodic boundary conditions $f(y+2\pi)=f(y)$ in the case that the background velocity field $u^{b}(y)=\cos(y)$. Here by eigenvalues, we refer to the values of $c$ (for fixed $k$) for which the Rayleigh equation has a square-integrable, periodic solution.

The modified version of Rayleigh’s equation above arises in the study of the linear stability of shear flows in the inviscid, incompressible shallow water equations in the limit of small Rossby number Ro. The quasi-geostrophic shallow water equation (QG) is given by

[TABLE]

where $\psi=\psi(x,y,t)$ is the stream function and $q=q(x,y,t)$ is the potential vorticity, related by $q=\Delta\psi-\psi/\text{Bu}$. Linearizing the QG equation around a shear background stream function $\psi^{b}=\psi^{b}(y)$, we obtain the linear partial differential equation for the perturbed stream function $\psi^{p}$

[TABLE]

Note that the background stream function determines a background velocity $u^{b}=-(\psi^{b})^{\prime}$. If $f(y)$ is a solution to Equation 1 with this value of $u^{b}$ for some pair $(k,c)$, then $\psi^{p}(x,y,t)=e^{ik(x-ct)}f(y)$ is a solution of the linearized QG equation. We derive the QG equation from the shallow water equation in the limit of small Ro in Appendix A.2 below. We derive the linearized QG equations and the Rayleigh equation in Appendix B.

For a given value of $k$, there will in general be countably (and often finitely many) values of $c$ for which Equation 1 will have a square-integrable, periodic solution (ie. for which $c$ is an eigenvalue). In this way the choice of background profile $u^{b}(y)$ determines a dispersion relation, ie. a relationship between (complex) frequencies $c$ and wave numbers $k$, which we can represent as a multi-valued function $c(k)$. Note that if $c$ is a value of $c(k)$, then so too is $\overline{c}$. A wave number $k$ is called unstable if one of the values of $c(k)$ is nonreal, and in this case an associated perturbed solution $\psi^{p}$ of the linear QG equation grows exponentially with growth rate $k\cdot\text{Im}(c(k))$.

For some very special background velocity profiles $u^{b}(y)$, the solutions of Equation 1 may be determined explicitly analytically and the dispersion relation thereby determined also. However, for the vast majority of profiles this is not the case. Instead, numerical methods of determining the dispersion relation $c(k)$ are required. One popular method is to replace the differential operators in Rayleigh’s equation with approximations in the form of finite-dimensional linear operators acting on a finite-dimensional vector space. To do so, we can replace the interval $[0,2\pi]$ with a finite grid, and the differential operators with difference operators on this grid wang2012ageostrophic menesguen2012ageostrophic gula2010instabilities . Alternatively, we can expand $f$ in terms of a orthonormal basis for $[0,2\pi]$ and take a finite truncation dolph1958application gallagher1962behaviour orszag1971accurate . Either way, this replaces Equation 1 with a simple eigenvalue problem on a finite-dimensional vector space, and we can imagine that as the accuracy of our approximation is increased that the scattering relations obtained by the various approximations will converge to the true scattering relation $c(k)$.

This presents us with a problem. We have to try to tell which of the eigenvalues of the various approximations are also approximations of the eigenvalues of Equation 1. This problem becomes even more apparent for the wide class of background velocity profiles $u^{b}(y)$ for which Equation 1 has finitely many eigenvalues for each fixed value of $k$. As the precision of our approximations increases, so too does the dimension of the linear system approximating Rayleigh’s equation, resulting in an ever increasing amount of eigenvalues. Even worse, we have no explicit estimates of the rate of convergence.

It is useful to rephrase this problem in the language of differential operators. Consider the Schrödinger operator $L_{c}$ which acts as an unbounded operator on (a dense subset of) the Hilbert space $\mathcal{H}=L^{2}([0,2\pi])$ of square-integrable functions on the interval $[0,2\pi]$ by

[TABLE]

for all $f$ in the domain $\mathcal{D}(L_{c})$ of $L_{c}$

[TABLE]

The spectrum of any unbounded operator $L$ on a Hilbert space $\mathcal{H}$ with domain $\mathcal{D}(L)$ is composed of three parts: a discrete component $\sigma_{d}(L)$, an continuous component $\sigma_{c}(L)$, and a singular component $\sigma_{s}(L)$. Points in each component of the spectrum are characterized as follows:

[TABLE]

where in the above $\mathcal{R}(L-\lambda I)$ denotes the range of $L-\lambda I$. With this in mind, the question of determining the eigenvalues of 1 for all values of $k$ is equivalent to determining for which values of $c$ the operator $L_{c}$ has a positive eigenvalue $k^{2}$ in its discrete spectrum.

In this paper, we consider the shear background profile $u^{b}(y)=\cos(y)$. This profile is complicated enough that the solution to Equation 1 cannot be obtained analytically. However, we will show that explicit and rigorous estimates for the dispersion relation above can be made. Our method for estimating the eigenvalues of Rayleigh’s equation 1 for a background cosine profile is based on relating the eigenvalues to the roots of a sequence of orthogonal polynomials for a certain measure defined on the real line $\mathbb{R}$. Roots of orthogonal polynomials satisfy an interlacing property which we apply to determine the number of eigenvalues and obtain monotonic sequences whose limits are the desired eigenvalues. Specifically we prove the following

Theorem 1.0.1**.**

Let $b_{0}=2z_{0}^{2}+z_{1}^{2}$, $b_{j}=z_{2j}^{2}+z_{2j+1}^{2}$ for $j\geq 1$, and $a_{j}=z_{2j+1}z_{2j+2}$ for all $j\geq 0$, where

[TABLE]

and let $p_{-1}(x)=0,p_{0}(x)=1$ and define $p_{n+1}(x)$ recursively for $n\geq 0$ by

[TABLE]

Take $u^{b}(y)=\cos(y)$ in Rayleigh’s equation 1. Then the following holds

1. (a)

Rayleigh’s equation 1 has complex eigenvalues if and only if $|k|<1$

Furthermore, assuming $|k|<1$, the following is true

1. (b)

For $n$ large enough, the polynomial $p_{n}(x)$ has a unique negative root $-r_{n}$. 2. (c)

The sequence of positive real number $r_{1},r_{2},\dots$ from (b) is monotone increasing and converges to a real number $r$. 3. (d)

The complex eigenvalues of Rayleigh’s equation 1 are given by $\pm i\sqrt{r}$.

The most significant part of Theorem 1.0.1 is that the roots of the polynomials $p_{n}(x)$ converge monotonically. Therefore for each $n$, we get a new, sharper lower bound for the growth rate of instabilities.

2 General Results on Rayleigh’s Equation

Before diving into some of the mathematical background on orthogonal polynomials used in this paper, it makes sense to recount some of the more basic results known for the Rayleigh equation. However, since Equation 1 is not exactly the Rayleigh equation, but a modified version, these results will change in some important ways.

Many general results on Rayleigh’s equation involve finding criteria for the existence of non-real eigenvlaues, and bounds for the growth rates of the associated unstable linear modes. We are dealing with a modified version of Rayleigh’s equation however, because of the presence of the Burger number factor, and this modifies some of the usual instability criteria in interesting ways. For example, a well-known criterion for instability is Rayleigh’s inflection point criterion, which says that for a smooth, shear profile to be linearly unstable, it must have an inflection point. We see in the next theorem, that we no longer need an inflection point if the Burger number is small enough.

Theorem 2.0.1** (Rayleigh’s Inflection Point Criterion).**

Let $u^{b}(y)$ be twice differentiable with continuous second derivative. Suppose that for fixed $k$, Equation 1 has a non-real eigenvalue $c$. Then there exists a point $y_{0}$ satisfying

[TABLE]

Proof.

Let $f(y)$ be an eigenvector for the eigenvalue $c$. Then multiplying Equation 1 by $f^{*}(y)$, integrating by parts, and taking the imaginary part of the resultant identity, we find

[TABLE]

Since $c_{i}\neq 0$, the statement of our theorem is follows from the intermediate value theorem. ∎

We show in this paper that $u^{b}(y)=\cos(y)$ is an unstable profile. Yet, by the above theorem $u^{b}(y)=\cos(y)+2$ is a stable profile for $\text{Bu}<1$. This is indicative of an important difference between Rayleigh’s original equation and Equation 1, namely that the choice of inertial reference frame matters. This is a consequence of the fact that the modified equation was derived in a rotating reference frame. We note that as the rate of rotation is decreased to [math], the Burger number increases to $\infty$ and Rayleigh’s equation takes it’s traditional form; in particular the choice of reference frame no longer matters.

Another well-known result is Howard’s semicircle theorem, which states that for a background velocity profile taking values in the finite interval $[u_{min},u_{max}]$, the eigenvalues of the equation corresponding to unstable modes occur in a circle of radius $(u_{max}-u_{min})/2$ centered at $(u_{max}+u_{min})/2$. However, it seems to be the case that Howard’s semicircle theorem as stated does not hold for the modified Rayleigh equation. This is again due to the fact that the choice of reference frame matters. Instead, we have a modified form, which turns out to be equivalent to Howard’s semicircle theorem in the limit $\text{Bu}\rightarrow\infty$.

Theorem 2.0.2** (Centered Howard’s Semicircle Theorem).**

Let $u^{b}(y)$ be twice differentiable with continuous second derivative, with $u^{b}(y)\in[-r,r]$. Suppose furthermore that for fixed $k$, Equation 1 has an eigenvalue $c$. Then $c$ lies in the complex plane within a circle of radius $r$ of the origin.

Proof.

Make the substitution $g(y)=f/(u^{b}-c)$. With this, Rayleigh’s equation becomes

[TABLE]

Multiplying both sides of this equation by $g^{*}(y)$ and integrating by parts, we obtain

[TABLE]

Taking real and imaginary parts and simplifying with $Q(y)=|g^{\prime}(y)|^{2}+(k^{2}+1/\text{Bu})|g(y)|^{2}$, we obtain

[TABLE]

[TABLE]

where we have written $c$ in terms of its real and imaginary components $c=c_{r}+ic_{i}$. This latter equation in particular says that

[TABLE]

In particular, $|c|^{2}\leq r^{2}$. ∎

3 Jacobi Matrices, Polynomials, and Measures

Next, we will provide a brief summary of the relation between orthogonal matrix polynomials and Jacobi matrices.

Definition 3.0.1**.**

A Jacobi matrix is an infinite or finite square, tri-diagonal matrix of the form

[TABLE]

for some sequences of complex numbers $a_{i},b_{i}\in\mathbb{C}$ with $a_{i}\neq 0$ for all $i$.

Definition 3.0.2**.**

Let $A$ be an $n\times n$ Jacobi matrix of the form 4. The associated sequence of polynomials is the sequence $p_{0}(x),p_{1}(x),\dots$ defined recursively by $p_{0}(x)=1$, $p_{1}(x)=a_{0}^{-1}(x-b_{0})$, and

[TABLE]

The roots of the associated sequence of polynomials describe the eigenvalues of a finite Jacobi matrix, as stated in the next proposition. This is proved in several places, including totik2005orthogonal .

Proposition 3.0.3**.**

Let $A$ be an $n\times n$ Jacobi matrix of the form 4, and let $p_{0}(x),\dots,p_{n}(x)$ be the associated sequence of polynomials. Then the eigenvalues of $A$ are the roots $\lambda_{1},\dots,\lambda_{n}$ of the polynomial $p_{n}(x)$. The corresponding eigenspaces are given by

[TABLE]

The spectrum $\sigma(A)$ of $A$ for infinite Jacobi matrices is more complicated, and consists of both a discrete part $\sigma_{d}(A)$, a continuous part $\sigma_{c}(A)$ and a singular part $\sigma_{s}(A)$. Often the singular part of the spectrum of $A$ is empty, as is the case when $A$ is essentially normal.

Proposition 3.0.4** (arlinskiui2006non ).**

Let $A$ be an infinite Jacobi matrix with bounded coefficients. Then $A$ defines a bounded linear operator on the Hilbert space $\ell^{2}(\mathbb{N})$ whose spectrum $\sigma(A)$ consists of limit points of $\bigcup_{n}\{\lambda:p_{n}(\lambda)=0\}$. The discrete component of the spectrum is

[TABLE]

There is a correspondence between Hermitian Jacobi matrices and probability measures on the real line. Under this correspondence, the support of the measure agrees with the spectrum of the Jacobi matrix. This result is often called Favard’s Theorem and is proved in favard1935polynomes but was also proved by others, including Stieltjes in stieltjes1894recherches .

Theorem 3.0.5** (Favard’s Theoremfavard1935polynomes stieltjes1894recherches ).**

Let $A$ be a Jacobi operator with real, bounded coefficients, and let $p_{0}(x),p_{1}(x),\dots$ be the associated sequence of polynomials. Then $\sigma(A)$ is a bounded subset of $\mathbb{R}$ and there exists a positive measure $\mu$ supported on $\sigma(A)$ satisfying

[TABLE]

and also

[TABLE]

A sequence of polynomials $p_{0}(x),p_{1}(x),\dots$ with $\deg(p_{i})=i$ for all $i$, satisfying the identity $\int p_{j}(x)p_{k}(x)d\mu(x)$ is called a sequence of orthogonal polynomials for the measure $\mu$.

The correspondence between Hermitian Jacobi matrices and the associated probability measures is clarified even further by considering the Stieltjes transform of the measure.

Proposition 3.0.6** (gesztesy1997m ).**

Let $A$ be a Jacobi operator with real, bounded coefficients. Consider the Stieltjes transform of $\mu(x)$

[TABLE]

Then $m(z)$ is a meromorphic function on $\mathbb{C}\backslash\sigma_{c}(A)$ for which the following is true

1. (a)

the set of poles of $m(z)$ is $\sigma_{d}(A)$ 2. (b)

$m(z)$* has the Laurent series expansion*

[TABLE]

If additionally $a_{j}>0$ for all $j$, then

1. (c)

if $a_{j}>0$ for all $j$, then $m(z)$ has the continued fraction expansion

[TABLE]

4 Orthogonal Polynomials, Root Interlacing, and Growth Rates

4.1 Orthogonal Polynomials

Definition 4.1.1**.**

Let $\mu$ be a real, positive measure on the real line. We say that $\mu$** has finite** moments if $\int_{\mathbb{R}}|x|^{n}d\mu(x)<\infty$ for all $n\geq 0$.

A positive measure with finite moments $\mu$ defines a real inner product on the vector space of polynomials $\mathbb{R}[x]$ via the formula

[TABLE]

By Gram-Schmidt orthogonalization, we may construct an orthogonal basis $p_{0}(x),p_{1}(x),\dots$ for $\mathbb{R}[x]$ such that $\deg(p_{n})=n$ for all $n\geq 0$.

Definition 4.1.2**.**

Let $\mu$ be a real, positive measure on the real line with finite moments. A sequence of orthogonal polynomials for $\mu$ is a sequence of polynomials $p_{0}(x),p_{1}(x),\dots$ satisfying $\deg(p_{n})=n$ for all $n\geq 0$ and $\int p_{m}(x)p_{n}(x)d\mu(x)=0$ for all pairs $m,n$ with $m\neq n$.

If $p_{0}(x),p_{1}(x),\dots$ and $q_{0}(x),q_{1}(x),\dots$ are two seqeuences of orthogonal polynomials for $\mu$, then there exists constants $c_{0},c_{1},\dots$ such that $p_{i}(x)=c_{i}q_{i}(x)$ for all $i$. Thus orthogonal polynomials are essentially unique for a given measure.

The converse of Favard’s theorem also holds. Given a probability measure and a sequence of (normalized) orthogonal polynomials, one may prove that the polynomials satisfy a $3$-term recursion relation, ie. they are the same as the associated polynomials of some Jacobi matrix. This is proven in many places, including totik2005orthogonal .

Proposition 4.1.3** (totik2005orthogonal ).**

Let $\mu$ be a real, positive measure on the real line with finite moments. Then there exists a Jordan matrix $J$ whose associated sequence of orthogonal polynomials is a sequence of orthogonal polynomials for $\mu$.

4.2 Root Interlacing

Suppose that $p_{0}(x),p_{1}(x),\dots$ are the orthogonal polynomials for some positive measure $\mu$ on the real line $\mathbb{R}$ with finite moments, and whose support has infinite cardinality. Then $p_{n}(x)$ has $n$ distinct roots for each integer $n$ and between any two roots of $p_{n+1}(x)$ there must lie a root of $p_{n}(x)$.

Proposition 4.2.1** (nevai1989orthogonal ).**

Let $\mu$ be a positive measure on the real line with finite moments, and suppose the cardinality of $\textnormal{supp}(\mu)$ is not finite. Let $p_{0}(x),p_{1}(x),\dots$ be a sequence of orthogonal polynomials for $\mu$. Then the following is true

1. (a)

$p_{n}(x)$* has $n$ distinct, real roots $r_{n,1}<r_{n,2}<\dots<r_{n,n}$ for all positive integers $n$* 2. (b)

the roots of the polynomials satisfy the interlacing property for all $n\geq 1$:

[TABLE]

As a corollary of this, we note that between any two roots of $p_{n}(x)$ there must exist an accumulation point of the spectrum of $J$.

Corollary 4.2.1.1**.**

Let $p_{0}(x),p_{1}(x),\dots$ be a sequence of orthogonal polynomials for a positive measure $\mu$ on the real line whose support has infinte cardinality. Let $r_{n,j}$ be defined as in the statement of the proposition for all $n>0$ and all $1\leq j\leq n$. Then for all $1\leq j<n$ and $n\geq 2$ we have that

[TABLE]

where here $A$ is the associated Hermitian Jacobi matrix and $\sigma(A)_{\ell}$ denotes the set of limit points of the spectrum $\sigma(A)$ of $A$.

Proof.

As a consequence of the interlacing property, the set

[TABLE]

has infinite cardinality. Since $[r_{n,j},r_{n,j+1}]$ is compact, it follows that $\bigcup_{n\geq 1}\{\lambda:p_{n}(\lambda)=0\}$ has a limit point in $[r_{n,j},r_{n,j+1}]$. Hence $[r_{n,j},r_{n,j+1}]\cap\sigma(A)$ is nonempty. This argument applied again to the successive root pairs in $[r_{n,j},r_{n,j+1}]$ actually shows that $[r_{n,j},r_{n,j+1}]\cap\sigma(A)$ has infinitely many points. Hence it has a limit point. ∎

4.3 Growth Rates

We will also require estimates for the growth rates of $p_{n}(x)$ for $x\notin\textnormal{supp}(\mu)$. The main idea is that if $x\notin\textnormal{supp}(\mu)$, then the magnitude of $p_{n}(x)$ grows exponentially in $n$. The following result was obtained by Brian Simanek based on results of Barry Simon simanek2011new simon2004orthogonal

Theorem 4.3.1** (simanek2011new ).**

Let $A$ be a Jacobi matrix of the form 4, and let $\mu(x)$ and $p_{0}(x),p_{1}(x),\dots$ be the associated measure and sequence of orthogonal polynomials. Suppose that $\mu$ has compact support on the real line, and moreover that

[TABLE]

Then for all $z\notin\textnormal{supp}(\mu)$

[TABLE]

If $A$, $\mu$, and $p_{0}(x),p_{1}(x),\dots$ satisfy the assumptions of Theorem (4.3.1), then the support of the absolutely continuous component of $\mu$ is contained in $[x_{0}-2r,x_{0}+2r]$. For real points $z$ outside this interval, and outside the support of $\mu$, the magnitude of $\frac{(z-x_{0})/r+\sqrt{(z-x_{0})^{2}/r^{2}-4}}{2}$ is greater than $1$, and therefore the magnitude of $p_{n}(z)$ grows exponentially fast in $n$ for large $n$.

5 Rayleigh’s Equation and Jacobi Matrices

Suppose that we wish to find square-integrable solutions of Rayleigh’s equation 1 for given $k$. Any such solution $f(y)$ has a Fourier expansion:

[TABLE]

and inserting this into Rayleigh’s equation along with the Fourier expansion of $u^{b}(y)$ and simplifying yeilds the following eigenvalue problem for $u^{b}(y)$:

[TABLE]

where here $*$ denotes the (discrete) convolution operator.

5.1 Instability of the Cosine Profile

We next consider specifically the case that $u^{b}(y)=\cos(y)$. In this case, the eigenvalue problem 5 becomes

[TABLE]

Setting

[TABLE]

we find

[TABLE]

Hence $[\dots,\hat{g}(-2),\hat{g}(-1),\hat{g}(0),\hat{g}(1),\dots]^{T}$ is an eigenvector with eigenvalue $c$ of the bi-infinite tri-diagonal matrix

[TABLE]

Here we have used the fact that $\hat{q}(j)=\hat{q}(-j)$, and that therefore $z_{-j}=z_{j-1}$ for $j<0$. Note that $B$ is not a Jacobi matrix, because it is bi-infinite.

We next show how to relate the eigenvalues of $B$ to the eigenvalues of an infinite Hermitian Jacobi matrix. The symmetry of $B$ implies that each eigenspace is invariant under the involution

[TABLE]

Therefore $B$ must have an eigenvector $\vec{v}$ with eigenvalue $c$ satisfying $\sigma(\vec{v})=\pm\vec{v}$, ie. each eigenvalue of $B$ must have either a $\sigma$-symmetric or $\sigma$-skew symmetric eigenvector. If $\sigma(\vec{v})=\vec{v}$, then $cv(0)=2z_{0}v(1)$ and therefore $[v(0),v(1),v(2),\dots]^{T}$ is an eigenvector of

[TABLE]

with eigenvalue $c$. If we square $\widetilde{B}$ and take the imaginary part, then we see that $[v(0),v(1),v(2),\dots]^{T}$ is an eigenvector of

[TABLE]

It follows from the above checkerboard pattern that $[v(1),v(3),v(5),\dots]^{T}$ is an eigenvector with eigenvalue $c^{2}$ of the infinite Hermitian Jacobi matrix

[TABLE]

The $\sigma$-symmetric eigenvectors are significant, because the unstable modes exhibit this symmetry.

Lemma 5.1.1**.**

Suppose that $c$ is a non-real eigenvalue of Equation 6. Then the associated eigenspace of $B$ is $1$-dimensional and consists of a single $\sigma$-symmetric vector.

Proof.

Suppose that $c$ is a non-real eigenvalue of Equation 6, and let $\vec{v}=[\dots,v(-1),v(0),v(1),\dots]$ be an eigenvector of $B$ with this eigenvalue. If $v(0)=0$, then $[v(1),v(2),v(3),\dots]$ is an eigenvector with eigenvalue $c$ of the infinite Hermitian Jacobi matrix

[TABLE]

However, this implies that $c$ is real, which is a contradiction. Therefore $v(0)\neq 0$. If the eigenspace with eigenvalue $c$ contains more than one linearly independent vector, then by taking an appropriate linear combination we can obtain a nonzero eigenvector with $v(0)=0$, which again leads to a contradiction. Therefore the eigenspace of $c$ must be one-dimensional. Since the eigenspace must contain a $\sigma$-symmetric or $\sigma$-skew symmetric vector, and $\sigma$-skew symmetric vectors satisfy $v(0)=0$, we also see that the eigenspace has a $\sigma$-symmetric eigenvector. ∎

Lemma 5.1.2**.**

Let $\mu$ be the measure associated with the Hermitian Jacobi matrix $A$. Then the support of the absolutely continuous component of $\mu$ is $[0,1]$.

Proof.

The support of the absolutely continuous component of $\mu$ is equal $\sigma_{c}(A)$. Furthermore $\sigma_{c}(A)=\sigma_{c}(\widetilde{B}^{2})$ and $\sigma_{c}(\widetilde{B}^{2})=\sigma_{c}(\widetilde{B})^{2}$. Moreover, $\widetilde{B}$ differs from the Chebyshev operator

[TABLE]

by a compact operator $K$, ie. $\widetilde{B}=C+K$. Hence $\sigma_{c}(\widetilde{B})=\sigma_{c}(C)=[-1,1]$, and it follows that $\sigma_{c}(A)=\sigma_{c}(\widetilde{B})^{2}=[0,1]$. ∎

Proposition 5.1.3**.**

Let $c\in\mathbb{C}$ be non-real. Then $c$ is an eigenvalue of Equation 1 if and only if $c$ is an element of the discrete spectrum of $A$.

Proof.

Suppose $c^{2}\in\sigma_{d}(A)$. Then there exists $v\in\ell^{2}(\mathbb{N})$ with $\vec{v}=[v(0),v(1),v(2),\dots]^{T}$ an eigenvector of $J$ with eigenvalue $c$. Define $v(-j)=v(j)$ and

[TABLE]

Then $[\dots,u(-2),u(-1),u(0),u(1),\dots]^{T}$ is an eigenvector of $B$. Moreover this extends $v$ to an element of $\ell^{2}(\mathbb{Z})$. It follows that $\widehat{f}(\ell)$ satisfies Equation 6 for $\widehat{f}(0)=0$ and for $\ell\neq 0$,

[TABLE]

Moreover, since $\widehat{f}(\ell)$ is the product of two functions in $\ell^{2}(\mathbb{Z})$ we have that $\widehat{f}(\ell)$ is in $\ell^{2}(\mathbb{Z})$.

To prove the converse, suppose that $c$ is an eigenvalue of Equation 1, and let $\widehat{f}(\ell)$ be the Fourier coefficients of the associated eigenfunction $f(y)$. Then $c^{2}$ is an eigenvalue of the Jacobi matrix $A$, whose eigenvector is $[v(1),v(3),v(5),\dots]^{T}$ for $v(j)=\widehat{f}(j)w(j)$ where

[TABLE]

However, the eigenspace of $A$ with eigenvalue $c^{2}$ is exactly $[p_{0}(c^{2}),p_{1}(c^{2}),p_{2}(c^{2}),\dots]$ where $p_{0}(x),p_{1}(x),p_{2}(x),\dots$ is the sequence of polynomials associated to $J$. This implies that for some constant $K$

[TABLE]

and therefore that $[p_{0}(c^{2})/w(0),p_{1}(c^{2})/w(1),\dots]\in\ell^{2}(\mathbb{N})$. The continuous support of the measure $\mu$ associated to $A$ is contained on the positive real axis. Therefore if $c^{2}\notin\sigma_{d}(A)$ then since $c^{2}$ is not a positive real number $c^{2}\notin\sigma(A)$. The growth rate of $p_{n}(x)$ for $x\notin\textnormal{supp}(\mu)$ is exponential in $n$ by Theorem 4.3.1, and since $w(n)$ has polynomial growth, this contradicts the possibility that $[p_{0}(c^{2})/w(0),p_{1}(c^{2})/w(1),\dots]$ is in $\ell^{2}(\mathbb{N})$. This completes the proof. ∎

Proposition 5.1.4**.**

Let $A$ be as above, and let $p_{0}(x),p_{1}(x),\dots$ be the sequence of polynomials associated to $A$, and let $r_{n1}<r_{n2}<\dots<r_{nn}$ be the roots of $p_{n}(x)$. If $|k|<1$, then for all $n$ large enough, the polynomial $p_{n}(x)$ has exactly one negative root $r_{n1}$. The spectrum of $A$ has exactly one negative value $r$, where $r$ is the limit of the monotone decreasing sequence $r_{11},r_{21},r_{31},\dots$.

Proof.

Fix $k$ with $|k|<1$. Since $A$ is Hermitian, its spectrum is a subset of $\mathbb{R}$. One may verify that $p_{n}(x)$ has a negative root for $n$ large enough. If $p_{n}(x)$ has more than one negative root for some $n$, then by Corollary 4.2.1.1, then $\sigma(A)$ will have a negative limit point, and thus $\sigma(A)$ will have infinitely many negative elements. Since the continuous part of the spectrum of $A$ is positive, this means that $\sigma_{d}(A)$ is infinite. This implies that Rayleigh’s equation has infinitely many eigenvalues $c$ for this value of $k$.

However, the number of eigenvalues of Rayleigh’s equation is finite, as can be seen by the fact that the eigenvalues of $A$ correspond to the poles of the Stieltjes transform $m(z)$ of the measure $\mu$ associated to $A$. Since the function $m(z)$ is meromorphic, its poles are discrete. Therefore by the centered Howard’s semicircle theorem, there are finitely many of them. Alternatively for $\text{Bu}=\infty$, one may use the result of Howard that the number of unstable modes is bounded by the number of inflection points of the background profile $u^{b}(y)$ howard1964number . Therefore $p_{n}(x)$ has at most one negative eigenvalue for each $n$. Thus $p_{n}(x)$ has exactly one negative root $r_{n,1}$ for $n$ large enough, and by root interlacing (Proposition 4.2.1) $r_{n,1}$ is monotone decreasing. This proves the proposition. ∎

We now turn to the proof of Theorem 1.0.1, as stated in the introduction. We have essentially proved it in the previous two propositions.

Proof of Theorem 1.0.1.

1. (a)

If $|k|\geq 1$, then $B$ is a Hermitian matrix, and therefore the eigenvalues are all real. Since the eigenvalues of $B$ determine the eigenvalues of Rayleigh’s equation, this proves (a). 2. (b)

This is a restatement of the conclusion of the previous proposition. 3. (c)

This is a restatement of the conclusion of the previous proposition. 4. (d)

This follows from Proposition 5.1.3.

∎

6 Numerical Results

In this section we will provide a numerical verification of the result of Theorem 1.0.1. We will also provide a demonstration of the change in the behavior of the growth rate curve as a function of the Burger number Bu. The numerical calculations were carried out in python using numpy for the root calculations.

7 Conclusions

In this paper we explore the linear stability of the QG shallow water equation for a shear cosine profile with periodic boundary conditions, and derive bounds for the growth rate of the instabilities. We relate the instabilities we see to the roots of a sequence of orthogonal polynomials, and our computational results verify our findings.

There are many unresolved questions that we would like to answer in the future, some of which we list here:

1. (a)

Is there an exact expression for the measure associated to our polynomials? 2. (b)

Rayleigh’s equation with the cosine profile $u^{b}(y)$ can be transformed into a Heun differential equation via the substitution $z=\cos(y)$, and therefore the associated Heun functions may be used as generating functions for our polynomials. Do our polynomials comprise a “nice” basis for the expansion of the associated Heun functions? 3. (c)

Is there an analytic expression for growth rates? Or else, can we establish estimates of the error in our approximation?

We would also like to relate the linear instability we calculated here to the linear instability of the cosine profile for the full shallow water equations.

Acknowledgement

The author is grateful for a summer research internship at Los Alamos National Laboratory, where part of the research for this paper was conducted.

Appendix A Derivation of the Quasi-Geostrophic Shallow Water Equation

A.1 Shallow Water Equations

In this paper, we consider the inviscid shallow water equations in a doubly periodic domain with a flat bottom. The nondimensional form of the shallow water equations is given by

[TABLE]

where $\vec{u}$ is the velocity field, $\eta$ is the free surface height, $\text{Ro}=U/(fL)$ is the Rossby number, $\text{Fr}=U/(gH)$ is the Froude number, $\text{Bu}=(\text{Ro}/\text{Fr})^{2}$ is the Burger number, $U$ is the characteristic magnitude of the velocity field, $H$ is the mean depth, $L$ is the characteristic length scale, $f$ is the Coriolis frequency, and $g$ is the (reduced) gravitational acceleration.

A.2 Shallow Water QG Equation

In mid-latitude regions of the ocean or atmosphere at length scales relevent to geophysical flows, the Rossby number is typically small $\text{Ro}\ll 1$. In this situation, solutions to the shallow water equations are predominantly geostrophically balanced, eg. $\vec{u}\times\hat{z}=\nabla\eta+\mathcal{O}(\text{Ro})$. With this in mind, we next derive a geostrophically balanced model which approximates the shallow water equations in the limit of small Rossby number.

Consider a perturbation expansion in terms of the Rossby number Ro to obtain a balanced model for the evolution. We expand the velocity and free surface height fields as

[TABLE]

[TABLE]

Inserting this back into the shallow water equations and comparing similar powers of Ro, we obtain two equations describing the leading balance for small Rossby number:

[TABLE]

as well as the following relations for all $i\geq 0$

[TABLE]

To understand Equations (10) and (11) better, we introduce the following decomposition of the velocity field into a geostrophically balanced and imbalanced component:

[TABLE]

Using this decomposition, Equations (10) simply say that $u^{0}_{1}=0$. In other words, to leading order in Ro the velocity field is geostrophically balanced.

Next, note that

[TABLE]

and therefore

[TABLE]

Equation 13 shows that the time evolution of $\vec{u}^{i}$ is determined by $\vec{u}^{i}$ and the time evolution of the ageostrophic $\vec{u}^{i}_{1}$ field.

We next look at the field $\vec{u}^{i}_{1}$ more closely, examining in particular its divergence and curl. Combining Equation 13 with the last of the four equations, we obtain the Helmholtz decomposition for $\vec{u}^{i+1}_{1}\times\widehat{z}$:

[TABLE]

where here $p^{i}$ is a pressure term given by

[TABLE]

where $\eta_{tt}^{-1}:=0$. Thus if we know $\eta^{i}$ and $\vec{u}^{i}$ for $i\leq n$, then we may obtain $\vec{u}^{n+1}_{1}$.

In the special case that $n=0$, Equation 10, allows us to determine $\eta^{0}$ from $\vec{u}^{0}$. Hence by combining Equation (14) and Equation (13), we obtain a closed equation for the time evolution of $\vec{u}^{0}$, which we refer to as the shallow water quasi-geostrophic equation

[TABLE]

Note that by taking the curl of the above equation, we obtain the usual potential vorticity form of the QG shallow water equation

[TABLE]

[TABLE]

Appendix B Shear Instabilities and Rayleigh’s Equation

The quasi-geostrophic shallow water (QG) equation describes the motion of a vertically homogeneous fluid in a rotating reference frame for fixed Burger number Bu in the limit of small Rossby number Ro. The QG equation is given by

[TABLE]

where here $q=q(x,y,t)$ is the potential vorticity, $\psi=\psi(x,y,t)$ is the stream function, and $J$ is the Jacobian

[TABLE]

The potential vorticity and the stream function are related to each other by

[TABLE]

We will consider the linear stability of solutions to the QG equation in a domain satisfying (normalized) periodic boundary conditions in the $y$ direction

[TABLE]

To obtain an equation for the linear stability of a given solution $\psi=\psi^{b}$ of 17, we consider a solution of 17 of the form $\psi=\psi^{b}+\psi^{p}$, where $\psi^{b}$ is a base state solution and $\psi^{p}$ is a perturbation. Inserting this back into the equation and ignoring quadratic terms in the perturbation, we obtain the linearized QG shallow water equation

[TABLE]

In the case of a shear instabilities, the background solution is homogeneous in one of the directions, which we take to be the $x$-direction. Then the background state is $\psi^{b}=-\int u^{b}(y)dy$ for some background velocity profile $u^{b}(y)$. Then the linear equation reduces to

[TABLE]

Rewriting this in terms of stream functions only, this says:

[TABLE]

Since the coefficients of the above differential equation are constant in $x$ and $t$, it makes sense to look for solutions of the form

[TABLE]

for some unknown function $f(y)$. Substituting this in, we obtain a modified form of Rayleigh’s Equation 1:

[TABLE]

This differs from the usual Rayleigh equation in the inclusion of the Bu term. As Bu increases however, this results in the usual form of Rayleigh’s equation.

Acknowledgement