(Evaluate block to execute LaTeX definitions) $ \newcommand{\ybar}{\overline{y}} $ $ \renewcommand{\d}[2]{\frac{d #1}{d #2}} $ $ \newcommand{\dd}[2]{\frac{d^2 #1}{d #2^2}} $ $ \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} $ $ \newcommand{\pdd}[2]{\frac{\partial^2 #1}{\partial #2^2}} $ $ \renewcommand{\b}{\beta} $ $ \newcommand{\m}{\mu} $ $ \renewcommand{\v}[1]{\mathbf{#1}} $ $ \newcommand{\N}{\mathcal{N}} $ $ \renewcommand{\l}{\lambda} $ $ \newcommand{\ol}[1]{\overline{#1}} $
This is the supplementary material for the paper Dynamics of zonal flows: failure of wave-kinetic theory, and new geometrical optics approximations. The article may also be found on arXiv here.
The dispersion relation for zonostrophic instability is shown in Figure 2 in the paper. We show how to derive the general form of the dispersion relation for zonostrophic instability, and then how to compute it numerically for a given set of parameters. Python code is included.
In each case, we reduce the general form of the dispersion relation to a simpler form by using a fluctuation spectrum that consists of a thin ring in wavevector space $\sim \delta(k - k_f)$ centered at $k_f$. This reduction was first considered for CE2 by Srinivasan and Young (2012); more details can be found in Parker (Ph.D. Thesis, 2014).
If this document is being read as a PDF or webpage, an interactive ipython notebook can be found here, or at a direct link (right click, save as). In the ipython notebook one can follow along the mathematics and run the Python code inside the notebook using the Jupyter software. To view it non-interactively within a browser without the Jupyter software, see here.
This document uses the geophysical convention for coordinates. To use the standard plasma coordinates, let $(x,y,\ybar,\b,U) \mapsto (-y, x, \ol{x}, \kappa, -U)$.
We start from Eq. (3.2). There is a homogeneous equilibrium, independent of $\ybar$, at $\N_H(k_x,k_y) = F(k_x,k_y) / 2\b\m$, $U=0$. Using the form \begin{equation} \N = \N_H + e^{iq\ybar} e^{\l t} \N_1(k_x, k_y), \end{equation} \begin{equation} U = e^{iq\ybar} e^{\l t} U_1, \end{equation} we linearize Eq. (3.2) for perturbations about the equilibrium. We obtain \begin{equation} \l \N_1 = iqk_x U_1 \pd{\N_H}{k_y} - \frac{2i\b q k_x k_y}{\ol{k}^4} \N_1 - 2\m \N_1 \end{equation} \begin{equation} (\l + \m) U_1 = iq \int \frac{d\v{k}}{(2\pi)^2} \frac{k_x k_y}{\ol{k}^4} \b \N_1. \end{equation}
The first equation is solved for $\N_1$ in terms of $U_1$ and then substituted into the second equation. This procedure yields a single nonlinear equation for the eigenvalues $\l$ corresponding to zonostrophic instability. One finds \begin{equation} \l + \m = -q^2 \int \frac{d\v{k}}{(2\pi)^2} \frac{k_x^2 k_y}{(\l + 2\m)\ol{k}^4 + 2i\b q k_x k_y} \b \pd{\N_H}{k_y} \end{equation}
When $q$ is large, $\l \sim q$.
We follow the same procedure, starting from Eq. (4.1). After linearizing about the homogeneous equilibrium $W_H(k_x,k_y) = F(k_x,k_y) / 2\m$, we find \begin{equation} \l W_1 = iqk_x U_1 \pd{}{k_y} \left[ \left( 1 - \frac{q^2}{\ol{k}^2} \right) W_H \right] - \frac{2i\b q k_x k_y}{\ol{k}^4} W_1 - 2\m W_1, \end{equation} \begin{equation} (\l + \m)U_1 = iq \int \frac{d\v{k}}{(2\pi)^2} \frac{k_x k_y}{\ol{k}^4} W_1. \end{equation}
Solving for $W_1$ and substituting into the second equation yields \begin{equation} \l + \m = -q^2 \int \frac{d\v{k}}{(2\pi)^2} \frac{k_x^2 k_y}{(\l + 2\m)\ol{k}^4 + 2i\b q k_x k_y} \pd{}{k_y} \left[ \left( 1 - \frac{q^2}{\ol{k}^2} \right) W_H \right] \end{equation}
We follow the same procedure, starting from Eq. (4.3). We linearize about the homogeneous equilibrium $\N_H(k_x,k_y) = F(k_x,k_y) / 2\m \b$, and find
\begin{equation} \l \N_1 = iqk_x U_1 \left( 1 - \frac{q^2}{\ol{k}^2} \right) \pd{\N_H}{k_y} - \frac{2i\b q k_x k_y}{\ol{k}^4} \N_1 - \frac{q^2 F}{\b^2} U_1 - 2 \m \N_1, \end{equation} \begin{equation} (\l + \m) U_1 = iq \int \frac{d\v{k}}{(2\pi)^2} \frac{k_x k_y}{\ol{k}^4} (\b \N_1 - U_1'' \N_H). \end{equation}
Under the assumption that $F$ and $\N_H$ are even in $k_x$ and in $k_y$, the integral over $U_1'' \N_H$ vanishes. Solving for $\N_1$ gives \begin{equation} \N_1 = \frac{iqk_x U_1 (1 - q^2/\ol{k}^2)}{\l + 2\m + 2i\b q k_x k_y / \ol{k}^4} \pd{\N_H}{k_y} - \frac{q^2 F U_1}{\b^2 (\l + 2\m + 2i\b q k_x k_y / \ol{k}^4)} \end{equation} Plugging this into the zonal flow equation yields the dispersion relation, \begin{equation} \l + \m = -q^2 \int \frac{d\v{k}}{(2\pi)^2} \left[ \frac{\b k_x^2 k_y (1 - q^2/\ol{k}^2)}{(\l+2\m) \ol{k}^4 + 2i\b q k_x k_y} \pd{\N_H}{k_y} + \frac{qF}{\b} \frac{i k_x k_y}{(\l + 2\m) \ol{k}^4 + 2i\b q k_x k_y} \right] \end{equation}
The last term here has no analog in the other models and its presence is a consequence of the invalidity of this model --- wave action is not conserved during this instability.
The same procedure can be applied to Eq. (2.3) with a little more algebraic complexity. See Srinivasan and Young (2012) or Parker (PhD Thesis, Section 3.2) for the derivation. Using Eqs. (3.25) and (3.26) of Parker (PhD Thesis), and taking $\ol{q}^2 = q^2$ (corresponding to $L_d^{-2} = 0$ for the zonal flows) and the limit of zero viscosity $\nu=0$, the dispersion relation is \begin{equation} \l + \m = -q \int \frac{d\v{k}}{(2\pi)^2} \frac{k_x^2 k_y}{(\l + 2\m) \ol{h}^2_1 \ol{h}^2_2 + 2i\b qk_xk_y } \left[ \left( 1 - \frac{q^2}{\ol{h}_1^2} \right) W_H(k_x, k_y + \tfrac{1}{2}q) - \left(1 - \frac{q^2}{\ol{h}_2^2} \right) W_H(k_x, k_y - \tfrac{1}{2}q) \right], \end{equation}
where $\ol{h}_{1,2}^2 = k_x^2 + (k_y \pm \tfrac{1}{2} q)^2 + L_d^{-2}$. After a Taylor expansion for small $q$, this equation reduces to the CE2-GO dispersion relation.
A thin-ring forcing isotropic in wavevector space is a convenient simplification. We use \begin{equation} F(k) = 4\pi \varepsilon k_f \delta(k-k_f). \end{equation} If $L_d^{-2} = 0$, then $\varepsilon$ is equal to the energy input into the system by the forcing. With this forcing, the dispersion relations above can be simplified by converting the integrals to polar coordinates. One is left with only one integral, the angle integral, that must be computed numerically.
For some of these computations we will integrate by parts, and we need the relation
\begin{equation} \pd{}{k_y} \frac{k_y}{(\l + 2\m) \ol{k}^4 + 2i\b qk_x k_y} = \frac{(\l + 2\m) \ol{k}^2 (\ol{k}^2 - 4k_y^2)}{[(\l + 2\m) \ol{k}^4 + 2i\b qk_x k_y]^2} \end{equation}
First we apply integration by parts in $k_y$ for the dispersion relation, then substitute $\N_H = F/2\m \b$. This yields
\begin{equation} \l + \m = q^2 \int \frac{d\v{k}}{(2\pi)^2} \frac{F}{2\m} \frac{(\l + 2\m) \ol{k}^2 k_x^2 (\ol{k}^2 - 4k_y^2)} {[(\l + 2\m) \ol{k}^4 + 2i\b qk_x k_y]^2} \end{equation}
For isotropic $F$, we use polar coordinates with $k_x = k\sin \phi$, $k_y = -k \cos \phi$. The dispersion relation becomes
\begin{equation} \l + \m = q^2 \int_0^\infty dk\, \frac{F(k)}{4\pi \m} k^2 (\l + 2\m)(1+m) \times \int_0^{2\pi} \frac{d\phi}{2\pi} \frac{\sin^2 \phi (1 + m - 4\cos^2 \phi)}{ [(\l + 2\m)k^2 (1+m)^2 - 2i\b q \cos\phi \sin\phi]^2} \end{equation}
where $m = (k L_d)^{-2}$. Substituting the thin-ring forcing $F = 4\pi \varepsilon k_f \delta(k-k_f)$, we obtain the final form of the dispersion relation,
\begin{equation} \l + \m = q^2 \frac{\varepsilon}{\m} k_f^4 (\l + 2\m) (1 + m_f) \int_0^{2\pi} \frac{d\phi}{2\pi} \frac{(1 + m_f - 4\cos^2\phi) \sin^2 \phi}{[(\l + 2\m)k_f^2 (1 + m_f)^2 - 2i\b q \cos\phi \sin \phi]^2} \end{equation}
with $m_f = (k_f L_d)^{-2}$.
Now, this is a nonlinear equation that we can solve for $\l$. For most choices of the spectrum W_H, an unstable perturbation with $\text{Re } \l > 0$, has an eigenvalue $\l$ that is purely real. It can be seen that if $\l$ is real, the imaginary part of the $\phi$ integral vanishes, so that one can work only with the real part. Working only with real $\l$ rather than complex $\l$ has the advantage that robust one-dimensional root finders can be used. For the integrand above, we have \begin{equation} \frac{A}{(B + iC)^2} = \frac{A(B^2 - C^2)}{(B^2 + C^2)^2} + \text{imag. part} \end{equation} where $A$, $B$, and $C$ are assumed real. If $\l$ is real, the imaginary part vanishes after integration, so we drop this term. We do indeed find real solutions for $\l$.
Below is Python code to solve the dispersion for zonostrophic instability within the Asymptotic WKE.
In [2]:
# Some python setup
from __future__ import division
import numpy as np
import scipy.optimize
import scipy.integrate
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
matplotlib.style.use('ggplot')
plt.rcParams['figure.figsize'] = (8.0, 8.0)
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 14}
plt.rc('font', **font)
In [26]:
def lambda_of_params(beta, mu, epsilon, Ldinvsq, kf, q):
"""Given parameters, calculate the eigenvalue lambda for the asymptotic WKE."""
funres = lambda lamb: disp_relation_residual(beta, mu, epsilon, Ldinvsq, kf, q, lamb)
domainmin = -2*mu
domainmax = 100
lamb = scipy.optimize.brentq(funres, domainmin, domainmax) #1D root finder
return lamb
def disp_relation_residual(beta, mu, epsilon, Ldinvsq, kf, q, lamb):
"""Compute the residual for the dispersion relation. lamb must be real."""
mf = Ldinvsq / (kf*kf)
lhs = lamb + mu
rhs = q*q*epsilon/mu * kf**4 * (lamb + 2*mu) * (1+mf) * polarint(beta, mu, epsilon, Ldinvsq, kf, q, lamb)
res = rhs - lhs
return res
def polarint(beta, mu, epsilon, Ldinvsq, kf, q, lamb):
"""Carry out the polar integral"""
fun = lambda phi: polarint_integrand(beta, mu, epsilon, Ldinvsq, kf, q, lamb, phi)
out, err = scipy.integrate.quad(fun, 0, 2*np.pi, epsabs=2e-6, epsrel=2e-6)
out = out / (2*np.pi)
return out
def polarint_integrand(beta, mu, epsilon, Ldinvsq, kf, q, lamb, phi):
""" As explained above, integrand is A/(B+iC)^2. Keep only the real part.
So use A(B^2 - C^2) / (B^2 + C^2)^2."""
c = np.cos(phi)
s = np.sin(phi)
mf = Ldinvsq / (kf*kf)
A = (1 + mf - 4*c*c) * s * s
B = (lamb + 2*mu) * kf**2 * (1+mf)**2
C = -2*beta*q*c*s
out = A*(B*B - C*C) / (B*B + C*C)**2
return out
# Now, let's go and calculate the dispersion relation lambda(q)
beta = 1
mu = 0.02
epsilon = 1
Ldinvsq = 1
kf = 1
q1 = np.logspace(-4, 0, 150)
q2 = np.linspace(1.02, 2.2, 150)
q = np.concatenate([q1, q2])
lamb = np.zeros_like(q)
for j in range(len(q)):
lamb[j] = lambda_of_params(beta, mu, epsilon, Ldinvsq, kf, q[j])
fig = plt.figure()
plt.plot(q, lamb)
plt.xlabel(r'$q$')
plt.ylabel(r'$\lambda$ (Asymptotic WKE)')
Out[26]:
The only difference for the general dispersion relation in CE2-GO as compared to the asymptotic WKE is the presence of a $(1 - q^2 / \ol{k}^2)$ term. For the same isotropic thin-ring forcing, the dispersion relation reduces to the same result except that in CE2-GO there is an additional factor of \begin{equation} 1 - \frac{q^2}{k_f^2 (1+m_f)} \end{equation} in front of the polar integral.
In [6]:
def CE2GO_lambda_of_params(beta, mu, epsilon, Ldinvsq, kf, q):
"""Given parameters, calculate the eigenvalue lambda."""
funres = lambda lamb: CE2GO_disp_relation_residual(beta, mu, epsilon, Ldinvsq, kf, q, lamb)
domainmin = -2*mu
domainmax = 100
lamb = scipy.optimize.brentq(funres, domainmin, domainmax) # 1D root finder
return lamb
def CE2GO_disp_relation_residual(beta, mu, epsilon, Ldinvsq, kf, q, lamb):
"""Compute the residual for the dispersion relation. lamb must be real."""
mf = Ldinvsq / (kf*kf)
lhs = lamb + mu
rhs = q*q*epsilon/mu * kf**4 * (lamb + 2*mu) * (1+mf) * (1-q*q/(kf*kf*(1+mf))) * CE2GO_polarint(beta, mu, epsilon, Ldinvsq, kf, q, lamb)
res = rhs - lhs
return res
def CE2GO_polarint(beta, mu, epsilon, Ldinvsq, kf, q, lamb):
"""Carry out the polar integral. Same as for asymptotic WKE"""
fun = lambda phi: CE2GO_polarint_integrand(beta, mu, epsilon, Ldinvsq, kf, q, lamb, phi)
out, err = scipy.integrate.quad(fun, 0, 2*np.pi, epsabs=2e-6, epsrel=2e-6)
out = out / (2*np.pi)
return out
def CE2GO_polarint_integrand(beta, mu, epsilon, Ldinvsq, kf, q, lamb, phi):
""" As explained above, integrand is A/(B+iC)^2. Keep only the real part.
So use A(B^2 - C^2) / (B^2 + C^2)^2."""
c = np.cos(phi)
s = np.sin(phi)
mf = Ldinvsq / (kf*kf)
A = (1 + mf - 4*c*c) * s * s
B = (lamb + 2*mu) * kf**2 * (1+mf)**2
C = -2*beta*q*c*s
out = A*(B*B - C*C) / (B*B + C*C)**2
return out
# Now, let's go and calculate the dispersion relation lambda(q)
beta = 1
mu = 0.02
epsilon = 1
Ldinvsq = 1
kf = 1
q1 = np.logspace(-4, 0, 150)
q2 = np.linspace(1.02, 2.2, 150)
q = np.concatenate([q1, q2])
lamb = np.zeros_like(q)
for j in range(len(q)):
lamb[j] = CE2GO_lambda_of_params(beta, mu, epsilon, Ldinvsq, kf, q[j])
fig = plt.figure()
plt.plot(q, lamb)
plt.xlabel(r'$q$')
plt.ylabel(r'$\lambda$ (CE2-GO)')
Out[6]:
For the WKE, the dispersion relation is calculated similarly (details and code are omitted here). However, a brief note is warranted on how the dispersion relation in Figure 2 is obtained. Because the WKE is invalid for calculating the dispersion relation, it is not wholly surprising to see strange behavior. When one calculates the dispersion relation about an equilibrium that balances forcing and dissipation, one finds that for some values of $q$, the eigenvalue $\l$ becomes complex, even though the correct answer is that $\l$ is real. The source of this can be traced to a linearization of the $F/(\b - U'')$ term, which gives $F U_1'' / \b^2$.
The WKE dispersion relation in Figure 2 is obtained by neglecting the $F U_1'' / \b^2$ term, which is justifiable through an alternate procedure. There are two ways of obtaining the equilibrium incoherent spectrum that we linearize about. The first way to realize the equilibrium is the route we have been using, which is a balance between external forcing and linear dissipation. An alternate route is to remove forcing and dissipation, in which case any homogeneous spectrum is an equilibrium. This alternate procedure yields effectively the same dispersion relation, although within the WKE linearization it has the effect of removing the $FU_1''/\b^2$ term in the linearization. This procedure yields a real $\l$ and is the one shown in Figure 2.
The main point here is solely to demonstrate quantitatively that the WKE is not correct, which is to be expected because one had to assume that $U$ varied slowly in time in order to derive the WKE.
The CE2 dispersion relation is obtained in a similar way (details and code omitted). Since there is no $k_y$ derivative, one does not use an integration by parts but rather a shift of the integration variable, after which the isotropic form for $W_H$ can be inserted. For details, see Srinivsan and Young (2012) or Parker (Ph.D. thesis, section 3.2.3).