2.1. Governing equations

In this paper we investigate the linear perturbation eigenmodes and eigenvalues of a swirling flow in an unbounded domain $\mathbb {R}^3$. We express the velocity and pressure eigenmodes in cylindrical coordinates $(r, \phi, z)$, as

(2.1a,b)\begin{equation} \boldsymbol{u}'= \tilde{\boldsymbol u}(r;m,\kappa) \, {\rm e}^{\mathrm{i}(m\phi + \kappa z) + \sigma t} , \quad p' = \tilde{p}(r;m,\kappa) \, {\rm e}^{\mathrm{i}(m\phi + \kappa z) + \sigma t} ,\end{equation}

where $m$ and $\kappa$ represent the azimuthal and axial wavenumbers of the eigenmode, respectively, and $\sigma$ denotes the complex growth (or decay) rate of the eigenmode. Here, $m \in \mathbb {Z}$ since the fields must be periodic in $\phi$ with a period of $2{\rm \pi}$, while $\kappa \in \mathbb {R}\setminus {0}$ since there are no restrictions on the axial wavelength $2{\rm \pi} /\kappa$. The real part of $\sigma$ represents the growth/decay rate, while the imaginary part represents its wave frequency in time. Although Khorrami, Malik & Ash (Reference Khorrami, Malik and Ash1989) formulated a more general problem, we employ a more specialised form of the steady, equilibrium, swirling flow $\bar {\boldsymbol {U}}$, i.e.

(2.2)\begin{equation} \bar{\boldsymbol{U}} (r) = \bar{U}_\phi (r) \hat{\boldsymbol{e}}_\phi + \bar{U}_z (r) \hat{\boldsymbol{e}}_z, \end{equation}

which is only $r$ dependent and has no radial velocity component $\bar {U}_r$. The unperturbed base flow profile we consider for a wake vortex model is Batchelor's similarity solution adapted by Lessen et al. (Reference Lessen, Singh and Paillet1974) with

(2.3a,b)\begin{equation} \frac{\bar{U}_\phi}{U_0} = \frac{1- {\rm e}^{-r^2 / R_0^2}}{r / R_0}, \quad \frac{\bar{U}_z}{U_0} = \frac{1}{q}\, {\rm e}^{-r^2/R_0^2}, \end{equation}

where $R_0$ and $U_0$ are the length and velocity scales defined in Lessen et al. (Reference Lessen, Singh and Paillet1974, p. 755), and $q \neq 0$ is a dimensionless swirl parameter. This flow is often called the $q$ vortex, which is steady, axisymmetric and analytically tractable as the far-field asymptotic solution under the viscous light-loading condition (see Saffman Reference Saffman1993, pp. 257–260). When the axial flow component vanishes, i.e. $1/q \rightarrow 0$, this flow is equivalent to the Lamb–Oseen vortex. A schematic of the geometry is shown in figure 1. The unperturbed vortex is oriented along the $z$ direction with a circulation over the entire plane $\varGamma \equiv 2{\rm \pi} R_0 U_0$. Here $R_0$ is referred to as the characteristic radius of the vortex. As for the vortex profile, we consider the azimuthal velocity $\bar {U}_\phi$, which is maximised at $r = 1.122 R_0$ (Lessen et al. Reference Lessen, Singh and Paillet1974).

Figure 1. Vortex with circulation $\varGamma$ of length scale $R_0$ and coordinate systems.

To establish governing equations, we assume the fluid has constant density $\rho$ and constant kinematic viscosity $\nu$. The total velocity $\boldsymbol {u} \equiv u_r \hat {\boldsymbol {e}}_r + u_\phi \hat {\boldsymbol {e}}_\phi + u_z \hat {\boldsymbol {e}}_z$ obeys

(2.4)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 , \end{gather}$$

(2.5)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} =- \left( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u} - \frac{1}{\rho} \boldsymbol{\nabla} p + \nu {\nabla}^2 \boldsymbol{u} =- \boldsymbol{\nabla} \varphi + \boldsymbol{u} \times \boldsymbol{\omega} + \nu {\nabla}^2 \boldsymbol{u}, \end{gather}$$

where the total pressure is $p$, the vorticity is $\boldsymbol {\omega } \equiv \boldsymbol {\nabla } \times \boldsymbol {u}$ and the total specific energy is $\varphi \equiv u^2 /2 + p/\rho$, where $u^2 \equiv \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {u}$. We non-dimensionalise the equations using $R_0$ as the unit of length and $R_0/U_0$ as the unit of time. After non-dimensionalising and linearising (2.4) and (2.5) about the unperturbed flow (indicated with overbars $\bar {*}$), we obtain the following equation for the perturbations (indicated with primes $*'$):

(2.6)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}' = 0, \end{gather}$$

(2.7)$$\begin{gather}\frac{\partial \boldsymbol{u}'}{\partial t} =- \boldsymbol{\nabla} \varphi' + \bar{\boldsymbol{U}}(r) \times \boldsymbol{\omega}' - \bar{\boldsymbol{\omega}}(r) \times \boldsymbol{u}' + \frac{1}{Re} {\nabla}^2 \boldsymbol{u}'. \end{gather}$$

Here the Reynolds number, denoted $Re$, is defined to be $U_0 R_0 / \nu$. Note that the non-dimensionalised $q$ vortex is

(2.8)\begin{equation} \bar{\boldsymbol{U}}(r) = \left(\frac{1- {\rm e}^{-r^2}}{r}\right) \hat{\boldsymbol{e}}_\phi + \left(\frac{1}{q}\, {\rm e}^{-r^2}\right) \hat{\boldsymbol{e}}_z .\end{equation}

The established governing equations are essentially the incompressible, linearised Navier–Stokes equations, which, in combination with the $q$ vortex, were also used in recent vortex stability analyses, such as Qiu et al. (Reference Qiu, Cheng, Xu, Xiang and Liu2021). By putting (2.1a,b) to (2.6) and (2.7), we obtain the equations that govern the perturbations

(2.9)$$\begin{gather} \boldsymbol{\nabla}_{m \kappa} \boldsymbol{\cdot} \tilde{\boldsymbol u} = 0, \end{gather}$$

(2.10)$$\begin{gather}\sigma \tilde{\boldsymbol u} =-\boldsymbol{\nabla}_{m \kappa} \tilde{\varphi} + \bar{\boldsymbol{U}} \times \tilde{\boldsymbol \omega} - \bar{\boldsymbol{\omega}} \times \tilde{\boldsymbol u} + \frac{1}{Re} {\nabla}_{m \kappa}^2 \tilde{\boldsymbol u}, \end{gather}$$

where $\sigma$ is a function of $m$ and $\kappa$ (i.e. it obeys the dispersion relationship), $\tilde {\boldsymbol \omega } \equiv \boldsymbol {\nabla }_{m \kappa } \times \tilde {\boldsymbol u}$ and $\tilde {\varphi } \equiv \bar {\boldsymbol {U}} \boldsymbol {\cdot } \tilde {\boldsymbol u} + \tilde {p}$. In the equations above, the subscript $(\ast )_{m\kappa }$ attached to the operators means that they act on modes of fixed azimuthal and axial wavenumbers $m$ and $\kappa$. Therefore, the differential operators $\partial / \partial \phi$ and $\partial / \partial z$ inside these operators are replaced with the simple multiplication operators $\mathrm {i} m$ and $\mathrm {i} \kappa$, respectively (see Appendix A).

2.2. Boundary and analyticity conditions

We require both the velocity and pressure fields to be analytic at $r = 0$ and to decay rapidly to 0 as $r \rightarrow \infty$. The conversion of these conditions to numerical boundary conditions can be found in previous works such as Batchelor & Gill (Reference Batchelor and Gill1962), Mayer & Powell (Reference Mayer and Powell1992) and Ash & Khorrami (Reference Ash and Khorrami1995). In this section we briefly discuss these conditions and how they will be treated in our method, where functions are treated as a truncated sum of the mapped Legendre functions.

The analyticity at the origin is equivalent to the pole condition that correctly removes the coordinate singularity (see Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988; Matsushima & Marcus Reference Matsushima and Marcus1995; Lopez, Marques & Shen Reference Lopez, Marques and Shen2002). The pole condition for a scalar function $f(r,\phi, z)$ to be analytic at $r=0$ is that it asymptotically behaves as a polynomial in $r$, with the degree dependent on the azimuthal wavenumber $m$ (see Matsushima & Marcus Reference Matsushima and Marcus1997, p. 323), that is,

(2.11)\begin{equation} f(r,\phi, z) = \sum_{m =-\infty}^{\infty} \, {\rm e}^{\mathrm{i}m\phi} r^{|m|} \left(\sum_{n = 0}^{\infty} {C_n(z;m) r^{2n}} \right)\text{ as }r\rightarrow 0,\end{equation}

for some set of functions, analytic in $z, C_n(z;m)$. Although the pole condition for velocity fields in polar or cylindrical coordinates is rather complicated because of the velocity coupling of $r$ and $\phi$ at the origin (see Matsushima & Marcus Reference Matsushima and Marcus1997, pp. 328–330), we use toroidal and poloidal streamfunctions, given in (2.12), instead of the primitive velocity components so that the analyticity can be determined by making these streamfunctions obey the requirements of scalars (see Appendix B).

On the other hand, the rapid decay condition as $r \rightarrow \infty$ is relevant to the feasibility of linear perturbations. Since a perturbation lasting even at radial infinity would require infinite kinetic energy, decay should be necessary to consider it physical (see our definition in § 1.3). The simplest description is $\tilde {\boldsymbol u}, \tilde {p} \rightarrow 0$ as $r \rightarrow \infty$ (Batchelor & Gill Reference Batchelor and Gill1962). Several numerical methods that require the domain truncation at large $r$ mimic this condition by imposing the homogeneous Dirichlet boundary condition for $\tilde {\boldsymbol u}$ and $\tilde {p}$ at the outer boundary of the radially truncated domain $r=r_\infty$. In other words, $\tilde {u}_r=\tilde {u}_\phi = \tilde {u}_z = \tilde {p} =0$ at $r=r_\infty$ (see Khorrami et al. Reference Khorrami, Malik and Ash1989; Khorrami Reference Khorrami1991). However, this approach involves two problems. First, it cannot preclude non-physical eigenmode solutions that do not decay properly but incidentally end up being zero at $r=r_\infty$ (i.e. wall-bounded modes). Such non-physical solutions may also appear with non-zero viscosity, triggering non-normalisable singularities at radial infinity, where more information can be found in Fabre et al. (Reference Fabre, Sipp and Jacquin2006, p. 268) or Mao & Sherwin (Reference Mao and Sherwin2011, pp. 17–21). Second, it does not explicitly take into account how rapidly the perturbation decays. Considering the velocity field, it must decay faster than algebraic decay rates of $O(r^{-1})$ for kinetic energy to be finite as $r \rightarrow \infty$ (cf. Bölle et al. Reference Bölle, Brion, Robinet, Sipp and Jacquin2021). Mathematically, the restriction is even more strict, requiring exponential or super-exponential decay rates (Ash & Khorrami Reference Ash and Khorrami1995). Our method is free from domain truncation and explicitly forces solutions to decay harmonically, i.e. $O ( r^{-|m|} \, {\rm e}^{\mathrm {i}m\phi } )$ as $r \rightarrow \infty$, due to the decaying nature of the basis functions.

By utilising the current method, it can be ensured that any scalar functions, represented by the sum of mapped Legendre functions that serve as Galerkin basis functions, comply with the aforementioned conditions. This is precisely how each basis function behaves as the radial distance approaches zero and approaches infinity. Therefore, an advantage of employing the mapped associated Legendre functions is that there is no need for additional treatment for numerical boundary conditions. For further information regarding the properties of the mapped Legendre functions, please refer to § 3.

2.3. Poloidal–toroidal decomposition

The governing equations (2.9) and (2.10), along with the correct boundary conditions and given values of $m$ and $\kappa$, formally constitute a set of four equations that make up a generalized eigenvalue problem in terms of $\tilde {p}$ (or $\tilde {\varphi }$) and the three components of $\tilde {\boldsymbol u} \equiv \tilde {u}_r \hat {\boldsymbol {e}}_r + \tilde {u}_\phi \hat {\boldsymbol {e}}_\phi + \tilde {u}_z \hat {\boldsymbol {e}}_z$, which are often referred to as primitive variables, with $\sigma$ as the eigenvalue. The formal expression of the eigenvalue problem can be found in Bölle et al. (Reference Bölle, Brion, Robinet, Sipp and Jacquin2021, p. 7). Some previous studies have taken additional steps to eliminate $\tilde {p}$ from the momentum equations or even reduce the problem in terms of, for example, only $\tilde {u}_\phi$ and $\tilde {u}_z$, resulting in the generalized eigenvalue problem form ${\boldsymbol{\mathsf{A}}}\boldsymbol {x} = \lambda {\boldsymbol{\mathsf{B}}} \boldsymbol {x}$ (Mayer & Powell Reference Mayer and Powell1992; Heaton & Peake Reference Heaton and Peake2007; Mao & Sherwin Reference Mao and Sherwin2011). However, such variable reduction inevitably increases the spatial order of the system and, consequently, requires a higher resolution for computation, which undermines the advantage of having a smaller number of state variables (Mayer & Powell Reference Mayer and Powell1992). To avoid this issue, we use a poloidal–toroidal decomposition of the velocity field to formulate the matrix eigenvalue problem while preserving the spatial order of the governing equations. Moreover, the use of the poloidal and toroidal streamfunctions is advantageous because the formulation results in the standard eigenvalue problem of the form ${\boldsymbol{\mathsf{A}}}\boldsymbol {x} = \lambda \boldsymbol {x}$.

To begin with, we apply the poloidal–toroidal decomposition to the governing equations of wake vortices that are linearised about the $q$ vortex. The basic formulation was performed by Matsushima & Marcus (Reference Matsushima and Marcus1997, p. 339), and we provide more details of its mathematical foundation in this section. Although the poloidal–toroidal decomposition of solenoidal vector fields is mainly discussed in spherical geometry (Chandrasekhar Reference Chandrasekhar1981, pp. 622–626), it can be employed in the cylindrical coordinate system while preserving some essential properties of the decomposition (Ivers Reference Ivers1989). When we select the unit vector in the $z$ direction $\hat {\boldsymbol {e}_z}$ as a reference vector, a solenoidal vector field $\boldsymbol {V} (r, \phi, z ) = V_r (r, \phi, z) \hat {\boldsymbol {e}_r} + V_\phi (r, \phi, z) \hat {\boldsymbol {e}_\phi } + V_z (r, \phi, z) \hat {\boldsymbol {e}_z}$ can be expressed as

(2.12)\begin{equation} \boldsymbol{V} = \boldsymbol{\nabla} \times \{ \psi (r, \phi, z) \hat{\boldsymbol{e}_z} \} + \boldsymbol{\nabla} \times [\boldsymbol{\nabla} \times \{ \chi (r, \phi, z) \hat{\boldsymbol{e}_z} \} ],\end{equation}

where $\psi$ and $\chi$ are the toroidal and poloidal streamfunctions of $\boldsymbol {V}$. Such a decomposition is feasible if $\boldsymbol {V}$ has zero spatial mean components in the radial and azimuthal directions over an infinite disk for all $z$ (cf. Jones Reference Jones2008). This zero-mean condition is satisfied in our study because our velocity fields are spatially periodic perturbations of the base flow. Ivers (Reference Ivers1989) concluded that the toroidal and poloidal fields are orthogonal over an infinite slab $a < z < b$ if $\psi$ and $\chi$ decay sufficiently rapidly as $r \rightarrow \infty$. The decay condition of $\psi$ and $\chi$ requires $\boldsymbol {V}$ to decay sufficiently rapidly to zero for large $r$.

In what follows, we find more rigorous statement for the decay condition of $\boldsymbol {V}$ as $r \rightarrow \infty$ where $\psi$ and $\chi$ are well defined. The $z$ component of (2.12) is

(2.13)\begin{equation} \frac{1}{r} \frac{\partial }{\partial r} \left( r \frac{\partial \chi}{ \partial r} \right) + \frac{1}{r^2} \frac{\partial^2 \chi}{\partial \phi^2} =- V_z . \end{equation}

Taking the curl of (2.12), we obtain

(2.14)\begin{equation} \boldsymbol{\nabla} \times \boldsymbol{V} = \boldsymbol{\nabla} \times \{(-\nabla^2 \chi) \hat{\boldsymbol{e}_z}\} + \boldsymbol{\nabla} \times \{\boldsymbol{\nabla} \times (\psi \hat{\boldsymbol{e}_z})\}, \end{equation}

with its $z$ component equal to

(2.15)\begin{equation} \frac{1}{r} \frac{\partial }{\partial r} \left( r \frac{\partial \psi}{ \partial r} \right) + \frac{1}{r^2} \frac{\partial^2 \psi}{\partial \phi^2} =- ( \boldsymbol{\nabla} \times \boldsymbol{V} )_z. \end{equation}

Solving (2.13) and (2.15), which are the two-dimensional Poisson equations, can yield the solution to $\psi$ and $\chi$. By ignoring the gauge freedom with respect to $z$, we can determine the solution using two-dimensional convolution as

(2.16)$$\begin{gather} \psi =- ( \boldsymbol{\nabla} \times \boldsymbol{V} )_z \ast \varPhi, \end{gather}$$

(2.17)$$\begin{gather}\chi =- V_z \ast \varPhi, \end{gather}$$

where $\varPhi$ is Green's function for the entire plane $\mathbb {R}^2$ equivalent to

(2.18)\begin{equation} \varPhi (r, \phi) = \frac{1}{2{\rm \pi}} \ln r. \end{equation}

In order for the convolutions in (2.16) and (2.17) to be meaningful everywhere, there exist $p_1 > 0, p_2 > 0$ and $p_3 >0$ such that

(2.19) \begin{equation} \left.\begin{array}{ll@{}} V_r \sim O ( r^{-1-p_1} ), &\\ V_\phi \sim O ( r^{-1-p_2} ), &\text{ as } r \rightarrow \infty,\\ V_z \sim O ( r^{-2-p_3} ), &\\ \end{array}\right\}\end{equation}

given that $\boldsymbol {V}$ decays algebraically. Otherwise, $\boldsymbol {V}$ may decay exponentially or super-exponentially. If $\boldsymbol {V}$ is referred to as a velocity field, it has finite total kinetic energy over the entire space since all components decay faster than $O(r^{-1})$ as $r \rightarrow \infty$. The finite kinetic energy condition is physically reasonable, especially when dealing with velocity fields representing small perturbations (cf. Bölle et al. Reference Bölle, Brion, Robinet, Sipp and Jacquin2021). On the other hand, Matsushima & Marcus (Reference Matsushima and Marcus1997) considered the case where $\psi$ and $\chi$ could be unbounded by including additional logarithmic terms in $\psi$ and $\chi$, providing a more comprehensive extension of the poloidal–toroidal decomposition to more general vector fields, including the mean axial components. However, in the present study, we choose $\boldsymbol {V}$ as a linear perturbation of no bulk movement and, therefore, the logarithmic terms do not need to be considered.

Suppressing the gauge freedom by adding restrictions that are independent of $z$ to $\psi$ and $\chi$, e.g.

(2.20)\begin{equation} \lim_{r\rightarrow\infty} \psi (r, \phi, z) = \lim_{r\rightarrow\infty} \chi (r, \phi, z)= 0,\end{equation}

we can define the following linear and invertible operator $\mathbb {P}: \mathcal {U} \rightarrow \mathcal {P}$ as

(2.21)\begin{equation} \mathbb{P}(\boldsymbol{V}) \equiv \begin{pmatrix} \psi (r, \phi, z) \\ \chi (r, \phi, z) \end{pmatrix}, \end{equation}

where $\mathcal {U}$ is the set of sufficiently rapidly decaying solenoidal vector fields from $\mathbb {R}^3$ to $\mathbb {R}^3$ ($\mathbb {C}^3$) that satisfy (2.19) and $\mathcal {P}$ is the set of functions from $\mathbb {R}^3$ to $\mathbb {R}^2$ ($\mathbb {C}^2$) that satisfy (2.20). Using Helmholtz's theorem, we may extensively define $\mathbb {P}$ on more generalized vector fields that are not solenoidal but their solenoidal portion can be decomposed toroidally and poloidally. If we expand the domain of $\mathbb {P}$, however, it should be kept in mind that the operator is no longer injective because for any $\boldsymbol {V} \in \mathcal {U}, \mathbb {P}(\boldsymbol {V}) = \mathbb {P}(\boldsymbol {V} + \boldsymbol {\nabla }v)$, where $v$ is an arbitrary scalar potential for a non-zero irrotational vector field. On the other hand, it is noted that $\mathbb {P} ( \nabla ^2 \boldsymbol {V} ) = \nabla ^2 \mathbb {P}(\boldsymbol {V}) \equiv ( \nabla ^2 \psi, \nabla ^2 \chi )$ for $\boldsymbol {V} \in \mathcal {U}$ because

(2.22)\begin{align} &\nabla^2 [ \boldsymbol{\nabla} \times \{ \psi (r, \phi, z) \hat{\boldsymbol{e}_z} \} + \boldsymbol{\nabla} \times [\boldsymbol{\nabla} \times \{ \chi (r, \phi, z) \hat{\boldsymbol{e}_z} \} ]] \nonumber\\ &\quad = \boldsymbol{\nabla} \times \{ \nabla^2 \psi (r, \phi, z) \hat{\boldsymbol{e}_z} \} + \boldsymbol{\nabla} \times [ \boldsymbol{\nabla} \times \{ \nabla^2 \chi (r, \phi, z) \hat{\boldsymbol{e}_z} \}]. \end{align}

Applying the operator $\mathbb {P}$ to both sides of (2.7), we obtain

(2.23)\begin{equation} \frac{\partial \mathbb{P}(\boldsymbol{u}')}{\partial t} = \mathbb{P} ( \bar{\boldsymbol{U}}(r) \times \boldsymbol{\omega}' ) - \mathbb{P} ( \bar{\boldsymbol{\omega}}(r) \times \boldsymbol{u}' ) + \frac{1}{Re} {\nabla}^2 \mathbb{P}(\boldsymbol{u}'),\end{equation}

because $\mathbb {P}(\boldsymbol {\nabla }\varphi ) = \mathbb {P}(\boldsymbol {0} + \boldsymbol {\nabla }\varphi ) = \mathbb {P}(\boldsymbol {0}) = \boldsymbol {0}$. Assuming $\boldsymbol {u}'$ to be solenoidal, $\boldsymbol {u}'$ automatically satisfies the continuity equation and can be determined from $\mathbb {P}(\boldsymbol {u}')$ by taking the inverse of it using (2.12).

Since we are interested in the perturbation velocity field as in (2.1a,b), we define two $r$-dependent scalar functions $\tilde {\psi }(r;m,\kappa )$ and $\tilde {\chi }(r;m,\kappa )$ such that

(2.24)\begin{equation} \mathbb{P} ( \tilde{\boldsymbol u}(r;m,\kappa) \, {\rm e}^{\mathrm{i}(m\phi + \kappa z)+\sigma t} ) = \begin{pmatrix} \tilde{\psi} (r;m,\kappa) \, {\rm e}^{\mathrm{i}(m\phi + \kappa z)+\sigma t} \\ \tilde{\chi} (r;m,\kappa) \, {\rm e}^{\mathrm{i}(m\phi + \kappa z)+\sigma t} \end{pmatrix} . \end{equation}

The fact that the poloidal and toroidal components in (2.24) preserve the exponential part can be verified by substituting the perturbation velocity field formula into $\boldsymbol {V}$ in (2.13) and (2.15). For convenience, we simplify the expression in (2.24) to

(2.25)\begin{equation} \mathbb{P}_{m \kappa}(\tilde{\boldsymbol u}(r;m,\kappa)) \equiv \begin{pmatrix} \tilde{\psi}(r;m,\kappa) \\ \tilde{\chi}(r;m,\kappa) \end{pmatrix} .\end{equation}

Finally, putting (2.1a,b) into (2.23) leads to the standard eigenvalue problem form in terms of $\mathbb {P}_{m \kappa }(\tilde {\boldsymbol u}(r;m,\kappa ))$,

(2.26)\begin{equation} \sigma [ \mathbb{P}_{m\kappa}(\tilde{\boldsymbol u}(r;m,\kappa)) ] = \mathcal{L}_{m\kappa}^{\nu} [ \mathbb{P}_{m\kappa}(\tilde{\boldsymbol u}(r;m,\kappa)) ], \end{equation}

where the linear operator $\mathcal {L}_{m\kappa }^{\nu }$ is defined as

(2.27)\begin{equation} \mathcal{L}_{m\kappa}^{\nu} [ \mathbb{P}_{m\kappa}(\tilde{\boldsymbol u}) ] \equiv \mathbb{P}_{m\kappa} ( \bar{\boldsymbol{U}}(r) \times \tilde{\boldsymbol \omega} ) - \mathbb{P}_{m\kappa} ( \bar{\boldsymbol{\omega}} (r) \times \tilde{\boldsymbol u} ) + \frac{1}{Re} {\nabla}^2_{m \kappa} \mathbb{P}_{m\kappa}(\tilde{\boldsymbol u}). \end{equation}

Excluding the viscous diffusion term, we additionally define the inviscid operator $\mathcal {L}_{m\kappa }^{0}$ as

(2.28)\begin{equation} \mathcal{L}_{m\kappa}^{0} [ \mathbb{P}_{m\kappa}(\tilde{\boldsymbol u}) ] \equiv \mathbb{P}_{m\kappa} ( \bar{\boldsymbol{U}}(r) \times \tilde{\boldsymbol \omega} ) - \mathbb{P}_{m\kappa} ( \bar{\boldsymbol{\omega}}(r) \times \tilde{\boldsymbol u} ) \end{equation}

for the inviscid linear analysis solving

(2.29)\begin{equation} \sigma [ \mathbb{P}_{m\kappa}(\tilde{\boldsymbol u}(r;m,\kappa)) ] = \mathcal{L}_{m\kappa}^{0} [ \mathbb{P}_{m\kappa}(\tilde{\boldsymbol u}(r;m,\kappa)) ]. \end{equation}

3.1. Mapped Legendre functions

Associated Legendre functions with algebraic mapping are used as basis functions to expand an arbitrary function over $0 \le r < \infty$, ultimately discretising the eigenvalue problems to be solved numerically. The expansion was first introduced by Matsushima & Marcus (Reference Matsushima and Marcus1997) and applied to three-dimensional vortex instability studies by Bristol et al. (Reference Bristol, Ortega, Marcus and Savaş2004) and Feys & Maslowe (Reference Feys and Maslowe2016). The algebraically mapped associated Legendre functions $P_{L_n}^m (r)$, or simply mapped Legendre functions, are equivalent to the mapping of the associate Legendre functions $P_n^m (\zeta )$ with order $m$ and degree $n$ defined on $-1 \le \zeta < 1$, where

(3.1)\begin{equation} \zeta \equiv \frac{r^2 - L^2}{r^2 + L^2} \Longleftrightarrow r = L \sqrt{\frac{1+\zeta}{1-\zeta}}. \end{equation}

An additional parameter $L>0$ is the map parameter, which can be arbitrarily set. However, when it is used for a spectral collocation method, a change in $L$ affects the spatial resolution of discretisation and the value should be carefully chosen to achieve fast convergence or eliminate spurious results. Matsushima & Marcus (Reference Matsushima and Marcus1997) showed that $P_{L_n}^m (r) \sim O ( r^{|m|} ) \text { as } r \rightarrow 0$ and $P_{L_n}^m (r) \sim O ( r^{-|m|} ) \text { as } r \rightarrow \infty$, which leads to the fact that any polar function $P_{L_n}^m (r) \, {\rm e}^{\mathrm {i} m \phi }$ behaves analytically at the origin (see Eisen, Heinrichs & Witsch Reference Eisen, Heinrichs and Witsch1991, pp. 243–244) and decays harmonically to zero at radial infinity. These asymptotic properties are suitable to apply the correct boundary conditions for the present problem.

Next, we prove that a set of some mapped Legendre functions can constitute a complete orthogonal basis of spectral space. Since the associate Legendre functions $P_n^m (\zeta )$ are the solutions to the associate Legendre equation

(3.2)\begin{equation} \frac{{\rm d}}{{\rm d} \zeta} \left[ \left( 1 - \zeta^2 \right) \frac{{\rm d} P_n^m}{{\rm d} \zeta} \right] + \left[ n(n+1) - \frac{m^2}{1-\zeta^2} \right]P_n^m (\zeta) = 0, \end{equation}

the mapped Legendre functions satisfy the following second-order differential equation:

(3.3)\begin{equation} \frac{{\rm d} }{{\rm d} r} \left[ r \frac{{\rm d} {P_{L_n}^m}}{{\rm d} r} \right] - \frac{m^2}{r}P_{L_n}^m (r) + \frac{4 n (n+1)L^2 r}{(L^2 +r^2)^2} P_{L_n}^{m} (r) = 0.\end{equation}

As (3.3) is the Sturm–Liouville equation with the weight function

(3.4)\begin{equation} w(r) \equiv \frac{4L^2 r}{(L^2 + r^2)^2}, \end{equation}

the mapped Legendre functions $P_{L_{|m|}}^m (r), P_{L_{|m|+1}}^m (r), P_{L_{|m|+2}}^m (r), \ldots$ form an orthogonal basis of the Hilbert space $L^2 ( \mathbb {R}^+, w(r)\, {\rm d} r )$. Thus, for two integers $n$ and $k$ larger than or equal to $|m|$,

(3.5)\begin{align} &\langle P_{L_n}^m, P_{L_k}^m \rangle = \int_{0}^{\infty} P_{L_n}^m(r) P_{L_k}^m (r) w(r) \, {\rm d} r \nonumber\\ &\quad = \int_{-1}^{1} P_{n}^m (\zeta) P_{k}^m (\zeta) \, {\rm d} \zeta = \frac{2(n+|m|)!}{(2n+1)(n-|m|)!} \delta_{nk}, \end{align}

where $\delta _{nk}$ denotes the Kronecker delta with respect to $n$ and $k$.

Considering a polar function $f_m (r) \, {\rm e}^{\mathrm {i} m \phi }$, where $f_m \in L^2 ( \mathbb {R}^+, w(r)\, {\rm d} r )$, it can be expanded by the mapped Legendre functions as

(3.6)\begin{equation} f_m (r) \, {\rm e}^{\mathrm{i} m \phi} = \sum_{n = |m|}^{\infty} f_n^m P_{L_n}^m (r) \, {\rm e}^{\mathrm{i} m \phi}, \end{equation}

and the coefficient $f_n^m$ can be calculated based on the orthogonality of the basis functions,

(3.7)\begin{align} f_n^m &= \frac{\langle \,f_m, P_{L_n}^m \rangle}{\langle P_{L_n}^m, P_{L_n}^m \rangle } = \frac{(2n+1) (n-|m|)!}{2(n +|m|)!} \int_{0}^{\infty} f_m (r) P_{L_n}^m (r) w(r)\, {\rm d} r \nonumber\\ &= \frac{(2n+1) (n-|m|)!}{2(n+|m|)!} \int_{-1}^{1} f_m \left( L \sqrt{\frac{1 + \zeta}{1 - \zeta}} \right) P_{n}^m (\zeta) \, {\rm d} \zeta. \end{align}

When we expand an analytic function on $0 \le r < \infty$ that vanishes at infinity, the expansion in (3.6) is especially suitable because they are able to serve as Galerkin basis functions. Even if we use the truncated series of (3.6), analyticity at the origin and vanishing behaviour at infinity remain valid.

3.2. Mapped Legendre spectral collocation method

In order to discretise the problem, we use a spectral collocation method using the mapped Legendre functions as basis functions. Given the azimuthal and axial wavenumbers $m$ and $\kappa$, we take a truncated basis set of first $M$ elements $\{ P_{L_{|m|}}^{m} , \ldots, P_{L_{|m|+M-1}}^{m} \}$ and expand $f_m(r) \, {\rm e}^{\mathrm {i}(m\phi + \kappa z)}$ as

(3.8)\begin{equation} f_m (r) \, {\rm e}^{\mathrm{i}(m \phi + \kappa z)} = \sum_{n=|m|}^{|m|+M-1} f_n^m P_{L_n}^{m} (r) \, {\rm e}^{\mathrm{i} (m\phi + \kappa z)},\end{equation}

so that the function is represented by $M$ discretised coefficients $(\,f_{|m|}^{m}, \ldots, f_{|m|+M-1}^{m} )$. The coefficients are numerically obtained by applying the Gauss–Legendre quadrature rule to (3.7). Let $\zeta _j$ and $\varpi _j$ be the $j$th root of the Legendre polynomial $P_N$ of degree $N$ in $(-1,1)$ with its quadrature weight defined as

(3.9)\begin{equation} \varpi_j = 2 (1-\zeta_j^2)^{-1} \left[ \left.\dfrac{{\rm d} P_N}{{\rm d} \zeta} \right|_{\zeta = \zeta_j} \right]^{-2}, \quad j = 1, \ldots, N, \end{equation}

and with radial collocation points $r_j$ determined from (3.1) as

(3.10)\begin{equation} r_j \equiv L\sqrt{\frac{1 + \zeta_j}{1-\zeta_j}}, \quad j = 1, \ldots, N, \end{equation}

which means that half of the collocation points are distributed in the inner high-resolution region $0 \le r < L$ whereas the other half are posed in the outer low-resolution region $r \ge L$ (Matsushima & Marcus Reference Matsushima and Marcus1997). In order to describe spatial resolution, we define the characteristic resolution parameter $\varDelta$ as

(3.11)\begin{equation} \varDelta (N,L) \equiv \frac{2L}{N},\end{equation}

which represents the mean spacing between the collocation points in $0 \le r < L$.

A quadrature algorithm presented by Press et al. (Reference Press, Teukolsky, Vetterling and Flannery2007, pp. 179–194) is implemented and all abscissas and weights are computed with an absolute precision error less than $10^{-15}$. The quadrature converts the integration formula to the weighted sum of the function values evaluated at the collocation points and, consequently, the integral of (3.7) finally becomes the discretised formula

(3.12)\begin{equation} f_n^m \simeq \frac{(2n+1) (n-|m|)!}{2(n+|m|)!} \sum_{j=1}^{N} \varpi_j f_m(r_j) P_n^m(\zeta_j).\end{equation}

It is convenient in practice to conceal the factorial coefficient term by defining the normalised mapped Legendre functions and coefficients as follows:

(3.13a,b)\begin{equation} \hat{P}_{L_n}^{m}(r) \equiv P_{L_n}^{m}(r) \sqrt{\frac{(2n+1) (n-|m|)!}{2(n+|m|)!}}, \quad \hat{f}_n^m \equiv f_n^m \sqrt{\frac{2(n+|m|)!}{(2n+1) (n-|m|)!}}. \end{equation}

Using these normalised terms, (3.12) can be expressed as

(3.14)\begin{equation} \hat{f}_n^m \simeq \sum_{j=1}^{N} \varpi_i f_m(r_j) \hat{P}_n^m(\zeta_j),\end{equation}

and, moreover, (3.8) at $r = r_j$ maintains the identical form

(3.15)\begin{equation} f_m (r_j) \, {\rm e}^{\mathrm{i}(m \phi + \kappa z)} = \sum_{n=|m|}^{|m|+M-1} \hat{f}_n^m \hat{P}_{L_n}^{m} (r_j) \, {\rm e}^{\mathrm{i} (m\phi + \kappa z)}. \end{equation}

As a preliminary step of the mapped Legendre spectral collocation method, we need to compute (1) the Gauss–Legendre abscissas $\zeta _i$, (2) weights $\varpi _i$, (3) radial collocation points $r_i$, and (4) normalised mapped Legendre functions evaluated at the collocation points $\hat {P}_{L_n}^m (r_i)$. The normalisation procedure may require temporary multiple-precision arithmetic to handle large function values and factorials if one uses $N$ larger than about 170. There have been several multi-precision arithmetic libraries available recently and we consider using the FM multiple-precision package (Smith Reference Smith2003). All essential computations ahead, however, can be performed under typical double-precision arithmetic.

It is noted that the number of abscissas (or collocation points) $N$ must be equal to or larger than the number of basis elements $M$ for the sake of a proper transform between physical space $( \,f_m(r_1) , \ldots, f_m(r_N) )$ and spectral space $(\, \hat {f}_{|m|}^{m} , \ldots, \hat {f}_{|m|+M-1}^{m} )$. On the other hand, due to the even and odd parity of the associate Legendre functions, taking even $N$ and $M$ can reduce the work by half in the transform procedure (Matsushima & Marcus Reference Matsushima and Marcus1997). Consequently, we set both $N$ and $M$ to be even and $N = M + 2$ in further analyses unless otherwise specified.

Finally, we discuss how to apply the mapped Legendre spectral collocation method to the present problem. Recalling (2.25) where $\mathbb {P}_{m \kappa } ( \tilde {\boldsymbol u} ) = (\tilde {\psi }, \tilde {\chi })$, we write

(3.16)$$\begin{gather} \tilde{\psi} (r;m,\kappa) \, {\rm e}^{\mathrm{i} (m\phi + \kappa z)} = \sum_{n=|m|}^{|m|+M-1} \tilde{\psi}_n^{m \kappa} \hat{P}_{L_n}^{m} (r) \, {\rm e}^{\mathrm{i} (m\phi + \kappa z)} , \end{gather}$$

(3.17)$$\begin{gather}\tilde{\chi} (r;m,\kappa) \, {\rm e}^{\mathrm{i} (m\phi + \kappa z)} = \sum_{n=|m|}^{|m|+M-1} \tilde{\chi}_n^{m \kappa} \hat{P}_{L_n}^{m} (r) \, {\rm e}^{\mathrm{i} (m\phi + \kappa z)} . \end{gather}$$

We point out that when $\tilde {\psi }$ is expressed in the partial sums above, it obeys the boundary conditions of an analytic scalar at the origin, i.e. as $r \rightarrow 0$,

(3.18)\begin{equation} \tilde{\psi}(r; m, \kappa) \rightarrow r^{|m|} \sum_{i =0}^{\infty} a_i^{m\kappa} r^{2 i},\end{equation}

where $a_0^{m\kappa }, a_1^{m\kappa }, \ldots$ are constants (see Eisen et al. Reference Eisen, Heinrichs and Witsch1991; Matsushima & Marcus Reference Matsushima and Marcus1995). Similar analyticity conditions are obeyed by $\tilde {\chi }(r; m, \kappa )$ and, therefore, the perturbation velocity field $\tilde {\boldsymbol u}(r) \, {\rm e}^{\mathrm {i} (m\phi + \kappa z)}$ is also analytic at the origin (see Appendix B). Due to the properties of the mapped Legendre functions, the perturbation vorticity also decays as $r \rightarrow \infty$ (Matsushima & Marcus Reference Matsushima and Marcus1997). As a consequence, $\mathbb {P}_{m \kappa } ( \tilde {\boldsymbol u} )$ can be uniquely represented by $2M$ spectral coefficients of $\tilde {\psi }_{|m|}^{m\kappa }, \ldots, \tilde {\psi }_{|m|+M-1}^{m\kappa }$ and $\tilde {\chi }_{|m|}^{m\kappa }, \ldots, \tilde {\chi }_{|m|+M-1}^{m\kappa }$. We may discretise the eigenvalue problem for viscous cases in (2.26) as

(3.19)\begin{equation} \sigma \begin{pmatrix} \tilde{\psi}_{|m|}^{m\kappa} \\ \vdots \\ \tilde{\psi}_{|m|+M-1}^{m\kappa} \\ \tilde{\chi}_{|m|}^{m\kappa} \\ \vdots \\ \tilde{\chi}_{|m|+M-1}^{m\kappa} \end{pmatrix} = {\boldsymbol{\mathsf{L}}}_{m\kappa}^{\nu} \begin{pmatrix} \tilde{\psi}_{|m|}^{m\kappa} \\ \vdots \\ \tilde{\psi}_{|m|+M-1}^{m\kappa} \\ \tilde{\chi}_{|m|}^{m\kappa} \\ \vdots \\ \tilde{\chi}_{|m|+M-1}^{m\kappa} \end{pmatrix},\end{equation}

where ${\boldsymbol{\mathsf{L}}}_{m\kappa }^{\nu }$ is a $2M \times 2M$ complex matrix representing the linear operator $\mathcal {L}_{m\kappa }^{\nu }$. In a similar sense, we can define ${\boldsymbol{\mathsf{L}}}_{m\kappa }^{0}$ representing $\mathcal {L}_{m\kappa }^{0}$ for the inviscid analysis and

(3.20)\begin{equation} {\boldsymbol{\mathsf{L}}}_{m\kappa}^{\nu} = {\boldsymbol{\mathsf{L}}}_{m\kappa}^{0} + Re^{-1} {\boldsymbol{\mathsf{H}}}. \end{equation}

Here ${\boldsymbol{\mathsf{H}}}$ is a matrix representation of the Laplacian $\nabla _{m\kappa }^2$ acting on the spectral coefficients $\tilde {\psi }_{|m|}^{m\kappa }, \ldots, \tilde {\psi }_{|m|+M-1}^{m\kappa }$ and $\tilde {\chi }_{|m|}^{m\kappa }, \ldots, \tilde {\chi }_{|m|+M-1}^{m\kappa }$, respectively. For a scalar function expanded by the mapped Legendre functions $a(r) = \sum _{n\ge |m|} a_n^m P_{L_n}^m (r)$, if we expand its Laplacian as $\nabla _{m\kappa }^2 a (r) = \sum _{n\ge |m|} b_n^m P_{L_n}^m (r)$, then the coefficients $a_n^m$ and $b_n^m$ constitute the relationship, for all $n \ge |m|$,

(3.21)\begin{align} b_n^m &=-\left[\frac{(n-|m|-1)(n-|m|)(n-2)(n-1)}{(2n-3)(2n-1)L^2} \right]a_{n-2}^{m} +\left[\frac{2n(n-|m|)(n-1)}{(2n-1)L^2} \right]a_{n-1}^{m} \nonumber\\ &\quad -\left[\frac{2n(n+1)(3n^2 + 3n - m^2 - 2)}{(2n-1)(2n+3)L^2} + \kappa^2 \right]a_{n}^{m} \nonumber\\ &\quad+\left[\frac{2(n+1)(n+|m|+1)(n+2)}{(2n+3)L^2} \right]a_{n+1}^{m} \nonumber\\ &\quad-\left[\frac{(n+|m|+1)(n+|m|+2)(n+2)(n+3)}{(2n+3)(2n+5)L^2} \right]a_{n+2}^{m}, \end{align}

under the assumption that $a_n^m \equiv 0$ if $n$ is less than $|m|$ (Matsushima & Marcus Reference Matsushima and Marcus1997, p. 344). Here ${\boldsymbol{\mathsf{H}}}$ can be formulated by (3.21).

The formulation of ${\boldsymbol{\mathsf{L}}}_{m\kappa }^{0}$ involves the vector products in physical space and is conducted using a pseudospectral approach based on the Gauss–Legendre quadrature rule. Reconstructing $\tilde {\boldsymbol u}$ from $\mathbb {P}_{m \kappa } ( \tilde {\boldsymbol u} )$ via (2.12), we evaluate the vector products $\bar {\boldsymbol {U}} \times \tilde {\boldsymbol \omega }$ and $\bar {\boldsymbol {\omega }} \times \tilde {\boldsymbol u}$ at $N$ radial collocation points and apply $\mathbb {P}_{m \kappa }$ again. As for the detailed algorithm including the numerical implementation of $\mathbb {P}_{m\kappa }$ as well as its inverse, refer to (69) and (70) in Matsushima & Marcus (Reference Matsushima and Marcus1997), providing the spectral coefficients of $\mathbb {P}_{m \kappa } (\bar {\boldsymbol {U}} \times \tilde {\boldsymbol \omega })$ and $\mathbb {P}_{m \kappa } (\bar {\boldsymbol {\omega }} \times \tilde {\boldsymbol u})$. Integration in these equations can be performed numerically by the Gauss–Legendre quadrature rule, as given in (3.14). Following this procedure, we can compute the $i$th column vector of ${\boldsymbol{\mathsf{L}}}_{m\kappa }^{0}$ by substituting the $i$th standard unit vector $\hat {\boldsymbol {e}}_i \in \mathbb {R}^{2M}$ for $( \tilde {\psi }_{|m|}^{m\kappa }, \ldots, \tilde {\psi }_{|m|+M-1}^{m\kappa }, \tilde {\chi }_{|m|}^{m\kappa }, \ldots, \tilde {\chi }_{|m|+M-1}^{m\kappa } )$.

A global eigenvalue problem solver with the QR algorithm for non-Hermitian matrices, based on the Lapack routine named Zgeev, is used to solve the discretised eigenvalue problem. The procedure of constructing a global matrix and finding all eigenvalues has been established in previous studies, such as Fabre et al. (Reference Fabre, Sipp and Jacquin2006, p. 241). However, as shown in (3.19), our formulation directly results in the standard eigenvalue problem rather than the generalized form. Thus, it is sufficient to construct only one matrix of dimension $2M \times 2M$, with a reduction in the number of state variables from 4 to 2.

3.3. Numerical parameters and their effects

The mapped Legendre spectral collocation method comprises of three adjustable numerical parameters: $M, N$ and $L$. The first two parameters are commonly used in most spectral collocation methods, while the last parameter is unique to our method. This section elaborates on the impact of each parameter on the numerical method's performance and provides guidelines on their selection.

3.3.3. Map parameter for Legendre functions $L$

The map parameter $L$ provides an additional level of computational freedom that distinguishes the present numerical method from others. We highlight three significant roles of this parameter, two of which are related to spatial resolution in physical space and the other to basis change in spectral space.

In physical space, when $M$ (and $N$) is fixed, a change in $L$ results in two anti-complementary effects with respect to spatial resolution, as shown in figure 2. When $L$ increases, the high-resolution region $0 \le r < L$, where half of the collocation points are clustered, expands, which has a positive effect. However, it negatively impacts the resolution, especially in the high-resolution region, where $\varDelta$ increases with $L$. Increasing $N = M+2$ may compensate for the loss in resolution. However, if $M$ is already at a practical limit due to the computing budget, expanding the high-resolution region by increasing $L$ should stop when $\varDelta$ remains satisfactorily small. The requirement for satisfaction should be specific to the eigenmodes to be resolved, which will be discussed in each analysis section later. Similar discussions can be made in the opposite direction when decreasing $L$.

Figure 2. Changes in distribution of the collocation points with respect to $L$ given $N = 52$. Some collocation points at large radii are omitted. The high-resolution region is $0\le r < L$, where half of the collocation points are clustered around the origin. As $L$ increases, the high-resolution region is expanded. However, the mean spacing $\varDelta$ grows simultaneously. Parameter $L$ should be chosen carefully to balance these anti-complementary effects.

In spectral space, changing $L$ entirely replaces the complete basis function set. For instance, when $L=A$ and $L=B$, the spectral method can be constructed on either of two different complete basis sets, i.e. $\{ P_{A_{|m|}}^{m} ,P_{A_{|m|+1}}^{m} \ldots \}$ or $\{ P_{B_{|m|}}^{m},P_{B_{|m|+1}}^{m} , \ldots \}$. Since orthogonality among the basis functions does not necessarily hold across the basis sets, an eigenmode found with $L=A$ can differ from that found with $L=B$. If $B$ differs from $A$ by an infinitesimal amount, our method makes it possible to find eigenmodes that continuously vary if they exist. This was thought to be hardly achievable via classic eigenvalue solvers due to discretisation (cf. Mao & Sherwin Reference Mao and Sherwin2011, p. 11). Once the numerical method's convergence is secured by sufficiently large $M$ and $N$, we explore such non-normal eigenmodes that vary continuously by fine tuning $L$.

3.4. Validation

To confirm the numerical validity of our method, we compared some eigenvalues from the discrete branch of the spectra with those previously calculated by Mayer & Powell (Reference Mayer and Powell1992). They also used a spectral collocation method but with Chebyshev polynomials as radial basis functions over an artificially truncated radial domain, rather than the mapped Legendre basis functions over an unbounded radial domain we use. For comparison, we linearly scaled the eigenvalues reported in Mayer & Powell (Reference Mayer and Powell1992) to match the $q$-vortex model used in our study because the azimuthal velocity component is scaled by $q$ in their study, whereas we adjust the axial velocity component.

We compared the most unstable eigenvalue calculations for the inviscid case $m=1, \kappa =0.5, q=-0.5$ (or, equivalently, $m=1, \kappa =-0.5, q=0.5$) and the viscous case $m=0, \kappa =0.5, q=1, Re=10^4$ in table 1. We conducted the calculations using three different numbers of basis elements $M$ (20, 40 and 80) and three different map parameters $L$ (8, 4 and 2). Our results show that the trend towards convergence is apparent as $M$ increases and $L$ decreases. As we discuss in terms of the characteristic resolution parameter $\varDelta$ defined in (3.11), both parameters influence the numerical resolution. Increasing $M$ leads to an increase in the number of radial collocation points $N$, while decreasing $L$ improves spatial resolution by filling the inner high-resolution region ($0 \le r < L$) with more collocation points (see figure 2). However, this comes at the expense of reducing the range of the high-resolution region and effectively shrinking the radial domain by placing the collocation point with the largest radius at $r_N = L \sqrt {(1+\zeta _N)/(1-\zeta _N)}$, which can lead to inaccuracies if any significant portion of the solution exists either in the outer low-resolution region or outside the effective limit. The convergence test of $\sigma ^{{\dagger} }_{viscous}$ with $M = 20$ in table 1 partially demonstrates this concern. When we compare the eigenvalues computed with $L=4$ and $L=2$, the latter shows no clear improvement in convergence compared with the former, despite having a smaller $L$. Even small $L$ causes the eigenvalue's real part to move further away from the reference value of Mayer & Powell (Reference Mayer and Powell1992). Therefore, we must keep in mind that blindly pursuing small $L$ does not guarantee better convergence, although using large $M$ is always favoured for numerical convergence.

Table 1. Comparison of the eigenvalues associated with the most unstable mode (indicated with a superscript ${\dagger}$) for the inviscid case with $m=1, \kappa = 0.5, q=-0.5$ and for the viscous case with $m=0, \kappa = 0.5, q=1, Re = 10^4$. The table illustrates how the values change when we alter the map parameter $L$ and the number of radial mapped Legendre basis functions $M$. The last row displays the values obtained by Mayer & Powell (Reference Mayer and Powell1992), who employed up to 200 radial Chebyshev basis functions. Their published eigenvalues were appropriately rescaled to fit the $q$-vortex model employed in our study. Our numerically computed eigenvalues tend towards a fixed point as we increase $M$ beyond 40. It should be noted that the size of the matrix eigenvalue problem system is $2M$ for our method and $3M$ for that of Mayer & Powell (Reference Mayer and Powell1992). Thus, even when using the same $M$, our method is expected to require $(2/3)^3$ less work than theirs.

The high-resolution range of the present method, represented by $L$, should not match the domain truncation radius in the method of Mayer & Powell (Reference Mayer and Powell1992). Adjusting the high-resolution range through $L$ has no impact on the unbounded nature of the domain and can be customised essentially. However, altering the domain truncation radius fundamentally harms the unbounded nature of the domain and must be set to its maximum computing limit. On the other hand, we achieve the same accuracy as Mayer & Powell (Reference Mayer and Powell1992) with roughly three times smaller $M$, which supports the numerical efficiency of our method. Presumably, our method is around ten times more computationally efficient in solving matrix eigenvalue problems that scale as $O(M^3)$. We believe this is mainly because their simple algebraic mapping of Chebyshev collocation points (see Ash & Khorrami Reference Ash and Khorrami1995, p. 357) clusters approximately one-third of the collocation points near the artificial outer radial boundary, where vortex motion is near zero and not important by assumption. Such collocation points do not significantly contribute to solving the problem, resulting in an inefficient use of numerical resources.

Note that the eigenmodes shown here are regular and have no singularities, as depicted in figure 3. Such regular eigenmodes are expanded by a finite number of radial basis elements that are already regular, and as shown in table 1, their numerical results converge exponentially with increasing $M$. However, singular eigenmodes can only be expressed exactly when an infinite sum of mapped Legendre functions is taken (see Gottlieb & Orszag Reference Gottlieb and Orszag1977). Nonetheless, as stated in the preliminary remarks, we are essentially interested in physical eigenmodes, i.e. those without singularities and computed numerically with a small spectral residual error. The current validation is strong enough to underpin this objective. Later in this paper (see § 6.1.4), we present some eigenmodes that have convincing signatures of viscous remnants after regularising the inviscid critical-layer singularities. These singularities become regularised but are still nearly singular regions of local rapid oscillations. We can find the value of $M$ at which these eigenmodes are spatially resolved, even if it typically goes beyond 80. Also note that in this respect, we only peripherally examine their inviscid counterparts with the critical-layer singularities using our numerical method (see § 5.1.2).

Figure 3. A comparison of our numerical calculation with that of Mayer & Powell (Reference Mayer and Powell1992). Shown is the radial velocity component of the most unstable eigenmode for the validation cases (a) $(m, \kappa, q, Re) = (1, 0.5, -0.5, \infty )$ and (b) $(m, \kappa, q, Re) = (0, 0.5, 1, 10^4)$, where the maximum of ${\rm Re}(\tilde {u}_r)$ is normalised to unity. Numerical parameters are $M=80$ and $L=2$. Note that Mayer & Powell (Reference Mayer and Powell1992) only plotted the real parts of the eigenmodes.

5.1. Numerical spectra and eigenmodes

In figure 6 we present the eigenvalues of two inviscid vortices: the Lamb–Oseen vortex with $(m, \kappa, q) = (1, 1.0, \infty )$ and the strong swirling Batchelor vortex with $(m, \kappa, q) = (2, 3.0, 4.0)$. By comparing these two vortices, we demonstrate their common properties and extract features that can be generalized to vortices with large $q$ and moderate $m$ and $\kappa$ of order unity, which are thought to be relevant for practical aeronautical applications (see Fabre & Jacquin Reference Fabre and Jacquin2004, pp. 258–259). To observe the effect of the numerical parameter $M$, we computed each vortex in four ways: with $M = 100, 200,300$ and $400$. Analytically, every eigenvalue is expected to lie on $\mathrm {i}\mathbb {R}$. The shaded area in each plot is the non-normal region of the spectra, indicating the frequency band that includes the analytic range of $\sigma _c^{0}$.

Figure 6. Numerical spectra computed with zero viscosity (a) for the Lamb–Oseen vortex ($q \rightarrow \infty$) in $(m, \kappa ) = (1, 1.0)$ and (b) for the strong swirling Batchelor vortex ($q = 4.0)$ in $(m, \kappa ) = (2, 3.0)$ with respect to $M=100,200,300$ and $400$. Here $L$ is fixed at $6.0$ and $N=M+2$. A shaded band in each plot indicates the non-normal region where $\sigma _c^{0}$ appears. The larger $M$ we use, the closer the numerical spectra is to their true shape (see figure 4a). However, with sufficiently large values of $M$ and with appropriately tuned values of $L$, the under-resolved can be corrected, making all eigenvalues lie on the imaginary axis; see figure 10.

Clearly, all these numerical spectra contain some eigenvalues that are incorrect (i.e. not on the imaginary axis). We can observe three families of numerical eigenvalues. A discrete family ($+$) corresponds to $\sigma _d^{0}$, where the eigenvalues are discrete and located outside the shaded area. An inviscid critical-layer family ($\bullet$) corresponds to $\sigma _c^{0}$. Its eigenvalues lie on the imaginary axis, are within the shaded area and the number of them increases as $M$ increases. Finally, a family of under-resolved eigenvalues ($\times$), which, had they been spatially well resolved, would have been eigenvalues belonging to $\sigma _c^{0}$ and lie on the imaginary axis. Instead, these eigenvalues lie off the imaginary axis and within the shaded area. These under-resolved eigenvalues are characterised by non-zero real parts with absolute values typically greater than $10^{-10}$ as a result of numerical discretisation errors. The eigenvalues form clouds of structures that are symmetric about the imaginary axis. The cloud structures are due to insufficient spatial resolution, and the absolute values of the real parts of the eigenvalues tend to increase as the value of $q$ decreases. As $M$ increases, the absolute values of the real parts of the eigenvalues tend to decrease, and the cloud of eigenvalues gets ‘squeezed’ to the imaginary axis, which is similar to the ‘squeeze’ observed by Mayer & Powell (Reference Mayer and Powell1992) when they increased the number of Chebyshev basis elements in their spectral method calculation.

5.1.1. Discrete eigenmodes

Although $\sigma _d^{0}$ and the discrete eigenmodes are not the main focus of this paper, it is worthwhile to confirm their convergence properties. Figure 6 shows that the discrete eigenmodes associated with eigenvalues away from the accumulation points (see Gallay & Smets Reference Gallay and Smets2020, pp. 14–16) (i.e. intersections of the imaginary axis with the lower boundary of the shaded regions in figure 6) are spatially resolved for $M \geq 100, L=6$ and $N=M+2$. For these values of $L, M$ and $N$, each eigenvalue approaches a fixed point as $M$ increases. The discrete eigenmodes are distinguishable from each other by their radial structures and, in particular, by the number of ‘wiggles’ (intervals between two neighbouring zeros) as a function of radius. Typically, the eigenmodes with eigenvalues farthest from the accumulation points have the fewest wiggles, as shown in figure 7. The discrete eigenmodes have an increasing number of wiggles as the eigenvalue approaches the accumulation point, forming a countably infinite, linearly independent set in the eigenspace of $\mathcal {L}_{m\kappa }^{0}$.

Figure 7. Radial velocity profiles of the inviscid discrete eigenmodes associated with three largest $|{\rm Im} (\sigma )|$ (a) for the Lamb–Oseen vortex ($q \rightarrow \infty$) in $(m, \kappa ) = (1, 1.0)$ and (b) for the strong swirling Batchelor vortex ($q = 4.0)$ in $(m, \kappa ) = (2, 3.0)$. The maximum of ${\rm Re}(r \tilde {u}_r)$ is normalised to unity. Here $M=400$ and $L=6.0$ are used. The number of ‘wiggles’ in and around the vortex core distinguishes each discrete eigenmode. Note that, for the eigenmodes that are neutrally stable, the phase of the eigenmodes can be chosen such that the radial velocity components are made to be either real or pure imaginary for all $r$. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,\infty )$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,\infty )$.

The eigenmodes with discrete eigenvalues and ${\rm Im}(\sigma ) / m > 0$ in figure 6 were referred to as ‘countergrade’ by Fabre et al. (Reference Fabre, Sipp and Jacquin2006). They appear to exist only for eigenmodes with specific values of $m$, including $m = \pm 1$ (see Gallay & Smets Reference Gallay and Smets2020). However, we remark that these eigenmodes are also legitimate solutions to the problem and can be spatially resolved using our numerical method, just like those shown in figure 7. They are also expected to be crucial for triad-resonant interactions among the eigenmodes and will be actively considered in further instability studies.

The numerically computed eigenmodes correspond to the eigenvectors of the $2M \times 2M$ matrix ${\boldsymbol{\mathsf{L}}}_{m\kappa }^{0}$, which implies that the maximum number of numerical eigenmodes that can be obtained is $2M$ under double-precision arithmetic. The number of discrete eigenmodes that our numerical solver can find increases with an increase in $M$. For instance, in the case of a strong swirling Batchelor vortex illustrated in figure 6(b), the number of discrete eigenmodes (i.e. in the $\sigma _d^{0}$ spectrum) is $4, 7, 9$ and $11$ for $M=100, 200, 300$ and $400$, respectively. This behaviour is expected because a finer spatial resolution is required to resolve more wiggles in the eigenmode structure. If $n$ wiggles exist in the vortex core region $(r \le 1.122)$, whose non-dimensionalised scale is of order unity, the necessary spatial resolution to resolve all the wiggles is $O(1/n)$. As $\varDelta = 2L/(M+2) \sim O(1/M)$ in our analysis, the proportionality of $n$ to $M$ is verified. The implication of this scaling is that the number of discrete eigenmodes accounts for only a small portion of the total number of numerical eigenmodes computed, and the vast majority are associated with the non-regular, continuous part of the spectrum, $\sigma _{c}^{0}$.

5.1.2. Inviscid critical-layer eigenmodes

We emphasise that our essential interest lies in eigenmodes with small, but non-zero viscosity. This ensures that the eigenmodes can be physical and do not have difficult-to-compute singularities. Nevertheless, it is still intriguing to compute the eigenmodes with $\nu \equiv 0$, which are numerically (not physically) regularisable by the spatial discretisation. By selecting a suitably large value of $M$ and an appropriate value for the mapping parameter $L$ (see § 5.2), we can resolve the spatial structure of the inviscid eigenmode outside the critical-layer singularity neighbourhood well. In addition, the numerical error in the eigenvalue, caused by the slow decay of spectral coefficients or the Gibbs phenomenon around the critical-layer singularity, can be kept adequately small and the eigenvalues correctly lie on the imaginary axis.

Figure 8 shows some critical-layer eigenmodes, which were numerically obtained with $M=400$. The real parts of the eigenvalues are zero, and the velocity components are either real or purely imaginary for all $r$, with a suitable phase choice. Typically, $r_c$ increases as $|\sigma |$ decreases along the critical-layer spectrum. The singular behaviour of abrupt slope change commonly occurs at the critical-layer singularity, as predicted analytically by (4.9). As stated in § 1.3, we cannot claim that they are perfectly resolved due to the presence of the singularity and the continuous nature of their associated spectrum. However, our focus is not on their exact convergence but rather on their well-behaved spatial structure outside the neighbourhood of the singularity, achieved by using a large $M$, along with purely imaginary eigenvalues that conform to analytic expectations. We use this information later to study the spatial correspondence of eigenmodes with non-zero viscosity to determine which viscous eigenmodes are of physical relevance.

Figure 8. Radial velocity profiles of three inviscid, critical-layer eigenmodes (a) for the Lamb–Oseen vortex ($q \rightarrow \infty$) in $(m, \kappa ) = (1, 1.0)$ and (b) for the strong swirling Batchelor vortex ($q = 4.0)$ in $(m, \kappa ) = (2, 3.0)$. The maximum of the real part of $r \tilde {u}_r$ is normalised to unity. Here $M=400, N=M+2$ and $L=6.0$ are used. For each eigenmode, the vertical dashed line indicates the critical-layer location $r_c$ determined by (4.9). Note that all of the radial components of the velocity can be made to be real valued for all $r$ by a proper choice of phase as they are neutrally stable. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,\infty )$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,\infty )$.

For $r < r_c$, the radial velocity components of the inviscid critical-layer eigenmodes oscillate in $r$, and the number of oscillations decreases as the value of $r_c$ increases (or equivalently, as $|\sigma |$ decreases). Consequently, when $r_c > r_{cc}$ for some value $r_{cc}$, there is no longer one full oscillation. In our numerical investigation we found that for the Lamb–Oseen vortex with $(m,\kappa ) = (1,1.0), r_{cc}$ equals $2.2$, which corresponds to $\sigma = -0.21\mathrm {i}$. We believe that our numerically found value of $r_{cc}$ approximately coincides with the theoretical threshold of $r=2.124$, at which the analytic solutions obtained by the Frobenius method change form regarding the roots of the indicial equation (see Gallay & Smets Reference Gallay and Smets2020, pp. 20 and 50). For $r > r_c$, the radial velocity components of the critical-layer eigenmodes are not oscillatory, and the amplitudes of $r\tilde {u}_r$ achieve the local maximum or minimum values close to $r=r_c$, before decreasing monotonically as rapidly as those of the discrete eigenmodes, as shown in figure 7.

5.1.3. Under-resolved eigenmodes

The under-resolved eigenmodes, which, if resolved, would be part of the spectrum with $\sigma _c^{0}$, have eigenvalues in the complex $\sigma$ plane on either side of the imaginary axis in the shaded region in figure 6. The eigenvalues come in pairs, with one unstable and one stable eigenmode. The reflection symmetry with respect to the imaginary axis is due to the fact that the analytic operator $\mathcal {L}_{m\kappa }^{0}$ is time reversible (cf. Bölle et al. Reference Bölle, Brion, Robinet, Sipp and Jacquin2021, p. 10). Therefore, the eigenmode $(\tilde {u}_{r}, \tilde {u}_{\phi }, \tilde {u}_{z}, \tilde {p})$ with eigenvalue $\sigma _s$ corresponds to the eigenmode $(\tilde {u}_{r}^{*}, -\tilde {u}_{\phi }^{*}, -\tilde {u}_{z}^{*}, -\tilde {p}^{*})$ with eigenvalue $-\sigma _s^*$.

Some examples of under-resolved eigenmodes are shown in figure 9. These eigenmodes are qualitatively incorrect because (1) unlike the eigenmodes in figures 7 and 8, there is no choice of phase that makes their radial components real for all $r$, and more importantly, because (2) we know that their eigenvalues should be purely imaginary when $q$ is sufficiently large, and they are not. However, these eigenmodes appear to exhibit no other distinguishing properties, except for the two properties listed above, from the inviscid critical-layer eigenmodes in figure 8. It should be noted that they have been called ‘spurious’ in previous numerical studies (see Mayer & Powell Reference Mayer and Powell1992; Heaton Reference Heaton2007), of which the usage was similar to our clarification given in § 1.3. However, instead of following convention, we propose naming these numerical eigenmodes ‘under-resolved’ eigenmodes of the continuous part of the inviscid spectrum. In this way, we put more emphasis on the fact that adjusting the numerical parameters can ‘correct’ these eigenmodes so that neither of the two key properties listed above apply.

Figure 9. Radial velocity profiles of two inviscid under-resolved eigenmodes whose eigenvalues are symmetric about the imaginary axis (a) for the Lamb–Oseen vortex ($q \rightarrow \infty$) in $(m, \kappa ) = (1, 1.0)$ and (b) for the strong swirling Batchelor vortex ($q = 4.0)$ in $(m, \kappa ) = (2, 3.0)$. The maximum of the real part of $r \tilde {u}_r$ is normalised to unity. Here $M=400$ and $L=6.0$ are used. For each eigenmode, an abrupt slope change occurs at the vertical dashed line at the critical-layer location $r =r_c$ (which is determined from (4.9) by ignoring the real part of the eigenvalue), indicating that they will become correct critical-layer eigenmodes given more resolution. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,\infty )$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,\infty )$.

By examining the spatial structure of the under-resolved eigenmodes, we can detect sudden changes in slope at the critical-layer singularity point at $r = r_c$. The value of $r_c$ is obtained by setting the imaginary part of either of the eigenvalues ${\rm Im}(\sigma _s)$ to $\sigma _c$ in (4.9). The break in slope confirms that the under-resolved eigenmodes originate from $\sigma _c^{0}$ and indicates that they have lost their neutrally stable property due to numerical errors at the critical-layer singularity.

Correcting the under-resolved eigenmodes is crucial, not only for correctly evaluating $\sigma _c^{0}$ but also for the following reasons. Despite their invalid origin, half of the under-resolved eigenmodes in ${\rm Re}(\sigma ) > 0$ erroneously suggest that the wake vortex is linearly unstable. In the future, we plan to use the computed velocity eigenmodes from the present numerical method to initialise an initial-value code that solves the full nonlinear equations of motion given by (2.6) and (2.7). Inappropriately computed eigenmodes that grow erroneously, rather than remain neutrally stable, are likely to corrupt these calculations.

5.2. Correction of the under-resolved eigenmodes

An intriguing question is whether the under-resolved eigenmodes tend towards something as $M$ increases and, furthermore, what is the potential outcome of such a convergence? In the beginning of this section it was argued that the real part of eigenvalues remains at zero (i.e. all eigenmodes are neutrally stable) when $q$ is sufficiently large. In figure 6 this can be observed as the ‘squeeze’ of the eigenvalue cluster towards the imaginary axis. However, we have also indicated that the imaginary part of the eigenvalues may not converge to a fixed point, instead continuing to evolve along the imaginary axis. Therefore, instead of concentrating on the convergence of individual under-resolved eigenmodes to a fixed point, it is more pragmatic to aim to ‘correct’ the set of eigenmodes as a whole, that is, to restore their neutrally stable nature. The ‘correction’ means that we comprehensively treat the entire set of eigenmodes as a single entity, which complies with the usage of this term in this section up to this point.

To ‘correct’ the under-resolved eigenmodes, we first consider increasing $M$ to its largest possible value within the available computing resources. However, increasing $M$ is generally undesirable because it always comes at a steep computational expense; the cost of finding the eigenmodes is proportional to $(2M)^3$. Instead, we may consider dealing with the mapping parameter $L$, where the novelty and usefulness of our method come from. Here $L$ controls the spatial resolution locally as a function of $r$. As seen from the resolution parameter $\varDelta$ in (3.11), $L$ controls the spatial resolution by providing more resolution near the radial origin (i.e. $0 \leq r < L$). It is important to note that changing or tuning $L$ does not affect the cost of computation.

For a fixed $M$, with $N=M+2$, figure 10 shows five numerical eigenvalue spectra for two prescribed cases with different values of $L$, varying from $1.0$ to $5.0$. Overall, decreasing $L$ brings the numerical spectra closer to the imaginary axis. In particular, some values of $L$ enable complete resolution of $\sigma _c^{0}$ on the imaginary axis, which cannot be achieved by increasing $M$ within a modest computing budget. However, decreasing $L$ does not always shrink the clouds of eigenvalues closer to the imaginary axis. We separate the numerically computed eigenvalues and eigenmodes into two categories: those with the critical layer $r=r_c$ located in the high-resolution region $0 \le r < L$, and those where $r_c$ is in the low-resolution region $r \ge L$. Figure 10 indicates this separation with horizontal dashed lines. For the former, the cloud structures vanish as $L$ decreases, and $\sigma _c^0$ is correctly resolved. In contrast, for the latter, the cloud structures persist or even recur if $L$ is too small, resulting from excessive concentration of collocation points solely around the centre. Once $\sigma _c^{0}$ is satisfactorily resolved, adjusting $L$ should stop to keep the portion of $\sigma _c^{0}$ resolved in the high-resolution region as large as possible. For instance, in figure 10 we propose setting $L$ between $3.0$ and $4.0$ for the Lamb–Oseen vortex case and between $1.0$ and $2.0$ for the Batchelor vortex case.

Figure 10. Numerical spectra computed at zero viscosity (a) for the Lamb–Oseen vortex ($q\rightarrow \infty$) in $(m,\kappa )=(1, 1.0)$ and (b) for the strong swirling Batchelor vortex ($q=4.0$) in $(m,\kappa )=(2, 3.0)$ with respect to $L= 5.0, 4.0, 3.0, 2.0$ and $1.0$. Here $M$ is fixed at $400$ and $N=M+2$. In each plot a shaded band indicates the non-normal region in which $\sigma _c^{0}$ appears, and a horizontal dashed line represents the threshold used to determine if the critical layer $r=r_c$ is located within the high-resolution region $0\le r < L$. It should be noted that there is a one-to-one correspondence between a critical-layer eigenvalue $\sigma$ and a critical-layer radius $r_c$, as seen in (4.9). Furthermore, $r_c$ approaches zero at the bottom of the shaded band, ${\rm Im}(\sigma ) = m + \kappa /q$, and monotonically increases towards infinity as $|\sigma |$ becomes smaller. By tuning $L$, under-resolved eigenmodes can be corrected without requiring additional computing resources. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,\infty )$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,\infty )$.

To provide a detailed explanation of what we have seen, we must revisit the differences in the way $M$ and $L$ operate in the current numerical method, as stated in § 3.3. One of the roles of $L$ is to serve as a tuning parameter for spatial resolution in physical space, whereas $M$ determines the number of basis elements used in spectral space. Increasing $M$ allows us to handle eigenmodes with more complex shapes, such as (nearly) singular functions, which often have more wiggling and are thus more numerically sensitive. Here $M$ has only an indirect effect on spatial resolution through $N$, which is required to be greater than or equal to $M$. On the other hand, the critical-layer singularity is essentially a phenomenon that occurs in physical space. Although using more spectral basis elements relates to improving spatial resolution because we set $N = M+2$, the main contribution to dealing with the critical-layer singularity with minimal errors comes from the latter. Therefore, it can be more effective to use $L$ to directly control resolution and suppress the emergence of under-resolved eigenmodes, rather than using $M$. It is worth noting that increasing $N$ to very large values while keeping $M$ constant can also reduce the number of under-resolved eigenvalues to some extent. This observation supports that high spatial resolution is crucial for suppressing under-resolved eigenmodes.

If one aims to correct the under-resolved eigenmodes and obtain $\sigma _c^{0}$ using the present numerical method, the following steps are suggested to properly set up the numerical parameters. Assuming that $M$ is already at the practical maximum due to finite computing budget, and $N$ follows $M+2$.

1. (i) Start with an arbitrarily chosen value of $L$ and gradually decrease it if under-resolved eigenmodes exist, until they vanish in the high-resolution region $0 \le r < L$. This step improves spatial resolution, helping to identify the critical-layer singularity with less numerical error despite the discretisation.

2. (ii) If there are no eigenvalue clouds around the imaginary axis, increase $L$ as long as they do not appear in the numerical spectra. This step expands the high-resolution region where the critical-layer singularity can be accurately treated.

5.3. Pairing in the inviscid critical-layer spectrum

In the numerical spectrum of $\sigma _c^{0}$, we observe that numerical eigenvalues tend to appear in pairs. This pairing phenomenon is demonstrated in the left panel of figure 11, which illustrates the Lamb–Oseen vortex case with $(m,\kappa ) = (1, 1.0)$ that we computed with $M=400$ and $L=3.0$ (see figure 10a). We argue that the pairing in our numerical results arises from a degeneracy resulting from the critical-layer singularity. We refer to Gallay & Smets (Reference Gallay and Smets2020, pp. 19–21), who used the Frobenius method to construct analytic solutions to the problem under the assumption of non-zero $m$ and $\kappa$ with $q \rightarrow \infty$. They showed that if a critical-layer singularity occurs at $r=r_c$, there exists a unique solution with scalar multiplication that is only non-zero on $(0, r_c)$ and another one that is only non-zero on $(r_c, \infty )$. Here we call the inner and outer solutions, respectively. For both the inner and outer solutions, the radial velocity components can be made real valued by an appropriate choice of phase, since their degenerate eigenvalue is purely imaginary. These two solutions are independent of each other, and their linear combination should be the general form of an inviscid critical-layer eigenmode that is singular at $r=r_c$.

Figure 11. (a) Numerical inviscid spectra with no under-resolved eigenmodes for the Lamb–Oseen vortex ($q\rightarrow \infty$) in $(m,\kappa )=(1, 1.0)$ along with a magnified part exhibiting the pairing phenomenon, and (b) four radial velocity profiles of the critical-layer eigenmodes from two neighbouring pairs, labelled as A1/2 and B1/2. Here, $M=400, L=3.0$ and $N=M+2$. Note the similarity in structure within each pair, and the change in the critical-layer location (marked by vertical dashed lines) by one collocation point between these neighbouring pairs. This pairing phenomenon stems from the singular degeneracy in $\sigma _c^{0}$. The linear combination of the pair constructs two independent solutions that are singular at the same critical-layer location and are nearly zero on either $(0,r_c)$ or $(r_c, \infty )$.

We can observe these analytic characteristics in our numerically computed pairs. In the right panel of figure 11 we present the $r \tilde {u}_r$ profiles of the critical-layer eigenmodes from two neighbouring pairs. In each pair, the velocity profiles have an abrupt change in slope across an interval between two collocation points, whose location matches the critical-layer singularity radius calculated by (4.9). The difference in $r_c$ among the neighbouring pairs corresponds to the collocation interval, indicating their continuous emergence. Furthermore, by linearly combining these paired eigenmodes, we can construct the inner and outer solutions as derived analytically, each of which is approximately zero on $(0,r_c)$ or $(r_c, \infty )$. Although their eigenvalues are slightly different, we believe that it is due to the numerical error resulting from the spatial discretisation, which slightly breaks the degeneracy. This error decreases with increasing $M$.

Strictly speaking, the discussion made here is limited to infinite $q$ because the analytic results found in Gallay & Smets (Reference Gallay and Smets2020) were verified in the Lamb–Oseen vortex case, and we can only compare this case. Nonetheless, we remark that we have numerically observed this same pairing phenomenon with finite $q$ (e.g. $q=4.0$). We conjecture that the pairing phenomenon exists for values of $q$ where the eigenmodes are all neutrally stable.

6.1. Numerical spectra and eigenmodes

The Lamb–Oseen vortex with $(m, \kappa, q) = (1, 1.0, \infty )$ and the strong swirling Batchelor vortex with $(m, \kappa, q) = (2, 3.0, 4.0)$ are analysed at $Re = 10^5$. As with the prior inviscid analysis, our aim is to identify common characteristics in the linear vortex dynamics with viscosity for large $q$ and moderate $m$ and $\kappa$. In the current analysis, however, we have a specific focus on the physical relevance of each eigensolution. To observe how the viscous numerical spectra converge, we compute four numerical spectra for each case, with different values of $M$ ranging from $100$ to $400$. The spectra are presented in figure 12, or supplementary movie 1 is available at https://doi.org/10.1017/jfm.2023.455. Based on these spectra, we classify the numerical eigenmodes into five families: unresolved ($-$), discrete ($+$), spurious ($\times$), potential ($\Box$) and viscous critical layer ($\bullet$). Note that with increasing numerical resolution, an eigenmode in the unresolved family will always evolve into an eigenmode in one of the other four families.

Figure 12. Numerical viscous spectra at $Re = 10^5$ (a) for the Lamb–Oseen vortex ($q\rightarrow \infty$) in $(m,\kappa )=(1,1.0)$ and (b) for the strong swirling Batchelor vortex ($q=4.0$) in $(m,\kappa ) = (2,3.0)$ with respect to $M = 100, 200, 300$ and $400$. Here $L$ is fixed at $2.0$ and $N=M+2$. Larger $M$ enables more portion of the spectra to be resolved. Near the right boundary of the potential spectrum there are two distinct branches of the viscous critical-layer spectrum. See supplementary movie 1 for animation. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,10^5)$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,10^5)$.

The eigenvalues in the discrete spectrum, $\sigma _d^{\nu }$, converge to fixed points with increasing $M$. These eigenvalues populate two distinct discrete branches shown near the bottom of the panels in figure 12, in addition to a few locations outside of these branches. The eigenmodes in the discrete spectrum were known previously and were studied by other authors. For example, according to Fabre et al. (Reference Fabre, Sipp and Jacquin2006), the lower branch was designated as the C-family, while the upper one was labelled the V-family.

With viscosity, none of the eigenmodes have critical-layer singularities as they are regularised, and no eigenvalues lie exactly along the imaginary axis. In the non-normal region, the spectrum of eigenmodes, $\sigma _p^{\nu }$, that we have labelled as potential, occupies an area in the complex $\sigma$ plane that stretches out towards ${\rm Re}(\sigma )\rightarrow -\infty$ as $M$ increases. However, unlike the region shown in the schematic in figure 4, there is a gap between the upper bound of this numerical spectrum and the real axis on which the free-stream spectrum $\sigma _f^{\nu }$ is located. The reason is that we force solutions to vanish at radial infinity due to the decaying nature of the spectral basis elements. Therefore, the gap should be considered a peculiar product of our numerical method that excludes solutions with decay rates that are too slow in $r$, such as those with velocity decaying slower than $O(r^{-1})$ in $r$ and $\phi$, or $O(r^{-2})$ in $z$ in the far field; see (2.19) in § 2.3.

The fact that the numerically computed eigenvalues in $\sigma _p^{\nu }$ shift towards the left side of the complex $\sigma$ plane as $M$ increases coincides with Mao & Sherwin (Reference Mao and Sherwin2011). Additionally, the number of potential eigenmodes increases with increasing $M$. Up to the largest value of $M$ that we have explored, the numerical eigenvalues in $\sigma _p^{\nu }$ tend to emerge randomly. This random scattering can be understood as the spectrum's extremely high sensitivity to numerical errors even in the order of machine precision (see § 6.4.1).

On the other hand, we observe a moving branch of numerical eigenvalues attached to the left end of $\sigma _p^{\nu }$. They also never converge with respect to $M$, and the values of their $|{\rm Re}(\sigma )|$ increase rapidly. We explicitly label them as spurious because of their absolutely irregular spatial characteristics, as shown later in § 6.1.2. Although they are not removable, we can pull them away by setting $M$ to a large value.

Last but not least, we report the presence of two new distinct branches of viscous critical-layer eigenvalues, which are seen on the right side of the area containing the potential eigenvalues. The two branches of these eigenvalues, $\sigma _{c}^{\nu }$, as depicted in figure 4, converge to distinct loci. We distinguish $\sigma _{c}^{\nu }$ from $\sigma _{p}^{\nu }$ because of their unique bifurcating shape. Furthermore, this is the only part of the spectra in the non-normal region that approach fixed points at finite $M$, along with the discrete spectrum. As the name suggests, we argue that their associated eigenmodes are not only non-spurious but also physical, since they are the true viscous remnants of the inviscid critical-layer eigenmodes, as explained with details in § 6.1.4.

6.1.1. Discrete eigenmodes

Figure 13 presents three viscous discrete eigenmodes with respect to each base flow, whose spatial structures are inherited from the inviscid discrete eigenmodes displayed in figure 7. They remain non-singular throughout the whole radial domain. Viscosity only marginally affects the spatial structures of these eigenmodes compared with their inviscid counterparts, making the velocity components have slightly non-zero imaginary parts due to viscous perturbation. The number of wiggles in the eigenmodes still determines their spatial characteristics. Moreover, those eigenmodes with more wiggles near $r=0$ are more stable over time, i.e. $|{\rm Re}(\sigma )|$ increases. This phenomenon is physically justifiable since the spatial gradient of velocity components becomes steeper when the spacing between the wiggles is reduced, and viscous diffusion should, therefore, be more intensive. These eigenmodes are physical as they are regular, well-resolved solutions to the linearised Navier–Stokes equations on the $q$ vortex. They are typically characterised by modest wiggles that are spatially resolved, have rapid monotonous decay in $r$ and clearly correspond to the inviscid discrete eigenmodes associated with $\sigma _{d}^{0}$.

Figure 13. Radial velocity profiles of the viscous discrete eigenmodes associated with three smallest ${\rm Im} (\sigma )$ (a) for the Lamb–Oseen vortex ($q \rightarrow \infty$) in $(m, \kappa ) = (1, 1.0)$ and (b) for the strong swirling Batchelor vortex ($q = 4.0)$ in $(m, \kappa ) = (2, 3.0)$. The maximum of ${\rm Re}(r \tilde {u}_r)$ is normalised to unity. Here $M=400$ and $L=2.0$ are used. Comparing with the inviscid counterparts in figure 7, we note that viscosity only marginally affects these eigenmodes. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,10^5)$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,10^5)$.

6.1.2. Spurious eigenmodes

Two numerically computed eigenmodes that represent the viscous spurious eigenmodes are shown in figure 14. We have not observed any signs of convergence up to $M=400$. These eigenmodes are not spatially resolved, as evidenced by irregularly fast oscillations that alternate at every collocation point. It is apparent that they are neither analytically nor physically meaningful. Therefore, we will not perform an in-depth analysis of them.

Figure 14. Radial velocity profiles of a representative viscous spurious eigenmode (a) for the Lamb–Oseen vortex ($q \rightarrow \infty$) in $(m, \kappa ) = (1, 1.0)$ and (b) for the strong swirling Batchelor vortex ($q = 4.0)$ in $(m, \kappa ) = (2, 3.0)$. The maximum of ${\rm Re}(r \tilde {u}_r)$ is normalised to unity. Here $M=400$ and $L=2.0$ are used. Non-trivial and irregularly fast oscillations with alternating sign at every collocation point, as shown in each inset for magnification, manifest that they are spurious. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,10^5)$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,10^5)$.

6.1.3. Potential eigenmodes

Next, we examine the numerical eigenmodes associated with $\sigma _p^{\nu }$, or the potential eigenmodes. If we look at the randomly scattered eigenvalues while increasing $M$, it is possible to observe common spatial characteristics of these eigenmodes that are spatially resolved with a sufficiently large value of $M$, unlike the spurious family mentioned above. Figure 15 presents the three representative potential eigenmodes for each vortex case, using $M=400$. We note that we have selected eigenmodes whose smallest wiggle is captured with more than two collocation points to ensure that we validly discuss their common spatial features. These eigenmodes are characterised by excessive wiggles, resulting in slow radial decay rates (cf. Mao & Sherwin Reference Mao and Sherwin2011). They exhibit generally faster decay rates in time (i.e. larger $|{\rm Re}(\sigma )|$) than the discrete ones, as more wiggles demand steeper spatial gradients vulnerable to viscous diffusion.

Figure 15. Radial velocity profiles of three viscous potential eigenmodes (a) for the Lamb–Oseen vortex ($q \rightarrow \infty$) in $(m, \kappa ) = (1, 1.0)$ and (b) for the strong swirling Batchelor vortex ($q = 4.0)$ in $(m, \kappa ) = (2, 3.0)$. The maximum of ${\rm Re}(r \tilde {u}_r)$ is normalised to unity with the use of $M=400$ and $L=2.0$. The first and middle two potential eigenmodes exhibit similar ${\rm Re}(\sigma )$, and their number of major oscillations is comparable. The middle and last two eigenmodes have similar ${\rm Im}(\sigma )$, and their major oscillatory positions are similar. Each vertical dashed line indicates the critical-layer location $r_c$, which is estimated by setting each ${\rm Im}(\sigma )$ to $\sigma _c$ in (4.9). Each inset within a dashed box reveals small-amplitude wiggles where $r\tilde {u}_{r} \sim O(10^{-5})$ that persist at large $r$ even when the amplitude seems to be nearly zero, indicating their slow radial decay rates. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,10^5)$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,10^5)$.

The potential eigenmodes have wiggles that are usually concentrated roughly near the inviscid critical-layer singularity locations, which is estimated by setting their ${\rm Im}(\sigma )$ to $\sigma _c$ in (4.9). In other words, as noted by Mao & Sherwin (Reference Mao and Sherwin2011), they take the form of ‘wavepackets’, whose major oscillatory components are localised both in physical and spectral spaces. The correspondence between these ‘wavepackets’ and the inviscid critical-layer singularity locations leads us to posit that the potential eigenmodes originate from the viscous regularisation of the critical layers. From a mathematical standpoint, the introduction of the viscous term serves only to ensure their regularisation and does not impose any restrictions on their appearance following regularisation, such as thickness and wave amplitude. This may explain why potential eigenmodes exhibit various wavepacket widths at different locations.

In figure 15 the first and second eigenmodes have similar decay rates in time, i.e. ${\rm Re}(\sigma _1) \simeq {\rm Re}(\sigma _2)$, which relates to the fact that they also have a similar number of wiggles at their major oscillatory positions. On the other hand, the second and third eigenmodes have similar wave frequencies, i.e. ${\rm Im}(\sigma _2) \simeq {\rm Im}(\sigma _3)$, which means that their major oscillatory locations are close. As the number of wiggles increases, $|{\rm Re}(\sigma )|$ becomes large, and the major oscillatory structure extends to a wide range in $r$. This extension likely contributes to the retardation of radial decay rates, as the wiggles remain at large radii in small scales (see the insets in figure 15).

There are several noteworthy factors that should be pointed out regarding the spatial characteristics of these eigenmodes. Although they appear physical, they make it difficult to believe that they are the true viscous remnants of the inviscid critical-layer eigenmodes. First, potential eigenmodes’ wavepackets can have varying widths even at the same $Re$, indicating the absence of a clear scaling relationship between wavepacket widths and the important physical parameter $Re$. Second, it is challenging to identify a clear spatial similarity to the inviscid critical-layer eigenmodes. The typical radial decaying behaviour of the viscous eigenmodes appears slow and oscillatory, as shown in figure 15, in contrast to the inviscid critical-layer eigenmodes that exhibit monotonically rapid radial decay (see figure 8). We postulate that the rapid radial decaying behaviour in $\sigma _c^{0}$ must be sustained for its true viscous remnants since the viscous regularisation effect should be highly localised around the critical-layer singularity. Therefore, a subsequent question should arise as to which other eigenmodes in the non-normal region can be considered the true viscous remnants of $\sigma _c^{0}$. As the name suggests, we claim that the viscous critical-layer eigenmodes associated with $\sigma _c^{\nu }$ offer the answer, which we set forth in the following.

6.1.4. Viscous critical-layer eigenmodes

Figure 16(a) shows two viscous critical-layer eigenmodes of a Lamb–Oseen vortex with values of ${\rm Im}(\sigma _{c}^{\nu })$ that are within 6 % of each other. Due to the closeness of their eigenvalues and their similar appearances, we believe that they evolved from a pair of degenerate inviscid critical-layer eigenmodes. The eigenvalue of the eigenmode in the upper row of figure 16(a) is in the left branch, while the lower row is in the right branch of $\sigma _c^{\nu }$ in figures 4 and 12. We believe that the regions of large amplitude oscillations shown in the middle column of figure 16(a) are the true remnants of the inviscid critical layers.

Figure 16. Two viscous critical-layer eigenmodes with nearly identical ${\rm Im}(\sigma )$. (a) Radial component of the velocity eigenmode of the Lamb–Oseen vortex ($q \rightarrow \infty$) with $(m, \kappa ) = (1, 1.0)$ and (b) of the Batchelor vortex ($q = 4.0)$ with $(m, \kappa ) = (2, 3.0)$. The maximum of ${\rm Re}(r \tilde {u}_r)$ is normalised to unity. Here $M=400$ and $L=2.0$ are used. Each vertical dashed line indicates the location of the viscous critical layer estimated by setting ${\rm Im}(\sigma )$ equal to $\sigma _c$ in (4.9). These locations are nearly equal to the centroid of the magnitude of $r \tilde {u}_r$. Due to the similarity of the shape of small-amplitude structures in the right and left columns, where $r \tilde {u}_r \sim O(10^{-5})$, to the inviscid critical-layer eigenmodes (compare them with the middle column panels in figures 8(a) and 8(b), respectively), we hypothesise that these nearly degenerate viscous critical-layer eigenmodes are the viscous analogues of the inviscid two-fold degenerate critical-layer eigenmodes. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,10^5)$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,10^5)$.

The central locations of the critical layers in the middle column of figure 16(a), which we define as the centroid of the magnitude of $r\tilde {u}_r$, are nearly equal to the inviscid critical-layer singularity locations $r_c$, as estimated by setting ${\rm Im}(\sigma _c^{\nu })$ to $\sigma _c$ in (4.9). An important qualitative difference from the potential eigenmodes in figure 15 and the critical-layer eigenmodes is that in the radial regions outside the large amplitude oscillations the viscous critical-layer eigenmodes decay monotonically, while the decay of the potential eigenmodes is highly oscillatory. Figure 16(b) shows two eigenmodes of the Batchelor vortex, which have similar eigenvalues (differing by only 6 %). Their properties are similar to those of the eigenmode of the Lamb–Oseen vortex.

These numerical eigenmodes and eigenvalues exhibit good convergence with increasing $M$ and are spatially resolved. For physical relevance, it is worthwhile to investigate their structures outside the remnant critical layers. By normalising the oscillation amplitude in the remnant critical layer to be of order unity, we can identify small-scale perturbation structures outside the critical layer of $O(10^{-5})$ or less. We note the similarity in shape of these small-scale perturbations to the inviscid critical-layer eigenmodes of similar ${\rm Im}(\sigma )$ (see the middle column of panels in figure 8), where each part in $(0, r_c)$ and $(r_c, \infty )$ appears to be a scalar multiple of each side of the inviscid solutions (see § 5.3). This is one indication that the viscous critical-layer eigenmodes are truly inherited from the inviscid critical-layer eigenmodes. Note that viscosity has a profound influence on the structure of the eigenmode at radial locations inside the remnant critical layer, where it locally regularises the critical layer's singularity, but viscosity has only a marginal impact on the eigenmode at radial locations outside the remnant critical layer. Therefore, we expect the inviscid critical-layer eigenmodes (in figure 8) and the viscous critical-layer eigenmodes (in figure 16) to look similar in the regions outside the critical layer.

As can be seen in the viscous spectra, the decay rates in time of the viscous critical-layer eigenmodes are comparable to those of the viscous discrete eigenmodes, indicating that they can last for a relatively long time against viscous diffusion. Moreover, when comparing an eigenmode in the left branch with another in the right branch of $\sigma _c^{\nu }$, no notable structural difference is observed between them. This observation is further supported by the fact that these eigenmodes lie in the non-normal region. A more detailed analysis of the viscous critical-layer eigenmodes is presented later in this paper, dealing with $L$ and including the viscous remnant critical layers conforming to the $Re^{-1/3}$ scaling law (see § 6.2) and the continuity of the viscous critical-layer spectrum $\sigma _c^{\nu }$ (see § 6.4).

If $\sigma _c^{\nu }$ is the truly regularised descendant of $\sigma _c^{0}$ with the correct critical-layer thickness, an important question that remains to be answered is how the spectrum of a single straight line bifurcates into two distinct branches. This bifurcation is physically meaningful because the separation of the branches, or, equivalently, the difference in ${\rm Re}(\sigma )$, is significantly larger compared with the extent of purely numerical error at the same level of $M$, such as the eigenvalue difference found in the pairing phenomenon in $\sigma _c^{0}$ (see figure 11). Recall that there exist numerous singular, degenerate eigenmodes associated with the same eigenvalue due to the critical-layer singularities in $\sigma _c^{0}$. We can infer that the viscous effect perturbs these two-fold degenerate singular eigenmodes and splits them into two regularised eigenmodes with marginally different eigenvalues. Hence, we expect that the emergence of $\sigma _c^{\nu }$ in two bifurcating curvy branches is not accidental but explicable by means of perturbation theory dealing with two-fold degeneracy (Sakurai & Napolitano Reference Sakurai and Napolitano2021, pp. 300–305).

It is worth discussing why $\sigma _c^{\nu }$ was not distinguished by previous researchers. When we compare our numerical method with that of Mao & Sherwin (Reference Mao and Sherwin2011) or Bölle et al. (Reference Bölle, Brion, Robinet, Sipp and Jacquin2021), we see that they truncated the radial domain at a large but finite $r$ and applied a homogeneous boundary condition there. In contrast, our method essentially involves the entire radial domain $0 \leq r < \infty$, and each basis function $P_{L_m}^n(r)$ obeys the boundary conditions that we want to apply. As a result, our truncated spectral sums, expanded by $P_{L_m}^n(r)$ as the basis elements of the Galerkin method, implicitly and exactly impose the boundary conditions on the solutions, regardless of the value of $M$ used. However, the boundary condition at $r \rightarrow \infty$ is only approximately satisfied by the others. Considering the sensitivity of the numerical spectra to numerical errors (see § 6.4.1), the truncation is likely to impede the numerical convergence of $\sigma _c^{\nu }$ because an approximate far-field radial boundary condition introduces errors. For instance, in the numerical spectrum plot provided by Mao & Sherwin (Reference Mao and Sherwin2011, p. 8), we can see faint traces of the two bifurcating branches at the location of $\sigma _c^{\nu }$ found in our results. Nonetheless, the results were substantially disturbed with respect to the radial truncation as well as the number of spectral elements, and the authors could not distinguish them from $\sigma _p^{\nu }$.

6.2. Optimal choice of $L$ to resolve the viscous critical layers

One of our goals is to accurately compute the viscous critical-layer modes. Clearly, we should use the largest $M$ (with $N \equiv M+2$) that our computational budget allows, which in this analysis is $M=400$. We are interested in finding all of the viscous critical-layer eigenmodes and eigenvalues, not just one, nor are we interested in finding them one at a time. Unlike previous studies that looked at individual eigemmodes and stretched the radial domain locally around the location of that eigenmode's critical layer to maximise the resolution there (e.g. Le Dizés & Lacaze Reference Le Dizés and Lacaze2005), our numerical method is designed for a fixed $Re, m$ and $\kappa$ to compute the entire radial domain for all of the eigenmodes, regardless of the locations of their critical layers, using the same radial collocation points.

Choosing a small value of $L$ is advantageous because the spatial resolution of our method is $\varDelta = L/(M+2)$ (see (3.11)), and we need to have $\varDelta$ smaller than the critical-layer thickness to resolve it. However, only half of the collocation points lie in the vast range between $L$ and infinity, so eigenmodes with critical-layer radii with $r_c > L$ will have few collocation points (if any) within their critical layers and, therefore, be spatially under-resolved. The optimal value of $L$, denoted $L_{opt}$, must be a ‘Goldilocks’ value: not too big or too small. Figure 17 demonstrates another reason why $L_{opt}$ is a ‘Goldilocks’ value. The figure displays the eigenvalues in the imaginary plane for $Re = 10^5$ with three different values of $L$. The left column of the figure represents a scenario with small $L$, where viscous critical layer eigenmodes with small $r_c$ (and, therefore, large $|\sigma _c^{\nu }|$; see figure 5) are spatially well resolved. However, eigenmodes with small values of $|\sigma _c^{\nu }|$ and large $r_c$ are not adequately resolved. The right column of the figure shows the case with a large $L$. In this case, only eigenmodes with small $|\sigma _c^{\nu }|$ and large $r_c$ are well resolved.

Figure 17. Changes of numerical viscous spectra (a) for the Lamb–Oseen vortex ($q\rightarrow \infty$) in $(m,\kappa )=(1,1.0)$ and (b) for the strong swirling Batchelor vortex ($q=4.0$) in $(m,\kappa ) = (2,3.0)$ with respect to three different $L$ values. Here $M$ is fixed at $400$ and $N=M+2$. If we aim to optimally resolve the critical-layer spectrum, we should appropriately tune $L$ to find a balance between (left) the expansion of the high-resolution region $0 \le r < L$, and (right) the deterioration of the overall resolution represented by $\varDelta \sim O(L)$. The middle one shows the optimal $L$, denoted $L_{opt}$, which minimises the emergence of the numerical potential spectrum. Thus, most numerical eigenvalues in the non-normal region belong to the viscous critical-layer eigenvalues. See supplementary movie 2 for animation. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,10^5)$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,10^5)$.

Nevertheless, figure 17 reveals another reason for the ‘Goldilocks’ behaviour. The panels in the left column exhibit a clear separation between the potential eigenvalues $\sigma _p^{\nu }$ and the two new branches of viscous critical-layer eigenvalues $\sigma _c^{\nu }$. As $L$ increases, the potential eigenvalues shift towards the right in the complex plane (middle column). When $L$ becomes sufficiently large, the potential eigenvalues become intertwined with those of the viscous critical-layer eigenmodes (right column), and the latter set of eigenmodes are no longer well-resolved spatially.

Upon detailed examination of the viscous critical-layer eigenmodes of the Lamb–Oseen vortex with $(m, \kappa ) = (1, 1.0)$ and the Batchelor vortex with $(q = 4.0)$ and $(m, \kappa ) = (2, 3.0)$, with $M = 400$ and $Re = 10^5$, we found that $\varDelta = L/(M+2)$ is just small enough to resolve the viscous critical-layer thicknesses when $L = 4.0$ and $2.5$, respectively. Figure 17 illustrates that these values of $L$ also represent the maximum values where the eigenvalues of the potential eigenmodes remain distinct from those of the viscous critical-layer eigenmodes. Thus, we believe that these values of $L$ are the ‘Goldilocks’ values: large enough to maximise the region $0 \le r_c < L$, providing a sufficient number of collocation points to resolve the eigenmodes and small enough that $\varDelta = L/(M+2)$ adequately resolves the critical-layer thicknesses. Our procedure for determining the optimal value $L_{opt}$ is similar to how we found the optimal $L$ for resolving the inviscid critical-layer eigenvalues $\sigma _c^{0}$ in $\S$ 5.2.

1. (i) Start with $L$ of order unity (i.e. the core radius of the unperturbed aircraft wake vortex), and increase $L$ to expand the high-resolution region $0 \le r_c < L$.

2. (ii) Stop increasing $L$ just before the spatial resolution is so poor that the $\sigma _p^{\nu }$ and $\sigma _c^{\nu }$ eigenvalues intertwine as shown in the middle panels of figure 17.

6.3. Use of $L_{opt}$ to find the scaling behaviour of the critical-layer thickness with $Re$

We hypothesise that the values of $L$ at which the potential and viscous critical-layer eigenvalues intermingle in figure 17 and where $L/(M+2)$ just barely resolves the critical-layer thickness are the same for all $Re$. We believe that the loss of numerical spatial resolution causes the two families of eigenvalues to become non-distinct from one another. To partially test this hypothesis, we calculate $L_{opt}$ using the two-step procedure mentioned above, using the data in figures 17 and 18. We then assume that $\varDelta _{opt} \equiv L_{opt}/(M+2)$ represents the critical-layer thickness. Plotting $\varDelta _{opt}$ as a function of $Re$ in figure 19 demonstrates that the critical-layer thickness (if our hypothesis is correct) scales approximately as $Re^{-1/3}$. This scaling agrees with previous analyses using asymptotic expansions (Maslowe Reference Maslowe1986; Le Dizès Reference Le Dizès2004).

Figure 18. Numerical viscous spectra with $L_{opt}$ at $Re = 10^4$ and $10^3$ (a) for the Lamb–Oseen vortex ($q\rightarrow \infty$) in $(m,\kappa )=(1,1.0)$ and (b) for the strong swirling Batchelor vortex ($q=4.0$) in $(m,\kappa ) = (2,3.0)$. Here $M$ is fixed at $400$ and $N=M+2$. Results are shown for (a) $(m, \kappa, q)=(1,1.0,\infty )$ and (b) $(m, \kappa, q)=(2,3.0,4.0)$.

Figure 19. The optimal numerical resolution $\varDelta _{opt} \equiv 2L_{opt}/(M+2)$, at fixed $M=400$, to resolve the critical-layer spectrum with respect to $Re$. The trend indicates $\varDelta _{opt} \propto Re^{-1/3}$. The presented cases of $Re = 10^3, 10^4$ and $10^5$ for each vortex can be found in figures 17 and 18.

6.4. Continuity in the viscous critical-layer spectrum

6.4.1. Pseudospectral analysis

Finding the pseudospectra of the viscous operator ${\boldsymbol{\mathsf{L}}}_{m\kappa }$, we can obtain evidence that the spectra $\sigma _p^{\nu }$ and $\sigma _c^{\nu }$ fill the continuous region in the complex $\sigma$ plane, as depicted in the schematic in figure 4. According to Mao & Sherwin (Reference Mao and Sherwin2011), the $\varepsilon$-pseudospectra around the potential and critical-layer eigenvalues seem to enclose the entire area when $\varepsilon$ is small, as shown in figure 20. In addition, we present the $\varepsilon$-pseudospectrum with $\varepsilon$ as small as $10^{-14}$, which is much smaller than the values used by Mao & Sherwin (Reference Mao and Sherwin2011) or Bölle et al. (Reference Bölle, Brion, Robinet, Sipp and Jacquin2021). Therefore, we believe that our observation provides strong empirical support for the continuity of the non-normal region that we have numerically resolved.

Figure 20. The $\varepsilon$-pseudospectrum bounds of $\varepsilon = 10^{-14}, 10^{-8}$ and $10^{-2}$ with respect to ${\boldsymbol{\mathsf{L}}}_{m\kappa }$ at $Re=10^4$ (a) for the Lamb–Oseen vortex ($q\rightarrow \infty$) in $(m,\kappa ) = (1,1.0)$ and (b) for the strong swirling Batchelor vortex ($q=4.0$) in $(m,\kappa ) = (2,3.0)$. To construct the matrix, we use $M=400$ and $N=M+2$. Here $L$ is optimally chosen. We can infer from their formation which part of the spectra is continuous and how big the maximum transient growth is. Results are shown for (a) $(m, \kappa, q, Re)=(1,1.0,\infty,10^4)$ and (b) $(m, \kappa, q, Re)=(2,3.0,4.0,10^4)$.

Furthermore, based on the alternative statement of the pseudospectra given by Trefethen & Embree (Reference Trefethen and Embree2005), any point in the $\varepsilon$-pseudospectra of ${\boldsymbol{\mathsf{L}}}_{m\kappa }$ can be on the spectrum of ${\boldsymbol{\mathsf{L}}}_{m\kappa } + {\boldsymbol{\mathsf{E}}}$ for some small disturbance ${\boldsymbol{\mathsf{E}}}$ where $\lVert {\boldsymbol{\mathsf{E}}} \rVert < \varepsilon$. Since $\varepsilon = 10^{-14}$ is almost comparable to the double-precision machine arithmetic used in modern computing, one possible explanation for the random scattering of the numerical eigenvalues in the numerical representation of $\sigma _{p}^{\nu }$ is that they are perturbed by machine-dependent precision errors serving as ${\boldsymbol{\mathsf{E}}}$.

As an aside, we observe that the $\varepsilon$-pseudospectrum of $\varepsilon = 10^{-2}$ protrudes into the right half-plane of the complex $\sigma$ plane, as shown in figure 20. It is well known that the supremum of the real parts of $\sigma \in \sigma _{\varepsilon }({\boldsymbol{\mathsf{L}}}_{m\kappa })$, denoted $\alpha _{\varepsilon }$ and referred to as the $\varepsilon$-pseudospectral abscissa (Trefethen & Embree Reference Trefethen and Embree2005), is relevant to the lower bound of the maximum transient growth of the stable system with an arbitrary initial state of $\boldsymbol {x} = \boldsymbol {x}_0$ where $\lVert \boldsymbol {x}_0 \rVert = 1$,

(6.1)\begin{equation} \frac{\partial \boldsymbol{x}}{\partial t} = {\boldsymbol{\mathsf{L}}}_{m\kappa} \boldsymbol{x}. \end{equation}

The supremum of $\alpha _{\varepsilon } / \varepsilon$ in $\varepsilon > 0$ determines the lower bound of the maximum transient growth of the system (Apkarian & Noll Reference Apkarian and Noll2020), or

(6.2)\begin{equation} \sup_{t \ge 0}\lVert {\rm e}^{{\boldsymbol{\mathsf{L}}}_{m\kappa}t} \rVert \ge \sup_{\varepsilon > 0} \frac{\alpha_{\varepsilon ({\boldsymbol{\mathsf{L}}}_{m\kappa})}}{\varepsilon} . \end{equation}

The fact that the $\varepsilon$-pseudospectral abscissa of $\varepsilon = 10^{-2}$ occurs in the frequency band coinciding with the critical-layer spectrum implies the significance of this spectrum in regards to the transient vortex growth, which needs more investigation in further studies.

6.4.2. Loci of the numerical spectra

One issue with pseudospectra is that they cannot provide non-normal eigenmodes corresponding to each eigenvalue point in the continuum. Instead, pseudomodes can be constructed in association with pseudospectra as an approximation of the eigenmodes, which were introduced and described by Trefethen & Embree (Reference Trefethen and Embree2005). Unfortunately, pseudomodes do not necessarily satisfy the exact governing equations and boundary conditions (see Mao & Sherwin Reference Mao and Sherwin2011, p. 11).

In our numerical method it is possible to find critical-layer eigenmodes whose spatial structures continuously vary by fine tuning $L$. Recalling the role of $L$ (see § 3), we know that it changes the entire $P_{L_m}^n(r)$ in the basis function set. If we replace $L$ and solve the eigenvalue problem again, we can expect the eigenmodes generated from a new $L$ not necessarily to be identical to the eigenmodes generated from an old $L$. Moreover, if this parametric change occurs in parts of the spectra where numerical convergence with respect to $M$ is ensured, including $\sigma _d^{\nu }$ and $\sigma _{c}^{\nu }$, the loci of them with respect to $L$ should genuinely reflect the analytic spectra.

Based on the idea described above, we create the loci of the numerical spectra with respect to $L$ for the Lamb–Oseen vortex case, where $(m, \kappa, q) = (1, 1.0, \infty )$ at $Re = 10^4$ in figure 21(a). To draw the loci, the viscous eigenvalue problem is solved multiple times with fine tuning $L$ from $8.3$ to $8.7$ with $M = 400$, where both $\sigma _d^{\nu }$ and $\sigma _c^{\nu }$ are found to be well resolved. The other parts of the spectra, including $\sigma _p^{\nu }$, are excluded due to no clear convergence with respect to $M$. That being said, we note that the loci of $\sigma _p^{\nu }$ with varying $L$ sweep over the shaded area depicted in the schematic in figure 4.

Figure 21. (a) Loci of the numerical spectra for the Lamb–Oseen vortex ($q\rightarrow \infty$) in $(m,\kappa )= (1,1.0)$ obtained by fine tuning $L$ from $8.3$ to $8.7$, where $Re = 10^4$, and (b) three viscous critical-layer eigenmodes that marginally vary, all of which are obtained from different $L$. Unlike the discrete spectrum that does not change with respect to $L$, the critical-layer spectrum is continuously filled by numerical eigenvalues associated with valid critical-layer eigenmodes. (a) Loci of the numerical spectra by fine tuning L and (b) numerical critical-layer eigenmodes varying marginally.

As for $\sigma _d^{\nu }$, its locus is completely invariant against changes in $L$. It makes sense because there is no chance to find an intermediate form of two discrete eigenmodes. The locus of $\sigma _d^{\nu }$ remaining discrete rather strengthens our method's robustness for any $L$. On the contrary, the locus of $\sigma _c^{\nu }$ is notably different from that of $\sigma _d^{\nu }$; as $L$ changes, the eigenvalue points on two branches of $\sigma _c^{\nu }$ also move and eventually fill in two distinct curves as in figure 4. In figure 21(b) it can be confirmed that the critical-layer eigenmodes with slightly different eigenvalues, having only a marginal difference in their spatial structures, are obtained from varying $L$. By comparing the two loci of $\sigma _d^{\nu }$ and $\sigma _c^{\nu }$, we can conclude the continuity of the critical-layer spectrum.

For reference, we report the polynomial fitting results up to sixth order of some loci of $\sigma _c^{\nu }$ among what we have explored. In the case $(m, \kappa, q) = (1, 1.0, \infty )$ at $Re = 10^4$, the left and right branches of $\sigma _c^{\nu }$ in the complex $\sigma$ plane are fitted as

(6.3)\begin{align} \sigma_r &=- (1.905 \times 10^1) \,{\cdot}\, \sigma_i^6 - (4.562 \times 10^1)\,{\cdot}\, \sigma_i^5 - (4.138 \times 10^1)\,{\cdot}\, \sigma_i^4 \nonumber\\ &\quad - (1.741 \times 10^1)\,{\cdot}\, \sigma_i^3 - (2.761 \times 10^0)\,{\cdot}\, \sigma_i^2 + (5.348 \times 10^{-1})\,{\cdot}\, \sigma_i \end{align}

and

(6.4)\begin{align} \sigma_r &=- (4.682 \times 10^0)\,{\cdot}\, \sigma_i^6 - (8.233 \times 10^0)\,{\cdot}\, \sigma_i^5 - (6.816 \times 10^0)\,{\cdot}\, \sigma_i^4 \nonumber\\ &\quad - (3.243 \times 10^0)\,{\cdot}\, \sigma_i^3 - (2.636 \times 10^{-1})\,{\cdot}\, \sigma_i^2 + (6.108 \times 10^{-1})\,{\cdot}\, \sigma_i, \end{align}

where $\sigma _r$ and $\sigma _i$ indicate the real and imaginary parts of $\sigma$, respectively. In the case $(m, \kappa, q) = (2, 3.0, 4.0)$ at $Re = 10^4$, the left and right branches of $\sigma _c^{\nu }$ are fitted as

(6.5)\begin{align} \sigma_r &=+ (1.071 \times 10^{-2})\,{\cdot}\, \sigma_i^6 + (2.553 \times 10^{-2})\,{\cdot}\, \sigma_i^5 - (9.022 \times 10^{-2})\,{\cdot}\, \sigma_i^4 \nonumber\\ &\quad - (3.417 \times 10^{-1})\,{\cdot}\, \sigma_i^3 - (1.906 \times 10^{-1})\,{\cdot}\, \sigma_i^2 + (4.764 \times 10^{-1})\,{\cdot}\, \sigma_i \end{align}

and

(6.6)\begin{align} \sigma_r &=+ (7.098 \times 10^{-2})\,{\cdot}\, \sigma_i^6 + (4.052 \times 10^{-1})\,{\cdot}\, \sigma_i^5 + (8.111 \times 10^{-1})\,{\cdot}\, \sigma_i^4 \nonumber\\ &\quad + (6.323 \times 10^{-1})\,{\cdot}\, \sigma_i^3 + (2.508 \times 10^{-1})\,{\cdot}\, \sigma_i^2 + (4.833 \times 10^{-1})\,{\cdot}\, \sigma_i. \end{align}

We will work on the analytic formulation of $\sigma _c^{\nu }$ to better understand the bifurcation in future studies. These fitting forms will be considered for comparison and validation.