1. Introduction
Modons, or dipolar vortices, are common coherent structures in the ocean and atmosphere. In the ocean, they are remarkably stable, existing over long time scales and transporting fluid over large distances (Nycander & Isichenko Reference Nycander and Isichenko1990). Modons are typically observed at the ocean surface, although they can extend over several kilometres into the deep ocean (Ni et al. Reference Ni, Zhai, Wang and Hughes2020). In the study of atmospheric dynamics, modons have been used to model various atmospheric processes such as atmospheric blocking (McWilliams Reference McWilliams1980) and Madden–Julian oscillation events (Rostami & Zeitlin Reference Rostami and Zeitlin2021).
Early studies of modon solutions to the two-dimensional Euler equation were carried out independently by Lamb and Chaplygin (Lamb Reference Lamb1932; Meleshko & van Heijst Reference Meleshko and van Heijst1994) and similar analytical solutions to the equivalent barotropic quasi-geostrophic problem were identified by Larichev & Reznik (Reference Larichev and Reznik1976). The long-term behaviour of these structures remains an active area of study, with much recent work devoted to understanding modon breakdown processes such as instabilities (Brion, Sipp & Jacquin Reference Brion, Sipp and Jacquin2014; Davies et al. Reference Davies, Sutyrin, Crowe and Berloff2023; Protas Reference Protas2024), energy loss through wave generation (Flierl & Haines Reference Flierl and Haines1994; Johnson & Crowe Reference Johnson and Crowe2021) and viscous decay (Flór, van Heijst & Delfos Reference Flór, van Heijst and Delfos1995; Nielsen & Rasmussen Reference Nielsen and Rasmussen1997).
Layered quasi-geostrophic (QG) models are commonly used in modelling the ocean and atmosphere. These models consist of two-dimensional layers coupled through a vertical pressure gradient and are generally valid on scales where the Earth's rotation is dominant. Determining modon solutions in such models allows for the study of coherent structures in a more realistic setting than simple two-dimensional dynamics. Additionally, QG solutions can be used to provide a balanced initial condition for the more realistic primitive equation models. Kizner et al. (Reference Kizner, Berson, Reznik and Sutyrin2003) thoroughly explore modon solutions in a two-layer model, however, calculating solutions in models with arbitrarily many layers would be difficult using the same analytical methods.
Here, we extend a method we developed for finding surface QG dipoles (Crowe & Johnson Reference Crowe and Johnson2023; Johnson & Crowe Reference Johnson and Crowe2023) to the layered QG system by analytically reducing the problem to a linear eigenvalue problem which can be solved numerically. This approach, and the techniques used, have wide-ranging applications across fluid mechanics and applied mathematics. We begin by presenting the model and solution method in §§ 2 and 3, before applying it to some examples in § 4. Our method is found to be consistent with previous one-layer and two-layer studies and effective at finding new solutions in a three-layer model. Finally, we discuss our method and results in § 5.
2. Problem set-up
Consider a dipolar vortex (or modon) of radius 
$a$ moving in the zonal (
$x$) direction with speed 
$U$. In the frame of the moving vortex, the 
$N$-layer QG equations are
\begin{equation} (\partial_t - U \partial_x) q_i + J(\psi_i,q_i) + \beta_i \partial_x \psi_i = 0,\quad i \in \{1,2,\ldots, N\}, \end{equation}
for streamfunction 
$\psi _i$ and potential vorticity anomaly 
$q_i$ in each layer. Here
\begin{align*}J(f,g) = \partial_x f \partial_y g - \partial_y f \partial_x g\end{align*}
denotes the Jacobian operator, 
$\beta _i$ denotes the background vorticity gradient in the meridional (
$y$) direction which may vary by layer to incorporate the effects of bottom topography. The potential vorticities in each layer are given by
$$\begin{gather} q_1 = \nabla^2 \psi_1 + R_1^{{-}2} (\psi_2 - \psi_1), \end{gather}$$
$$\begin{gather}q_i = \nabla^2 \psi_i + R_i^{{-}2} (\psi_{i-1}-2\psi_i+\psi_{i+1}),\quad i \in \{2,\ldots N-1\}, \end{gather}$$
$$\begin{gather}q_N = \nabla^2 \psi_N + R_N^{{-}2} (\psi_{N-1}-\psi_N), \end{gather}$$
where 
$R_i = \sqrt {g'H_i}/f$ is the Rossby radius of deformation in each layer. Here, 
$g'$ denotes the buoyancy difference between consecutive layers and is taken to be the same for all layers for simplicity. Layer depth is denoted by 
$H_i$ and 
$f$ denotes the Coriolis parameter. The effects of a background flow can be easily incorporated into this model but will not be considered here.
Equation (2.1) can be written as
\begin{equation} \partial_t q_i + J(\psi_i+Uy,q_i+\beta_i y) = 0, \quad i \in \{1,2,\ldots, N\}, \end{equation}
and we now seek steady, nonlinear modon solutions by letting 
$\partial _t = 0$. Equation (2.5) may be written as
\begin{equation} q_i+\beta_i y = F_i(\psi_i+Uy), \quad i \in \{1,2,\ldots, N\} \end{equation}
where 
$F_i$ is an arbitrary, piece-wise differentiable function of a single variable. To proceed, we take 
$F_i$ to be a piece-wise linear function
\begin{equation} F_i(z) = \begin{cases} (\mathcal{K}_i/a^2)z, & x^2+y^2 < a^2,\\ (\beta_i/U) z, & x^2 + y^2 > a^2, \end{cases} \end{equation}
where 
$a$ denotes the radius of the vortex. The form of the 
$F_i$ outside the vortex is set by the requirement that vorticity and streamfunction perturbations decay towards infinity. Inside the vortex, we may either choose 
$\mathcal {K}_i$ to match the exterior solution, i.e. 
$\mathcal {K}_i = \beta _i a^2/U$, or take 
$\mathcal {K}_i = -K_i^2$, where 
$K_i$ is an eigenvalue which will be determined by continuity of 
$\psi _i$ across the vortex boundary. These cases will be respectively referred to as ‘passive’ and ‘active’ vortex regions. An active region consists of an isolated region of vorticity whereas the vorticity in a passive region is determined only by the advection of the background vorticity, 
$\beta _i y$, by the vortex flow field. Considering a single thin active upper layer above an arbitrary number of passive layers gives the 
$\beta$-plane extension of the surface quasi-geostrophic model of Johnson & Crowe (Reference Johnson and Crowe2023) with the 
$f$-plane model recovered by setting 
$\beta = 0$.
3. Solution method
Using the form of 
$F_i$ from (2.7), equation (2.6) can be written as
$$\begin{gather} (\nabla^2-R_1^{{-}2})\psi_1 + R_1^{{-}2}\psi_2 + \beta_1 y = (\mathcal{K}_1/a^2) (\psi_1 + Uy), \end{gather}$$
$$\begin{gather}(\nabla^2 - 2R_i^{{-}2})\psi_i + R_i^{{-}2}(\psi_{i-1}+\psi_{i+1}) +\beta_i y = (\mathcal{K}_i/a^2) (\psi_i + Uy), \end{gather}$$
$$\begin{gather}(\nabla^2 - R_N^{{-}2})\psi_N + R_N^{{-}2}\psi_{N-1} + \beta_N y = (\mathcal{K}_N/a^2) (\psi_N + Uy), \end{gather}$$
inside the vortex (
$x^2+y^2< a^2$) and
$$\begin{gather} (\nabla^2-R_1^{{-}2})\psi_1 + R_1^{{-}2}\psi_2 = (\beta_1/U) \psi_1, \end{gather}$$
$$\begin{gather}(\nabla^2 - 2R_i^{{-}2})\psi_i + R_i^{{-}2}(\psi_{i-1}+\psi_{i+1})= (\beta_i/U) \psi_i, \end{gather}$$
$$\begin{gather}(\nabla^2 - R_N^{{-}2})\psi_N + R_N^{{-}2}\psi_{N-1} = (\beta_N/U) \psi_N, \end{gather}$$
outside the vortex (
$x^2+y^2>a^2$) for 
$i \in \{2,\ldots N-1\}$.
3.1. Solution using Hankel transforms
To proceed, we move to plane polar coordinates, 
$(x,y) = (r\cos \theta,r\sin \theta )$, and represent 
$\psi _i$ using the Hankel transform
\begin{equation} \psi_i = Ua\sin\theta \int_0^\infty \hat{\psi}_i(\xi) \operatorname{J}_1(s\xi)\xi \,\textrm{d} \xi, \end{equation}
where 
$s = r/a$ and 
$\operatorname {J}_1$ denotes the Bessel function of the first kind of order 
$1$. The angular dependence in (3.7) is chosen to match that of a typical dipolar vortex. Substituting (3.7) into (3.1) to (3.6) gives
$$\begin{gather} \int_0^\infty [(\xi^2+\lambda_1^2 +\mathcal{K}_1)\hat{\psi}_1 - \lambda_1^2 \hat{\psi}_2] \operatorname{J}_1(s\xi) \xi \,\textrm{d} \xi = (\mu_1-\mathcal{K}_1) s, \end{gather}$$
$$\begin{gather}\int_0^\infty [(\xi^2+2\lambda_i^2 +\mathcal{K}_i)\hat{\psi}_i - \lambda_i^2 (\hat{\psi}_{i+1} + \hat{\psi}_{i-1})] \operatorname{J}_1(s\xi)\xi \,\textrm{d} \xi = (\mu_i-\mathcal{K}_i) s, \end{gather}$$
$$\begin{gather}\int_0^\infty [(\xi^2+\lambda_N^2 +\mathcal{K}_N)\hat{\psi}_N - \lambda_N^2 \hat{\psi}_{N-1}] \operatorname{J}_1(s\xi) \xi \,\textrm{d} \xi = (\mu_N-\mathcal{K}_N) s, \end{gather}$$
inside the vortex (
$s<1$) and
$$\begin{gather} \int_0^\infty [(\xi^2+\lambda_1^2 +\mu_1)\hat{\psi}_1 - \lambda_1^2 \hat{\psi}_2] \operatorname{J}_1(s\xi) \xi \,\textrm{d} \xi = 0, \end{gather}$$
$$\begin{gather}\int_0^\infty [(\xi^2+2\lambda_i^2 +\mu_i)\hat{\psi}_i - \lambda_i^2 (\hat{\psi}_{i+1} + \hat{\psi}_{i-1})] \operatorname{J}_1(s\xi) \xi \,\textrm{d} \xi = 0, \end{gather}$$
$$\begin{gather}\int_0^\infty [(\xi^2+\lambda_N^2 +\mu_N)\hat{\psi}_N - \lambda_N^2 \hat{\psi}_{N-1}] \operatorname{J}_1(s\xi) \xi \,\textrm{d} \xi = 0, \end{gather}$$
outside the vortex (
$s>1$) for 
$i \in \{2,\ldots N-1\}$. Here, 
$\mu _i = \beta _i a^2/U$ and 
$\lambda _i = a/R_i$, respectively, denote the non-dimensional background vorticity gradient and the ratio of the vortex radius to the Rossby radius in each layer.
We now define the matrix function
\begin{equation} \boldsymbol{\mathsf{K}}(\xi) = \begin{pmatrix} \xi^2+\lambda_1^2 & -\lambda_1^2 & 0 & \cdots & 0 & 0\\ -\lambda_2^2 & \xi^2+2\lambda_2^2 & -\lambda_2^2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -\lambda_N^2 & \xi^2+\lambda_N^2 \end{pmatrix}, \end{equation}and diagonal matrix
\begin{equation} \boldsymbol{\mathsf{D}}(\boldsymbol{c}) = \begin{pmatrix} c_1 & 0 & \cdots & 0 \\ 0 & c_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & c_N \end{pmatrix}, \end{equation}
for some vector 
$\boldsymbol {c} = (c_1, c_2, \ldots c_N)^{\rm T}$. We also introduce the vectors 
$\boldsymbol {\mu } = (\mu _1, \mu _2, \ldots, \mu _N)^{\rm T}$, 
$\boldsymbol {\mathcal {K}} = (\mathcal {K}_1, \mathcal {K}_2, \ldots, \mathcal {K}_N)^{\rm T}$ and 
$\hat {\boldsymbol {\psi }} = (\hat {\psi }_1, \hat {\psi }_2, \ldots, \hat {\psi }_N)^{\rm T}$ so (3.8) to (3.13) can be written concisely as
$$\begin{gather} \int_0^\infty [\boldsymbol{\mathsf{K}}(\xi)+\boldsymbol{\mathsf{D}}(\boldsymbol{\mathcal{K}})] \hat{\boldsymbol{\psi}}(\xi) \operatorname{J}_1(s\xi)\xi \,\textrm{d} \xi = (\boldsymbol{\mu} - \boldsymbol{\mathcal{K}})s,\quad s < 1, \end{gather}$$
$$\begin{gather}\int_0^\infty [\boldsymbol{\mathsf{K}}(\xi)+\boldsymbol{\mathsf{D}} (\boldsymbol{\mu})]\hat{\boldsymbol{\psi}}(\xi) \operatorname{J}_1(s\xi) \xi \,\textrm{d} \xi = 0,\quad s > 1. \end{gather}$$To proceed, we follow the approach of Johnson & Crowe (Reference Johnson and Crowe2023) and let
\begin{equation} \boldsymbol{A}(\xi) = [\boldsymbol{\mathsf{K}}(\xi)+ \boldsymbol{\mathsf{D}}(\boldsymbol{\mu})]\hat{\boldsymbol{\psi}}(\xi)\xi, \end{equation}
where 
$\boldsymbol {A}$ is expanded in terms of Bessel functions as
\begin{equation} \boldsymbol{A}(\xi) = \sum_{j = 0}^\infty \boldsymbol{a}_j \operatorname{J}_{2j+2}(\xi), \end{equation}
where the 
$\boldsymbol {a}_j$ are vectors of expansion coefficients. This expansion allows us to exploit the integral relation
\begin{equation} \int_0^\infty \operatorname{J}_{2j+2}(\xi) \operatorname{J}_1(s\xi) \,\textrm{d} \xi = \begin{cases} \operatorname{R}_j(s) & \textrm{for}\ s < 1,\\ 0 & \textrm{for}\ s > 1, \end{cases} \end{equation}and its inverse result
\begin{equation} \operatorname{J}_{2k+2}(\xi)/\xi = \int_0^1 s \operatorname{R}_k(s) \operatorname{J}_1(s\xi) \,\textrm{d} s, \end{equation}where
\begin{equation} \operatorname{R}_n(s) = ({-}1)^n s \operatorname{P}^{(0,1)}_n(2s^2-1), \end{equation}
denotes the Zernike radial function, a set of polynomials of degree 
$2n+1$ (Born & Wolf Reference Born and Wolf2019) which are orthogonal over 
$s\in [0,1]$ with weight 
$s$, and 
$\operatorname {P}_n^{(\alpha,\beta )}$ denotes the Jacobi polynomial. From (3.20) for 
$s>1$, we observe that (3.17) is automatically satisfied by our expansion of 
$\boldsymbol {A}(\xi )$ in (3.19) for all choices of 
$\boldsymbol {a}_j$.
Equation (3.16) may now be written as
\begin{equation} \sum_{j = 0}^\infty\left[\int_0^\infty [\boldsymbol{\mathsf{K}}(\xi)+ \boldsymbol{\mathsf{D}}(\boldsymbol{\mathcal{K}})] [\boldsymbol{\mathsf{K}}(\xi)+\boldsymbol{\mathsf{D}}(\boldsymbol{\mu})]^{{-}1} \operatorname{J}_{2j+2}(\xi) \operatorname{J}_1(s\xi) \,\textrm{d} \xi\right]\boldsymbol{a}_j = (\boldsymbol{\mu} - \boldsymbol{\mathcal{K}})s, \end{equation}
for 
$s < 1$. We now multiply (3.23) by 
$s\operatorname {R}_k(s)$ and integrate over 
$s\in [0,1]$ to get
\begin{equation} \sum_{j = 0}^\infty \left[\int_0^\infty [\boldsymbol{\mathsf{K}}(\xi)+\boldsymbol{\mathsf{D}}(\boldsymbol{\mathcal{K}})] [\boldsymbol{\mathsf{K}}(\xi)+\boldsymbol{\mathsf{D}}(\boldsymbol{\mu})]^{{-}1} \xi^{{-}1} \operatorname{J}_{2j+2}(\xi) \operatorname{J}_{2k+2}(\xi) \,\textrm{d} \xi \right]\boldsymbol{a}_j = \frac{\delta_{k0}}{4}(\boldsymbol{\mu} - \boldsymbol{\mathcal{K}}), \end{equation}
where 
$k \in \{0,1,2,\ldots \}$. Defining
$$\begin{gather} \boldsymbol{\mathsf{A}}_{kj} = \int_0^\infty \boldsymbol{\mathsf{K}}(\xi) [\boldsymbol{\mathsf{K}}(\xi) +\boldsymbol{\mathsf{D}}(\boldsymbol{\mu})]^{{-}1} \xi^{{-}1} \operatorname{J}_{2j+2}(\xi) \operatorname{J}_{2k+2}(\xi) \,\textrm{d} \xi, \end{gather}$$
$$\begin{gather}\boldsymbol{\mathsf{B}}_{kj} = \int_0^\infty [\boldsymbol{\mathsf{K}}(\xi) + \boldsymbol{\mathsf{D}}(\boldsymbol{\mu})]^{{-}1} \xi^{{-}1} \operatorname{J}_{2j+2}(\xi) \operatorname{J}_{2k+2}(\xi) \,\textrm{d} \xi, \end{gather}$$
$$\begin{gather}c_k = \frac{\delta_{k0}}{4}, \end{gather}$$allows us to write (3.24) as
\begin{equation} \sum_{j=0}^\infty [\boldsymbol{\mathsf{A}}_{kj} + \boldsymbol{\mathsf{D}}(\boldsymbol{\mathcal{K}}) \boldsymbol{\mathsf{B}}_{kj}] \boldsymbol{a}_j = (\boldsymbol{\mu}-\boldsymbol{\mathcal{K}})c_k. \end{equation}
We may now truncate this sum after a finite number of terms, say 
$M$, so 
$j,k \in \{0,1,2,\ldots, M-1\}$. The resulting problem is an 
$NM \times NM$ inhomogeneous, eigenvalue problem with up to 
$N$ eigenvalues. Up to 
$N$ additional equations are needed to solve this system and may be obtained by the requirement that the vortex boundary is a streamline in the vortex frame, 
$\psi _i +Uy = 0$ on 
$s = 1$. This condition gives 
$N$ equations in the form of a single equation for the coefficient vectors 
$\boldsymbol {a}_j$
\begin{align} \int_0^\infty \boldsymbol{A}(\xi) \operatorname{J}_1(\xi)\,\xi \,\textrm{d} \xi = \sum_{j = 0}^\infty \boldsymbol{a}_j \int_0^\infty \operatorname{J}_{2j+2}(\xi) \operatorname{J}_1(\xi) \,\textrm{d} \xi = 0 \quad \implies \quad \sum_{j = 0}^\infty ({-}1)^j \boldsymbol{a}_j = 0, \end{align}
which may similarly be truncated after 
$M$ terms.
We note that modon solutions do not exist if there is a singularity of 
$[\boldsymbol{\mathsf{K}}(\xi ) +\boldsymbol{\mathsf{D}}(\boldsymbol {\mu })]^{-1}$ for 
$\xi \in (0,\infty )$ as the integrals 
$\boldsymbol{\mathsf{A}}_{kj}$ and 
$\boldsymbol{\mathsf{B}}_{jk}$ do not converge. These singularities corresponds to the excitation of linear waves. In the two-layer case, the condition for the existence of singularities is equivalent to the conditions discussed in Kizner et al. (Reference Kizner, Berson, Reznik and Sutyrin2003).
3.2. The truncated eigenvalue problem
Our truncated eigenvalue problem from (3.28) and (3.29) may be written as
\begin{equation} \left[\boldsymbol{\mathsf{A}} + \sum_{i = 1}^N \mathcal{K}_i \boldsymbol{\mathsf{B}}_i\right]\boldsymbol{a} = \sum_{i=1}^N (\mu_i-\mathcal{K}_i)\boldsymbol{c}_i, \quad \textrm{s.t.}\quad \boldsymbol{d}_i^{\rm T}\boldsymbol{a} = 0 \quad \textrm{for}\ i\in\{1,\ldots,N\}, \end{equation}for
\begin{align} \boldsymbol{\mathsf{A}} = \begin{pmatrix} \boldsymbol{\mathsf{A}}_{0,0} & \cdots & \boldsymbol{\mathsf{A}}_{0,M-1} \\ \vdots & \ddots & \vdots \\ \boldsymbol{\mathsf{A}}_{M-1,0} & \cdots & \cdots \boldsymbol{\mathsf{A}}_{M-1,M-1} \end{pmatrix},\quad \boldsymbol{\mathsf{B}} = \begin{pmatrix} \boldsymbol{\mathsf{B}}_{0,0} & \cdots & \boldsymbol{\mathsf{B}}_{0,M-1} \\ \vdots & \ddots & \vdots \\ \boldsymbol{\mathsf{B}}_{M-1,0} & \cdots & \cdots \boldsymbol{\mathsf{B}}_{M-1,M-1} \end{pmatrix}, \end{align}
where 
$\boldsymbol {a} = (\boldsymbol {a}_0^{\rm T},\boldsymbol {a}_1^{\rm T},\ldots,\boldsymbol {a}_{M-1}^{\rm T})^{\rm T}$, 
$\boldsymbol {c} = (\boldsymbol {1},\boldsymbol {0},\ldots,\boldsymbol {0})^{\rm T}$ and 
$\boldsymbol {d} = (\boldsymbol {1},-\boldsymbol {1},\ldots,(-1)^{M-1} \boldsymbol {1})^{\rm T}$ where 
$\boldsymbol {1}$ and 
$\boldsymbol {0}$ denote length 
$N$ row vectors of ones and zeros, respectively. The indexed quantities 
$\boldsymbol{\mathsf{B}}_i$ and 
$\boldsymbol {c}_i$ are given by 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol {c}$, respectively, where all rows with a row index, 
$n \neq i \pmod {N}$, are set to zero.
Equation (3.30) is a multi-parameter eigenvalue problem (Atkinson Reference Atkinson1972) and may be solved using a range of methods (see Appendix A). For active layers, we must solve for the eigenvalue 
$\mathcal {K}_i = -K_i^2$ while for passive layers we let 
$\mathcal {K}_i = \mu _i$ so the term 
$\mathcal {K}_i\boldsymbol{\mathsf{B}}_i$ may be merged into 
$\boldsymbol{\mathsf{A}}$ and the right-hand side term 
$(\mu _i-\mathcal {K}_i)\boldsymbol {c}_i$ vanishes. Additional conditions of the form 
$\boldsymbol {d}_i^{\rm T}\boldsymbol {a} = 0$ are required for active layers but not for passive layers where continuity is automatically enforced by the fact that the components of 
$\boldsymbol {a}$ corresponding to a passive layer are zero. Mathematically, this is equivalent to the requirement to include an additional equation each time a new (unknown) eigenvalue is introduced. There will be a discrete spectrum of eigenvalues for each layer, corresponding to different modes in the radial direction. However, modon studies typically focus on the first-order radial mode which has the smallest value of 
$K_i$ and is generally the most stable. For multi-layer models, modon solutions may have different radial mode numbers in each layer.
The error in 
$K$ arising from the finite truncation decreases exponentially with increasing 
$M$. The supplementary material ‘Scaling_and_accuracy.pdf’ available at https://doi.org/10.1017/jfm.2024.619 discusses the accuracy and computational complexity of our Matlab implementation.
3.3. Determining streamfunctions and potential vortices
Once the 
$K_i$ and 
$\boldsymbol {a}_j$ are calculated, we may determine the streamfunctions, 
$\psi _i$, by inverting (3.18) for 
$\hat {\boldsymbol {\psi }}_i(\xi )\xi$ and evaluating the Hankel transform in (3.7) directly. Alternatively, using (3.18), (3.19) and (3.20) we can show that
\begin{equation} \boldsymbol{\varDelta}_N(\boldsymbol{\beta})\boldsymbol{\psi} ={-}\frac{U}{a} \sin\theta \sum_{j=0}^{M-1} \boldsymbol{a}_j \operatorname{R}_j(r/a), \end{equation}
where 
$\boldsymbol {\psi } = (\psi _1, \psi _2,\ldots, \psi _N)^{\rm T}$ and
\begin{align}
\boldsymbol{\varDelta}_N(\boldsymbol{\beta}) =
\begin{pmatrix} \nabla^2 -
\dfrac{1}{R_1^2}-\dfrac{\beta_1}{U} & \dfrac{1}{R_1^2} & 0
& \cdots & 0 & 0\\ \dfrac{1}{R_2^2} & \nabla^2 -
\dfrac{2}{R_2^2}-\dfrac{\beta_2}{U} & \dfrac{1}{R_2^2} &
\cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots
& \vdots \\ 0 & 0 & 0 & \cdots & \dfrac{1}{R_N^2} &
\nabla^2 - \dfrac{1}{R_N^2} - \dfrac{\beta_N}{U}
\end{pmatrix}, \end{align}
where 
$\boldsymbol {\beta } = (\beta _1, \beta _2,\ldots,\beta _N)^{\rm T}$. Equation (3.32) may be easily inverted using discrete Fourier transforms so this method is often easier numerically than inverting the Hankel transforms. From (2.2) to (2.4), potential vorticity anomalies may be similarly calculated using 
$\boldsymbol {q} = \varDelta _N(\boldsymbol {0}) \boldsymbol {\psi }$ where 
$\boldsymbol {q} = (q_1,q_2,\ldots,q_N)^{\rm T}$.
4. Examples
We now present three examples to demonstrate that this method gives consistent results with previous studies and can be used to find new modon solutions in layered QG models. A Matlab script to solve for the general N-layer problem is included as supplementary material and scales well for models with up to 
$N\sim O(100)$ layers.
4.1. The Larichev and Reznik dipole
Consider a one-layer model and choose parameters 
$(U,a,R_1,\beta _1) = 1$ and 
$\boldsymbol {\mathcal {K}}_1 = -K_1^2$. Modon solutions may be found using the method of Larichev & Reznik (Reference Larichev and Reznik1976), which combines a Bessel function solution with root finding to determine 
$K_1$. Using our approach, we construct a linear problem of the form
\begin{equation} [\boldsymbol{\mathsf{A}} - K_1^2\boldsymbol{\mathsf{B}}]\boldsymbol{a} = (\mu_1+K_1^2)\boldsymbol{c},\quad \textrm{s.t.} \quad \boldsymbol{d}^{\rm T}\boldsymbol{a} = 0, \end{equation}
using (3.28) and (3.29). This system may be solved using the method outlined in Appendix A. We determine a value of 
$K = 4.10787\cdots$ which matches the value calculated via root finding to seven significant figures using only 
$M = 7$ terms. Plots of the streamfunction 
$\psi _1$ and potential vorticity anomaly 
$q_1$ are given in figure 1. Higher-order radial modes and solutions for different parameters, including the case of 
$R_1 = \infty$, may be easily calculated using the same method.
Figure 1. A one-layer modon solution with 
$(U,a,R_1,\beta _1) = (1,1,1,1)$. We show (a) 
$\psi _1$ and (b) 
$q_1$. This solution corresponds to an eigenvalue of 
$K_1 \approx 4.108$.
4.2. Beta-plane baroclinic topographic modons
Kizner et al. (Reference Kizner, Berson, Reznik and Sutyrin2003) provide an extensive discussion of baroclinic beta-plane topographic modons using a two-layer QG model. Here, we present a two-layer solution with 
$(U,a,R_1,R_2,\beta _1,\beta _2) = (1,1,1,1,0,1)$ and 
$(\mathcal {K}_1,\mathcal {K}_2) = -(K_1^2,K_2^2)$ corresponding to two active layers. Using our approach, this problem is straightforwardly reduced to solving the two-parameter eigenvalue problem
\begin{equation} [\boldsymbol{\mathsf{A}} - K_1^2 \boldsymbol{\mathsf{B}}_1 - K_2^2 \boldsymbol{\mathsf{B}}_2]\boldsymbol{a} = (\mu_1+K_1^2)\boldsymbol{c}_1 + (\mu_2+K_2^2)\boldsymbol{c}_2,\quad\textrm{s.t.}\quad \boldsymbol{d}_1^{\rm T}\boldsymbol{a} = \boldsymbol{d}_2^{\rm T}\boldsymbol{a} = 0, \end{equation}
which may be solved using the methods outlined in Appendix A. Our two-layer modon solutions are plotted in figure 2 and show qualitative agreement with the plots of Kizner et al. (Reference Kizner, Berson, Reznik and Sutyrin2003). Solutions with one passive layer – referred to as ‘modons with one interior domain’ by Kizner et al. (Reference Kizner, Berson, Reznik and Sutyrin2003) – may be obtained by solving the same problem using 
$K_2^2 = -\mu _2$ and neglecting the condition 
$\boldsymbol {d}_2^{\rm T}\boldsymbol {a} = 0$. This demonstrates a major advantage of our approach; different vortex regimes can be considered using the same problem set-up. By contrast, when using the Bessel function approach of Kizner et al. (Reference Kizner, Berson, Reznik and Sutyrin2003), care is needed to ensure that the correct Bessel function is used depending on the sign of 
$\mathcal {K}_i$. As such, switching between active and passive layers is less straightforward as it requires both a solution with a different functional form, and new matching conditions across the vortex boundary.
Figure 2. A two-layer modon solution with 
$(U,a,R_1,R_2,\beta _1,\beta _2) = (1,1,1,1,0,1)$. We show (a) 
$\psi _1$, (b) 
$q_1$, (c) 
$\psi _2$ and (d) 
$q_2$. This solution has eigenvalues 
$(K_1,K_2) \approx (3.800, 3.950)$.
4.3. A mid-depth topographic modon
To demonstrate the flexibility of our approach, we consider a mid-depth vortex in a three-layer model. We assume that the top and bottom layers are passive and consider an active middle layer. A topographic slope with gradient in the 
$y$ direction is modelled by taking 
$\beta _1 = 0, \beta _2 = 0$ in the top two layers and 
$\beta _3 = 1$ in the bottom layer. For simplicity we take 
$R_i = 1$ for 
$i \in \{1,2,3\}$, 
$a = 1$ and 
$U = 1$. This problem reduces to the eigenvalue problem
\begin{equation} [\boldsymbol{\mathsf{A}} +\mu_1 \boldsymbol{\mathsf{B}}_1 - K_2^2 \boldsymbol{\mathsf{B}}_2 + \mu_3\boldsymbol{\mathsf{B}}_3]\boldsymbol{a} = (\mu_2+K_2^2)\boldsymbol{c}_2, \quad\textrm{s.t.}\quad \boldsymbol{d}_2^{\rm T}\boldsymbol{a} = 0, \end{equation}
which can be straightforwardly solved using standard methods (see Appendix A). The solution consists of an isolated region of strong vorticity in the middle layer, with velocity fields in all three layers. Weak vorticity is generated in the bottom layer due to the advection of fluid in the up-slope (
$y$) direction while no vorticity anomaly exists in the surface layer. Figure 3 shows 
$\psi _i$ for 
$i \in \{1,2,3\}$ and 
$q_2$ for this solution, the fields 
$q_1 = 0$ and 
$q_3 = \psi _3$ are not shown.
Figure 3. A three-layer mid-depth vortex solution with 
$(U,a,R_1,R_2,R_3,\beta _1,\beta _2,\beta _3) = (1,1,1,1,1,0,0,1)$. We show (a) 
$\psi _1$, (b) 
$\psi _2$, (c) 
$q_2$ and (d) 
$\psi _3$. In this case, 
$q_1 = 0$ and 
$q_3 = \psi _3$. Layers 1 and 3 are passive and layer 2 has eigenvalue 
$K_2 \approx 4.1835$.
5. Discussion and conclusions
We have presented a semi-analytical method for finding modon, or dipolar vortex, solutions in a layered QG model. These solutions enable the study of layered modons in idealised QG models while also providing a balanced initial condition for the study of three-dimensional mesoscale dipoles in a primitive equation model. While modons are steady solutions to the QG system, they are unlikely to be steady over long times in a more realistic model. As such, the effects of dissipation, wave generation and model parametrisations on layered modons is an area for future research.
The layered QG system is, in principle, solvable using combinations of Bessel functions inside and outside the vortex region. However, the matching conditions required to determine the coefficients become complicated for multiple layers. Further, these matching conditions must be solved numerically with solutions depending on the vortex parameters in a complicated, nonlinear way. As such, any additional insight obtained from having a closed form solution is limited. By contrast, our procedural method can be straightforwardly used to convert a model with arbitrarily many layers into a multi-parameter eigenvalue problem which can be solved numerically.
Additional advantages of our method include the ability to solve for any combination of active and passive layers without having to recalculate any terms in the linear system or change the form of the solution. Although we have not examined this in detail, we expect conditions for the existence of steady modon solutions – corresponding to a lack of linear wave excitation – to follow from examining the singularities of the quantity 
$[\boldsymbol{\mathsf{K}}+\boldsymbol{\mathsf{D}}(\boldsymbol {\mu })]^{-1}$. Finally, cases which would have to be treated separately using the standard analytical approach, such as an infinite Rossby radius, require no special treatment here.
We suggest that methods combining Hankel transforms and Zernike radial functions may be relevant to a range of physical problems where important dynamics occurs in an isolated circular region. These methods are not restricted to variants of the Laplace and Helmholtz equations, as shown for the Dirichlet-to-Neumann operators of Crowe & Johnson (Reference Crowe and Johnson2023) and Johnson & Crowe (Reference Johnson and Crowe2023). Additional effects – such as variations in density difference between layers – could be easily included into our model. It may also be possible to include the effects of layer-dependent background flow or extend the method to a three-dimensional QG model.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2024.619.
Declaration of interests
The authors report no conflict of interest.
Code availability
Matlab scripts implementing this method for the N-layer case and demonstrating the three examples considered are included as supplementary material.
Appendix A. Notes on inhomogeneous eigenvalue problems
Consider the problem
\begin{equation} (\boldsymbol{\mathsf{A}} - \lambda \boldsymbol{\mathsf{B}}) \boldsymbol{a} = \boldsymbol{c}_0 + \lambda \boldsymbol{c}_1, \quad \textrm{s.t.} \quad \boldsymbol{d}^{\rm T}\boldsymbol{a} = 0. \end{equation}
Assuming that 
$\boldsymbol{\mathsf{A}}$ is non-singular, we can left multiply (A1) by 
$\boldsymbol {d}^{\rm T}\boldsymbol{\mathsf{A}}^{-1}$ and multiply the resulting scalar by 
$\boldsymbol {c}_0+\lambda \boldsymbol {c}_1$ to relate 
$\boldsymbol {a}$ and 
$\boldsymbol {c}_0+\lambda \boldsymbol {c}_1$. This relation allows us to eliminate the inhomogeneity from (A1) giving the quadratic homogeneous problem
\begin{align} & [(\boldsymbol{d}^{\rm T} \boldsymbol{\mathsf{A}}^{{-}1} \boldsymbol{c}_0) \boldsymbol{\mathsf{A}} + \lambda ((\boldsymbol{d}^{\rm T}\boldsymbol{\mathsf{A}}^{{-}1} \boldsymbol{c}_1)\boldsymbol{\mathsf{A}} - (\boldsymbol{d}^{\rm T}\boldsymbol{\mathsf{A}}^{{-}1}\boldsymbol{c}_0)\boldsymbol{\mathsf{B}} - (\boldsymbol{c}_0 \boldsymbol{d}^{\rm T})\boldsymbol{\mathsf{A}}^{{-}1}\boldsymbol{\mathsf{B}}) \nonumber\\ &\quad -\lambda^2 ((\boldsymbol{d}^{\rm T}\boldsymbol{\mathsf{A}}^{{-}1} \boldsymbol{c}_1)\boldsymbol{\mathsf{B}} + (\boldsymbol{c}_1\boldsymbol{d}^{\rm T})\boldsymbol{\mathsf{A}}^{{-}1}\boldsymbol{\mathsf{B}})] \boldsymbol{a} = 0. \end{align}
Equation (A2) may be put into canonical form and solved as a linear eigenvalue problem for 
$\lambda$. The corresponding eigenvector, 
$\boldsymbol {a}$, should be directly calculated from (A1).
Consider now the general problem
\begin{equation} \left[\boldsymbol{\mathsf{A}} - \sum_{i = 1}^N \lambda_i \boldsymbol{\mathsf{B}}_i\right] \boldsymbol{a} = \boldsymbol{c}_0 + \sum_{i=1}^N \lambda_i \boldsymbol{c}_i, \quad \textrm{s.t.}\quad \boldsymbol{d}_j^{\rm T}\boldsymbol{a} = 0 \quad \textrm{for}\ j\in\{1,\ldots,N\}, \end{equation}
where 
$\boldsymbol{\mathsf{A}}$ and the 
$\boldsymbol{\mathsf{B}}_i$ are 
$M\times M$ matrices and 
$M>N$. We can begin by extending the 
$\boldsymbol {d}_j$ to a basis of 
$\mathrm {I\!R}^M$ by introducing 
$M-N$ orthonormal vectors 
$\boldsymbol {e}_k$ for 
$k\in \{1,\ldots,M-N\}$ such that each 
$\boldsymbol {e}_k$ is orthogonal to the 
$\boldsymbol {d}_j$. Therefore, the vector 
$\boldsymbol {a}$ lies in the 
$M-N$-dimensional space spanned by the 
$\boldsymbol {e}_k$ and can be expanded in terms of the 
$\boldsymbol {e}_k$ with coefficients 
$a_k'$. We now define the vector valued function
\begin{equation} \boldsymbol{f}(\boldsymbol{x}) = \left[\boldsymbol{\mathsf{A}} - \sum_{i = 1}^N \lambda_i \boldsymbol{\mathsf{B}}_i\right] \left[\sum_{k = 1}^{M-N} a'_k \boldsymbol{e}_k\right] - \boldsymbol{c}_0 - \sum_{i=1}^N \lambda_i \boldsymbol{c}_i, \end{equation}
for 
$\boldsymbol {x} = [\lambda _1,\ldots,\lambda _N,a_1',\ldots,a_{M-N}']^{\rm T}$ so solving for our eigenvalues and eigenvectors is equivalent to finding roots of 
$\boldsymbol {f}(\boldsymbol {x}) = \boldsymbol {0}$. Gradients of 
$\boldsymbol {f}$ with respect to 
$\boldsymbol {x}$ can easily be determined enabling standard gradient-based root finding techniques. Alternatively, numerical techniques developed specifically for multi-parameter problems (Atkinson Reference Atkinson1972; Muhic̆ & Plestenjak Reference Muhic̆ and Plestenjak2010) may prove effective, however, these rely on converting the problem to a large single-parameter problem which scales poorly with 
$N$ and 
$M$. Scalable, linear algebra approaches to this problem remain an area of active research.
Open access
